aboutsummaryrefslogtreecommitdiff
path: root/src/pipeline.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/pipeline.cpp')
-rw-r--r--src/pipeline.cpp503
1 files changed, 503 insertions, 0 deletions
diff --git a/src/pipeline.cpp b/src/pipeline.cpp
new file mode 100644
index 00000000..76fc58a5
--- /dev/null
+++ b/src/pipeline.cpp
@@ -0,0 +1,503 @@
+/*******************************************************************************
+
+ Transformation pipeline manager
+
+ Thomas Knudsen, 2016-05-20/2016-11-20
+
+********************************************************************************
+
+ Geodetic transformations are typically organized in a number of
+ steps. For example, a datum shift could be carried out through
+ these steps:
+
+ 1. Convert (latitude, longitude, ellipsoidal height) to
+ 3D geocentric cartesian coordinates (X, Y, Z)
+ 2. Transform the (X, Y, Z) coordinates to the new datum, using a
+ 7 parameter Helmert transformation.
+ 3. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height)
+
+ If the height system used is orthometric, rather than ellipsoidal,
+ another step is needed at each end of the process:
+
+ 1. Add the local geoid undulation (N) to the orthometric height
+ to obtain the ellipsoidal (i.e. geometric) height.
+ 2. Convert (latitude, longitude, ellipsoidal height) to
+ 3D geocentric cartesian coordinates (X, Y, Z)
+ 3. Transform the (X, Y, Z) coordinates to the new datum, using a
+ 7 parameter Helmert transformation.
+ 4. Convert (X, Y, Z) back to (latitude, longitude, ellipsoidal height)
+ 5. Subtract the local geoid undulation (N) from the ellipsoidal height
+ to obtain the orthometric height.
+
+ Additional steps can be added for e.g. change of vertical datum, so the
+ list can grow fairly long. None of the steps are, however, particularly
+ complex, and data flow is strictly from top to bottom.
+
+ Hence, in principle, the first example above could be implemented using
+ Unix pipelines:
+
+ cat my_coordinates | geographic_to_xyz | helmert | xyz_to_geographic > my_transformed_coordinates
+
+ in the grand tradition of Software Tools [1].
+
+ The proj pipeline driver implements a similar concept: Stringing together
+ a number of steps, feeding the output of one step to the input of the next.
+
+ It is a very powerful concept, that increases the range of relevance of the
+ proj.4 system substantially. It is, however, not a particularly intrusive
+ addition to the PROJ.4 code base: The implementation is by and large completed
+ by adding an extra projection called "pipeline" (i.e. this file), which
+ handles all business, and a small amount of added functionality in the
+ pj_init code, implementing support for multilevel, embedded pipelines.
+
+ Syntactically, the pipeline system introduces the "+step" keyword (which
+ indicates the start of each transformation step), and reintroduces the +inv
+ keyword (indicating that a given transformation step should run in reverse, i.e.
+ forward, when the pipeline is executed in inverse direction, and vice versa).
+
+ Hence, the first transformation example above, can be implemented as:
+
+ +proj=pipeline +step proj=cart +step proj=helmert <ARGS> +step proj=cart +inv
+
+ Where <ARGS> indicate the Helmert arguments: 3 translations (+x=..., +y=...,
+ +z=...), 3 rotations (+rx=..., +ry=..., +rz=...) and a scale factor (+s=...).
+ Following geodetic conventions, the rotations are given in arcseconds,
+ and the scale factor is given as parts-per-million.
+
+ [1] B. W. Kernighan & P. J. Plauger: Software tools.
+ Reading, Massachusetts, Addison-Wesley, 1976, 338 pp.
+
+********************************************************************************
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-05-20
+
+********************************************************************************
+* 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"),
+* 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.
+*
+********************************************************************************/
+
+#define PJ_LIB__
+
+#include <errno.h>
+#include <math.h>
+#include <stddef.h>
+#include <string.h>
+
+#include "geodesic.h"
+#include "proj.h"
+#include "proj_internal.h"
+#include "projects.h"
+
+PROJ_HEAD(pipeline, "Transformation pipeline manager");
+
+/* Projection specific elements for the PJ object */
+namespace { // anonymous namespace
+struct pj_opaque {
+ int steps;
+ char **argv;
+ char **current_argv;
+ PJ **pipeline;
+};
+} // anonymous namespace
+
+
+
+static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P);
+static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P);
+static XYZ pipeline_forward_3d (LPZ lpz, PJ *P);
+static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P);
+static XY pipeline_forward (LP lp, PJ *P);
+static LP pipeline_reverse (XY xy, PJ *P);
+
+
+
+
+static PJ_COORD pipeline_forward_4d (PJ_COORD point, PJ *P) {
+ int i, first_step, last_step;
+
+ first_step = 1;
+ last_step = static_cast<struct pj_opaque*>(P->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);
+
+ return point;
+}
+
+
+static PJ_COORD pipeline_reverse_4d (PJ_COORD point, PJ *P) {
+ int i, first_step, last_step;
+
+ first_step = static_cast<struct pj_opaque*>(P->opaque)->steps;
+ last_step = 0;
+
+ for (i = first_step; i != last_step; i--)
+ point = proj_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_INV, point);
+
+ return point;
+}
+
+
+
+
+static XYZ pipeline_forward_3d (LPZ lpz, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ int i;
+ point.lpz = lpz;
+
+ for (i = 1; i <= static_cast<struct pj_opaque*>(P->opaque)->steps; i++)
+ point = pj_approx_3D_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_FWD, point);
+
+ return point.xyz;
+}
+
+
+static LPZ pipeline_reverse_3d (XYZ xyz, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ int i;
+ point.xyz = xyz;
+
+ for (i = static_cast<struct pj_opaque*>(P->opaque)->steps; i > 0 ; i--)
+ point = pj_approx_3D_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_INV, point);
+
+ return point.lpz;
+}
+
+
+
+
+static XY pipeline_forward (LP lp, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ int i;
+ point.lp = lp;
+
+ for (i = 1; i <= static_cast<struct pj_opaque*>(P->opaque)->steps; i++)
+ point = pj_approx_2D_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_FWD, point);
+
+ return point.xy;
+}
+
+
+static LP pipeline_reverse (XY xy, PJ *P) {
+ PJ_COORD point = {{0,0,0,0}};
+ int i;
+ point.xy = xy;
+ for (i = static_cast<struct pj_opaque*>(P->opaque)->steps; i > 0 ; i--)
+ point = pj_approx_2D_trans (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i], PJ_INV, point);
+
+ return point.lp;
+}
+
+
+
+
+static PJ *destructor (PJ *P, int errlev) {
+ int i;
+ if (nullptr==P)
+ return nullptr;
+
+ if (nullptr==P->opaque)
+ return pj_default_destructor (P, errlev);
+
+ /* Deallocate each pipeine step, then pipeline array */
+ if (nullptr!=static_cast<struct pj_opaque*>(P->opaque)->pipeline)
+ for (i = 0; i < static_cast<struct pj_opaque*>(P->opaque)->steps; i++)
+ proj_destroy (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i+1]);
+ pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->pipeline);
+
+ pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->argv);
+ pj_dealloc (static_cast<struct pj_opaque*>(P->opaque)->current_argv);
+
+ return pj_default_destructor(P, errlev);
+}
+
+
+static PJ *pj_create_pipeline (PJ *P, size_t steps) {
+
+ /* 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)
+ return nullptr;
+
+ static_cast<struct pj_opaque*>(P->opaque)->steps = (int)steps;
+
+ return P;
+}
+
+
+
+
+/* count the number of args in pipeline definition, and mark all args as used */
+static size_t argc_params (paralist *params) {
+ size_t argc = 0;
+ for (; params != nullptr; params = params->next) {
+ argc++;
+ params->used = 1;
+ }
+ return ++argc; /* one extra for the sentinel */
+}
+
+/* Sentinel for argument list */
+static const char *argv_sentinel = "step";
+
+/* turn paralist into argc/argv style argument list */
+static char **argv_params (paralist *params, size_t argc) {
+ char **argv;
+ size_t i = 0;
+ argv = static_cast<char**>(pj_calloc (argc, sizeof (char *)));
+ if (nullptr==argv)
+ return nullptr;
+ for (; params != nullptr; params = params->next)
+ argv[i++] = params->param;
+ argv[i++] = const_cast<char*>(argv_sentinel);
+ 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. */
+/* Unfortunately we can't use the normal ellipsoid setter in pj_init, since */
+/* it adds a +ellps parameter to the global args if nothing else is specified*/
+/* This is problematic since that ellipsoid spec is then passed on to the */
+/* pipeline children. This is rarely what we want, so here we implement our */
+/* own logic instead. If an ellipsoid is set in the global args, it is used */
+/* as the pipeline ellipsoid. Otherwise we use WGS84 parameters as default. */
+/* At last we calculate the rest of the ellipsoid parameters and */
+/* re-initialize P->geod. */
+static void set_ellipsoid(PJ *P) {
+ paralist *cur, *attachment;
+ int err = proj_errno_reset (P);
+
+ /* Break the linked list after the global args */
+ attachment = nullptr;
+ for (cur = P->params; cur != nullptr; cur = cur->next)
+ /* cur->next will always be non 0 given argv_sentinel presence, */
+ /* but this is far from being obvious for a static analyzer */
+ if (cur->next != nullptr && strcmp(argv_sentinel, cur->next->param) == 0) {
+ attachment = cur->next;
+ cur->next = nullptr;
+ break;
+ }
+
+ /* Check if there's any ellipsoid specification in the global params. */
+ /* If not, use WGS84 as default */
+ if (0 != pj_ellipsoid (P)) {
+ P->a = 6378137.0;
+ P->es = .00669438002290341575;
+
+ /* reset an "unerror": In this special use case, the errno is */
+ /* not an error signal, but just a reply from pj_ellipsoid, */
+ /* telling us that "No - there was no ellipsoid definition in */
+ /* the PJ you provided". */
+ proj_errno_reset (P);
+ }
+ P->a_orig = P->a;
+ P->es_orig = P->es;
+
+ pj_calc_ellipsoid_params (P, P->a, P->es);
+
+ geod_init(P->geod, P->a, (1 - sqrt (1 - P->es)));
+
+ /* Re-attach the dangling list */
+ /* Note: cur will always be non 0 given argv_sentinel presence, */
+ /* but this is far from being obvious for a static analyzer */
+ if( cur != nullptr )
+ cur->next = attachment;
+ proj_errno_restore (P, err);
+}
+
+
+
+
+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->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 = static_cast<struct pj_opaque*>(pj_calloc (1, sizeof(struct pj_opaque)));
+ if (nullptr==P->opaque)
+ return destructor(P, ENOMEM);
+
+ argc = (int)argc_params (P->params);
+ static_cast<struct pj_opaque*>(P->opaque)->argv = argv = argv_params (P->params, argc);
+ if (nullptr==argv)
+ return destructor (P, ENOMEM);
+
+ static_cast<struct pj_opaque*>(P->opaque)->current_argv = current_argv = static_cast<char**>(pj_calloc (argc, sizeof (char *)));
+ if (nullptr==current_argv)
+ return destructor (P, ENOMEM);
+
+ /* Do some syntactical sanity checking */
+ for (i = 0; i < argc; i++) {
+ if (0==strcmp (argv_sentinel, argv[i])) {
+ if (-1==i_pipeline) {
+ proj_log_error (P, "Pipeline: +step before +proj=pipeline");
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE);
+ }
+ if (0==nsteps)
+ i_first_step = i;
+ nsteps++;
+ continue;
+ }
+
+ 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");
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: nested pipelines */
+ }
+ i_pipeline = i;
+ }
+ }
+ nsteps--; /* Last instance of +step is just a sentinel */
+ static_cast<struct pj_opaque*>(P->opaque)->steps = nsteps;
+
+ if (-1==i_pipeline)
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: no pipeline def */
+
+ if (0==nsteps)
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: no pipeline def */
+
+ /* Make room for the pipeline and execution indicators */
+ if (nullptr==pj_create_pipeline (P, nsteps))
+ return destructor (P, ENOMEM);
+
+ 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++) {
+ int j;
+ int current_argc = 0;
+ int err;
+ PJ *next_step = nullptr;
+
+ /* Build a set of setup args for the current step */
+ proj_log_trace (P, "Pipeline: Building arg list for step no. %d", i);
+
+ /* First add the step specific args */
+ for (j = i_current_step + 1; 0 != strcmp ("step", argv[j]); j++)
+ current_argv[current_argc++] = argv[j];
+
+ i_current_step = j;
+
+ /* Then add the global args */
+ for (j = i_pipeline + 1; 0 != strcmp ("step", argv[j]); j++)
+ current_argv[current_argc++] = argv[j];
+
+ proj_log_trace (P, "Pipeline: init - %s, %d", current_argv[0], current_argc);
+ for (j = 1; j < current_argc; j++)
+ proj_log_trace (P, " %s", current_argv[j]);
+
+ err = proj_errno_reset (P);
+
+ next_step = proj_create_argv (P->ctx, current_argc, current_argv);
+ proj_log_trace (P, "Pipeline: Step %d (%s) at %p", i, current_argv[0], next_step);
+
+ if (nullptr==next_step) {
+ /* The step init failed, but possibly without setting errno. If so, we say "malformed" */
+ int err_to_report = proj_errno(P);
+ if (0==err_to_report)
+ err_to_report = PJD_ERR_MALFORMED_PIPELINE;
+ proj_log_error (P, "Pipeline: Bad step definition: %s (%s)", current_argv[0], pj_strerrno (err_to_report));
+ return destructor (P, err_to_report); /* ERROR: bad pipeline def */
+ }
+
+ proj_errno_restore (P, err);
+
+ /* Is this step inverted? */
+ for (j = 0; j < current_argc; j++)
+ 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;
+ }
+
+ static_cast<struct pj_opaque*>(P->opaque)->pipeline[i+1] = next_step;
+
+ proj_log_trace (P, "Pipeline at [%p]: step at [%p] (%s) done", P, next_step, current_argv[0]);
+ }
+
+ /* Require a forward path through the pipeline */
+ for (i = 1; i <= nsteps; i++) {
+ PJ *Q = static_cast<struct pj_opaque*>(P->opaque)->pipeline[i];
+ if ( ( Q->inverted && (Q->inv || Q->inv3d || Q->fwd4d) ) ||
+ (!Q->inverted && (Q->fwd || Q->fwd3d || Q->fwd4d) ) ) {
+ continue;
+ } else {
+ proj_log_error (P, "Pipeline: A forward operation couldn't be constructed");
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE);
+ }
+ }
+
+ /* determine if an inverse operation is possible */
+ for (i = 1; i <= nsteps; i++) {
+ PJ *Q = static_cast<struct pj_opaque*>(P->opaque)->pipeline[i];
+ if ( pj_has_inverse(Q) ) {
+ continue;
+ } else {
+ P->inv = nullptr;
+ P->inv3d = nullptr;
+ P->inv4d = nullptr;
+ break;
+ }
+ }
+
+ /* Check that output units from step i are compatible with expected units in step i+1 */
+ for (i = 1; i < nsteps; i++) {
+ enum pj_io_units unit_returned = pj_right (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i]);
+ enum pj_io_units unit_expected = pj_left (static_cast<struct pj_opaque*>(P->opaque)->pipeline[i+1]);
+
+ if ( unit_returned == PJ_IO_UNITS_WHATEVER || unit_expected == PJ_IO_UNITS_WHATEVER )
+ continue;
+ if ( unit_returned != unit_expected ) {
+ proj_log_error (P, "Pipeline: Mismatched units between step %d and %d", i, i+1);
+ return destructor (P, PJD_ERR_MALFORMED_PIPELINE);
+ }
+ }
+
+ proj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps);
+
+ /* Determine forward input (= reverse output) data type */
+ P->left = pj_left (static_cast<struct pj_opaque*>(P->opaque)->pipeline[1]);
+
+ /* Now, correspondingly determine forward output (= reverse input) data type */
+ P->right = pj_right (static_cast<struct pj_opaque*>(P->opaque)->pipeline[nsteps]);
+ return P;
+}
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-15 21:44:38 +0000