From 7b52d75ff14096d05bdc06a1cc1de2ee6997da22 Mon Sep 17 00:00:00 2001 From: zvezdochiot Date: Thu, 25 Nov 2021 22:30:15 +0300 Subject: [PATCH] 1.0.3: ellipsoid --- .gitignore | 1 + Makefile | 1 + README.md | 4 +- examples/GRS80.txt | 11 --- examples/testpoints_blh.txt | 11 +++ man/man1/helmblhtoxyz.1 | 26 ++++-- man/man1/helmdiff3d.1 | 4 +- man/man1/helmert3d.1 | 4 +- man/man1/helmparms3d.1 | 4 +- src/helmblhtoxyz.c | 164 +++++++++++++++++++++++++++--------- src/helmdiff3d.c | 15 +++- src/helmert3d.c | 17 ++-- src/helmert3d.h | 2 +- src/helmparms3d.c | 34 +++++--- 14 files changed, 212 insertions(+), 86 deletions(-) create mode 100644 examples/testpoints_blh.txt diff --git a/.gitignore b/.gitignore index d6696da..bef4294 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ libsvdm.a helmert3d helmparms3d helmdiff3d +helmblhtoxyz diff --git a/Makefile b/Makefile index c6b30d5..3ff8ff1 100644 --- a/Makefile +++ b/Makefile @@ -56,6 +56,7 @@ install: $(INSTALL) -m 0644 LICENSE $(DOCS) $(INSTALL) -m 0644 $(EXAMPS)/testpoints_src.txt $(DOCS)/$(EXAMPS) $(INSTALL) -m 0644 $(EXAMPS)/testpoints_dest.txt $(DOCS)/$(EXAMPS) + $(INSTALL) -m 0644 $(EXAMPS)/testpoints_blh.txt $(DOCS)/$(EXAMPS) $(INSTALL) -m 0644 $(EXAMPS)/GRS80.txt $(DOCS)/$(EXAMPS) $(INSTALL) -m 0644 $(MANS)/helmparms3d.1 $(PREFIX)/$(MANS) $(INSTALL) -m 0644 $(MANS)/helmert3d.1 $(PREFIX)/$(MANS) diff --git a/README.md b/README.md index f311bae..7e64eb4 100644 --- a/README.md +++ b/README.md @@ -91,10 +91,10 @@ Helmert parameter file format: XYZ data file format: ``` - x[1] y[1] z[1] + X[1] Y[1] Z[1] .. .. .. .. .. .. - x[n] y[n] z[n] + X[n] Y[n] Z[n] ``` ---- diff --git a/examples/GRS80.txt b/examples/GRS80.txt index 3b0801c..e0dff87 100644 --- a/examples/GRS80.txt +++ b/examples/GRS80.txt @@ -1,12 +1 @@ GRS80 6378137.0 6356752.3141 -56.32258375 21.22652722 0 -56.4935 21.783 0 -56.61155556 22.93708333 0 -56.48455556 23.55005556 0 -56.48191667 24.31558333 0 -56.59155556 25.15272222 0 -56.25533333 25.54025 0 -56.01013889 26.22444444 0 -55.69546942 26.61753336 0 -56.94747222 24.10877528 0 -56.94747556 24.1087875 0 diff --git a/examples/testpoints_blh.txt b/examples/testpoints_blh.txt new file mode 100644 index 0000000..aa76c57 --- /dev/null +++ b/examples/testpoints_blh.txt @@ -0,0 +1,11 @@ +56.32258375 21.22652722 0 +56.4935 21.783 0 +56.61155556 22.93708333 0 +56.48455556 23.55005556 0 +56.48191667 24.31558333 0 +56.59155556 25.15272222 0 +56.25533333 25.54025 0 +56.01013889 26.22444444 0 +55.69546942 26.61753336 0 +56.94747222 24.10877528 0 +56.94747556 24.1087875 0 diff --git a/man/man1/helmblhtoxyz.1 b/man/man1/helmblhtoxyz.1 index 30e9435..7d5ee6d 100644 --- a/man/man1/helmblhtoxyz.1 +++ b/man/man1/helmblhtoxyz.1 @@ -1,19 +1,30 @@ -.TH "helmblhtoxyz" 1 1.0.2 "25 Nov 2021" "User Manual" +.TH "helmblhtoxyz" 1 1.0.3 "26 Nov 2021" "User Manual" .SH NAME helmblhtoxyz .SH DESCRIPTION -3D calculate from B-L-H to X-Y-Z for Helmert +3D calculate from B-L-H to X-Y-Z for Helmert and back .SH SYNOPSIS -helmblhtoxyz blh_src_infilename [xyz_diff_outfilename] +helmblhtoxyz commands blh|xyz_src_infilename ellipsoid_src_infilename [blh|xyz_diff_outfilename] + +.SH COMMANDS +.TP +xyz +convert B-L-H to X-Y-Z +.TP +blh +convert X-Y-Z to B-L-H .SH FILES .TP -blh data file format: +ellipsoid data file format: name a b +.TP +blh data file format: + B[1] L[1] H[1] .. .. .. .. .. .. @@ -24,12 +35,14 @@ xyz data file format: X[1] Y[1] Z[1] .. .. .. .. .. .. - H[n] Y[n] Z[n] + X[n] Y[n] Z[n] .SH EXAMPLE export EXAMP=/usr/share/doc/helmert3d/examples/ - helmblhtoxyz $EXAMP/GRS80.txt XYH.txt + helmblhtoxyz xyz $EXAMP/testpoints_blh.txt $EXAMP/GRS80.txt xyh.txt + helmblhtoxyz blh xyh.txt $EXAMP/GRS80.txt blh.txt + helmdiff3d examples/testpoints_blh.txt blh.txt .SH COPYRIGHT helmert3d: @@ -57,4 +70,3 @@ Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. .SH CONTACT Website: https://github.com/dr-ni/helmert3d - diff --git a/man/man1/helmdiff3d.1 b/man/man1/helmdiff3d.1 index f5b579f..c821169 100644 --- a/man/man1/helmdiff3d.1 +++ b/man/man1/helmdiff3d.1 @@ -1,4 +1,4 @@ -.TH "helmdiff3d" 1 1.0.2 "25 Nov 2021" "User Manual" +.TH "helmdiff3d" 1 1.0.3 "26 Nov 2021" "User Manual" .SH NAME helmdiff3d @@ -7,7 +7,7 @@ helmdiff3d 3D Helmert Diff-Test Tool .SH SYNOPSIS -helmdiff3d [xyz_src_infilename] [xyz_dest_infilename] [xyz_diff_outfilename] +helmdiff3d xyz_src_infilename xyz_dest_infilename [xyz_diff_outfilename] .SH FILES .TP diff --git a/man/man1/helmert3d.1 b/man/man1/helmert3d.1 index 26e26a4..13a3e68 100644 --- a/man/man1/helmert3d.1 +++ b/man/man1/helmert3d.1 @@ -1,4 +1,4 @@ -.TH "helmert3d" 1 1.0.2 "25 Nov 2021" "User Manual" +.TH "helmert3d" 1 1.0.3 "26 Nov 2021" "User Manual" .SH NAME helmert3d @@ -17,7 +17,7 @@ where s = scale factor .SH SYNOPSIS -helmert3d [xyz_src_infilename] [param_infilename] [xyz_outfilename] +helmert3d xyz_src_infilename param_infilename [xyz_outfilename] .SH FILES .TP diff --git a/man/man1/helmparms3d.1 b/man/man1/helmparms3d.1 index 88508d7..7f1565d 100644 --- a/man/man1/helmparms3d.1 +++ b/man/man1/helmparms3d.1 @@ -1,4 +1,4 @@ -.TH "helmparms3d" 1 1.0.2 "25 Nov 2021" "User Manual" +.TH "helmparms3d" 1 1.0.3 "26 Nov 2021" "User Manual" .SH NAME helmparms3d @@ -17,7 +17,7 @@ where s = scale factor .SH SYNOPSIS -helmparms3d [xyz_src_infilename] [xyz_dest_infilename] [parms_outfilename] +helmparms3d xyz_src_infilename xyz_dest_infilename [parms_outfilename] .SH FILES .TP diff --git a/src/helmblhtoxyz.c b/src/helmblhtoxyz.c index e2c4b15..6ad70b1 100644 --- a/src/helmblhtoxyz.c +++ b/src/helmblhtoxyz.c @@ -50,87 +50,173 @@ static size_t get_m_size(char *filename) int main(int argc, char* argv[]) { FILE *ifile; + FILE *cfile; FILE *ofile; char *ifilename = ""; - char *ofilename = "out.xyz"; - char ibuf[256], name[128]; - double a = 6378137.0, b = 6356752.3142; - double bc, bs, lc, ls, ta, tb, nb, lt, bt, ht, xt, yt, zt; - size_t s = 0, i = 0; + char *cfilename = ""; + char *ofilename = "out.txt"; + char *command = "xyz"; + char ibuf[256], cbuf[256], name[128]; + double a = 6378137.0, b = 6356752.3142, ro, roi; + double a2, b2, bc, bs, lc, ls, ta, tb, nb, lt, bt, ht, xt, yt, zt; + double x2, y2, z2, r2, r, e2, et2, kf, kg, kc, ks, ks2, kp, kq, r0, re, re2, ku, kv, kz, z0; + + + size_t l = 0; int stat = 0; - fprintf(stdout,"\n*******************************\n"); - fprintf(stdout, "* helmblhtoxyz v%s *\n",VERS); - fprintf(stdout, "* (c) U. Niethammer 2020 *\n"); - fprintf(stdout, "*******************************\n"); + fprintf(stdout, "\n*******************************\n"); + fprintf(stdout, "* blh2xyz v%s *\n",VERS); + fprintf(stdout, "* (c) U. Niethammer 2020 *\n"); + fprintf(stdout, "*******************************\n"); - if(argc < 2) + if(argc < 4) { - fprintf(stdout,"\nSyntax: %s blh_src_infilename [xyz_dest_infilename]\n\n",argv[0]); + fprintf(stdout,"\nSyntax: %s commands blh|xyz_src_infilename ellipsoid_src_infilename [blh|xyz_dest_infilename]\n\n",argv[0]); + fprintf(stdout,"commands:\n"); + fprintf(stdout," xyz - convert B-L-H to X-Y-Z\n"); + fprintf(stdout," blh - convert X-Y-Z to B-L-H\n\n"); + fprintf(stdout,"ellipsoid file format:\n"); + fprintf(stdout," Name a b\n\n"); fprintf(stdout,"blh data file format:\n"); - fprintf(stdout," Name a b\n B[1] L[1] H[1]\n .. .. ..\n .. .. ..\n B[n] L[n] H[n]\n\n"); + fprintf(stdout," B[1] L[1] H[1]\n .. .. ..\n .. .. ..\n B[n] L[n] H[n]\n\n"); fprintf(stdout,"xyz data file format:\n"); fprintf(stdout," X[1] Y[1] Z[1]\n .. .. ..\n .. .. ..\n X[n] Y[n] Z[n]\n\n"); exit(EXIT_FAILURE); } + command = argv[1]; + if (!strcmp(command,"xyz")) + { + fprintf(stdout,"B-L-H -> X-Y-Z...\n"); + } + else if (!strcmp(command,"blh")) + { + fprintf(stdout,"X-Y-Z -> B-H-L...\n"); + } + else + { + fprintf(stderr,"Unknow command %s\n",command); + exit(EXIT_FAILURE); + } fprintf(stdout,"Reading points...\n"); - ifilename = argv[1]; - s = get_m_size(ifilename) - 1; - fprintf(stdout,"Found %lu points\n",(unsigned long)s); + ifilename = argv[2]; + l = get_m_size(ifilename); + fprintf(stdout,"Found %lu points\n",(unsigned long)l); ifile = fopen( ifilename, "r"); if(ifile == NULL) { fprintf(stderr,"Error opening %s\n",ifilename); exit(EXIT_FAILURE); } - if(argc > 2) + cfilename = argv[3]; + l = get_m_size(cfilename); + fprintf(stdout,"Found ellipsoid\n"); + cfile = fopen( cfilename, "r"); + if(cfile == NULL) { - ofilename = argv[2]; + fprintf(stderr,"Error opening %s\n",cfilename); + exit(EXIT_FAILURE); + } + if(argc > 4) + { + ofilename = argv[4]; + ofile = fopen( ofilename, "w"); } - ofile = fopen( ofilename, "w"); + else + ofile = stdout; + if(ofile == NULL) { fprintf(stderr,"Error writing %s\n",ofilename); exit(EXIT_FAILURE); } - fprintf(stdout,"Starting calculate...\n"); - while(fgets( ibuf, 128, ifile)!=NULL) + fprintf(stdout,"Starting calculation...\n"); + fgets( cbuf, 256, cfile); + stat = sscanf( cbuf, "%s %lf %lf", name, &a, &b); + if(stat != 3) { - if (i == 0) + fprintf(stderr,"Error wrong data format in %s\n",cfilename); + exit(EXIT_FAILURE); + } + fprintf(stdout,"Ellipsoid %s, a = %lf, b = %lf\n", name, a, b); + a2 = a * a; + b2 = b * b; + ro = M_PI / 180.0; + roi = 180.0 / M_PI; + if (!strcmp(command,"xyz")) + { + while(fgets( ibuf, 256, ifile)!=NULL) { - stat = sscanf( ibuf, "%s %lf %lf", &name, &a, &b); + stat=sscanf( ibuf, "%lf %lf %lf", &bt, <, &ht); if(stat != 3) { fprintf(stderr,"Error wrong data format in %s\n",ifilename); exit(EXIT_FAILURE); } - fprintf(stdout,"Ellipsoid %s, a = %lf, b = %lf\n", name, a, b); + bc = cos(ro * bt); + bs = sin(ro * bt); + lc = cos(ro * lt); + ls = sin(ro * lt); + ta = a * bc; + tb = b * bs; + nb = a2 / sqrt(ta * ta + tb * tb); + xt = (nb + ht) * bc * lc; + yt = (nb + ht) * bc * ls; + zt = (b2 / a2 * nb + ht) * bs; + fprintf(ofile,"%lf %lf %lf\n", xt , yt , zt); } - else + } + else if (!strcmp(command,"blh")) + { + while(fgets( ibuf, 256, ifile)!=NULL) { - stat=sscanf( ibuf, "%lf %lf %lf", &bt, <, &ht); + stat=sscanf( ibuf, "%lf %lf %lf", &xt, &yt, &zt); if(stat != 3) { fprintf(stderr,"Error wrong data format in %s\n",ifilename); exit(EXIT_FAILURE); } - bc = cos(M_PI / 180.0 * bt); - bs = sin(M_PI / 180.0 * bt); - lc = cos(M_PI / 180.0 * lt); - ls = sin(M_PI / 180.0 * lt); - ta = a * bc; - tb = b * bs; - nb = (a * a) / sqrt(ta * ta + tb * tb); - xt = (nb + ht) * bc * lc; - yt = (nb + ht) * bc * ls; - zt = (b * b / a / a * nb + ht) * bs; - fprintf(ofile,"%lf %lf %lf\n", xt , yt , zt); + x2 = xt * xt; + y2 = yt * yt; + r2 = x2 + y2; + r = sqrt(r2); + e2 = 1.0 - b2 / a2; + et2 = (a2 - b2) / b2; + z2 = zt * zt; + kf = 54.0 * b2 * z2; + kg = r2 + (1 - e2) * z2 - e2 * (a2 - b2); + kc = e2 * e2 * kf * r2 / (kg * kg * kg); + ks = cbrt(1.0 + kc + sqrt(kc * kc + 2.0 * kc)); + ks2 = ks + 1.0 / ks + 1.0; + kp = kf / (3.0 * ks2 * ks2 * kg * kg); + kq = sqrt(1.0 + 2.0 * e2 * e2 * kp); + r0 = -kp * e2 * r / (1.0 + kq) + sqrt(0.5 * a2 *(1.0 + 1.0 / kq) - kp * (1.0 - e2) * z2 / (kq * (1.0 + kq)) - 0.5 * kp * r2); + re = r - e2 * r0; + re2 = re * re; + ku = sqrt(re2 + z2); + kv = sqrt(re2 + (1.0 - e2) * z2); + kz = b2 / (a * kv); + z0 = kz * zt; + ht = ku * (1.0 - kz); + bt = roi * atan((zt + et2 * z0) / r); + lt = roi * atan2(yt, xt); + fprintf(ofile,"%3.8lf %3.8lf %lf\n", bt , lt , ht); } - i++; } - fprintf(stdout,"...done\nResults written to %s\n", ofilename); + else + { + fprintf(stderr,"Unknow command %s\n",command); + exit(EXIT_FAILURE); + } + if(argc > 4) + { + fprintf(stdout,"\nResults written to %s\n", ofilename); + (void)fclose(ofile); + } (void)fclose(ifile); - (void)fclose(ofile); + (void)fclose(cfile); + + fprintf(stdout,"...done\n"); return(0); } diff --git a/src/helmdiff3d.c b/src/helmdiff3d.c index 4be2110..3732b13 100644 --- a/src/helmdiff3d.c +++ b/src/helmdiff3d.c @@ -69,7 +69,7 @@ int main(int argc, char* argv[]) if(argc < 3) { - fprintf(stdout,"\nSyntax: %s [xyz_src_infilename] [xyz_dest_infilename] [xyz_diff_outfilename]\n\n",argv[0]); + fprintf(stdout,"\nSyntax: %s xyz_src_infilename xyz_dest_infilename [xyz_diff_outfilename]\n\n",argv[0]); fprintf(stdout,"xyz data file format:\n"); fprintf(stdout," X[1] Y[1] Z[1]\n .. .. ..\n .. .. ..\n X[n] Y[n] Z[n]\n\n"); exit(EXIT_FAILURE); @@ -96,8 +96,10 @@ int main(int argc, char* argv[]) if(argc > 3) { ofilename = argv[3]; + ofile = fopen( ofilename, "w"); } - ofile = fopen( ofilename, "w"); + else + ofile = stdout; if(ofile == NULL) { fprintf(stderr,"Error writing %s\n",ofilename); @@ -124,9 +126,14 @@ int main(int argc, char* argv[]) zout=zs-zc; fprintf(ofile,"%lf %lf %lf\n", xout , yout , zout); } - fprintf(stdout,"...done\nResults written to %s\n", ofilename); + if(argc > 3) + { + fprintf(stdout,"\nResults written to %s\n", ofilename); + (void)fclose(ofile); + } (void)fclose(ifile); (void)fclose(cfile); - (void)fclose(ofile); + + fprintf(stdout,"...done\n"); return(0); } diff --git a/src/helmert3d.c b/src/helmert3d.c index ada3139..d98fa55 100644 --- a/src/helmert3d.c +++ b/src/helmert3d.c @@ -73,7 +73,7 @@ int main(int argc, char* argv[]) if(argc < 3) { - fprintf(stdout,"\nSyntax: %s [xyz_src_infilename] [helmert_param_infilename] [xyz_dest_outfilename]\n\n",argv[0]); + fprintf(stdout,"\nSyntax: %s xyz_src_infilename helmert_param_infilename [xyz_dest_outfilename]\n\n",argv[0]); fprintf(stdout,"helmert parameter file format:\n"); fprintf(stdout," r11 r12 r13\n r21 r22 r23\n r31 r32 r33\n tx ty tz\n s\n\n"); fprintf(stdout,"xyz data file format:\n"); @@ -106,8 +106,10 @@ int main(int argc, char* argv[]) if(argc > 3) { ofilename = argv[3]; + ofile = fopen( ofilename, "w"); } - ofile = fopen( ofilename, "w"); + else + ofile = stdout; if(ofile == NULL) { fprintf(stderr,"Error writing %s\n",ofilename); @@ -192,8 +194,6 @@ int main(int argc, char* argv[]) fprintf(stdout,"%lf %lf %lf\n", tx , ty , tz); fprintf(stdout,"%lf\n", m); - fprintf(stdout,"...done\n"); - fprintf(stdout,"Starting transformation...\n"); while(fgets( buf, 128, ifile)!=NULL) { @@ -208,9 +208,14 @@ int main(int argc, char* argv[]) zout=tz+m*(x*r31+y*r32+z*r33); fprintf(ofile,"%lf %lf %lf\n", xout , yout , zout); } - fprintf(stdout,"...done\nResults written to %s\n", ofilename); + if(argc > 3) + { + fprintf(stdout,"\nResults written to %s\n", ofilename); + (void)fclose(ofile); + } (void)fclose(ifile); (void)fclose(parmfile); - (void)fclose(ofile); + + fprintf(stdout,"...done\n"); return(0); } diff --git a/src/helmert3d.h b/src/helmert3d.h index dba2188..18ba40d 100644 --- a/src/helmert3d.h +++ b/src/helmert3d.h @@ -22,7 +22,7 @@ #include #include #include -#define VERS "1.0.2" +#define VERS "1.0.3" #define DEBUG 0 #endif // HELMERT3D_H diff --git a/src/helmparms3d.c b/src/helmparms3d.c index 5f5d19e..ba8b148 100644 --- a/src/helmparms3d.c +++ b/src/helmparms3d.c @@ -301,9 +301,9 @@ int main(int argc, char* argv[]) fprintf(stdout, "* helmparms3d v%s *\n",VERS); fprintf(stdout, "* (c) U. Niethammer 2020 *\n"); fprintf(stdout, "*******************************\n"); - if(argc < 4) + if(argc < 3) { - fprintf(stdout,"\nSyntax: %s [xyz_src_infilename] [xyz_dest_infilename] [helmert_param_outfilename]\n\n",argv[0]); + fprintf(stdout,"\nSyntax: %s xyz_src_infilename xyz_dest_infilename [helmert_param_outfilename]\n\n",argv[0]); fprintf(stderr,"helmert parameter file format:\n"); fprintf(stderr," r11 r12 r13\n r21 r22 r23\n r31 r32 r33\n tx ty tz\n s\n\n"); fprintf(stderr,"xyz data file format:\n"); @@ -312,7 +312,6 @@ int main(int argc, char* argv[]) } src_pts_name = argv[1]; dest_pts_name = argv[2]; - out_param_name = argv[3]; m = get_m_size(src_pts_name); m2 = get_m_size(dest_pts_name); @@ -462,7 +461,7 @@ int main(int argc, char* argv[]) #endif scal=trace1/trace2; - ppm=scal-1.0; + ppm=(scal-1.0)*1000000.0; #if DEBUG fprintf(stdout,"\nscal = %10.10lf\nscal = %10.10lf ppm\n\n",scal, ppm); #endif @@ -510,31 +509,45 @@ int main(int argc, char* argv[]) print_vector(stdout, 3, T_vec); #endif - outfile = fopen(out_param_name, "w"); + if(argc > 3) + { + out_param_name = argv[3]; + outfile = fopen(out_param_name, "w"); + } + else + outfile = stdout; if(outfile == NULL) { fprintf(stderr,"Error writing %s\n",out_param_name); exit(EXIT_FAILURE); } + init_matrix(m,m,src_mat_T); transpose_matrix(m, m, R_mat, src_mat_T); - print_matrix(outfile, n, n, src_mat_T); fprintf(stdout,"R =\n"); print_matrix(stdout, n, n, src_mat_T); fprintf(stdout,"\n"); - print_vector(outfile, 3, T_vec); fprintf(stdout,"T =\n"); print_vector(stdout, 3, T_vec); fprintf(stdout,"\n"); - fprintf(outfile, "%10.10lf\n", scal); fprintf(stdout,"s = %10.10lf (= %10.10lf ppm)\n\n",scal, ppm); - (void)fclose(outfile); - + if(argc < 4) + { + fprintf(stdout,"Results matrix:\n"); + } + print_matrix(outfile, n, n, src_mat_T); + print_vector(outfile, 3, T_vec); + fprintf(outfile, "%10.10lf\n", scal); + if(argc > 3) + { + fprintf(stdout,"Results written to %s\n", out_param_name); + (void)fclose(outfile); + } freevector(D_vec); freevector(T_vec); freevector(one_vec); @@ -552,6 +565,7 @@ int main(int argc, char* argv[]) freematrix(m, C_mat_interm); freematrix(m, src_mat_T); freematrix(m, D_mat_interm); + fprintf(stdout,"...done\n"); return(0); }