[WIP] Pipeline coordinate stack

OSGeo · Feb 6, 2019 · b56c1f0 · b56c1f0
1 parent ffbaaa8
commit b56c1f0
Show file tree

Hide file tree

Showing 2 changed files with 110 additions and 10 deletions.
diff --git a/src/pipeline.cpp b/src/pipeline.cpp
@@ -100,6 +100,7 @@ Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
 #include <math.h>
 #include <stddef.h>
 #include <string.h>
+#include <stack>
 
 #include "geodesic.h"
 #include "proj.h"
@@ -115,6 +116,8 @@ struct pj_opaque {
     char **argv;
     char **current_argv;
     PJ **pipeline;
+    bool *push[4];
+    bool *pop[4];
 };
 } // anonymous namespace
 
@@ -132,26 +135,59 @@ static PJ_LP     pipeline_reverse (PJ_XY xy, PJ *P);
 
 static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) {
     int i, first_step, last_step;
+    pj_opaque *opaque = static_cast<struct pj_opaque*>(P->opaque);
+    std::stack<double> stack[4];
 
     first_step = 1;
-    last_step  = static_cast<struct pj_opaque*>(P->opaque)->steps + 1;
+    last_step  = opaque->steps + 1;
 
-    for (i = first_step;  i != last_step;  i++)
-        point = proj_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_FWD, point);
+    for (i = first_step;  i != last_step;  i++) {
 
+        for (int j=0; j<4; j++) {
+            if (opaque->push[j][i-1]) {
+                stack[j].push(point.v[j]);
+            }
+        }
+
+        point = proj_trans (opaque->pipeline[i], PJ_FWD, point);
+
+        for (int j=0; j<4; j++) {
+            if (opaque->pop[j][i-1] && stack[j].size() > 0)  {
+                point.v[j] = stack[j].top();
+                stack[j].pop();
+            }
+        }
+    }
     return point;
 }
 
 
 static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
     int i, first_step, last_step;
+    pj_opaque *opaque = static_cast<struct pj_opaque*>(P->opaque);
+    std::stack<double> stack[4];
 
     first_step = static_cast<struct pj_opaque*>(P->opaque)->steps;
     last_step  =  0;
 
-    for (i = first_step;  i != last_step;  i--)
+    for (i = first_step;  i != last_step;  i--) {
+        /* pop first, which in reverse is push */
+        for (int j=0; j<4; j++) {
+            if (opaque->pop[j][i-1]) {
+                stack[j].push(point.v[j]);
+            }
+        }
+
         point = proj_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_INV, point);
 
+        /* then push, which in reverse is pop */
+        for (int j=0; j<4; j++) {
+            if (opaque->push[j][i-1] && stack[j].size() > 0)  {
+                point.v[j] = stack[j].top();
+                stack[j].pop();
+            }
+        }
+    }
     return point;
 }
 
@@ -231,13 +267,22 @@ static PJ *destructor (PJ *P, int errlev) {
 
 
 static PJ *pj_create_pipeline (PJ *P, size_t steps) {
-
+    pj_opaque *opaque = static_cast<struct pj_opaque*>(P->opaque);
     /* Room for the pipeline: An array of PJ * with room for sentinels at both ends */
-    static_cast<struct pj_opaque*>(P->opaque)->pipeline = static_cast<PJ**>(pj_calloc (steps + 2, sizeof(PJ *)));
-    if (nullptr==static_cast<struct pj_opaque*>(P->opaque)->pipeline)
+    opaque->pipeline = static_cast<PJ**>(pj_calloc (steps + 2, sizeof(PJ *)));
+    if (nullptr==opaque->pipeline)
         return nullptr;
 
-    static_cast<struct pj_opaque*>(P->opaque)->steps = (int)steps;
+    opaque->steps = (int)steps;
+
+    for (int i=0; i<4; i++) {
+        opaque->push[i] = static_cast<bool*>(pj_calloc(steps, sizeof(bool)));
+        opaque->pop[i]  = static_cast<bool*>(pj_calloc(steps, sizeof(bool)));
+        for (int j=0; j<(int)steps; j++) {
+            opaque->push[i][j] = false;
+            opaque->pop[i][j] = false;
+        }
+    }
 
     return P;
 }
@@ -427,8 +472,10 @@ PJ *OPERATION(pipeline,0) {
     if (nullptr==pj_create_pipeline (P, nsteps))
         return destructor (P, ENOMEM);
 
+    pj_opaque *opaque = static_cast<struct pj_opaque*>(P->opaque);
     set_ellipsoid(P);
 
+
     /* Now loop over all steps, building a new set of arguments for each init */
     i_current_step = i_first_step;
     for (i = 0;  i < nsteps;  i++) {
@@ -470,12 +517,23 @@ PJ *OPERATION(pipeline,0) {
 
         proj_errno_restore (P, err);
 
-        /* Is this step inverted? */
-        for (j = 0;  j < current_argc; j++)
+        /* take care of pipeline specific parameters */
+        for (j = 0;  j < current_argc; j++) {
+            /* Is this step inverted? */
             if (0==strcmp("inv", current_argv[j])) {
                 /* if +inv exists in both global and local args the forward operation should be used */
                 next_step->inverted = next_step->inverted == 0 ? 1 : 0;
             }
+            /* do we need to push or pop something on the coordinate stacks? */
+            if (strcmp("push=1", current_argv[j])==0) opaque->push[0][i] = true;
+            if (strcmp("push=2", current_argv[j])==0) opaque->push[1][i] = true;
+            if (strcmp("push=3", current_argv[j])==0) opaque->push[2][i] = true;
+            if (strcmp("push=4", current_argv[j])==0) opaque->push[3][i] = true;
+            if (strcmp("pop=1",  current_argv[j])==0) opaque->pop[0][i]  = true;
+            if (strcmp("pop=2",  current_argv[j])==0) opaque->pop[1][i]  = true;
+            if (strcmp("pop=3",  current_argv[j])==0) opaque->pop[2][i]  = true;
+            if (strcmp("pop=4",  current_argv[j])==0) opaque->pop[3][i]  = true;
+        }
 
         static_cast<struct pj_opaque*>(P->opaque)->pipeline[i+1] = next_step;
 

diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
@@ -294,6 +294,48 @@ accept      12          56
 expect      1335.8339   7522.963
 -------------------------------------------------------------------------------
 
+-------------------------------------------------------------------------------
+Test Pipeline Coordinate Stack
+-------------------------------------------------------------------------------
+operation  +proj=pipeline
+            +step +proj=utm +zone=32 +push=1
+            +step +proj=utm +zone=33 +pop=1 +inv
+
+accept      12  56  0   2020
+expect      12  56  0   2020
+roundtrip   10
+
+operation  +proj=pipeline
+            +step +proj=latlon # dummy step
+            +step +proj=utm +zone=32 +push=1
+            +step +proj=utm +zone=33 +pop=1 +inv
+            +step +proj=affine # dummy step
+
+accept      12  56  0   2020
+expect      12  56  0   2020
+roundtrip   10
+
+operation   +proj=pipeline
+            +step +inv +proj=eqearth ++push=2
+            +step +proj=laea +pop=2
+
+accept      900000      6000000 0   2020
+expect      896633.0226 6000000 0   2020
+
+# Datum shift in cartesian space but keeping the height
+# (simulates a datum-shift with affin since ISO19111 code
+# currrently obfuscates proj-strings using cart/helmert/invcart)
+operation   +proj=pipeline +ellps=GRS80
+            +step +proj=cart +push=3
+            +step +proj=affine +xoff=1000 +yoff=2000 +xoff=3000
+            +step +proj=cart +pop=3 +inv
+
+accept      12              56              0 2020
+expect      12.0280112877   55.9896187413   0 2020
+roundtrip   1
+
+
+-------------------------------------------------------------------------------
 
 -------------------------------------------------------------------------------
 Test bugfix of https://github.com/OSGeo/proj.4/issues/1002