Merge pull request #949 from kbevers/return-nans-quickly

Return NaN coordinate immediately when receiving NaN input

12 files changed, 116 insertions, 43 deletions

diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 4e0eda32..3c811a18 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -4,14 +4,15 @@
  */
 
 #include <errno.h>
-#include <float.h>
 #include <math.h>
+#include <float.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <string.h>
 
 #define PJ_LIB__
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "proj.h"
 #include "projects.h"
 
@@ -79,7 +80,7 @@ static void hex_iso(struct hex *h) {
 static void hexbin2(double width, double x, double y, int *i, int *j) {
     double z, rx, ry, rz;
     double abs_dx, abs_dy, abs_dz;
-    int ix, iy, iz, s;
+    long ix, iy, iz, s;
     struct hex h;
 
     x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
@@ -92,11 +93,11 @@ static void hexbin2(double width, double x, double y, int *i, int *j) {
     z = -x - y;
 
     rx = floor(x + 0.5);
-    ix = (int)rx;
+    ix = lround(rx);
     ry = floor(y + 0.5);
-    iy = (int)ry;
+    iy = lround(ry);
     rz = floor(z + 0.5);
-    iz = (int)rz;
+    iz = lround(rz);
 
     s = ix + iy + iz;
 
@@ -666,7 +667,7 @@ static int isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt,
     double          hexwidth;
     double          sidelength; /* in hexes */
     int             d, i;
-    int             maxcoord;
+    long            maxcoord;
     struct hex      h;
 
     /* This is the number of hexes from apex to base of a triangle */
@@ -678,7 +679,7 @@ static int isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt,
     /* TODO I think sidelength is always x.5, so
      * (int)sidelength * 2 + 1 might be just as good
      */
-    maxcoord = (int) (sidelength * 2.0 + 0.5);
+    maxcoord = lround((sidelength * 2.0));
 
     v = *pt;
     hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
@@ -741,7 +742,7 @@ static int isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt,
                      struct isea_pt *di) {
     struct isea_pt  v;
     double          hexwidth;
-    int             sidelength; /* in hexes */
+    long            sidelength; /* in hexes */
     struct hex      h;
 
     if (g->aperture == 3 && g->resolution % 2 != 0) {
@@ -749,7 +750,7 @@ static int isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt,
     }
     /* todo might want to do this as an iterated loop */
     if (g->aperture >0) {
-        sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+        sidelength = lround(pow(g->aperture, g->resolution / 2.0));
     } else {
         sidelength = g->resolution;
     }
@@ -821,6 +822,7 @@ static int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
 }
 
 /* q2di to seqnum */
+
 static int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
     int             sidelength;
     int             sn, height;
@@ -831,20 +833,20 @@ static int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
         return g->serial;
     }
     /* hexes in a quad */
-    hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
+    hexes = lround(pow(g->aperture, g->resolution));
     if (quad == 11) {
         g->serial = 1 + 10 * hexes + 1;
         return g->serial;
     }
     if (g->aperture == 3 && g->resolution % 2 == 1) {
-        height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
+        height = lround(floor((pow(g->aperture, (g->resolution - 1) / 2.0))));
         sn = ((int) di->x) * height;
         sn += ((int) di->y) / height;
         sn += (quad - 1) * hexes;
         sn += 2;
     } else {
-        sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
-        sn = (int) ((quad - 1) * hexes + sidelength * di->x + di->y + 2);
+        sidelength = lround((pow(g->aperture, g->resolution / 2.0)));
+        sn = lround(floor(((quad - 1) * hexes + sidelength * di->x + di->y + 2)));
     }
 
     g->serial = sn;
@@ -860,8 +862,8 @@ static int isea_hex(struct isea_dgg *g, int tri,
                     struct isea_pt *pt, struct isea_pt *hex) {
     struct isea_pt v;
 #ifdef FIXME
-    int sidelength;
-    int d, i, x, y;
+    long sidelength;
+    long d, i, x, y;
 #endif
     int quad;
 
@@ -872,12 +874,12 @@ static int isea_hex(struct isea_dgg *g, int tri,
 
     return 1;
 #ifdef FIXME
-    d = (int)v.x;
-    i = (int)v.y;
+    d = lround(floor(v.x));
+    i = lround(floor(v.y));
 
     /* Aperture 3 odd resolutions */
     if (g->aperture == 3 && g->resolution % 2 != 0) {
-        int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
+        long offset = lround((pow(3.0, g->resolution - 1) + 0.5));
 
         d += offset * ((g->quad-1) % 5);
         i += offset * ((g->quad-1) % 5);
@@ -901,7 +903,7 @@ static int isea_hex(struct isea_dgg *g, int tri,
     }
 
     /* aperture 3 even resolutions and aperture 4 */
-    sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+    sidelength = lround((pow(g->aperture, g->resolution / 2.0)));
     if (g->quad == 0) {
         hex->x = 0;
         hex->y = sidelength;
diff --git a/src/PJ_robin.c b/src/PJ_robin.c
index a1d5ad4d..ca9e5fc7 100644
--- a/src/PJ_robin.c
+++ b/src/PJ_robin.c
@@ -78,12 +78,12 @@ static const struct COEFS Y[] = {
 
 static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
     XY xy = {0.0,0.0};
-    int i;
+    long i;
     double dphi;
     (void) P;
 
     dphi = fabs(lp.phi);
-    i = isnan(lp.phi) ? -1 : (int)floor(dphi * C1);
+    i = isnan(lp.phi) ? -1 : lround(floor(dphi * C1));
     if( i < 0 ){
         proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
         return xy;
@@ -100,7 +100,7 @@ static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
 
 static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
     LP lp = {0.0,0.0};
-    int i;
+    long i;
     double t, t1;
     struct COEFS T;
     int iters;
@@ -118,7 +118,7 @@ static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
         }
     } else { /* general problem */
         /* in Y space, reduce to table interval */
-        i = isnan(lp.phi) ? -1 : (int)floor(lp.phi * NODES);
+        i = isnan(lp.phi) ? -1 : lround(floor(lp.phi * NODES));
         if( i < 0 || i >= NODES ) {
             proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
             return lp;
diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index 5c25a16b..a3051ad7 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -71,6 +71,7 @@ Last update: 2017-05-16
 #include <time.h>
 
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 PROJ_HEAD(unitconvert, "Unit conversion");
@@ -156,14 +157,14 @@ static double decimalyear_to_mjd(double decimalyear) {
 /***********************************************************************
     Epoch of modified julian date is 1858-11-16 00:00
 ************************************************************************/
-    int year;
+    long year;
     double fractional_year;
     double mjd;
 
     if( decimalyear < -10000 || decimalyear > 10000 )
         return 0;
 
-    year = (int)floor(decimalyear);
+    year = lround(floor(decimalyear));
     fractional_year = decimalyear - year;
     mjd = (year - 1859)*365 + 14 + 31;
     mjd += fractional_year*days_in_year(year);
@@ -232,9 +233,9 @@ static double yyyymmdd_to_mjd(double yyyymmdd) {
     Date given in YYYY-MM-DD format.
 ************************************************************************/
 
-    int year  = (int)floor( yyyymmdd / 10000 );
-    int month = (int)floor((yyyymmdd - year*10000) / 100);
-    int day   = (int)floor( yyyymmdd - year*10000 - month*100);
+    long year  = lround(floor( yyyymmdd / 10000 ));
+    long month = lround(floor((yyyymmdd - year*10000) / 100));
+    long day   = lround(floor( yyyymmdd - year*10000 - month*100));
     double mjd = daynumber_in_year(year, month, day);
 
     for (year -= 1; year > 1858; year--)
diff --git a/src/gie.c b/src/gie.c
index 7b829c7f..d5c5ee8b 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -1943,6 +1943,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, NAN, NAN, 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, NAN, NAN, 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;
 }
 
diff --git a/src/nad_intr.c b/src/nad_intr.c
index 80e428ed..1f9d1e0c 100644
--- a/src/nad_intr.c
+++ b/src/nad_intr.c
@@ -14,9 +14,9 @@ nad_intr(LP t, struct CTABLE *ct) {
 	int in;
 
 	t.lam /= ct->del.lam;
-	indx.lam = isnan(t.lam) ? 0 : (int)floor(t.lam);
+	indx.lam = isnan(t.lam) ? 0 : (pj_int32)lround(floor(t.lam));
 	t.phi /= ct->del.phi;
-	indx.phi = isnan(t.phi) ? 0 : (int)floor(t.phi);
+	indx.phi = isnan(t.phi) ? 0 : (pj_int32)lround(floor(t.phi));
 
 	frct.lam = t.lam - indx.lam;
 	frct.phi = t.phi - indx.phi;
diff --git a/src/pj_apply_vgridshift.c b/src/pj_apply_vgridshift.c
index ca932674..c1344951 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -48,8 +48,8 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
     int  itable = 0;
     double value = HUGE_VAL;
     double grid_x, grid_y;
-    int    grid_ix, grid_iy;
-    int    grid_ix2, grid_iy2;
+    long   grid_ix, grid_iy;
+    long   grid_ix2, grid_iy2;
     float  *cvs;
     /* do not deal with NaN coordinates */
     /* cppcheck-suppress duplicateExpression */
@@ -108,8 +108,8 @@ static double read_vgrid_value( PJ *defn, LP input, int *gridlist_count_p, PJ_GR
         /* Interpolation a location within the grid */
         grid_x = (input.lam - ct->ll.lam) / ct->del.lam;
         grid_y = (input.phi - ct->ll.phi) / ct->del.phi;
-        grid_ix = (int) floor(grid_x);
-        grid_iy = (int) floor(grid_y);
+        grid_ix = lround(floor(grid_x));
+        grid_iy = lround(floor(grid_y));
         grid_x -= grid_ix;
         grid_y -= grid_iy;
 
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index fe04f733..b5f1b36e 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -32,6 +32,7 @@
 #include <math.h>
 
 #include "proj_internal.h"
+#include "proj_math.h"
 #include "projects.h"
 
 #define INPUT_UNITS  P->left
@@ -39,7 +40,7 @@
 
 
 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 */
@@ -193,7 +194,7 @@ XY pj_fwd(LP lp, PJ *P) {
 
     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 */
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 72e06fe0..ca149674 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -31,13 +31,14 @@
 #include <math.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 ();
     }
diff --git a/src/pj_math.c b/src/pj_math.c
index 751c2750..540ab9eb 100644
--- a/src/pj_math.c
+++ b/src/pj_math.c
@@ -27,6 +27,21 @@
 
 #include "proj_math.h"
 
+/* pj_isnan is used in gie.c which means that is has to */
+/* be exported in the Windows DLL and therefore needs   */
+/* to be declared even though we have isnan() on the    */
+/* system.                                              */
+
+#ifdef HAVE_C99_MATH
+int pj_isnan (double x);
+#endif
+
+/* Returns 0 if not a NaN and non-zero if val is a NaN */
+int pj_isnan (double x) {
+    /* cppcheck-suppress duplicateExpression */
+    return x != x;
+}
+
 #if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
 
 /* Compute hypotenuse */
@@ -61,11 +76,33 @@ double pj_asinh(double x) {
     return x > 0 ? y : (x < 0 ? -y : x);
 }
 
-/* Returns 0 if not a NaN and non-zero if val is a NaN */
-int pj_isnan (double x) {
-    /* cppcheck-suppress duplicateExpression */
-    return x != x;
+double pj_round(double x) {
+  /* The handling of corner cases is copied from boost; see
+   *   https://github.com/boostorg/math/pull/8
+   * with improvements to return -0.0 when appropriate */
+  double t;
+  if (x == 0)
+    return x;               /* Retain sign of 0 */
+  else if (0 < x && x <  0.5)
+    return +0.0;
+  else if (0 > x && x > -0.5)
+    return -0.0;
+  else if (x > 0) {
+    t = ceil(x);
+    return 0.5 < t - x ? t - 1 : t;
+  } else {                  /* Includes NaN */
+    t = floor(x);
+    return 0.5 < x - t ? t + 1 : t;
+  }
 }
 
+long pj_lround(double x) {
+  /* Default value for overflow + NaN + (x == LONG_MIN) */
+  long r = LONG_MIN;
+  x = round(x);
+  if (fabs(x) < -(double)LONG_MIN) /* Assume (double)LONG_MIN is exact */
+    r = (int)x;
+  return r;
+}
 
 #endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */
diff --git a/src/proj.def b/src/proj.def
index d6c84dc3..9b0a95d8 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -155,3 +155,5 @@ EXPORTS
 
     proj_geod                      @140
     proj_context_errno             @141
+
+    pj_isnan                       @142
diff --git a/src/proj_etmerc.c b/src/proj_etmerc.c
index 06da604d..457edcb8 100644
--- a/src/proj_etmerc.c
+++ b/src/proj_etmerc.c
@@ -310,7 +310,7 @@ PJ *PROJECTION(etmerc) {
 
 
 PJ *PROJECTION(utm) {
-    int zone;
+    long zone;
     struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
     if (0==Q)
         return pj_default_destructor (P, ENOMEM);
@@ -337,7 +337,7 @@ PJ *PROJECTION(utm) {
     }
     else /* nearest central meridian input */
     {
-        zone = (int)(floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI));
+        zone = lround((floor ((adjlon (P->lam0) + M_PI) * 30. / M_PI)));
         if (zone < 0)
             zone = 0;
         else if (zone >= 60)
diff --git a/src/proj_math.h b/src/proj_math.h
index 84c9b0da..f6dd7e09 100644
--- a/src/proj_math.h
+++ b/src/proj_math.h
@@ -29,6 +29,7 @@
 #define PROJ_MATH_H
 
 #include <math.h>
+#include <limits.h>
 
 #ifdef __cplusplus
 extern "C" {
@@ -36,17 +37,32 @@ extern "C" {
 
 #if !(defined(HAVE_C99_MATH) && HAVE_C99_MATH)
 
+#ifndef NAN
+#ifdef _WIN32
+#define NAN sqrt(-1.0)
+#else
+#define NAN 0.0/0.0
+#endif
+#endif
+
 double  pj_hypot(double x, double y);
 double  pj_log1p(double x);
 double  pj_asinh(double x);
+double  pj_round(double x);
+long    pj_lround(double x);
 int     pj_isnan(double x);
 
 #define hypot   pj_hypot
 #define log1p   pj_log1p
 #define asinh   pj_asinh
+#define round   pj_round
+#define lround  pj_lround
+
+
 #ifdef isnan
 #undef isnan
 #endif
+
 #define isnan   pj_isnan
 
 #endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */