From 31597922fd1d61639ab238fc2100c725beb4de96 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Wed, 23 Mar 2022 20:16:33 +0200 Subject: [PATCH 01/10] Add inversible horner real polynomials This enables the user to create horner polynomial transformations in the real number space without the need to supply "inv" parameters: "+inv_origin", "+inv_u" and "+inv_v". We can get the inverse transformation by only using the "fwd" parameters and solving a system of equations iteratively. The system is usually solved in a very small number of iterations. A new optional parameter "+inv_tolerance" is added that can make the iterative process stop earlier if a certain coordinate accuracy is met. For example "+inv_tolerance=0.1" will make the iterative process stop when the calculated coordinates differ less than 0.1 meters, if the units of the transformation destination is meters. --- src/transformations/horner.cpp | 86 ++++++++++++++++++++++++++++++---- test/unit/gie_self_tests.cpp | 74 +++++++++++++++++++++++++++++ 2 files changed, 152 insertions(+), 8 deletions(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 7c8ad19230..8826fb4863 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -99,6 +99,8 @@ struct horner { int order; /* maximum degree of polynomium */ int coefs; /* number of coefficients for each polynomium */ double range; /* radius of the region of validity */ + int has_only_real_fwd; /* in case of real coefficients, has only fwd, inverse is done by gauss-newton iteration */ + double inverse_tolerance; /* in the units of the destination coords, specifies when to stop iterating if has_only_read_fwd and direction is reverse */ double *fwd_u; /* coefficients for the forward transformations */ double *fwd_v; /* i.e. latitude/longitude to northing/easting */ @@ -244,6 +246,7 @@ summing the tiny high order elements first. sz = horner_number_of_coefficients(transformation->order); range = transformation->range; + const bool iterative_inverse = direction == PJ_INV && transformation->has_only_real_fwd; if (direction==PJ_FWD) { /* forward */ tcx = transformation->fwd_u + sz; @@ -251,17 +254,73 @@ summing the tiny high order elements first. e = position.u - transformation->fwd_origin->u; n = position.v - transformation->fwd_origin->v; } else { /* inverse */ - tcx = transformation->inv_u + sz; - tcy = transformation->inv_v + sz; - e = position.u - transformation->inv_origin->u; - n = position.v - transformation->inv_origin->v; + if (!iterative_inverse) { + tcx = transformation->inv_u + sz; + tcy = transformation->inv_v + sz; + e = position.u - transformation->inv_origin->u; + n = position.v - transformation->inv_origin->v; + } else { + // in this case fwd_origin needs to be added in the end + e = position.u; + n = position.v; + } } if ((fabs(n) > range) || (fabs(e) > range)) { proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN); return uv_error; + } else if (iterative_inverse) { + /* + * solve iteratively + * + * | E | | u00 | | u01 + u02*x + ... ' u10 + u11*x + u20*y + ... | | x | + * | | = | | + |-------------------------- ' --------------------------| | | + * | N | | v00 | | v10 + v11*y + v20*x + ... ' v01 + v02*y + ... | | y | + * + * | x | | Ma ' Mb |-1 | E-u00 | + * | | = |-------- | | | + * | y | | Mc ' Md | | N-v00 | + */ + const int order = transformation->order; + const double tol = transformation->inverse_tolerance; + double de = e - transformation->fwd_u[0]; + double dn = n - transformation->fwd_v[0]; + double x0 = 0.0; + double y0 = 0.0; + int loops = 32; // usually converges really fast (1-2 loops) + while (loops-- > 0) { + double Ma = 0.0; + double Mb = 0.0; + double Mc = 0.0; + double Md = 0.0; + tcx = transformation->fwd_u; + tcy = transformation->fwd_v; + for (int i = 0; i <= order; ++i) { + for (int j = 0; j <= (order-i); ++j) { + if (i == 0 && j == 0) { + // do nothing + } else if (i == 0) { + Ma += (*tcx) * pow(x0, j-1); + Md += (*tcy) * pow(y0, j-1); + } else { + Mb += (*tcx) * pow(y0, i-1) * pow(x0, j); + Mc += (*tcy) * pow(x0, i-1) * pow(y0, j); + } + tcx++; + tcy++; + } + } + double idet = 1.0 / (Ma*Md - Mb*Mc); + double x = idet * (Md*de - Mb*dn); + double y = idet * (Ma*dn - Mc*de); + bool ok = (fabs(x-x0) < tol) && (fabs(y-y0) < tol); + x0 = x; + y0 = y; + if (ok) break; + } + position.u = x0 + transformation->fwd_origin->u; + position.v = y0 + transformation->fwd_origin->v; } - /* The melody of this block is straight out of the great Engsager/Poder songbook */ else { int g = transformation->order; @@ -443,6 +502,7 @@ static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { PJ *PROJECTION(horner) { /*********************************************************************/ int degree = 0, n, complex_polynomia = 0; + int has_only_real_fwd = 0; HORNER *Q; P->fwd4d = horner_forward_4d; P->inv4d = horner_reverse_4d; @@ -474,6 +534,14 @@ PJ *PROJECTION(horner) { return horner_freeup (P, PROJ_ERR_OTHER /*ENOMEM*/); P->opaque = Q; + if (!complex_polynomia) { + has_only_real_fwd = + !pj_param_exists(P->params, "inv_u") && + !pj_param_exists(P->params, "inv_v") && + !pj_param_exists(P->params, "inv_origin"); + } + Q->has_only_real_fwd = has_only_real_fwd; + if (complex_polynomia) { /* Westings and/or southings? */ Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0; @@ -506,12 +574,12 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (0==parse_coefs (P, Q->inv_u, "inv_u", n)) + if (!has_only_real_fwd && 0==parse_coefs (P, Q->inv_u, "inv_u", n)) { proj_log_error (P, _("missing inv_u")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (0==parse_coefs (P, Q->inv_v, "inv_v", n)) + if (!has_only_real_fwd && 0==parse_coefs (P, Q->inv_v, "inv_v", n)) { proj_log_error (P, _("missing inv_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); @@ -523,13 +591,15 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) + if (!has_only_real_fwd && 0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) { proj_log_error (P, _("missing inv_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } if (0==parse_coefs (P, &Q->range, "range", 1)) Q->range = 500000; + if (0==parse_coefs (P, &Q->inverse_tolerance, "inv_tolerance", 1)) + Q->inverse_tolerance = 0.001; return P; } diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index 5c14d74709..793f548d3c 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -763,6 +763,80 @@ TEST(gie, horner_selftest) { proj_destroy(P); } +static const char tc32_utm32_fwd_only[] = { + " +proj=horner" + " +ellps=intl" + " +range=10000000" + " +fwd_origin=877605.269066,6125810.306769" + " +deg=4" + " +fwd_v=6.1258112678e+06,9.9999971567e-01,1.5372750011e-10,5.9300860915e-" + "15,2.2609497633e-19,4.3188227445e-05,2.8225130416e-10,7.8740007114e-16,-1." + "7453997279e-19,1.6877465415e-10,-1.1234649773e-14,-1.7042333358e-18,-7." + "9303467953e-15,-5.2906832535e-19,3.9984284847e-19" + " +fwd_u=8.7760574982e+05,9.9999752475e-01,2.8817299305e-10,5.5641310680e-" + "15,-1.5544700949e-18,-4.1357045890e-05,4.2106213519e-11,2.8525551629e-14,-" + "1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1." + "2429019719e-15,5.3886155968e-19,-1.0167505000e-18" }; + +static const char hatt_to_ggrs[] = { + " +proj=horner" + " +ellps=bessel" + " +fwd_origin=0.0, 0.0" + " +deg=2" + " +range=10000000" + " +fwd_u=370552.68, 0.9997155, -1.08e-09, 0.0175123, 2.04e-09, 1.63e-09" + " +fwd_v=4511927.23, 0.9996979, 5.60e-10, -0.0174755, -1.65e-09, -6.50e-10" }; + +TEST(gie, horner_only_fwd_selftest) { + + { + PJ *P = proj_create(PJ_DEFAULT_CTX, tc32_utm32_fwd_only); + ASSERT_TRUE(P != nullptr); + + PJ_COORD a = proj_coord(0, 0, 0, 0); + a.uv.v = 6125305.4245; + a.uv.u = 878354.8539; + + /* Check roundtrip precision for 1 iteration each way, starting in forward + * direction */ + double dist = proj_roundtrip(P, PJ_FWD, 1, &a); + EXPECT_LE(dist, 0.01); + + proj_destroy(P); + } + + { + PJ_COORD a; + a = proj_coord(0, 0, 0, 0); + a.xy.x = -10157.950; + a.xy.y = -21121.093; + PJ_COORD c; + c = proj_coord(0, 0, 0, 0); + c.enu.e = 360028.794; + c.enu.n = 4490989.862; + + PJ *P = proj_create(PJ_DEFAULT_CTX, hatt_to_ggrs); + ASSERT_TRUE(P != nullptr); + + /* Forward projection */ + PJ_COORD b = proj_trans(P, PJ_FWD, a); + double dist = proj_xy_dist(b, c); + EXPECT_LE(dist, 0.001); + + /* Inverse projection */ + b = proj_trans(P, PJ_INV, c); + dist = proj_xy_dist(b, a); + EXPECT_LE(dist, 0.001); + + /* Check roundtrip precision for 1 iteration each way, starting in forward + * direction */ + dist = proj_roundtrip(P, PJ_FWD, 1, &a); + EXPECT_LE(dist, 0.01); + + proj_destroy(P); + } +} + // --------------------------------------------------------------------------- TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) { From f14feeb5dde53f7a108591215a28c7a5a67d7730 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Thu, 24 Mar 2022 19:04:49 +0200 Subject: [PATCH 02/10] Fix windows compilation error This is due to a false positive warning. Also add const to some variables that do not change. --- src/transformations/horner.cpp | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 8826fb4863..8824567d6b 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -222,7 +222,6 @@ summing the tiny high order elements first. /* These variable names follow the Engsager/Poder implementation */ int sz; /* Number of coefficients per polynomial */ - double *tcx, *tcy; /* Coefficient pointers */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ double n, e; PJ_UV uv_error; @@ -249,14 +248,10 @@ summing the tiny high order elements first. const bool iterative_inverse = direction == PJ_INV && transformation->has_only_real_fwd; if (direction==PJ_FWD) { /* forward */ - tcx = transformation->fwd_u + sz; - tcy = transformation->fwd_v + sz; e = position.u - transformation->fwd_origin->u; n = position.v - transformation->fwd_origin->v; } else { /* inverse */ if (!iterative_inverse) { - tcx = transformation->inv_u + sz; - tcy = transformation->inv_v + sz; e = position.u - transformation->inv_origin->u; n = position.v - transformation->inv_origin->v; } else { @@ -283,8 +278,8 @@ summing the tiny high order elements first. */ const int order = transformation->order; const double tol = transformation->inverse_tolerance; - double de = e - transformation->fwd_u[0]; - double dn = n - transformation->fwd_v[0]; + const double de = e - transformation->fwd_u[0]; + const double dn = n - transformation->fwd_v[0]; double x0 = 0.0; double y0 = 0.0; int loops = 32; // usually converges really fast (1-2 loops) @@ -293,8 +288,9 @@ summing the tiny high order elements first. double Mb = 0.0; double Mc = 0.0; double Md = 0.0; - tcx = transformation->fwd_u; - tcy = transformation->fwd_v; + /* Coefficient pointers */ + double *tcx = transformation->fwd_u; + double *tcy = transformation->fwd_v; for (int i = 0; i <= order; ++i) { for (int j = 0; j <= (order-i); ++j) { if (i == 0 && j == 0) { @@ -326,6 +322,9 @@ summing the tiny high order elements first. int g = transformation->order; int r = g, c; double u, v, N, E; + /* Coefficient pointers */ + double *tcx = direction == PJ_FWD ? (transformation->fwd_u + sz) : (transformation->inv_u + sz); + double *tcy = direction == PJ_FWD ? (transformation->fwd_v + sz) : (transformation->inv_v + sz); /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */ N = *--tcy; From 12c3377b04a4d73cdf1f4caedf28addfd8bd741f Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Thu, 24 Mar 2022 19:18:02 +0200 Subject: [PATCH 03/10] Handle some inversible horner corner cases In the case of inversible horner, return error values if the iterations required for convergence are too high. --- src/transformations/horner.cpp | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 8824567d6b..892e5e0e88 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -283,7 +283,8 @@ summing the tiny high order elements first. double x0 = 0.0; double y0 = 0.0; int loops = 32; // usually converges really fast (1-2 loops) - while (loops-- > 0) { + bool converged = false; + while (loops-- > 0 && !converged) { double Ma = 0.0; double Mb = 0.0; double Mc = 0.0; @@ -309,13 +310,19 @@ summing the tiny high order elements first. double idet = 1.0 / (Ma*Md - Mb*Mc); double x = idet * (Md*de - Mb*dn); double y = idet * (Ma*dn - Mc*de); - bool ok = (fabs(x-x0) < tol) && (fabs(y-y0) < tol); + converged = (fabs(x-x0) < tol) && (fabs(y-y0) < tol); x0 = x; y0 = y; - if (ok) break; } - position.u = x0 + transformation->fwd_origin->u; - position.v = y0 + transformation->fwd_origin->v; + // if loops have been exhausted and we have not converged yet, + // we are never going to converge + if (!converged) { + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM); + return uv_error; + } else { + position.u = x0 + transformation->fwd_origin->u; + position.v = y0 + transformation->fwd_origin->v; + } } /* The melody of this block is straight out of the great Engsager/Poder songbook */ else { From 04baa726c004aff061e2ed3355078635db299234 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Fri, 25 Mar 2022 17:29:47 +0200 Subject: [PATCH 04/10] Add documentation for iterative real polynomial inversion --- .../operations/transformations/horner.rst | 69 +++++++++++++++++-- 1 file changed, 65 insertions(+), 4 deletions(-) diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst index ad09e8abd0..6636a2c82c 100644 --- a/docs/source/operations/transformations/horner.rst +++ b/docs/source/operations/transformations/horner.rst @@ -5,6 +5,7 @@ Horner polynomial evaluation ================================================================================ .. versionadded:: 5.0.0 +.. versionadded:: 9.1.0 Iterative real polynormal inversion +-----------------+-------------------------------------------------------------------+ | **Alias** | horner | @@ -61,7 +62,59 @@ The final coordinates are determined as Y_{out} = Y_{in} + \Delta Y The inverse transform is the same as the above but requires a different set of -coefficients. +coefficients. If only the forward set of coefficients and origin is known the inverse transform can +be done by iteratively solving a system of equations. By writing :eq:`real_poly` as: + +.. math:: + :label: real_poly_iterative_inversion + + \begin{align} + \begin{bmatrix} + \Delta X \\ + \Delta Y \\ + \end{bmatrix} = + \begin{bmatrix} + u_{0,0} \\ + v_{0,0} \\ + \end{bmatrix} + + \begin{bmatrix} + u_{0,1} + u_{0,2} U + ... & u_{1,0} + u_{1,1} U + u_{2,0} V + ... \\ + v_{1,0} + v_{1,1} V + v_{2,0} U + ... & v_{0,1} + v_{0,2} V \\ + \end{bmatrix} + \begin{bmatrix} + U \\ + V \\ + \end{bmatrix} = + \begin{bmatrix} + u_{0,0} \\ + v_{0,0} \\ + \end{bmatrix} + + \begin{bmatrix} + MA & MB \\ + MC & MD \\ + \end{bmatrix} + \begin{bmatrix} + U \\ + V \\ + \end{bmatrix} + \end{align} + + \begin{align} + \begin{bmatrix} + U \\ + V \\ + \end{bmatrix} = + \begin{bmatrix} + MA & MB \\ + MC & MD \\ + \end{bmatrix}^-1 + \begin{bmatrix} + \Delta X - u_{0,0} \\ + \Delta Y - v_{0,0} \\ + \end{bmatrix} + \end{align} + +We can iteratively solve with initial values of :math:`U = 0` and :math:`V = 0` and find :math:`U` and :math:`V`. Evaluation of the complex polynomials are defined by the following equations: @@ -133,7 +186,7 @@ describing real and complex polynomials can't be mixed. .. option:: +inv_origin= - Coordinate of origin for the inverse mapping + Coordinate of origin for the inverse mapping. Not required for iterative inversion. Real polynomials .............................................................................. @@ -160,12 +213,12 @@ of the polynomial: .. option:: +inv_u= Coefficients for the inverse transformation i.e. latitude to northing - as described in :eq:`real_poly`. + as described in :eq:`real_poly`. Not required for iterative inversion. .. option:: +inv_v= Coefficients for the inverse transformation i.e. longitude to easting - as described in :eq:`real_poly`. + as described in :eq:`real_poly`. Not required for iterative inversion. @@ -206,6 +259,14 @@ Optional Express longitude as westing. Only applies for complex polynomials. +.. option:: +inv_tolerance= + + Only applies to real polynomials and iterative inversion. + The procedure converges to the correct results with each step. + Iteration stops when the result differs from the previous calculated + result by less than . + should be the same units as :math:`U` and :math:`V` of :eq:`UV` + Defaults to 0.001 meters. Further reading ################################################################################ From 3dbe7017481a06e1c88eabe878e7d821f9bdb3da Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Sat, 26 Mar 2022 15:31:16 +0200 Subject: [PATCH 05/10] Address docs review comments Move +inv_u, +inv_v to optional parameters in horner docs Make some consistency changes for matrices and default values --- .../operations/transformations/horner.rst | 54 +++++++++++-------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst index 6636a2c82c..8836dce18b 100644 --- a/docs/source/operations/transformations/horner.rst +++ b/docs/source/operations/transformations/horner.rst @@ -5,7 +5,7 @@ Horner polynomial evaluation ================================================================================ .. versionadded:: 5.0.0 -.. versionadded:: 9.1.0 Iterative real polynormal inversion +.. versionchanged:: 9.1.0 Iterative real polynormal inversion +-----------------+-------------------------------------------------------------------+ | **Alias** | horner | @@ -66,9 +66,6 @@ coefficients. If only the forward set of coefficients and origin is known the in be done by iteratively solving a system of equations. By writing :eq:`real_poly` as: .. math:: - :label: real_poly_iterative_inversion - - \begin{align} \begin{bmatrix} \Delta X \\ \Delta Y \\ @@ -84,6 +81,12 @@ be done by iteratively solving a system of equations. By writing :eq:`real_poly` \begin{bmatrix} U \\ V \\ + \end{bmatrix} \\ + +.. math:: + \begin{bmatrix} + \Delta X \\ + \Delta Y \\ \end{bmatrix} = \begin{bmatrix} u_{0,0} \\ @@ -96,10 +99,9 @@ be done by iteratively solving a system of equations. By writing :eq:`real_poly` \begin{bmatrix} U \\ V \\ - \end{bmatrix} - \end{align} + \end{bmatrix} \\ - \begin{align} +.. math:: \begin{bmatrix} U \\ V \\ @@ -107,12 +109,11 @@ be done by iteratively solving a system of equations. By writing :eq:`real_poly` \begin{bmatrix} MA & MB \\ MC & MD \\ - \end{bmatrix}^-1 + \end{bmatrix}^{-1} \begin{bmatrix} \Delta X - u_{0,0} \\ \Delta Y - v_{0,0} \\ \end{bmatrix} - \end{align} We can iteratively solve with initial values of :math:`U = 0` and :math:`V = 0` and find :math:`U` and :math:`V`. @@ -210,18 +211,6 @@ of the polynomial: Coefficients for the forward transformation i.e. longitude to easting as described in :eq:`real_poly`. -.. option:: +inv_u= - - Coefficients for the inverse transformation i.e. latitude to northing - as described in :eq:`real_poly`. Not required for iterative inversion. - -.. option:: +inv_v= - - Coefficients for the inverse transformation i.e. longitude to easting - as described in :eq:`real_poly`. Not required for iterative inversion. - - - Complex polynomials .............................................................................. @@ -247,6 +236,24 @@ of the polynomial: Optional ------------------------------------------------------------------------------- +.. option:: +inv_u= + + .. versionchanged:: 9.1.0 + + Coefficients for the inverse transformation i.e. latitude to northing + as described in :eq:`real_poly`. Only applies for real polynomials. + Without this option iterative real polynomial evaluation is used for + the inverse tranformation. + +.. option:: +inv_v= + + .. versionchanged:: 9.1.0 + + Coefficients for the inverse transformation i.e. longitude to easting + as described in :eq:`real_poly`. Only applies for real polynomials. + Without this option iterative real polynomial evaluation is used for + the inverse tranformation. + .. option:: +range= Radius of the region of validity. @@ -261,12 +268,15 @@ Optional .. option:: +inv_tolerance= + .. versionadded:: 9.1.0 + Only applies to real polynomials and iterative inversion. The procedure converges to the correct results with each step. Iteration stops when the result differs from the previous calculated result by less than . should be the same units as :math:`U` and :math:`V` of :eq:`UV` - Defaults to 0.001 meters. + + *Defaults to 0.001.* Further reading ################################################################################ From 98e4a89947ec616b25dff5c4ee9098907f145a99 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Sat, 26 Mar 2022 18:55:45 +0200 Subject: [PATCH 06/10] Add inversible horner complex polynomials This enables the user to create horner polynomial transformations in the complex number space without the need to supply "inv" parameters: "+inv_origin" and "+inv_c". Also added a unit test and respective docs --- .../operations/transformations/horner.rst | 36 ++++--- src/transformations/horner.cpp | 96 ++++++++++++++----- test/unit/gie_self_tests.cpp | 37 +++++++ 3 files changed, 133 insertions(+), 36 deletions(-) diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst index 8836dce18b..9eb8d39c8e 100644 --- a/docs/source/operations/transformations/horner.rst +++ b/docs/source/operations/transformations/horner.rst @@ -5,7 +5,7 @@ Horner polynomial evaluation ================================================================================ .. versionadded:: 5.0.0 -.. versionchanged:: 9.1.0 Iterative real polynormal inversion +.. versionchanged:: 9.1.0 Iterative polynormal inversion +-----------------+-------------------------------------------------------------------+ | **Alias** | horner | @@ -126,7 +126,7 @@ Evaluation of the complex polynomials are defined by the following equations: Where :math:`n` is the degree of the polynomial. :math:`U` and :math:`V` are defined as in :eq:`UV` and the resulting coordinates are again determined -by :eq:`xy_out`. +by :eq:`xy_out`. Complex polynomials can be solved iteratively similar to real polynomials. Examples ################################################################################ @@ -185,10 +185,6 @@ describing real and complex polynomials can't be mixed. Coordinate of origin for the forward mapping -.. option:: +inv_origin= - - Coordinate of origin for the inverse mapping. Not required for iterative inversion. - Real polynomials .............................................................................. @@ -228,21 +224,24 @@ of the polynomial: Coefficients for the complex forward transformation as described in :eq:`complex_poly`. -.. option:: +inv_c= - - Coefficients for the complex inverse transformation - as described in :eq:`complex_poly`. - Optional ------------------------------------------------------------------------------- +.. option:: +inv_origin= + + .. versionchanged:: 9.1.0 + + Coordinate of origin for the inverse mapping. + Without this option iterative polynomial evaluation is used for + the inverse tranformation. + .. option:: +inv_u= .. versionchanged:: 9.1.0 Coefficients for the inverse transformation i.e. latitude to northing as described in :eq:`real_poly`. Only applies for real polynomials. - Without this option iterative real polynomial evaluation is used for + Without this option iterative polynomial evaluation is used for the inverse tranformation. .. option:: +inv_v= @@ -251,7 +250,16 @@ Optional Coefficients for the inverse transformation i.e. longitude to easting as described in :eq:`real_poly`. Only applies for real polynomials. - Without this option iterative real polynomial evaluation is used for + Without this option iterative polynomial evaluation is used for + the inverse tranformation. + +.. option:: +inv_c= + + .. versionchanged:: 9.1.0 + + Coefficients for the complex inverse transformation + as described in :eq:`complex_poly`. Only applies for complex polynomials. + Without this option iterative polynomial evaluation is used for the inverse tranformation. .. option:: +range= @@ -270,7 +278,7 @@ Optional .. versionadded:: 9.1.0 - Only applies to real polynomials and iterative inversion. + Only applies to cases of iterative inversion. The procedure converges to the correct results with each step. Iteration stops when the result differs from the previous calculated result by less than . diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 892e5e0e88..0c5219006a 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -82,6 +82,7 @@ #include #include #include +#include #include "proj.h" #include "proj_internal.h" @@ -99,8 +100,8 @@ struct horner { int order; /* maximum degree of polynomium */ int coefs; /* number of coefficients for each polynomium */ double range; /* radius of the region of validity */ - int has_only_real_fwd; /* in case of real coefficients, has only fwd, inverse is done by gauss-newton iteration */ - double inverse_tolerance; /* in the units of the destination coords, specifies when to stop iterating if has_only_read_fwd and direction is reverse */ + int has_only_fwd; /* inv parameters are not specified, inverse is done by gauss-newton iteration */ + double inverse_tolerance; /* in the units of the destination coords, specifies when to stop iterating if has_only_fwd and direction is reverse */ double *fwd_u; /* coefficients for the forward transformations */ double *fwd_v; /* i.e. latitude/longitude to northing/easting */ @@ -245,7 +246,7 @@ summing the tiny high order elements first. sz = horner_number_of_coefficients(transformation->order); range = transformation->range; - const bool iterative_inverse = direction == PJ_INV && transformation->has_only_real_fwd; + const bool iterative_inverse = direction == PJ_INV && transformation->has_only_fwd; if (direction==PJ_FWD) { /* forward */ e = position.u - transformation->fwd_origin->u; @@ -384,7 +385,6 @@ polynomial evaluation engine. /* These variable names follow the Engsager/Poder implementation */ int sz; /* Number of coefficients */ - double *c, *cb; /* Coefficient pointers */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ double n, e, w, N, E; PJ_UV uv_error; @@ -408,9 +408,9 @@ polynomial evaluation engine. sz = 2*transformation->order + 2; range = transformation->range; + const bool iterative_inverse = direction == PJ_INV && transformation->has_only_fwd; + if (direction==PJ_FWD) { /* forward */ - cb = transformation->fwd_c; - c = cb + sz; e = position.u - transformation->fwd_origin->u; n = position.v - transformation->fwd_origin->v; if (transformation->uneg) @@ -418,14 +418,18 @@ polynomial evaluation engine. if (transformation->vneg) n = -n; } else { /* inverse */ - cb = transformation->inv_c; - c = cb + sz; - e = position.u - transformation->inv_origin->u; - n = position.v - transformation->inv_origin->v; - if (transformation->uneg) - e = -e; - if (transformation->vneg) - n = -n; + if (!iterative_inverse) { + e = position.u - transformation->inv_origin->u; + n = position.v - transformation->inv_origin->v; + if (transformation->uneg) + e = -e; + if (transformation->vneg) + n = -n; + } else { + // in this case fwd_origin and any existing flipping needs to be added in the end + e = position.u; + n = position.v; + } } if ((fabs(n) > range) || (fabs(e) > range)) { @@ -433,6 +437,51 @@ polynomial evaluation engine. return uv_error; } + if (iterative_inverse) { + const int order = transformation->order; + const double tol = transformation->inverse_tolerance; + const std::complex dZ(n-transformation->fwd_c[0], e-transformation->fwd_c[1]); + std::complex w0(0.0, 0.0); + int loops = 32; // usually converges really fast (1-2 loops) + bool converged = false; + while (loops-- > 0 && !converged) { + double *cb = transformation->fwd_c + 2; // coefficient pointers after c0 + std::complex det(0.0, 0.0); + for (int i = 1; i <= order; ++i) { + double cbn = *cb++; + double cbe = *cb++; + if (i == 1) { + det += std::complex(cbn, cbe); + } else { + det += std::complex(cbn, cbe) * std::pow(w0, double(i-1)); + } + } + std::complex w1 = dZ / det; + converged = (fabs(w1.real()-w0.real()) < tol) && (fabs(w1.imag()-w0.imag()) < tol); + w0 = w1; + } + // if loops have been exhausted and we have not converged yet, + // we are never going to converge + if (!converged) { + proj_errno_set(P, PROJ_ERR_COORD_TRANSFM); + position = uv_error; + } else { + E = w0.imag(); + N = w0.real(); + if (transformation->uneg) + E = -E; + if (transformation->vneg) + N = -N; + position.u = E + transformation->fwd_origin->u; + position.v = N + transformation->fwd_origin->v; + } + return position; + } + + // coefficient pointers + double *cb = direction == PJ_FWD ? transformation->fwd_c : transformation->inv_c; + double *c = cb + sz; + /* Everything's set up properly - now do the actual polynomium evaluation */ E = *--c; N = *--c; @@ -441,7 +490,6 @@ polynomial evaluation engine. N = n*N - e*E + *--c; E = w; } - position.u = E; position.v = N; return position; @@ -508,7 +556,7 @@ static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { PJ *PROJECTION(horner) { /*********************************************************************/ int degree = 0, n, complex_polynomia = 0; - int has_only_real_fwd = 0; + int has_only_fwd = 0; HORNER *Q; P->fwd4d = horner_forward_4d; P->inv4d = horner_reverse_4d; @@ -541,12 +589,16 @@ PJ *PROJECTION(horner) { P->opaque = Q; if (!complex_polynomia) { - has_only_real_fwd = + has_only_fwd = !pj_param_exists(P->params, "inv_u") && !pj_param_exists(P->params, "inv_v") && !pj_param_exists(P->params, "inv_origin"); + } else { + has_only_fwd = + !pj_param_exists(P->params, "inv_c") && + !pj_param_exists(P->params, "inv_origin"); } - Q->has_only_real_fwd = has_only_real_fwd; + Q->has_only_fwd = has_only_fwd; if (complex_polynomia) { /* Westings and/or southings? */ @@ -559,7 +611,7 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_c")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (0==parse_coefs (P, Q->inv_c, "inv_c", n)) + if (!has_only_fwd && 0==parse_coefs (P, Q->inv_c, "inv_c", n)) { proj_log_error (P, _("missing inv_c")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); @@ -580,12 +632,12 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_real_fwd && 0==parse_coefs (P, Q->inv_u, "inv_u", n)) + if (!has_only_fwd && 0==parse_coefs (P, Q->inv_u, "inv_u", n)) { proj_log_error (P, _("missing inv_u")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_real_fwd && 0==parse_coefs (P, Q->inv_v, "inv_v", n)) + if (!has_only_fwd && 0==parse_coefs (P, Q->inv_v, "inv_v", n)) { proj_log_error (P, _("missing inv_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); @@ -597,7 +649,7 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_real_fwd && 0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) + if (!has_only_fwd && 0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) { proj_log_error (P, _("missing inv_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index 793f548d3c..e91464d5e6 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -778,6 +778,16 @@ static const char tc32_utm32_fwd_only[] = { "1.9107771273e-18,3.3615590093e-10,2.4380247154e-14,-2.0241230315e-18,1." "2429019719e-15,5.3886155968e-19,-1.0167505000e-18" }; +static const char sb_utm32_fwd_only[] = { + " +proj=horner" + " +ellps=intl" + " +range=10000000" + " +fwd_origin=4.94690026817276e+05,6.13342113183056e+06" + " +deg=3" + " +fwd_c=6.13258562111350e+06,6.19480105709997e+05,9.99378966275206e-01,-2." + "82153291753490e-02,-2.27089979140026e-10,-1.77019590701470e-09,1." + "08522286274070e-14,2.11430298751604e-15" }; + static const char hatt_to_ggrs[] = { " +proj=horner" " +ellps=bessel" @@ -835,6 +845,33 @@ TEST(gie, horner_only_fwd_selftest) { proj_destroy(P); } + + { + PJ *P = proj_create(PJ_DEFAULT_CTX, sb_utm32_fwd_only); + ASSERT_TRUE(P != nullptr); + + PJ_COORD a = proj_coord(0, 0, 0, 0); + PJ_COORD b = proj_coord(0, 0, 0, 0); + PJ_COORD c = proj_coord(0, 0, 0, 0); + a.uv.v = 6130821.2945; + a.uv.u = 495136.8544; + c.uv.v = 6130000.0000; + c.uv.u = 620000.0000; + + /* Forward projection */ + b = proj_trans(P, PJ_FWD, a); + double dist = proj_xy_dist(b, c); + EXPECT_LE(dist, 0.001); + + /* Inverse projection */ + b = proj_trans(P, PJ_INV, c); + dist = proj_xy_dist(b, a); + EXPECT_LE(dist, 0.001); + + /* Check roundtrip precision for 1 iteration each way */ + dist = proj_roundtrip(P, PJ_FWD, 1, &a); + EXPECT_LE(dist, 0.01); + } } // --------------------------------------------------------------------------- From 16aab79011b45749f196ecedb8b88c6bc0dfd237 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Sat, 26 Mar 2022 20:29:58 +0200 Subject: [PATCH 07/10] Optimize iterative inversion of complex horner Use a similar algorithm as the normal path. This enables us to not use std::pow at all. --- src/transformations/horner.cpp | 32 ++++++++++++++++---------------- 1 file changed, 16 insertions(+), 16 deletions(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 0c5219006a..fd9ba1bfca 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -386,7 +386,7 @@ polynomial evaluation engine. /* These variable names follow the Engsager/Poder implementation */ int sz; /* Number of coefficients */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ - double n, e, w, N, E; + double n, e; PJ_UV uv_error; uv_error.u = uv_error.v = HUGE_VAL; @@ -438,24 +438,23 @@ polynomial evaluation engine. } if (iterative_inverse) { - const int order = transformation->order; const double tol = transformation->inverse_tolerance; const std::complex dZ(n-transformation->fwd_c[0], e-transformation->fwd_c[1]); std::complex w0(0.0, 0.0); int loops = 32; // usually converges really fast (1-2 loops) bool converged = false; while (loops-- > 0 && !converged) { - double *cb = transformation->fwd_c + 2; // coefficient pointers after c0 - std::complex det(0.0, 0.0); - for (int i = 1; i <= order; ++i) { - double cbn = *cb++; - double cbe = *cb++; - if (i == 1) { - det += std::complex(cbn, cbe); - } else { - det += std::complex(cbn, cbe) * std::pow(w0, double(i-1)); - } + // coefficient pointers from back to front until the first complex pair (fwd_c0+i*fwd_c1) + double *cb = transformation->fwd_c + 2; + double *c = cb + sz; + double E = *--c; + double N = *--c; + while (c > cb) { + double w = w0.real()*E + w0.imag()*N + *--c; + N = w0.real()*N - w0.imag()*E + *--c; + E = w; } + std::complex det(N, E); std::complex w1 = dZ / det; converged = (fabs(w1.real()-w0.real()) < tol) && (fabs(w1.imag()-w0.imag()) < tol); w0 = w1; @@ -466,8 +465,8 @@ polynomial evaluation engine. proj_errno_set(P, PROJ_ERR_COORD_TRANSFM); position = uv_error; } else { - E = w0.imag(); - N = w0.real(); + double E = w0.imag(); + double N = w0.real(); if (transformation->uneg) E = -E; if (transformation->vneg) @@ -483,8 +482,9 @@ polynomial evaluation engine. double *c = cb + sz; /* Everything's set up properly - now do the actual polynomium evaluation */ - E = *--c; - N = *--c; + double E = *--c; + double N = *--c; + double w; while (c > cb) { w = n*E + e*N + *--c; N = n*N - e*E + *--c; From 1b556a869075d153fab13f84dd39b5226d1d51c4 Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Sat, 26 Mar 2022 23:07:52 +0200 Subject: [PATCH 08/10] Fix sanitizer errors --- src/transformations/horner.cpp | 3 ++- test/unit/gie_self_tests.cpp | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index fd9ba1bfca..4d7fda9f48 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -445,8 +445,9 @@ polynomial evaluation engine. bool converged = false; while (loops-- > 0 && !converged) { // coefficient pointers from back to front until the first complex pair (fwd_c0+i*fwd_c1) - double *cb = transformation->fwd_c + 2; + double *cb = transformation->fwd_c; double *c = cb + sz; + cb += 2; double E = *--c; double N = *--c; while (c > cb) { diff --git a/test/unit/gie_self_tests.cpp b/test/unit/gie_self_tests.cpp index e91464d5e6..c8933ea21a 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -871,6 +871,8 @@ TEST(gie, horner_only_fwd_selftest) { /* Check roundtrip precision for 1 iteration each way */ dist = proj_roundtrip(P, PJ_FWD, 1, &a); EXPECT_LE(dist, 0.01); + + proj_destroy(P); } } From d75519bda97cd0b2cdccce325e9c21f70ae566ff Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Thu, 31 Mar 2022 22:38:15 +0300 Subject: [PATCH 09/10] Apply review comments - repurpose has_only_fwd variable to has_inv - optimize polynomial inversion - refactor horner evaluation code to reusable functions --- src/transformations/horner.cpp | 190 +++++++++++++++++---------------- 1 file changed, 99 insertions(+), 91 deletions(-) diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 4d7fda9f48..19b657c8f5 100644 --- a/src/transformations/horner.cpp +++ b/src/transformations/horner.cpp @@ -100,8 +100,9 @@ struct horner { int order; /* maximum degree of polynomium */ int coefs; /* number of coefficients for each polynomium */ double range; /* radius of the region of validity */ - int has_only_fwd; /* inv parameters are not specified, inverse is done by gauss-newton iteration */ - double inverse_tolerance; /* in the units of the destination coords, specifies when to stop iterating if has_only_fwd and direction is reverse */ + bool has_inv; /* inv parameters are specified */ + double inverse_tolerance; /* in the units of the destination coords, + specifies when to stop iterating if !has_inv and direction is reverse */ double *fwd_u; /* coefficients for the forward transformations */ double *fwd_v; /* i.e. latitude/longitude to northing/easting */ @@ -180,8 +181,65 @@ static HORNER *horner_alloc (size_t order, int complex_polynomia) { return nullptr; } +inline static PJ_UV double_real_horner_eval(int order, const double *cx, const double *cy, PJ_UV en, int order_offset = 0) +{ + /* + The melody of this block is straight out of the great Engsager/Poder songbook. + For numerical stability, the summation is carried out backwards, + summing the tiny high order elements first. + Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout + */ + const double n = en.v; + const double e = en.u; + const int sz = horner_number_of_coefficients(order); /* Number of coefficients per polynomial */ + cx += sz; + cy += sz; + double N = *--cy; + double E = *--cx; + for (int r = order; r > order_offset; r--) { + double u = *--cy; + double v = *--cx; + for (int c = order; c >= r; c--) { + u = n*u + *--cy; + v = e*v + *--cx; + } + N = e*N + u; + E = n*E + v; + } + return { E, N }; +} +inline static double single_real_horner_eval(int order, const double *cx, double x, int order_offset = 0) +{ + const int sz = order + 1; /* Number of coefficients per polynomial */ + cx += sz; + double u = *--cx; + for (int r = order; r > order_offset; r--) { + u = x*u + *--cx; + } + return u; +} + +inline static PJ_UV complex_horner_eval(int order, const double *c, PJ_UV en, int order_offset = 0) +{ + // the coefficients are ordered like this: + // (Cn0+i*Ce0, Cn1+i*Ce1, ...) + const int sz = 2*order + 2; // number of coefficients + const double e = en.u; + const double n = en.v; + const double *cbeg = c + order_offset*2; + c += sz; + double E = *--c; + double N = *--c; + double w; + while (c > cbeg) { + w = n*E + e*N + *--c; + N = n*N - e*E + *--c; + E = w; + } + return { E, N }; +} /**********************************************************************/ static PJ_UV horner_func (PJ* P, const HORNER *transformation, PJ_DIRECTION direction, PJ_UV position) { @@ -216,13 +274,9 @@ P = sum (i = [0 : order]) sum (j = [0 : order - i]) pow(par_1, i) * pow(par_2, j) * coef(index(order, i, j)) -For numerical stability, the summation is carried out backwards, -summing the tiny high order elements first. - ***********************************************************************/ /* These variable names follow the Engsager/Poder implementation */ - int sz; /* Number of coefficients per polynomial */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ double n, e; PJ_UV uv_error; @@ -243,10 +297,9 @@ summing the tiny high order elements first. } /* Prepare for double Horner */ - sz = horner_number_of_coefficients(transformation->order); range = transformation->range; - const bool iterative_inverse = direction == PJ_INV && transformation->has_only_fwd; + const bool iterative_inverse = direction == PJ_INV && !transformation->has_inv; if (direction==PJ_FWD) { /* forward */ e = position.u - transformation->fwd_origin->u; @@ -290,23 +343,17 @@ summing the tiny high order elements first. double Mb = 0.0; double Mc = 0.0; double Md = 0.0; - /* Coefficient pointers */ - double *tcx = transformation->fwd_u; - double *tcy = transformation->fwd_v; - for (int i = 0; i <= order; ++i) { - for (int j = 0; j <= (order-i); ++j) { - if (i == 0 && j == 0) { - // do nothing - } else if (i == 0) { - Ma += (*tcx) * pow(x0, j-1); - Md += (*tcy) * pow(y0, j-1); - } else { - Mb += (*tcx) * pow(y0, i-1) * pow(x0, j); - Mc += (*tcy) * pow(x0, i-1) * pow(y0, j); - } - tcx++; - tcy++; - } + { + const double *tcx = transformation->fwd_u; + const double *tcy = transformation->fwd_v; + PJ_UV x0y0 = { x0, y0 }; + // sum the i > 0 coefficients + PJ_UV Mbc = double_real_horner_eval(order, tcx, tcy, x0y0, 1); + Mb = Mbc.u; + Mc = Mbc.v; + // sum the i = 0, j > 0 coefficients + Ma = single_real_horner_eval(order, tcx, x0, 1); + Md = single_real_horner_eval(order, tcy, y0, 1); } double idet = 1.0 / (Ma*Md - Mb*Mc); double x = idet * (Md*de - Mb*dn); @@ -325,31 +372,11 @@ summing the tiny high order elements first. position.v = y0 + transformation->fwd_origin->v; } } - /* The melody of this block is straight out of the great Engsager/Poder songbook */ else { - int g = transformation->order; - int r = g, c; - double u, v, N, E; - /* Coefficient pointers */ - double *tcx = direction == PJ_FWD ? (transformation->fwd_u + sz) : (transformation->inv_u + sz); - double *tcy = direction == PJ_FWD ? (transformation->fwd_v + sz) : (transformation->inv_v + sz); - - /* Double Horner's scheme: N = n*Cy*e -> yout, E = e*Cx*n -> xout */ - N = *--tcy; - E = *--tcx; - for (; r > 0; r--) { - u = *--tcy; - v = *--tcx; - for (c = g; c >= r; c--) { - u = n*u + *--tcy; - v = e*v + *--tcx; - } - N = e*N + u; - E = n*E + v; - } - - position.u = E; - position.v = N; + const double *tcx = direction == PJ_FWD ? transformation->fwd_u : transformation->inv_u; + const double *tcy = direction == PJ_FWD ? transformation->fwd_v : transformation->inv_v; + PJ_UV en = { e, n }; + position = double_real_horner_eval(transformation->order, tcx, tcy, en); } return position; @@ -384,7 +411,6 @@ polynomial evaluation engine. ***********************************************************************/ /* These variable names follow the Engsager/Poder implementation */ - int sz; /* Number of coefficients */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ double n, e; PJ_UV uv_error; @@ -405,10 +431,9 @@ polynomial evaluation engine. } /* Prepare for double Horner */ - sz = 2*transformation->order + 2; range = transformation->range; - const bool iterative_inverse = direction == PJ_INV && transformation->has_only_fwd; + const bool iterative_inverse = direction == PJ_INV && !transformation->has_inv; if (direction==PJ_FWD) { /* forward */ e = position.u - transformation->fwd_origin->u; @@ -438,24 +463,18 @@ polynomial evaluation engine. } if (iterative_inverse) { + // complex real part corresponds to Northing, imag part to Easting const double tol = transformation->inverse_tolerance; const std::complex dZ(n-transformation->fwd_c[0], e-transformation->fwd_c[1]); - std::complex w0(0.0, 0.0); + std::complex w0(0.0, 0.0); int loops = 32; // usually converges really fast (1-2 loops) bool converged = false; while (loops-- > 0 && !converged) { - // coefficient pointers from back to front until the first complex pair (fwd_c0+i*fwd_c1) - double *cb = transformation->fwd_c; - double *c = cb + sz; - cb += 2; - double E = *--c; - double N = *--c; - while (c > cb) { - double w = w0.real()*E + w0.imag()*N + *--c; - N = w0.real()*N - w0.imag()*E + *--c; - E = w; - } - std::complex det(N, E); + // sum coefficient pointers from back to front until the first complex pair (fwd_c0+i*fwd_c1) + const double *c = transformation->fwd_c; + PJ_UV en = { w0.imag(), w0.real() }; + en = complex_horner_eval(transformation->order, c, en, 1); + std::complex det(en.v, en.u); std::complex w1 = dZ / det; converged = (fabs(w1.real()-w0.real()) < tol) && (fabs(w1.imag()-w0.imag()) < tol); w0 = w1; @@ -480,19 +499,8 @@ polynomial evaluation engine. // coefficient pointers double *cb = direction == PJ_FWD ? transformation->fwd_c : transformation->inv_c; - double *c = cb + sz; - - /* Everything's set up properly - now do the actual polynomium evaluation */ - double E = *--c; - double N = *--c; - double w; - while (c > cb) { - w = n*E + e*N + *--c; - N = n*N - e*E + *--c; - E = w; - } - position.u = E; - position.v = N; + PJ_UV en = { e, n }; + position = complex_horner_eval(transformation->order, cb, en); return position; } @@ -557,7 +565,7 @@ static int parse_coefs (PJ *P, double *coefs, const char *param, int ncoefs) { PJ *PROJECTION(horner) { /*********************************************************************/ int degree = 0, n, complex_polynomia = 0; - int has_only_fwd = 0; + bool has_inv = false; HORNER *Q; P->fwd4d = horner_forward_4d; P->inv4d = horner_reverse_4d; @@ -590,16 +598,16 @@ PJ *PROJECTION(horner) { P->opaque = Q; if (!complex_polynomia) { - has_only_fwd = - !pj_param_exists(P->params, "inv_u") && - !pj_param_exists(P->params, "inv_v") && - !pj_param_exists(P->params, "inv_origin"); + has_inv = + pj_param_exists(P->params, "inv_u") || + pj_param_exists(P->params, "inv_v") || + pj_param_exists(P->params, "inv_origin"); } else { - has_only_fwd = - !pj_param_exists(P->params, "inv_c") && - !pj_param_exists(P->params, "inv_origin"); + has_inv = + pj_param_exists(P->params, "inv_c") || + pj_param_exists(P->params, "inv_origin"); } - Q->has_only_fwd = has_only_fwd; + Q->has_inv = has_inv; if (complex_polynomia) { /* Westings and/or southings? */ @@ -612,7 +620,7 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_c")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_fwd && 0==parse_coefs (P, Q->inv_c, "inv_c", n)) + if (has_inv && 0==parse_coefs (P, Q->inv_c, "inv_c", n)) { proj_log_error (P, _("missing inv_c")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); @@ -633,12 +641,12 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_fwd && 0==parse_coefs (P, Q->inv_u, "inv_u", n)) + if (has_inv && 0==parse_coefs (P, Q->inv_u, "inv_u", n)) { proj_log_error (P, _("missing inv_u")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_fwd && 0==parse_coefs (P, Q->inv_v, "inv_v", n)) + if (has_inv && 0==parse_coefs (P, Q->inv_v, "inv_v", n)) { proj_log_error (P, _("missing inv_v")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); @@ -650,7 +658,7 @@ PJ *PROJECTION(horner) { proj_log_error (P, _("missing fwd_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); } - if (!has_only_fwd && 0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) + if (has_inv && 0==parse_coefs (P, (double *)(Q->inv_origin), "inv_origin", 2)) { proj_log_error (P, _("missing inv_origin")); return horner_freeup (P, PROJ_ERR_INVALID_OP_MISSING_ARG); From a02f426e9a95cd885b2681b4b8524e21de264f2a Mon Sep 17 00:00:00 2001 From: georgeouzou Date: Sat, 16 Apr 2022 13:51:40 +0300 Subject: [PATCH 10/10] Add doc note on reversible polynomials Cite IOGP Publication on reversible polynomials and the difference with our general solution --- docs/source/operations/transformations/horner.rst | 3 +++ docs/source/references.bib | 11 +++++++++++ 2 files changed, 14 insertions(+) diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst index 9eb8d39c8e..b35b2b367e 100644 --- a/docs/source/operations/transformations/horner.rst +++ b/docs/source/operations/transformations/horner.rst @@ -117,6 +117,9 @@ be done by iteratively solving a system of equations. By writing :eq:`real_poly` We can iteratively solve with initial values of :math:`U = 0` and :math:`V = 0` and find :math:`U` and :math:`V`. +.. note:: + This iterative inverse transformation is a more general solution to *reversible polynormials of degree n* as presented in :cite:`IOGP2019`. These can provide a satisfactory solution in a single step when certain conditions are met. + Evaluation of the complex polynomials are defined by the following equations: .. math:: diff --git a/docs/source/references.bib b/docs/source/references.bib index cda584026b..3584a979a9 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -169,6 +169,17 @@ @TechReport{IOGP2018 Url = {https://www.iogp.org/bookstore/product/coordinate-conversions-and-transformation-including-formulas/} } +@TechReport{IOGP2019, + Title = {Geomatics Guidance Note 7, part 2: Coordinate Conversions \& Transformations including Formulas}, + Author = {IOGP}, + Institution = {International Association For Oil And Gas Producers}, + Year = {2019}, + Number = {373-7-2}, + Type = {IOGP Publication}, + + Url = {https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf} +} + @TechReport{ISO19111, Title = {{Geographic information -- Referencing by coordinates}}, Author = {ISO},