Introduce preparation/finalization steps in fwd/inv subsystem, supporting arbitrary dimensionality in test code

* Call trans func of same dimensionality as input in gie

* Refactor prep/fin code for pj_fwd/pj_inv 2D,3D,4D

* Remove prime meridian handling from pj_transform (now handled in pj_fwd_prepare/pj_inv_finalize)

* Introduce prep/fin skips, mostly in support of axisswap and pipeline drivers

* Refactor fwd/inv subsystem

* pj_transform: Let pj_fwd/inv handle scaling

* Let pj_fwd/inv3d fall back to 2D eventually

1 files changed, 150 insertions, 41 deletions

diff --git a/src/pj_inv.c b/src/pj_inv.c
index 68a5595b..e1fb05f6 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -1,60 +1,169 @@
-/* general inverse projection */
-#define PJ_LIB__
-#include <proj.h>
-#include <projects.h>
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Inverse operation invocation
+ * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2018-01-02
+ *           Based on material from Gerald Evenden (original pj_inv)
+ *           and Piyush Agram (original pj_inv3d)
+ *
+ ******************************************************************************
+ * Copyright (c) 2000, Frank Warmerdam
+ * Copyright (c) 2018, Thomas Knudsen / SDFE
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
 #include <errno.h>
-# define EPS 1.0e-12
 
-/* inverse projection entry */
-LP pj_inv(XY xy, PJ *P) {
-    LP lp;
-    LP err;
-    int last_errno;
+#include "proj_internal.h"
+#include "projects.h"
 
-    /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
-    err.lam = err.phi = HUGE_VAL;
 
-    if (0==P->inv)
-        return err;
 
-    /* can't do as much preliminary checking as with forward */
-    if (xy.x == HUGE_VAL || xy.y == HUGE_VAL) {
-        pj_ctx_set_errno( P->ctx, -15);
-        return err;
+PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
+    if (coo.xyz.x == HUGE_VAL) {
+        proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
+        return proj_coord_error ();
     }
 
-    last_errno = proj_errno_reset (P);
-
     /* de-scale and de-offset */
-    xy.x = (xy.x * P->to_meter - P->x0);
-    xy.y = (xy.y * P->to_meter - P->y0);
+    coo.xyz.x = (coo.xyz.x * P->to_meter  - P->x0);
+    coo.xyz.y = (coo.xyz.y * P->to_meter  - P->y0);
+    coo.xyz.z = (coo.xyz.z * P->vto_meter - P->z0);
 
-    /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis */
-    /* Multiplying by ra, rather than dividing by a because the CALCOFI projection      */
-    /* stomps on a and hence depends on this */
+    /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
+    /* Multiplying by ra, rather than dividing by a because the CALCOFI projection       */
+    /* stomps on a and hence (apparently) depends on this to roundtrip correctly         */
+    /* (CALCOFI avoids further scaling by stomping - a better solution must be possible) */
     if (P->right==PJ_IO_UNITS_CLASSIC) {
-        xy.x *= P->ra;
-        xy.y *= P->ra;
+        coo.xyz.x *= P->ra;
+        coo.xyz.y *= P->ra;
     }
 
-    /* Do inverse transformation */
-    lp = (*P->inv) (xy, P);
-    if (P->ctx->last_errno)
-        return err;
+    return coo;
+}
+
+
+
+PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
+    if (coo.xyz.x == HUGE_VAL) {
+        proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
+        return proj_coord_error ();
+    }
 
     if (P->left==PJ_IO_UNITS_RADIANS) {
-        /* reduce from del lp.lam */
-        lp.lam += P->lam0;
+        /* Distance from central meridian, taking system zero meridian into account */
+        coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
 
          /* adjust longitude to central meridian */
-        if (!P->over)
-            lp.lam = adjlon(lp.lam);
+        if (0==P->over)
+            coo.lpz.lam = adjlon(coo.lpz.lam);
+
+        /* If input latitude was geocentrical, convert back to geocentrical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_FWD, coo);
+    }
+
+    return coo;
+}
+
+
+
+LP pj_inv(XY xy, PJ *P) {
+    PJ_COORD coo = {{0,0,0,0}};
+    coo.xy = xy;
+
+    if (!P->skip_inv_prepare)
+        coo = pj_inv_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().lp;
+
+    /* Do the transformation, using the lowest dimensional transformer available */
+    if (P->inv)
+        coo.lp = P->inv(coo.xy, P);
+    else if (P->inv3d)
+        coo.lpz = P->inv3d (coo.xyz, P);
+    else if (P->inv4d)
+        coo = P->inv4d (coo, P);
+    else {
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ().lp;
+    }
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().lp;
+
+    if (!P->skip_inv_finalize)
+        coo = pj_inv_finalize (P, coo);
+    return coo.lp;
+}
+
+
+
+LPZ pj_inv3d (XYZ xyz, PJ *P) {
+    PJ_COORD coo = {{0,0,0,0}};
+    coo.xyz = xyz;
+
+    if (!P->skip_inv_prepare)
+        coo = pj_inv_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().lpz;
+
+    /* Do the transformation, using the lowest dimensional transformer feasible */
+    if (P->inv3d)
+        coo.lpz = P->inv3d (coo.xyz, P);
+    else if (P->inv4d)
+        coo = P->inv4d (coo, P);
+    else if (P->inv)
+        coo.lp = P->inv (coo.xy, P);
+    else {
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ().lpz;
+    }
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().lpz;
+
+    if (!P->skip_inv_finalize)
+        coo = pj_inv_finalize (P, coo);
+    return coo.lpz;
+}
+
+
+
+PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
+    if (!P->skip_inv_prepare)
+        coo = pj_inv_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ();
 
-        /* This may be redundant and never used */
-        if (P->geoc && fabs(fabs(lp.phi)-M_HALFPI) > EPS)
-            lp.phi = atan(P->one_es * tan(lp.phi));
+    /* Call the highest dimensional converter available */
+    if (P->inv4d)
+        coo = P->inv4d (coo, P);
+    else if (P->inv3d)
+        coo.lpz = P->inv3d (coo.xyz, P);
+    else if (P->inv)
+        coo.lp = P->inv (coo.xy, P);
+    else {
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ();
     }
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ();
 
-    proj_errno_restore (P, last_errno);
-    return lp;
-}
\ No newline at end of file
+    if (!P->skip_inv_finalize)
+        coo = pj_inv_finalize (P, coo);
+    return coo;
+}