Pipeline plus api - in continuation of the Genereic Coordinates pull request (#445)

This commit reflects continued work on the new rationalized API with
geodetic extensions  (see rationale in proj.h). It also reflects the
parallel work on the transformation pipeline reimplementation, by
introducing the PJ_cart cartesian/ellipsoidal converter.

See example/pj_obs_api_test.c for demo example code 

1 files changed, 248 insertions, 0 deletions

diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
new file mode 100644
index 00000000..ccec8a64
--- /dev/null
+++ b/src/pj_obs_api.c
@@ -0,0 +1,248 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Implement a (currently minimalistic) proj API based primarily
+ *           on the PJ_OBS generic geodetic data type.
+ *
+ *           proj thread contexts have not seen widespread use, so one of the
+ *           intentions with this new API is to make them less visible on the
+ *           API surface: Contexts do not have a life by themselves, they are
+ *           visible only through their associated PJs, and the number of
+ *           functions supporting them is limited.
+ *
+ * Author:   Thomas Knudsen,  thokn@sdfe.dk,  2016-06-09/2016-11-06
+ *
+ ******************************************************************************
+ * Copyright (c) 2016, 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_OBS_C
+#include <proj.h>
+#include <projects.h>
+#include <float.h>
+#include <math.h>
+
+
+/* Used as return value in case of errors */
+const PJ_OBS pj_obs_error = {
+    /* Cannot use HUGE_VAL here: MSVC misimplements HUGE_VAL as something that is not compile time constant */
+    {{DBL_MAX,DBL_MAX,DBL_MAX,DBL_MAX}},
+    {{DBL_MAX,DBL_MAX,DBL_MAX}},
+    0, 0
+};
+
+/* Used for zero-initializing new objects */
+const PJ_OBS pj_obs_null = {
+    {{0, 0, 0, 0}},
+    {{0, 0, 0}},
+    0, 0
+};
+
+/* Magic object signaling proj system shutdown mode to routines taking a PJ * arg */
+const PJ *pj_shutdown = (PJ *) &pj_shutdown;
+
+/* Euclidean distance between two 2D coordinates stored in PJ_OBSs */
+double pj_obs_dist_2d (PJ_OBS a, PJ_OBS b) {
+    double *A = a.coo.v, *B = b.coo.v;
+    return hypot (A[0] - B[0],  A[1] - B[1]);
+}
+
+/* Euclidean distance between two 3D coordinates stored in PJ_OBSs */
+double pj_obs_dist_3d (PJ_OBS a, PJ_OBS b) {
+    double *A = a.coo.v, *B = b.coo.v;
+    return hypot (hypot (A[0] - B[0],  A[1] - B[1]),  A[2] - B[2]);
+}
+
+
+PJ_OBS pj_fwdobs (PJ_OBS obs, PJ *P) {
+    if (0!=P->fwd3d) {
+        obs.coo.xyz  =  pj_fwd3d (obs.coo.lpz, P);
+        return obs;
+    }
+    if (0!=P->fwd) {
+        obs.coo.xy  =  pj_fwd (obs.coo.lp, P);
+        return obs;
+    }
+    pj_error_set (P, EINVAL);
+    return pj_obs_error;
+}
+
+
+PJ_OBS pj_invobs (PJ_OBS obs, PJ *P) {
+    if (0!=P->inv3d) {
+        obs.coo.lpz  =  pj_inv3d (obs.coo.xyz, P);
+        return obs;
+    }
+    if (0!=P->inv) {
+        obs.coo.lp  =  pj_inv (obs.coo.xy, P);
+        return obs;
+    }
+    pj_error_set (P, EINVAL);
+    return pj_obs_error;
+}
+
+
+/* Apply the transformation P to the observation obs */
+PJ_OBS pj_apply (PJ *P, enum pj_direction direction, PJ_OBS obs) {
+    if (0==P)
+        return obs;
+
+    switch (direction) {
+        case PJ_FWD:
+            return pj_fwdobs (obs, P);
+        case PJ_INV:
+            return  pj_invobs (obs, P);
+        case PJ_IDENT:
+            return obs;
+        default:
+            break;
+    }
+
+    pj_error_set (P, EINVAL);
+    return pj_obs_error;
+}
+
+
+/* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
+double pj_roundtrip (PJ *P, enum pj_direction direction, int n, PJ_OBS obs) {
+    int i;
+    PJ_OBS o, u;
+
+    if (0==P)
+        return HUGE_VAL;
+
+    if (n < 1) {
+        pj_error_set (P, EINVAL);
+        return HUGE_VAL;
+    }
+
+    o.coo = obs.coo;
+
+    for (i = 0;  i < n;  i++) {
+        switch (direction) {
+            case PJ_FWD:
+                u.coo.xyz  =  pj_fwd3d (o.coo.lpz, P);
+                o.coo.lpz  =  pj_inv3d (u.coo.xyz, P);
+                break;
+            case PJ_INV:
+                u.coo.lpz  =  pj_inv3d (o.coo.xyz, P);
+                o.coo.xyz  =  pj_fwd3d (u.coo.lpz, P);
+                break;
+            default:
+                pj_error_set (P, EINVAL);
+                return HUGE_VAL;
+        }
+    }
+
+    return pj_obs_dist_3d (o, obs);
+}
+
+
+PJ *pj_create (const char *definition) {
+    return pj_init_plus (definition);
+}
+
+
+PJ *pj_create_argv (int argc, char **argv) {
+    return pj_init (argc, argv);
+}
+
+
+/* From here: Minimum viable support for contexts. The first four functions   */
+/* relate to error reporting, debugging, and logging, hence being generically */
+/* useful. The remaining is a compact implementation of the more low level    */
+/* proj_api.h thread contexts, which may or may not be useful  */
+
+int pj_error (PJ *P) {
+    return pj_ctx_get_errno (pj_get_ctx(P));
+}
+
+
+void pj_error_set (PJ *P, int err) {
+    pj_ctx_set_errno (pj_get_ctx(P), err);
+}
+
+
+/* Set debug level 0-3. Higher number means more debug info. 0 turns it off */
+void pj_debug_set (PJ *P, enum pj_debug_level debuglevel) {
+    PJ_CONTEXT *ctx;
+    if (0==P)
+        ctx = pj_get_default_ctx();
+    else
+        ctx = pj_get_ctx (P);
+    ctx->debug_level = debuglevel;
+}
+
+
+/* Put a new logging function into P's context. The opaque object app_data is passed as first arg at each call to the logger */
+void pj_log_set (PJ *P, void *app_data, void (*log)(void *, int, const char *)) {
+    PJ_CONTEXT *ctx = pj_get_ctx (P);
+    ctx->app_data = app_data;
+    if (0!=log)
+        ctx->logger = log;
+}
+
+
+/* Move P to a new context - initially a copy of the default context */
+int pj_context_renew (PJ *P) {
+    PJ_CONTEXT *ctx = pj_ctx_alloc ();
+    if (0==ctx) {
+        pj_error_set (P, ENOMEM);
+        return 1;
+    }
+
+    pj_set_ctx (P, ctx);
+    return 0;
+}
+
+/* Move daughter to mother's context */
+void pj_context_inherit (PJ *mother, PJ *daughter) {
+    pj_set_ctx (daughter, pj_get_ctx (mother));
+}
+
+
+void pj_context_free (const PJ *P) {
+    PJ_CONTEXT *ctx;
+
+    if (0==P)
+        return;
+
+    /* During shutdown it should be OK to free the default context - but this needs more work */
+    if (pj_shutdown==P) {
+        /* The obvious solution "pj_ctx_free (pj_get_default_ctx ())" fails */
+        return;
+    }
+
+    ctx = pj_get_ctx ((PJ *) P);
+
+    /* Otherwise, trying to free the default context is a no-op */
+    if (pj_get_default_ctx ()==ctx)
+        return;
+
+    /* The common (?) case: free the context and move the PJ to the default context */
+    pj_ctx_free (ctx);
+    pj_set_ctx ((PJ *) P, pj_get_default_ctx ());
+}
+
+
+/* Minimum support for fileapi - which may never have been used... */
+void pj_fileapi_set (PJ *P, void *fileapi) {
+    PJ_CONTEXT *ctx = pj_get_ctx (P);
+    ctx->fileapi = (projFileAPI *) fileapi;
+}