diff --git a/docs/source/operations/transformations/horner.rst b/docs/source/operations/transformations/horner.rst index ad09e8abd0..b35b2b367e 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 +.. versionchanged:: 9.1.0 Iterative polynormal inversion +-----------------+-------------------------------------------------------------------+ | **Alias** | horner | @@ -61,7 +62,63 @@ 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:: + \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} \\ + +.. math:: + \begin{bmatrix} + \Delta X \\ + \Delta Y \\ + \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} \\ + +.. math:: + \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} + +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: @@ -72,7 +129,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 ################################################################################ @@ -131,10 +188,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 - Real polynomials .............................................................................. @@ -157,18 +210,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`. - -.. option:: +inv_v= - - Coefficients for the inverse transformation i.e. longitude to easting - as described in :eq:`real_poly`. - - - Complex polynomials .............................................................................. @@ -186,13 +227,43 @@ of the polynomial: Coefficients for the complex forward 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 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 polynomial evaluation is used for + the inverse tranformation. + .. option:: +inv_c= - Coefficients for the complex inverse transformation - as described in :eq:`complex_poly`. + .. versionchanged:: 9.1.0 -Optional -------------------------------------------------------------------------------- + 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= @@ -206,6 +277,17 @@ Optional Express longitude as westing. Only applies for complex polynomials. +.. option:: +inv_tolerance= + + .. versionadded:: 9.1.0 + + 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 . + should be the same units as :math:`U` and :math:`V` of :eq:`UV` + + *Defaults to 0.001.* Further reading ################################################################################ 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}, diff --git a/src/transformations/horner.cpp b/src/transformations/horner.cpp index 7c8ad19230..19b657c8f5 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,6 +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 */ + 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 */ @@ -177,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) { @@ -213,14 +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 *tcx, *tcy; /* Coefficient pointers */ double range; /* Equivalent to the gen_pol's FLOATLIMIT constant */ double n, e; PJ_UV uv_error; @@ -241,49 +297,86 @@ 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_inv; 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 */ - 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) { + 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; - } - - /* 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; - - /* 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; + } 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; + 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) + bool converged = false; + while (loops-- > 0 && !converged) { + double Ma = 0.0; + double Mb = 0.0; + double Mc = 0.0; + double Md = 0.0; + { + 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); } - N = e*N + u; - E = n*E + v; + double idet = 1.0 / (Ma*Md - Mb*Mc); + double x = idet * (Md*de - Mb*dn); + double y = idet * (Ma*dn - Mc*de); + converged = (fabs(x-x0) < tol) && (fabs(y-y0) < tol); + x0 = x; + y0 = y; + } + // 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; } - - position.u = E; - position.v = N; + } + else { + 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; @@ -318,10 +411,8 @@ 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; + double n, e; PJ_UV uv_error; uv_error.u = uv_error.v = HUGE_VAL; @@ -340,12 +431,11 @@ polynomial evaluation engine. } /* Prepare for double Horner */ - sz = 2*transformation->order + 2; range = transformation->range; + const bool iterative_inverse = direction == PJ_INV && !transformation->has_inv; + 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) @@ -353,14 +443,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)) { @@ -368,17 +462,45 @@ polynomial evaluation engine. return uv_error; } - /* Everything's set up properly - now do the actual polynomium evaluation */ - E = *--c; - N = *--c; - while (c > cb) { - w = n*E + e*N + *--c; - N = n*N - e*E + *--c; - E = w; + 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); + int loops = 32; // usually converges really fast (1-2 loops) + bool converged = false; + while (loops-- > 0 && !converged) { + // 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; + } + // 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 { + double E = w0.imag(); + double 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; } - position.u = E; - position.v = N; + // coefficient pointers + double *cb = direction == PJ_FWD ? transformation->fwd_c : transformation->inv_c; + PJ_UV en = { e, n }; + position = complex_horner_eval(transformation->order, cb, en); return position; } @@ -443,6 +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; + bool has_inv = false; HORNER *Q; P->fwd4d = horner_forward_4d; P->inv4d = horner_reverse_4d; @@ -474,6 +597,18 @@ PJ *PROJECTION(horner) { return horner_freeup (P, PROJ_ERR_OTHER /*ENOMEM*/); P->opaque = Q; + if (!complex_polynomia) { + 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_inv = + pj_param_exists(P->params, "inv_c") || + pj_param_exists(P->params, "inv_origin"); + } + Q->has_inv = has_inv; + if (complex_polynomia) { /* Westings and/or southings? */ Q->uneg = pj_param_exists (P->params, "uneg") ? 1 : 0; @@ -485,7 +620,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_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); @@ -506,12 +641,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_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 (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); @@ -523,13 +658,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_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); } 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..c8933ea21a 100644 --- a/test/unit/gie_self_tests.cpp +++ b/test/unit/gie_self_tests.cpp @@ -763,6 +763,119 @@ 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 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" + " +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); + } + + { + 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); + + proj_destroy(P); + } +} + // --------------------------------------------------------------------------- TEST(gie, proj_create_crs_to_crs_PULKOVO42_ETRS89) {