Horner and helmert (#456)

Introducing the Horner polynomial evaluator also introduces the need for
very long +init:tag arguments (a n'th order 2D polynomium has
(n+1)(n+2)/2 coefficients, and n is typically in the range 5-10, i.e. up
to around 60 coefficients for each polynomium, and there are 4 polynomia
in a complete back/forward transformation set).

Hence, in this commit, along with the first part of the Horner code, the
code for reading +init files has been modified in a (for all practical
purposes) backwards compatible way, by making it possible to introduce
line continuations by escaping line breaks, i.e. preceding them with a
backslash.

An escaped line break works (as it would in TeX), by skipping all
following whitespace, including interspersed #-comments. This simple
extension makes it possible to create very long initialization elements
without losing track of the structure (cf. s45b.pol and pj_init_test.c
in the examples-directory for a demo).

The s45b.pol file was created by hand-editing the output of the software
doing the original constrained adjustment for the polynomial
coefficients. The simple adding of the “skip following whitespace and
comments” feature has made it possible to retain almost all metadata
from the source material.

This is considered very important, since 1) For the lack of a prior
common file format for geodetic polynomial coefficients, there is a good
chance that this will become THE standard, at least for the time being,
and 2) Without the metadata represented, it will be very hard for a
human to debug code involving a slightly misrepresented polynomium.

Due to the current architecture of the pj_init.c code (mostly around the
fill_buffer() function), it is next to impossible to implement the line
continuation functionality in full generality. Hence, it has been
necessary to limit this format extension to files smaller than 64 kB.

* Correction of spherical HEALpix test case

The first HEALpix test case in nad/testvarious is clearly intended to
invoke the spherical form of HEALpix.

It does, however, specify the spheroid using the +a=1 size parameter,
without specifying any shape parameter.

But since +no_defs is not specified either, a shape parameter is picked
up from the nad/proj_def.dat file (where ellps=WGS84 is given in the
<general> section).

It appears that this has not happened before I updated the pj_init code to support projection
pipelines (see below). I do, however, believe that the present behaviour is the correct one,
and rather than retrohacking the pj_init code, to (incorrectly, I
believe) reproduce the prior behaviour, I have corrected the test case
invocation in nad/testvarious to specify the spheroid using the +R=1
size parameter (which was already used in the following test case).

* Repair scaling of projections stomping on value of semimajor axis

* Workaround MSVC HUGE_VAL misimplementation.

The "return const err object" idiom (i.e.  const <type> err =
{HUGE_VAL,...}; ...  if (bad) return err) is problematic to implement
due to MSVC's misimplementation of HUGE_VAL as a non-const.

Hence, we need to run-time initialize these. In the pj_inv functions,
this was mistakenly done to the wrong object.

For pj_fwdobs/invobs and the remaining part of the obs-based API, this
is now worked around by providing functions returning a run time
HUGE_VAL initialized PJ_OBS or PJ_COO resp.

Obnoxious, but given MSVC's market penetration there is really not much
else we can do.

1 files changed, 33 insertions, 21 deletions

diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index ef86ccd3..27fe6ec9 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -114,6 +114,8 @@ struct pj_opaque {
     int *reverse_step;
     int *omit_forward;
     int *omit_inverse;
+    char **argv;
+    char **current_argv;
     PJ_OBS stack[PIPELINE_STACK_SIZE];
     PJ **pipeline;
 };
@@ -175,13 +177,17 @@ static LP     pipeline_reverse (XY xyz, PJ *P);
 
 
 static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
-    int i, first_step, last_step, incr;
+    int i, first_step, last_step;
 
     first_step = 1;
     last_step  = P->opaque->steps + 1;
-    incr  = 1;
 
-    for (i = first_step; i != last_step; i += incr) {
+    for (i = first_step;  i != last_step;  i++) {
+        pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %20.20s",
+            i-first_step, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+            P->opaque->pipeline[i]->descr
+        );
+
         if (P->opaque->omit_forward[i])
             continue;
         if (P->opaque->reverse_step[i])
@@ -191,19 +197,23 @@ static PJ_OBS pipeline_forward_obs (PJ_OBS point, PJ *P) {
         if (P->opaque->depth < PIPELINE_STACK_SIZE)
             P->opaque->stack[P->opaque->depth++] = point;
     }
+    pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
 
     P->opaque->depth = 0;    /* Clear the stack */
     return point;
 }
 
+
 static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
-    int i, first_step, last_step, incr;
+    int i, first_step, last_step;
 
     first_step = P->opaque->steps;
     last_step  =  0;
-    incr  = -1;
-
-    for (i = first_step; i != last_step; i += incr) {
+    for (i = first_step;  i != last_step;  i--) {
+        pj_log_trace (P, "In[%2.2d]: %20.15g %20.15g %20.15g - %.4f %.4f",
+            i - 1, point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z,
+            P->opaque->pipeline[i]->a, P->opaque->pipeline[i]->rf
+        );
         if (P->opaque->omit_inverse[i])
             continue;
         if (P->opaque->reverse_step[i])
@@ -213,6 +223,7 @@ static PJ_OBS pipeline_reverse_obs (PJ_OBS point, PJ *P) {
         if (P->opaque->depth < PIPELINE_STACK_SIZE)
             P->opaque->stack[P->opaque->depth++] = point;
     }
+    pj_log_trace (P, "Out[ ]: %20.15g %20.15g %20.15g", point.coo.xyz.x, point.coo.xyz.y, point.coo.xyz.z);
 
     P->opaque->depth = 0;    /* Clear the stack */
     return point;
@@ -264,7 +275,7 @@ static void *pipeline_freeup (PJ *P, int errlev) {         /* Destructor */
     if (0==P)
         return 0;
 
-    pj_error_set (P, errlev);
+    pj_err_level (P, errlev);
 
     if (0==P->opaque)
         return pj_dealloc (P);
@@ -275,6 +286,8 @@ static void *pipeline_freeup (PJ *P, int errlev) {         /* Destructor */
     pj_dealloc (P->opaque->reverse_step);
     pj_dealloc (P->opaque->omit_forward);
     pj_dealloc (P->opaque->omit_inverse);
+    pj_dealloc (P->opaque->argv);
+    pj_dealloc (P->opaque->current_argv);
     pj_dealloc (P->opaque->pipeline);
 
     pj_dealloc (P->opaque);
@@ -313,8 +326,9 @@ static PJ *pj_create_pipeline (PJ *P, size_t steps) {
     return P;
 }
 
+
 /* count the number of args in pipeline definition */
-size_t argc_params (paralist *params) {
+static size_t argc_params (paralist *params) {
     size_t argc = 0;
     for (; params != 0; params = params->next)
         argc++;
@@ -322,18 +336,18 @@ size_t argc_params (paralist *params) {
 }
 
 /* Sentinel for argument list */
-static char argv_sentinel[5] = "step";
+static char argv_sentinel[] = "step";
 
-/* turn paralist int argc/argv style argument list */
-char **argv_params (paralist *params) {
+/* turn paralist into argc/argv style argument list */
+static char **argv_params (paralist *params, size_t argc) {
     char **argv;
-    size_t argc = 0;
-    argv = pj_calloc (argc_params (params), sizeof (char *));
+    size_t i = 0;
+    argv = pj_calloc (argc, sizeof (char *));
     if (0==argv)
         return 0;
     for (; params != 0; params = params->next)
-        argv[argc++] = params->param;
-    argv[argc++] = argv_sentinel;
+        argv[i++] = params->param;
+    argv[i++] = argv_sentinel;
     return argv;
 }
 
@@ -356,17 +370,15 @@ PJ *PROJECTION(pipeline) {
         return 0;
 
     argc = argc_params (P->params);
-    argv = argv_params (P->params);
+    P->opaque->argv = argv = argv_params (P->params, argc);
     if (0==argv)
         return pipeline_freeup (P, ENOMEM);
 
-    /* The elements of current_argv are not used - we just use argv_params */
-    /* as allocator for a "large enough" container needed later            */
-    current_argv = argv_params (P->params);
+    P->opaque->current_argv = current_argv = pj_calloc (argc, sizeof (char *));
     if (0==current_argv)
         return pipeline_freeup (P, ENOMEM);
 
-    /* Do some syntactic sanity checking */
+    /* Do some syntactical sanity checking */
     for (i = 0;  i < argc;  i++) {
         if (0==strcmp ("step", argv[i])) {
             if (-1==i_pipeline) {