Replace int typecasts with calls to lround to avoid bad conversions on NaN input. Added test to check for those cases.

7 files changed, 59 insertions, 42 deletions

diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 3e941c77..14aec843 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -2,12 +2,18 @@
  * This code was entirely written by Nathan Wagner
  * and is in the public domain.
  */
+#define PJ_LIB__
 
+#include <errno.h>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 #include <float.h>
 
+#include "proj.h"
+#include "projects.h"
+#include "proj_math.h"
+
 #ifndef M_PI
 #  define M_PI 3.14159265358979323846
 #endif
@@ -59,7 +65,7 @@ ISEA_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 */
@@ -72,11 +78,11 @@ 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;
 
@@ -707,7 +713,7 @@ isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_p
     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 */
@@ -719,7 +725,7 @@ isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_p
     /* 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 + 0.5));
 
     v = *pt;
     hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
@@ -783,7 +789,7 @@ 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) {
@@ -791,7 +797,7 @@ isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
     }
     /* 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) + 0.5));
     } else {
         sidelength = g->resolution;
     }
@@ -866,29 +872,29 @@ int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
 /* q2di to seqnum */
 ISEA_STATIC
 int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
-    int             sidelength;
-    int             sn, height;
-    int             hexes;
+    long            sidelength;
+    long            sn, height;
+    long            hexes;
 
     if (quad == 0) {
         g->serial = 1;
         return g->serial;
     }
     /* hexes in a quad */
-    hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
+    hexes = lround((pow(g->aperture, g->resolution) + 0.5));
     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((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) + 0.5));
+        sn = lround(((quad - 1) * hexes + sidelength * di->x + di->y + 2));
     }
 
     g->serial = sn;
@@ -905,8 +911,8 @@ 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;
 
@@ -917,12 +923,12 @@ int isea_hex(struct isea_dgg *g, int tri,
 
     return 1;
 #ifdef FIXME
-    d = (int)v.x;
-    i = (int)v.y;
+    d = lround(v.x);
+    i = lround(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);
@@ -946,7 +952,7 @@ 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) + 0.5));
     if (g->quad == 0) {
         hex->x = 0;
         hex->y = sidelength;
@@ -1016,10 +1022,6 @@ isea_forward(struct isea_dgg *g, struct isea_geo *in)
  * Proj 4 integration code follows
  */
 
-#define PJ_LIB__
-#include <errno.h>
-#include "proj.h"
-#include "projects.h"
 
 PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
 
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 a076b76c..53bf5e81 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -66,6 +66,8 @@ Last update: 2017-05-16
 #define PJ_LIB__
 #include <time.h>
 #include <errno.h>
+
+#include "proj_math.h"
 #include "proj_internal.h"
 #include "projects.h"
 
@@ -152,14 +154,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);
@@ -228,9 +230,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 b81e6bdc..f9f76452 100644
--- a/src/gie.c
+++ b/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, 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 e7106b25..b5b35640 100644
--- a/src/pj_apply_vgridshift.c
+++ b/src/pj_apply_vgridshift.c
@@ -37,8 +37,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 */
@@ -97,8 +97,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/proj_etmerc.c b/src/proj_etmerc.c
index 9a89e2b0..a4088ad9 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)