Return NaN coordinate immediately when receiving NaN input

Ideally NaNs shouldn't get special treatment but since double to int
conversion is undefined behavior we are forced to treat NaNs as
something special. Therefor we return a NaN coordinate when ever we
encounter a NaN in one of the dimensions of an input coordinate. The
error level is not changed since this all math operations on NaNs should
return a NaN.

This policy of returning NaNs immediately should also eliminate a bunch
of potential double to int conversion bugs.

src/gie.c

@@ -1942,6 +1942,19 @@ static int cart_selftest (void) {
     if (P->f != 1.0/298.257223563) return 125;
     proj_destroy(P);
 
+    /* Test that pj_fwd* and pj_inv* returns NaNs when receiving NaN input */
+    P = proj_create(PJ_DEFAULT_CTX, "+proj=merc");
+    if (0==P) return 0;
+    a = proj_coord_nan();
+    a = proj_trans(P, PJ_FWD, a);
+    if ( !( isnan(a.v[0]) && isnan(a.v[1]) && isnan(a.v[2]) && isnan(a.v[3]) ) )
+        return 126;
+    a = proj_coord_nan();
+    a = proj_trans(P, PJ_INV, a);
+    if ( !( isnan(a.v[0]) && isnan(a.v[1]) && isnan(a.v[2]) && isnan(a.v[3]) ) )
+        return 127;
+    proj_destroy(P);
+
     return 0;
 }
 

src/pj_fwd.c

@@ -30,14 +30,15 @@
 #include <errno.h>
 
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #define INPUT_UNITS  P->left
 #define OUTPUT_UNITS P->right
 
 
 static PJ_COORD fwd_prepare (PJ *P, PJ_COORD coo) {
-    if (HUGE_VAL==coo.v[0])
+    if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1] || HUGE_VAL==coo.v[2])
         return proj_coord_error ();
 
     /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
@@ -189,9 +190,12 @@ XY pj_fwd(LP lp, PJ *P) {
 
     last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]))
+        return proj_coord_nan ().xy;
+
     if (!P->skip_fwd_prepare)
         coo = fwd_prepare (P, coo);
-    if (HUGE_VAL==coo.v[0])
+    if (HUGE_VAL==coo.v[0] || HUGE_VAL==coo.v[1])
         return proj_coord_error ().xy;
 
     /* Do the transformation, using the lowest dimensional transformer available */
@@ -223,6 +227,9 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
 
     last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]) || isnan(coo.v[2]))
+        return proj_coord_nan ().xyz;
+
     if (!P->skip_fwd_prepare)
         coo = fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
@@ -253,6 +260,9 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P) {
     int last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]) || isnan(coo.v[2]) || isnan(coo.v[3]))
+        return proj_coord_nan ();
+
     if (!P->skip_fwd_prepare)
         coo = fwd_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])

src/pj_internal.c

@@ -35,6 +35,7 @@
 
 #include "geodesic.h"
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 
@@ -60,6 +61,12 @@ PJ_COORD proj_coord_error (void) {
     return c;
 }
 
+/* Similar to proj_coord_error but produces a NaN coordinate instead */
+PJ_COORD proj_coord_nan (void) {
+    PJ_COORD c;
+    c.v[0] = c.v[1] = c.v[2] = c.v[3] = NAN;
+    return c;
+}
 
 
 /**************************************************************************************/

src/pj_inv.c

@@ -30,13 +30,14 @@
 #include <errno.h>
 
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #define INPUT_UNITS  P->right
 #define OUTPUT_UNITS P->left
 
 static PJ_COORD inv_prepare (PJ *P, PJ_COORD coo) {
-    if (coo.xyz.x == HUGE_VAL) {
+    if (coo.v[0] == HUGE_VAL || coo.v[1] == HUGE_VAL || coo.v[2] == HUGE_VAL) {
         proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
         return proj_coord_error ();
     }
@@ -186,6 +187,9 @@ LP pj_inv(XY xy, PJ *P) {
 
     last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]))
+        return proj_coord_nan ().lp;
+
     if (!P->skip_inv_prepare)
         coo = inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
@@ -220,6 +224,9 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) {
 
     last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]) || isnan(coo.v[2]))
+        return proj_coord_nan ().lpz;
+
     if (!P->skip_inv_prepare)
         coo = inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])
@@ -250,6 +257,9 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) {
 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P) {
     int last_errno = proj_errno_reset(P);
 
+    if (isnan(coo.v[0]) || isnan(coo.v[1]) || isnan(coo.v[2]) || isnan(coo.v[3]))
+        return proj_coord_nan ();
+
     if (!P->skip_inv_prepare)
         coo = inv_prepare (P, coo);
     if (HUGE_VAL==coo.v[0])

src/proj_internal.h

@@ -78,6 +78,7 @@ enum pj_io_units pj_left (PJ *P);
 enum pj_io_units pj_right (PJ *P);
 
 PJ_COORD proj_coord_error (void);
+PJ_COORD proj_coord_nan (void);
 
 void proj_context_errno_set (PJ_CONTEXT *ctx, int err);
 void proj_context_set (PJ *P, PJ_CONTEXT *ctx);

src/proj_math.h

@@ -35,6 +35,9 @@ extern "C" {
 
 #if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
 
+#ifndef NAN
+#define NAN sqrt(-1)
+#endif
 
 double  pj_hypot(double x, double y);
 double  pj_log1p(double x);