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

1 files changed, 27 insertions, 25 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";