Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add invertible Horner polynomials #3133

Merged
merged 10 commits into from
Apr 19, 2022
Merged

Add invertible Horner polynomials #3133

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
3159792
Add inversible horner real polynomials
georgeouzou Mar 23, 2022
f14feeb
Fix windows compilation error
georgeouzou Mar 24, 2022
12c3377
Handle some inversible horner corner cases
georgeouzou Mar 24, 2022
04baa72
Add documentation for iterative real polynomial inversion
georgeouzou Mar 25, 2022
3dbe701
Address docs review comments
georgeouzou Mar 26, 2022
98e4a89
Add inversible horner complex polynomials
georgeouzou Mar 26, 2022
16aab79
Optimize iterative inversion of complex horner
georgeouzou Mar 26, 2022
1b556a8
Fix sanitizer errors
georgeouzou Mar 26, 2022
d75519b
Apply review comments
georgeouzou Mar 31, 2022
a02f426
Add doc note on reversible polynomials
georgeouzou Apr 16, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 85 additions & 14 deletions docs/source/operations/transformations/horner.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@ Horner polynomial evaluation
================================================================================

.. versionadded:: 5.0.0
.. versionchanged:: 9.1.0 Iterative real polynormal inversion

+-----------------+-------------------------------------------------------------------+
| **Alias** | horner |
Expand Down Expand Up @@ -61,7 +62,60 @@ 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`.

Evaluation of the complex polynomials are defined by the following equations:

Expand Down Expand Up @@ -133,7 +187,7 @@ describing real and complex polynomials can't be mixed.

.. option:: +inv_origin=<northing,easting>

Coordinate of origin for the inverse mapping
Coordinate of origin for the inverse mapping. Not required for iterative inversion.
georgeouzou marked this conversation as resolved.
Show resolved Hide resolved

Real polynomials
..............................................................................
Expand All @@ -157,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=<u_11,u_12,...,u_ij,..,u_mn>

Coefficients for the inverse transformation i.e. latitude to northing
as described in :eq:`real_poly`.

.. option:: +inv_v=<v_11,v_12,...,v_ij,..,v_mn>

Coefficients for the inverse transformation i.e. longitude to easting
as described in :eq:`real_poly`.



Complex polynomials
..............................................................................

Expand All @@ -194,6 +236,24 @@ of the polynomial:
Optional
-------------------------------------------------------------------------------

.. option:: +inv_u=<u_11,u_12,...,u_ij,..,u_mn>

.. 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=<v_11,v_12,...,v_ij,..,v_mn>

.. 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=<value>

Radius of the region of validity.
Expand All @@ -206,6 +266,17 @@ Optional

Express longitude as westing. Only applies for complex polynomials.

.. option:: +inv_tolerance=<value>

georgeouzou marked this conversation as resolved.
Show resolved Hide resolved
.. 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 <value>.
<value> should be the same units as :math:`U` and :math:`V` of :eq:`UV`

*Defaults to 0.001.*

Further reading
################################################################################
Expand Down
98 changes: 87 additions & 11 deletions src/transformations/horner.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand Down Expand Up @@ -220,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;
Expand All @@ -244,29 +245,93 @@ 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;
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;
} 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;
/* 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++;
}
georgeouzou marked this conversation as resolved.
Show resolved Hide resolved
}
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;
}
}

/* 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;
Expand Down Expand Up @@ -443,6 +508,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;
Expand Down Expand Up @@ -474,6 +540,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;
Expand Down Expand Up @@ -506,12 +580,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);
Expand All @@ -523,13 +597,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;
}
74 changes: 74 additions & 0 deletions test/unit/gie_self_tests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down