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

20 files changed, 539 insertions, 380 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 096cc672..e7e53394 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -70,7 +70,7 @@ libproj_la_SOURCES = \
 	aasincos.c adjlon.c bch2bps.c bchgen.c \
 	biveval.c dmstor.c mk_cheby.c pj_auth.c \
 	pj_deriv.c pj_ell_set.c pj_ellps.c pj_errno.c \
-	pj_factors.c pj_fwd.c pj_init.c pj_inv.c pj_fwd3d.c pj_inv3d.c\
+	pj_factors.c pj_fwd.c pj_init.c pj_inv.c \
 	pj_list.c pj_malloc.c pj_mlfn.c pj_msfn.c proj_mdist.c \
 	pj_open_lib.c pj_param.c pj_phi2.c pj_pr_list.c \
 	pj_qsfn.c pj_strerrno.c \
diff --git a/src/PJ_axisswap.c b/src/PJ_axisswap.c
index f8f17380..aee3c56e 100644
--- a/src/PJ_axisswap.c
+++ b/src/PJ_axisswap.c
@@ -173,6 +173,16 @@ PJ *CONVERSION(axisswap,0) {
         return pj_default_destructor(P, PJD_ERR_MISSING_ARGS);
     }
 
+    /* Preparation and finalization steps are skipped, since the raison   */
+    /* d'etre of axisswap is to bring input coordinates in line with the  */
+    /* the internally expected order (ENU), such that handling of offsets */
+    /* etc. can be done correctly in a later step of a pipeline */
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+
+
     /* fill axis list with indices from 4-7 to simplify duplicate search further down */
     for (i=0; i<4; i++)
         Q->axis[i] = i+4;
diff --git a/src/PJ_latlong.c b/src/PJ_latlong.c
index 5919023a..b1909954 100644
--- a/src/PJ_latlong.c
+++ b/src/PJ_latlong.c
@@ -30,7 +30,6 @@
 /* very loosely based upon DMA code by Bradford W. Drew */
 #define PJ_LIB__
 #include "proj_internal.h"
-#include <proj.h>
 #include "projects.h"
 
 PROJ_HEAD(lonlat, "Lat/long (Geodetic)")  "\n\t";
@@ -39,87 +38,87 @@ PROJ_HEAD(latlong, "Lat/long (Geodetic alias)")  "\n\t";
 PROJ_HEAD(longlat, "Lat/long (Geodetic alias)")  "\n\t";
 
 
- static XY forward(LP lp, PJ *P) {
+ static XY latlong_forward(LP lp, PJ *P) {
     XY xy = {0.0,0.0};
-    xy.x = lp.lam / P->a;
-    xy.y = lp.phi / P->a;
+    (void) P;
+    xy.x = lp.lam;
+    xy.y = lp.phi;
     return xy;
 }
 
 
-static LP inverse(XY xy, PJ *P) {
+static LP latlong_inverse(XY xy, PJ *P) {
     LP lp = {0.0,0.0};
-    lp.phi = xy.y * P->a;
-    lp.lam = xy.x * P->a;
+    (void) P;
+    lp.phi = xy.y;
+    lp.lam = xy.x;
     return lp;
 }
 
-static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
+
+ static XYZ latlong_forward_3d (LPZ lpz, PJ *P) {
+    XYZ xyz = {0,0,0};
+    (void) P;
+    xyz.x = lpz.lam;
+    xyz.y = lpz.phi;
+    xyz.z = lpz.z;
+    return xyz;
+}
+
+
+static LPZ latlong_inverse_3d (XYZ xyz, PJ *P) {
+    LPZ lpz = {0,0,0};
+    (void) P;
+    lpz.lam = xyz.x;
+    lpz.phi = xyz.y;
+    lpz.z   = xyz.z;
+    return lpz;
+}
+
+static PJ_COORD latlong_forward_4d (PJ_COORD obs, PJ *P) {
     (void) P;
     return obs;
 }
 
-static PJ_COORD inverse_4d(PJ_COORD obs, PJ *P) {
+
+static PJ_COORD latlong_inverse_4d (PJ_COORD obs, PJ *P) {
     (void) P;
     return obs;
 }
 
-PJ *PROJECTION(latlong) {
+
+
+static PJ *latlong_setup (PJ *P) {
     P->is_latlong = 1;
-    P->x0 = 0.0;
-    P->y0 = 0.0;
-    P->inv = inverse;
-    P->fwd = forward;
-    P->inv4d = inverse_4d;
-    P->fwd4d = forward_4d;
+    P->x0 = 0;
+    P->y0 = 0;
+    P->inv = latlong_inverse;
+    P->fwd = latlong_forward;
+    P->inv3d = latlong_inverse_3d;
+    P->fwd3d = latlong_forward_3d;
+    P->inv4d = latlong_inverse_4d;
+    P->fwd4d = latlong_forward_4d;
     P->left = PJ_IO_UNITS_RADIANS;
     P->right = PJ_IO_UNITS_RADIANS;
-
     return P;
 }
 
+PJ *PROJECTION(latlong) {
+    return latlong_setup (P);
+}
 
-PJ *PROJECTION(longlat) {
-    P->is_latlong = 1;
-    P->x0 = 0.0;
-    P->y0 = 0.0;
-    P->inv = inverse;
-    P->fwd = forward;
-    P->inv4d = inverse_4d;
-    P->fwd4d = forward_4d;
-    P->left = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
 
-    return P;
+PJ *PROJECTION(longlat) {
+    return latlong_setup (P);
 }
 
 
 PJ *PROJECTION(latlon) {
-    P->is_latlong = 1;
-    P->x0 = 0.0;
-    P->y0 = 0.0;
-    P->inv = inverse;
-    P->fwd = forward;
-    P->inv4d = inverse_4d;
-    P->fwd4d = forward_4d;
-    P->left = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
-
-    return P;
+    return latlong_setup (P);
 }
 
 
 PJ *PROJECTION(lonlat) {
-    P->is_latlong = 1;
-    P->x0 = 0.0;
-    P->y0 = 0.0;
-    P->inv = inverse;
-    P->fwd = forward;
-    P->inv4d = inverse_4d;
-    P->fwd4d = forward_4d;
-    P->left = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
-
-    return P;
+    return latlong_setup (P);
 }
 
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index 4ce4bd4d..43211fe7 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -170,11 +170,6 @@ PJ *PROJECTION(ob_tran) {
     P->opaque = Q;
     P->destructor = destructor;
 
-#if 0
-    if (0 != P->es)
-        return destructor(P, PJD_ERR_ELLIPSOIDAL_UNSUPPORTED);
-#endif
-
     /* get name of projection to be translated */
     if (!(name = pj_param(P->ctx, P->params, "so_proj").s))
         return destructor(P, PJD_ERR_NO_ROTATION_PROJ);
@@ -235,6 +230,12 @@ PJ *PROJECTION(ob_tran) {
         P->inv = Q->link->inv ? t_inverse : 0;
     }
 
+    /* Support some rather speculative test cases, where the rotated projection */
+    /* is actually latlong. We do not want scaling in that case... */
+    if (Q->link->right==PJ_IO_UNITS_RADIANS)
+        P->right = PJ_IO_UNITS_METERS;
+
+
     return P;
 }
 
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index 2b012193..af9c5394 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -72,7 +72,7 @@
 Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 
 ********************************************************************************
-* Copyright (c) 2016, Thomas Knudsen / SDFE
+* Copyright (c) 2016, 2017, 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"),
@@ -96,6 +96,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 
 #define PJ_LIB__
 #include <geodesic.h>
+#include <proj.h>
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -107,7 +108,6 @@ PROJ_HEAD(pipeline,         "Transformation pipeline manager");
 struct pj_opaque {
     int reversible;
     int steps;
-    int verbose;
     char **argv;
     char **current_argv;
     PJ **pipeline;
@@ -124,52 +124,6 @@ static LP     pipeline_reverse (XY xyz, PJ *P);
 
 
 
-/********************************************************************
-
-                  ISOMORPHIC TRANSFORMATIONS
-
-*********************************************************************
-
-    In this context, an isomorphic transformation is a proj PJ
-    object returning the same kind of coordinates that it
-    receives, i.e. a transformation from angular to angular or
-    linear to linear coordinates.
-
-    The degrees-to-radians operation is an example of the former,
-    while the latter is typical for most of the datum shift
-    operations used in geodesy, e.g. the Helmert 7-parameter shift.
-
-    Isomorphic transformations trips up the pj_inv/pj_fwd
-    functions, which try to check input sanity and scale output to
-    internal proj units, under the assumption that if input is of
-    angular type, then output must be of linear (or vice versa).
-
-    Hence, to avoid having pj_inv/pj_fwd stomping on output (or
-    choking on input), we need a way to tell them that they should
-    accept whatever is being handed to them.
-
-    The P->left and P->right elements indicate isomorphism.
-
-    For classic proj style projections, P->left has the value
-    PJ_IO_UNITS_RADIANS, while P->right has the value
-    PJ_IO_UNITS_CLASSIC, indicating that the forward driver expects
-    angular input coordinates, and provides linear output coordinates,
-    scaled by the P->a semimajor axis length.
-
-    Newer projections may set P->left and P->right to either
-    PJ_IO_UNITS_METERS, PJ_IO_UNITS_RADIANS or PJ_IO_UNITS_ANY,
-    to indicate their I/O style.
-
-    For the forward driver, left indicates input coordinate
-    type, while right indicates output coordinate type.
-
-    For the inverse driver, left indicates output coordinate
-    type, while right indicates input coordinate type.
-
-*********************************************************************/
-
-
-
 
 static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) {
     int i, first_step, last_step;
@@ -184,7 +138,6 @@ static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) {
 }
 
 
-
 static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
     int i, first_step, last_step;
 
@@ -198,37 +151,59 @@ static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
 }
 
 
-/* Delegate the work to pipeline_forward_4d() */
+
+
 static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
+    int i;
     point.lpz = lpz;
-    point = pipeline_forward_4d (point, P);
+
+    for (i = 1;  i <= P->opaque->steps;  i++)
+        point = pj_approx_3D_trans (P->opaque->pipeline[i], 1, point);
+
     return point.xyz;
 }
 
-/* Delegate the work to pipeline_reverse_4d() */
+
 static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
+    int i;
     point.xyz = xyz;
-    point = pipeline_reverse_4d (point, P);
+
+    for (i = P->opaque->steps;  i > 0 ;  i--)
+        point = pj_approx_3D_trans (P->opaque->pipeline[i], -1, point);
+
     return point.lpz;
 }
 
+
+
+
 static XY pipeline_forward (LP lp, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
+    int i;
     point.lp = lp;
-    point = pipeline_forward_4d (point, P);
+
+    for (i = 1;  i <= P->opaque->steps;  i++)
+        point = pj_approx_2D_trans (P->opaque->pipeline[i], 1, point);
+
     return point.xy;
 }
 
+
 static LP pipeline_reverse (XY xy, PJ *P) {
     PJ_COORD point = {{0,0,0,0}};
+    int i;
     point.xy = xy;
-    point = pipeline_reverse_4d (point, P);
+    for (i = P->opaque->steps;  i > 0 ;  i--)
+        point = pj_approx_2D_trans (P->opaque->pipeline[i], -1, point);
+
     return point.lp;
 }
 
 
+
+
 static void *destructor (PJ *P, int errlev) {
     int i;
     if (0==P)
@@ -240,8 +215,7 @@ static void *destructor (PJ *P, int errlev) {
     /* Deallocate each pipeine step, then pipeline array */
     if (0!=P->opaque->pipeline)
         for (i = 0;  i < P->opaque->steps; i++)
-            if (0!=P->opaque->pipeline[i+1])
-                P->opaque->pipeline[i+1]->destructor (P->opaque->pipeline[i+1], errlev);
+            proj_destroy (P->opaque->pipeline[i+1]);
     pj_dealloc (P->opaque->pipeline);
 
     pj_dealloc (P->opaque->argv);
@@ -264,6 +238,8 @@ static PJ *pj_create_pipeline (PJ *P, size_t steps) {
 }
 
 
+
+
 /* count the number of args in pipeline definition */
 static size_t argc_params (paralist *params) {
     size_t argc = 0;
@@ -288,6 +264,9 @@ static char **argv_params (paralist *params, size_t argc) {
     return argv;
 }
 
+
+
+
 /* Being the special operator that the pipeline is, we have to handle the    */
 /* ellipsoid differently than usual. In general, the pipeline operation does */
 /* not need an ellipsoid, but in some cases it is beneficial nonetheless.    */
@@ -335,18 +314,29 @@ static void set_ellipsoid(PJ *P) {
 }
 
 
+
+
 PJ *OPERATION(pipeline,0) {
     int i, nsteps = 0, argc;
     int i_pipeline = -1, i_first_step = -1, i_current_step;
     char **argv, **current_argv;
 
-    P->fwd4d  = pipeline_forward_4d;
-    P->inv4d  = pipeline_reverse_4d;
+    P->fwd4d  =  pipeline_forward_4d;
+    P->inv4d  =  pipeline_reverse_4d;
     P->fwd3d  =  pipeline_forward_3d;
     P->inv3d  =  pipeline_reverse_3d;
     P->fwd    =  pipeline_forward;
     P->inv    =  pipeline_reverse;
     P->destructor  =  destructor;
+    P->is_pipeline =  1;
+
+    /* Currently, the pipeline driver is a raw bit mover, enabling other operations */
+    /* to collaborate efficiently. All prep/fin stuff is done at the step levels. */
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+
 
     P->opaque = pj_calloc (1, sizeof(struct pj_opaque));
     if (0==P->opaque)
@@ -376,7 +366,7 @@ PJ *OPERATION(pipeline,0) {
 
         if (0==strcmp ("proj=pipeline", argv[i])) {
             if (-1 != i_pipeline) {
-                proj_log_error (P, "Pipeline: Nesting only allowed when child pipelines are wrapped in +init's");
+                proj_log_error (P, "Pipeline: Nesting only allowed when child pipelines are wrapped in '+init's");
                 return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: nested pipelines */
             }
             i_pipeline = i;
diff --git a/src/gie.c b/src/gie.c
index 8a1d356e..c6440a2d 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -184,6 +184,7 @@ typedef struct {
     int total_ok, total_ko;
     int grand_ok, grand_ko;
     size_t operation_lineno;
+    size_t dimensions_given, dimensions_given_at_last_accept;
     double tolerance;
     const char *curr_file;
     FILE *fout;
@@ -581,10 +582,13 @@ Attempt to interpret args as a PJ_COORD.
     const char *endp, *prev = args;
     PJ_COORD a = proj_coord (0,0,0,0);
 
-    for (i = 0; i < 4; i++) {
+    for (i = 0, T.dimensions_given = 0;   i < 4;   i++, T.dimensions_given++) {
         double d = proj_strtod (prev,  (char **) &endp);
+
+        /* Break out if there were no more numerals */
         if (prev==endp)
             return i > 1? a: proj_coord_error ();
+
         a.v[i] = d;
         prev = endp;
     }
@@ -601,6 +605,7 @@ Read ("ACCEPT") a 2, 3, or 4 dimensional input coordinate.
     T.a = parse_coord (args);
     if (T.verbosity > 3)
         printf ("#  %s\n", args);
+    T.dimensions_given_at_last_accept = T.dimensions_given;
     return 0;
 }
 
@@ -693,6 +698,19 @@ static int expect_failure_with_errno_message (int expected, int got) {
 }
 
 
+/* For test purposes, we want to call a transformation of the same */
+/* dimensionality as the number of dimensions given in accept */
+static PJ_COORD expect_trans_n_dim (PJ_COORD ci) {
+    if (4==T.dimensions_given_at_last_accept)
+        return proj_trans (T.P, T.dir, ci);
+
+    if (3==T.dimensions_given_at_last_accept)
+        return pj_approx_3D_trans (T.P, T.dir, ci);
+
+    return pj_approx_2D_trans (T.P, T.dir, ci);
+}
+
+
 /*****************************************************************************/
 static int expect (const char *args) {
 /*****************************************************************************
@@ -733,7 +751,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
 
         /* Try to carry out the operation - and expect failure */
         ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
-        co = proj_trans (T.P, T.dir, ci);
+        co = expect_trans_n_dim (ci);
 
         /* Failed to fail? - that's a failure */
         if (co.xyz.x!=HUGE_VAL)
@@ -775,7 +793,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         printf ("ACCEPTS  %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
 
     /* angular output from proj_trans comes in radians */
-    co = proj_trans (T.P, T.dir, ci);
+    co = expect_trans_n_dim (ci);
     T.b = proj_angular_output (T.P, T.dir)? todeg_coord (co): co;
     if (T.verbosity > 3)
         printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 88d88a97..f1337a54 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -185,7 +185,6 @@ SET(SRC_LIBPROJ_CORE
         pj_errno.c
         pj_factors.c
         pj_fwd.c
-        pj_fwd3d.c
         pj_gauss.c
         pj_gc_reader.c
         pj_geocent.c
@@ -196,7 +195,6 @@ SET(SRC_LIBPROJ_CORE
         pj_init.c
         pj_initcache.c
         pj_inv.c
-        pj_inv3d.c
         pj_list.c
         pj_list.h
         pj_log.c
diff --git a/src/makefile.vc b/src/makefile.vc
index c9ebd24c..ee89beb3 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -45,7 +45,6 @@ support = \
 	biveval.obj dmstor.obj mk_cheby.obj pj_auth.obj \
 	pj_deriv.obj pj_ell_set.obj pj_ellps.obj pj_errno.obj \
 	pj_factors.obj pj_fwd.obj pj_init.obj pj_inv.obj \
-    pj_fwd3d.obj pj_inv3d.obj \
 	pj_list.obj pj_malloc.obj pj_mlfn.obj pj_msfn.obj \
 	pj_open_lib.obj pj_param.obj pj_phi2.obj pj_pr_list.obj \
 	pj_qsfn.obj pj_strerrno.obj pj_tsfn.obj pj_units.obj \
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index e010f6ec..a1025aac 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -1,58 +1,176 @@
-/* general forward projection */
-#define PJ_LIB__
-#include <proj.h>
-#include <projects.h>
-# define EPS 1.0e-12
-	XY /* forward projection entry */
-pj_fwd(LP lp, PJ *P) {
-	XY xy;
-    XY err;
-    double t;
-    int last_errno;
-
-    /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
-    err.x = err.y = HUGE_VAL;
-
-    if (0==P->fwd)
-        return err;
-    last_errno = proj_errno_reset (P);
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Forward operation invocation
+ * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2018-01-02
+ *           Based on material from Gerald Evenden (original pj_fwd)
+ *           and Piyush Agram (original pj_fwd3d)
+ *
+ ******************************************************************************
+ * 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>
+
+#include "proj_internal.h"
+#include "projects.h"
+
+PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ();
 
     /* Check validity of angular input coordinates */
     if (P->left==PJ_IO_UNITS_RADIANS) {
+        double t;
+        LP lp = coo.lp;
 
-        /* check for forward and latitude or longitude overange */
-        t = fabs(lp.phi)-M_HALFPI;
-        if (t > EPS || fabs(lp.lam) > 10.) {
-            pj_ctx_set_errno( P->ctx, -14);
-            return err;
+        /* check for latitude or longitude over-range */
+        t = (lp.phi < 0  ?  -lp.phi  :  lp.phi) - M_HALFPI;
+        if (t > PJ_EPS_LAT  ||  lp.lam > 10  ||  lp.lam < -10) {
+            proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+            return proj_coord_error ();
         }
 
         /* Clamp latitude to -90..90 degree range */
-        if (fabs(t) <= EPS)
-            lp.phi = lp.phi < 0. ? -M_HALFPI : M_HALFPI;
-        else if (P->geoc)   /* Maybe redundant and never used. */
-            lp.phi = atan(P->rone_es * tan(lp.phi));
-        lp.lam -= P->lam0;    /* compute del lp.lam */
-        if (!P->over)
-            lp.lam = adjlon(lp.lam); /* adjust del longitude */
+        if (lp.phi > M_HALFPI)
+            lp.phi = M_HALFPI;
+        if (lp.phi < -M_HALFPI)
+            lp.phi = -M_HALFPI;
+
+        /* If input latitude is geocentrical, convert to geographical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_INV, coo);
+
+        /* Distance from central meridian, taking system zero meridian into account */
+        lp.lam = (lp.lam - P->from_greenwich) - P->lam0;
+
+        /* Ensure longitude is in the -pi:pi range */
+        if (0==P->over)
+            lp.lam = adjlon(lp.lam);
+
+        coo.lp = lp;
     }
 
-    /* Do the transformation */
-    xy = (*P->fwd)(lp, P);
-    if ( proj_errno (P) )
-        return err;
+    return coo;
+}
+
+
+PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
 
     /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
     if (P->right==PJ_IO_UNITS_CLASSIC) {
-        xy.x *= P->a;
-        xy.y *= P->a;
+        coo.xy.x *= P->a;
+        coo.xy.y *= P->a;
     }
 
     /* Handle false eastings/northings and non-metric linear units */
-    xy.x = P->fr_meter * (xy.x + P->x0);
-    xy.y = P->fr_meter * (xy.y + P->y0);
-    /* z is not scaled since this is handled by vto_meter outside */
+    coo.xyz.x = P->fr_meter  * (coo.xyz.x + P->x0);
+    coo.xyz.y = P->fr_meter  * (coo.xyz.y + P->y0);
+    coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+
+    return coo;
+}
+
+
+
+XY pj_fwd(LP lp, PJ *P) {
+    PJ_COORD coo = {{0,0,0,0}};
+    coo.lp = lp;
+
+    if (!P->skip_fwd_prepare)
+        coo = pj_fwd_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().xy;
+
+    /* Do the transformation, using the lowest dimensional transformer available */
+    if (P->fwd)
+        coo.xy = P->fwd(coo.lp, P);
+    else if (P->fwd3d)
+        coo.xyz = P->fwd3d (coo.lpz, P);
+    else if (P->fwd4d)
+        coo = P->fwd4d (coo, P);
+    else {
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ().xy;
+    }
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().xy;
+
+    if (!P->skip_fwd_finalize)
+        coo = pj_fwd_finalize (P, coo);
+    return coo.xy;
+}
+
+
+
+XYZ pj_fwd3d(LPZ lpz, PJ *P) {
+    PJ_COORD coo = {{0,0,0,0}};
+    coo.lpz = lpz;
+
+    if (!P->skip_fwd_prepare)
+        coo = pj_fwd_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().xyz;
+
+    /* Do the transformation, using the lowest dimensional transformer feasible */
+    if (P->fwd3d)
+        coo.xyz = P->fwd3d(coo.lpz, P);
+    else if (P->fwd4d)
+        coo = P->fwd4d (coo, P);
+    else if (P->fwd)
+        coo.xy = P->fwd (coo.lp, P);
+    else {
+        proj_errno_set (P, EINVAL);
+        return proj_coord_error ().xyz;
+    }
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ().xyz;
+
+    if (!P->skip_fwd_finalize)
+        coo = pj_fwd_finalize (P, coo);
+    return coo.xyz;
+}
+
+
+
+PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
+    if (!P->skip_fwd_prepare)
+        coo = pj_fwd_prepare (P, coo);
+    if (HUGE_VAL==coo.v[0])
+        return proj_coord_error ();
+
+    /* Call the highest dimensional converter available */
+    if (P->fwd4d)
+        coo = P->fwd4d (coo, P);
+    else if (P->fwd3d)
+        coo.xyz  =  P->fwd3d (coo.lpz, P);
+    else if (P->fwd)
+        coo.xy  =  P->fwd (coo.lp, 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 xy;
+    if (!P->skip_fwd_finalize)
+        coo = pj_fwd_finalize (P, coo);
+    return coo;
 }
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
deleted file mode 100644
index d141178f..00000000
--- a/src/pj_fwd3d.c
+++ /dev/null
@@ -1,61 +0,0 @@
-#define PJ_LIB__
-#include <proj.h>
-#include <projects.h>
-#include <errno.h>
-# define EPS 1.0e-12
-
-/* 3D Forward transformation */
-
-XYZ pj_fwd3d(LPZ lpz, PJ *P) {
-    XYZ xyz;
-    XYZ err;
-    double t;
-    int last_errno;
-
-    /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
-    err.x = err.y = err.z = HUGE_VAL;
-
-    if (0==P->fwd3d)
-        return err;
-
-    last_errno = proj_errno_reset(P);
-
-    /* Check validity of angular input coordinates */
-    if (P->left==PJ_IO_UNITS_RADIANS) {
-
-        /* check for forward and latitude or longitude overange */
-        t = fabs(lpz.phi)-M_HALFPI;
-        if (t > EPS || fabs(lpz.lam) > 10.) {
-            pj_ctx_set_errno( P->ctx, -14);
-            return err;
-        }
-
-        /* Clamp latitude to -90..90 degree range */
-        if (fabs(t) <= EPS)
-            lpz.phi = lpz.phi < 0. ? -M_HALFPI : M_HALFPI;
-        else if (P->geoc)   /* Maybe redundant and never used. */
-            lpz.phi = atan(P->rone_es * tan(lpz.phi));
-        lpz.lam -= P->lam0;    /* compute del lp.lam */
-        if (!P->over)
-            lpz.lam = adjlon(lpz.lam); /* adjust del longitude */
-    }
-
-    /* Do the transformation */
-    xyz = (*P->fwd3d)(lpz, P);
-    if ( P->ctx->last_errno )
-        return err;
-
-    /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
-        xyz.x *= P->a;
-        xyz.y *= P->a;
-    }
-
-    /* Handle false eastings/northings and non-metric linear units */
-    xyz.x = P->fr_meter * (xyz.x + P->x0);
-    xyz.y = P->fr_meter * (xyz.y + P->y0);
-    /* z is not scaled since this is handled by vto_meter outside */
-
-    proj_errno_restore (P, last_errno);
-    return xyz;
-}
diff --git a/src/pj_internal.c b/src/pj_internal.c
index 7bfd192b..88fb2b62 100644
--- a/src/pj_internal.c
+++ b/src/pj_internal.c
@@ -8,7 +8,7 @@
  * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2017-07-05
  *
  ******************************************************************************
- * Copyright (c) 2016, 2017, Thomas Knudsen/SDFE
+ * Copyright (c) 2016, 2017, 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"),
@@ -28,16 +28,15 @@
  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
-#define PJ_INTERNAL_C
-#include "proj_internal.h"
-#include "projects.h"
-#include <geodesic.h>
-
 #include <ctype.h>
 #include <stddef.h>
 #include <stdarg.h>
 #include <errno.h>
 
+#include <geodesic.h>
+#include "proj_internal.h"
+#include "projects.h"
+
 
 enum pj_io_units pj_left (PJ *P) {
     enum pj_io_units u = P->inverted? P->right: P->left;
@@ -62,32 +61,59 @@ PJ_COORD proj_coord_error (void) {
 }
 
 
-PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
-    if (0!=P->fwd4d)
-        return P->fwd4d (coo, P);
-    if (0!=P->fwd3d) {
-        coo.xyz  =  pj_fwd3d (coo.lpz, P);
-        return coo;
-    }
-    if (0!=P->fwd) {
-        coo.xy  =  pj_fwd (coo.lp, P);
+
+/**************************************************************************************/
+PJ_COORD pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/***************************************************************************************
+Behave mostly as proj_trans, but attempt to use 2D interfaces only.
+Used in gie.c, to enforce testing 2D code, and by PJ_pipeline.c to implement
+chained calls starting out with a call to its 2D interface.
+***************************************************************************************/
+    if (0==P)
         return coo;
+    if (P->inverted)
+        direction = -direction;
+    switch (direction) {
+        case PJ_FWD:
+            coo.xy = pj_fwd (coo.lp, P);
+            return coo;
+        case PJ_INV:
+            coo.lp = pj_inv (coo.xy, P);
+            return coo;
+        case PJ_IDENT:
+            return coo;
+        default:
+            break;
     }
     proj_errno_set (P, EINVAL);
     return proj_coord_error ();
 }
 
 
-PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
-    if (0!=P->inv4d)
-        return P->inv4d (coo, P);
-    if (0!=P->inv3d) {
-        coo.lpz  =  pj_inv3d (coo.xyz, P);
-        return coo;
-    }
-    if (0!=P->inv) {
-        coo.lp  =  pj_inv (coo.xy, P);
+/**************************************************************************************/
+PJ_COORD pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/***************************************************************************************
+Companion to pj_approx_2D_trans.
+
+Behave mostly as proj_trans, but attempt to use 3D interfaces only.
+Used in gie.c, to enforce testing 3D code, and by PJ_pipeline.c to implement
+chained calls starting out with a call to its 3D interface.
+***************************************************************************************/
+    if (0==P)
         return coo;
+    if (P->inverted)
+        direction = -direction;
+    switch (direction) {
+        case PJ_FWD:
+            coo.xyz = pj_fwd3d (coo.lpz, P);
+            return coo;
+        case PJ_INV:
+            coo.lpz = pj_inv3d (coo.xyz, P);
+            return coo;
+        case PJ_IDENT:
+            return coo;
+        default:
+            break;
     }
     proj_errno_set (P, EINVAL);
     return proj_coord_error ();
@@ -96,7 +122,6 @@ PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
 
 
 
-
 /* Move P to a new context - or to the default context if 0 is specified */
 void proj_context_set (PJ *P, PJ_CONTEXT *ctx) {
     if (0==ctx)
@@ -105,6 +130,7 @@ void proj_context_set (PJ *P, PJ_CONTEXT *ctx) {
     return;
 }
 
+
 void proj_context_inherit (PJ *parent, PJ *child) {
     if (0==parent)
         pj_set_ctx (child, pj_get_default_ctx());
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;
+}
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
deleted file mode 100644
index 53e39a76..00000000
--- a/src/pj_inv3d.c
+++ /dev/null
@@ -1,60 +0,0 @@
-#define PJ_LIB__
-#include <proj.h>
-#include <projects.h>
-#include <errno.h>
-# define EPS 1.0e-12
-
-/* 3D inverse transformation */
-
-LPZ pj_inv3d (XYZ xyz, PJ *P) {
-    LPZ lpz;
-    LPZ err;
-    int last_errno;
-
-    /* cannot const-initialize this due to MSVC's broken (non const) HUGE_VAL */
-    err.lam = err.phi = err.z = HUGE_VAL;
-
-    if (0==P->inv3d)
-        return err;
-
-    /* can't do as much preliminary checking as with forward */
-    if (xyz.x == HUGE_VAL || xyz.y == HUGE_VAL || xyz.z == HUGE_VAL ) {
-        pj_ctx_set_errno( P->ctx, -15);
-        return err;
-    }
-
-    last_errno = proj_errno_reset (P);
-
-    /* de-scale and de-offset */
-    /* z is not de-scaled since that is handled by vto_meter before we get here */
-    xyz.x = (xyz.x * P->to_meter - P->x0);
-    xyz.y = (xyz.y * P->to_meter - P->y0);
-    /* 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 */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
-        xyz.x *= P->ra;
-        xyz.y *= P->ra;
-    }
-
-    /* Do inverse transformation */
-    lpz = (*P->inv3d) (xyz, P);
-    if (P->ctx->last_errno)
-        return err;
-
-    if (P->left==PJ_IO_UNITS_RADIANS) {
-        /* reduce from del lp.lam */
-        lpz.lam += P->lam0;
-
-         /* adjust longitude to central meridian */
-        if (!P->over)
-            lpz.lam = adjlon(lpz.lam);
-
-        /* This may be redundant and never used */
-        if (P->geoc && fabs(fabs(lpz.phi)-M_HALFPI) > EPS)
-            lpz.phi = atan(P->one_es * tan(lpz.phi));
-    }
-
-    proj_errno_restore (P, last_errno);
-    return lpz;
-}
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 21861331..c72df2e7 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -117,15 +117,6 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
     }
 
 /* -------------------------------------------------------------------- */
-/*      Transform Z to meters if it isn't already.                      */
-/* -------------------------------------------------------------------- */
-    if( srcdefn->vto_meter != 1.0 && z != NULL )
-    {
-        for( i = 0; i < point_count; i++ )
-            z[point_offset*i] *= srcdefn->vto_meter;
-    }
-
-/* -------------------------------------------------------------------- */
 /*      Transform geocentric source coordinates to lat/long.            */
 /* -------------------------------------------------------------------- */
     if( srcdefn->is_geocent )
@@ -145,6 +136,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
                 {
                     x[point_offset*i] *= srcdefn->to_meter;
                     y[point_offset*i] *= srcdefn->to_meter;
+                    z[point_offset*i] *= srcdefn->to_meter;
                 }
             }
         }
@@ -253,11 +245,12 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
             }
         }
     }
+
 /* -------------------------------------------------------------------- */
 /*      But if they are already lat long, adjust for the prime          */
 /*      meridian if there is one in effect.                             */
 /* -------------------------------------------------------------------- */
-    if( srcdefn->from_greenwich != 0.0 )
+    if ((srcdefn->is_geocent || srcdefn->is_latlong) && ( srcdefn->from_greenwich != 0.0 ))
     {
         for( i = 0; i < point_count; i++ )
         {
@@ -308,15 +301,14 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
 /*      But if they are staying lat long, adjust for the prime          */
 /*      meridian if there is one in effect.                             */
 /* -------------------------------------------------------------------- */
-    if( dstdefn->from_greenwich != 0.0 )
+    if ((dstdefn->is_geocent || dstdefn->is_latlong) && ( dstdefn->from_greenwich != 0.0 ))
     {
         for( i = 0; i < point_count; i++ )
         {
             if( x[point_offset*i] != HUGE_VAL )
                 x[point_offset*i] -= dstdefn->from_greenwich;
         }
-    }
-
+}
 
 /* -------------------------------------------------------------------- */
 /*      Transform destination latlong to geocentric if required.        */
@@ -340,6 +332,7 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
                 {
                     x[point_offset*i] *= dstdefn->fr_meter;
                     y[point_offset*i] *= dstdefn->fr_meter;
+                    z[point_offset*i] *= srcdefn->fr_meter;
                 }
             }
         }
@@ -454,15 +447,6 @@ int pj_transform( PJ *srcdefn, PJ *dstdefn, long point_count, int point_offset,
     }
 
 /* -------------------------------------------------------------------- */
-/*      Transform Z from meters if needed.                              */
-/* -------------------------------------------------------------------- */
-    if( dstdefn->vto_meter != 1.0 && z != NULL )
-    {
-        for( i = 0; i < point_count; i++ )
-            z[point_offset*i] *= dstdefn->vfr_meter;
-    }
-
-/* -------------------------------------------------------------------- */
 /*      Transform normalized axes into unusual output coordinate axis   */
 /*      orientation if needed.                                          */
 /* -------------------------------------------------------------------- */
diff --git a/src/proj.def b/src/proj.def
index 51b79558..f758477d 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -151,3 +151,6 @@ EXPORTS
     pj_calc_ellipsoid_params       @135
     pj_chomp                       @136
     pj_shrink                      @137
+
+    pj_approx_2D_trans             @138
+    pj_approx_3D_trans             @139
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index a0550727..4a416084 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -132,8 +132,15 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
 
 
 
-/* Apply the transformation P to the coordinate coo */
+/**************************************************************************************/
 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/***************************************************************************************
+Apply the transformation P to the coordinate coo, preferring the 4D interfaces if
+available.
+
+See also pj_approx_2D_trans and pj_approx_3D_trans in pj_internal.c, which work
+similarly, but prefers the 2D resp. 3D interfaces if available.
+***************************************************************************************/
     if (0==P)
         return coo;
     if (P->inverted)
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 13a1cac6..3dd04e62 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -57,6 +57,9 @@ extern "C" {
 #define PJ_TORAD(deg)  ((deg)*M_PI/180.0)
 #endif
 
+/* Maximum latitudinal overshoot accepted */
+#define PJ_EPS_LAT 1e-12
+
 
 /* This enum is also conditionally defined in projects.h - but we need it here */
 /* for the pj_left/right prototypes, and enums cannot be forward declared      */
@@ -79,6 +82,14 @@ void proj_context_inherit (PJ *parent, PJ *child);
 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P);
 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
 
+PJ_COORD pj_fwd_prepare  (PJ *P, PJ_COORD coo);
+PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo);
+PJ_COORD pj_inv_prepare  (PJ *P, PJ_COORD coo);
+PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo);
+PJ_COORD pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
+PJ_COORD pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
+
+
 /* Grid functionality */
 int             proj_vgrid_init(PJ *P, const char *grids);
 int             proj_hgrid_init(PJ *P, const char *grids);
diff --git a/src/projects.h b/src/projects.h
index 63ea44da..209d6f39 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -337,7 +337,12 @@ struct PJconsts {
     int  geoc;                      /* Geocentric latitude flag */
     int  is_latlong;                /* proj=latlong ... not really a projection at all */
     int  is_geocent;                /* proj=geocent ... not really a projection at all */
+    int  is_pipeline;               /* 1 if PJ represents a pipeline */
     int  need_ellps;                /* 0 for operations that are purely cartesian */
+    int  skip_fwd_prepare;
+    int  skip_fwd_finalize;
+    int  skip_inv_prepare;
+    int  skip_inv_finalize;
 
     enum pj_io_units left;          /* Flags for input/output coordinate types */
     enum pj_io_units right;
@@ -349,7 +354,7 @@ struct PJconsts {
 
     **************************************************************************************/
 
-    double  lam0, phi0;                /* central longitude, latitude */
+    double  lam0, phi0;                /* central meridian, parallel */
     double  x0, y0, z0, t0;            /* false easting and northing (and height and time) */
 
 
@@ -590,7 +595,7 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
 
 
 #ifdef PJ_LIB__
-#define PROJ_HEAD(id, name) static const char des_##id [] = name
+#define PROJ_HEAD(name, desc) static const char des_##name [] = desc
 
 #define OPERATION(name, NEED_ELLPS)                          \
                                                              \
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index a9acf4f9..3877b58a 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -509,7 +509,9 @@ Central Conic
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=pipeline +step +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +a=6390000 +x_0=330000 +y_0=-350000 +step +proj=axisswap +order=1,-2,3,4
+operation +proj=pipeline +R=6390000
+    +step +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +x_0=330000 +y_0=-350000
+	+step +proj=axisswap +order=1,-2
 -------------------------------------------------------------------------------
 tolerance 0.00010 mm
 accept 24 55
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 2d558f88..099934a7 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -18,13 +18,13 @@ operation +proj=ob_tran  +o_proj=moll  +R=6378137.0  +o_lon_p=0  +o_lat_p=0  +lo
 tolerance 1 mm
 direction inverse
 
-accept    300000 400000 0 0
-expect    -42.7562158333  85.5911341667  0  0
+accept    300000 400000
+expect    -42.7562158333  85.5911341667
 
 direction forward
 
-accept    10 20 0 0
-expect    -1384841.18787 7581707.88240 0 0
+accept    10 20
+expect    -1384841.18787 7581707.88240
 -------------------------------------------------------------------------------