From 78f6a82a2ed9a3ad98f49e73c85a2d14aa9ee93c Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Mon, 7 May 2018 22:29:53 +0200 Subject: Add NAN definition to proj_math.h --- src/proj_math.h | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/proj_math.h b/src/proj_math.h index 1bb0a94f..1e0c6217 100644 --- a/src/proj_math.h +++ b/src/proj_math.h @@ -36,15 +36,30 @@ 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 + +#ifndef isnan #define isnan pj_isnan +#endif #endif /* !(defined(HAVE_C99_MATH) && HAVE_C99_MATH) */ -- cgit v1.2.3 From 987603a4ab2b925363fdcb6bf99fa89658d67613 Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sun, 6 May 2018 17:38:52 +0300 Subject: Add round() and lround() to proj_math.h --- src/pj_math.c | 28 ++++++++++++++++++++++++++++ src/proj_math.h | 1 + 2 files changed, 29 insertions(+) diff --git a/src/pj_math.c b/src/pj_math.c index 751c2750..97d4ae56 100644 --- a/src/pj_math.c +++ b/src/pj_math.c @@ -67,5 +67,33 @@ int pj_isnan (double x) { 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_math.h b/src/proj_math.h index 1e0c6217..ca3ab401 100644 --- a/src/proj_math.h +++ b/src/proj_math.h @@ -29,6 +29,7 @@ #define PROJ_MATH_H #include +#include #ifdef __cplusplus extern "C" { -- cgit v1.2.3 From 8fef2126f1c9fa17b79e6669f4457c299c0e33bf Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sun, 6 May 2018 17:40:44 +0300 Subject: Handle HUGE_VAL input better in fwd and inv functions --- src/pj_fwd.c | 5 +++-- src/pj_inv.c | 3 ++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/pj_fwd.c b/src/pj_fwd.c index 347f8334..4b6cfc58 100644 --- a/src/pj_fwd.c +++ b/src/pj_fwd.c @@ -30,6 +30,7 @@ #include #include "proj_internal.h" +#include "proj_math.h" #include "projects.h" #define INPUT_UNITS P->left @@ -37,7 +38,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 */ @@ -191,7 +192,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 327fef9f..25d416c0 100644 --- a/src/pj_inv.c +++ b/src/pj_inv.c @@ -30,13 +30,14 @@ #include #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 (); } -- cgit v1.2.3 From 58cbb9fe4f89b9febd780f7bdcfa4c2bb74a723e Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Sun, 6 May 2018 17:43:26 +0300 Subject: Replace int typecasts with calls to lround to avoid bad conversions on NaN input. Added test to check for those cases. --- src/PJ_isea.c | 52 ++++++++++++++++++++++++----------------------- src/PJ_robin.c | 8 ++++---- src/PJ_unitconvert.c | 12 ++++++----- src/gie.c | 13 ++++++++++++ src/nad_intr.c | 4 ++-- src/pj_apply_vgridshift.c | 8 ++++---- src/proj_etmerc.c | 4 ++-- 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 #include #include #include #include +#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 -#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 #include + +#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) -- cgit v1.2.3 From 6db260ff7857fcec63e5988ed43ecebc358a811d Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Wed, 23 May 2018 15:02:28 +0200 Subject: Add pj_isnan to proj.def to avoid linking errors for some MSVC versions --- src/pj_math.c | 21 +++++++++++++++------ src/proj.def | 2 ++ 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/src/pj_math.c b/src/pj_math.c index 97d4ae56..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,12 +76,6 @@ 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 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 -- cgit v1.2.3 From 8795737e29efe59526565d757bc560d9d19163cb Mon Sep 17 00:00:00 2001 From: Kristian Evers Date: Thu, 24 May 2018 14:56:27 +0200 Subject: Handle double to int typecasts in ISEA better Originally the code was doing double to int conversions like y = (int)(x + 0.5) which results in rounding when typecasting. In an earlier attempt to avoid buffer overflows in integer typecasts this was changed to y = lround(x + 0.5) which doesn't give the origial results. We fix that here by instead doing y = lround(x) It is safe to so as long as x is positive. --- src/PJ_isea.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/PJ_isea.c b/src/PJ_isea.c index c1ac2f00..3c811a18 100644 --- a/src/PJ_isea.c +++ b/src/PJ_isea.c @@ -679,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 = lround((sidelength * 2.0 + 0.5)); + maxcoord = lround((sidelength * 2.0)); v = *pt; hexbin2(hexwidth, v.x, v.y, &h.x, &h.y); @@ -750,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 = lround((pow(g->aperture, g->resolution / 2.0) + 0.5)); + sidelength = lround(pow(g->aperture, g->resolution / 2.0)); } else { sidelength = g->resolution; } @@ -833,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 = lround((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 = lround((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 = lround((pow(g->aperture, g->resolution / 2.0) + 0.5)); - sn = lround(((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; @@ -874,8 +874,8 @@ static int isea_hex(struct isea_dgg *g, int tri, return 1; #ifdef FIXME - d = lround(v.x); - i = lround(v.y); + d = lround(floor(v.x)); + i = lround(floor(v.y)); /* Aperture 3 odd resolutions */ if (g->aperture == 3 && g->resolution % 2 != 0) { @@ -903,7 +903,7 @@ static int isea_hex(struct isea_dgg *g, int tri, } /* aperture 3 even resolutions and aperture 4 */ - sidelength = lround((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; -- cgit v1.2.3