Skip to content

Commit

Permalink
Replace int typecasts with calls to lround to avoid bad conversions o…
Browse files Browse the repository at this point in the history 
…n NaN input. Added test to check for those cases.
  • Loading branch information
@kbevers
kbevers committed May 8, 2018
1 parent 3d1fa57 commit d9f0d6c
Show file tree
Hide file tree
Showing 7 changed files with 59 additions and 42 deletions.
52 changes: 27 additions & 25 deletions src/PJ_isea.c
View file
Expand Up @@ -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
Expand Down Expand Up @@ -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 */
Expand All @@ -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;

Expand Down Expand Up @@ -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 */
Expand All @@ -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);
Expand Down Expand Up @@ -783,15 +789,15 @@ 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) {
return isea_dddi_ap3odd(g, quad, 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;
}
Expand Down Expand Up @@ -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;
Expand All @@ -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;

Expand All @@ -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);
Expand All @@ -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;
Expand Down Expand Up @@ -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";

Expand Down
8 changes: 4 additions & 4 deletions src/PJ_robin.c
View file
Expand Up @@ -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;
Expand All @@ -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;
Expand All @@ -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;
Expand Down
12 changes: 7 additions & 5 deletions src/PJ_unitconvert.c
View file
Expand Up @@ -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"

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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--)
Expand Down
13 changes: 13 additions & 0 deletions src/gie.c
View file
Expand Up @@ -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;
}

Expand Down
4 changes: 2 additions & 2 deletions src/nad_intr.c
View file
Expand Up @@ -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;
Expand Down
8 changes: 4 additions & 4 deletions src/pj_apply_vgridshift.c
View file
Expand Up @@ -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 */
Expand Down Expand Up @@ -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;

Expand Down
4 changes: 2 additions & 2 deletions src/proj_etmerc.c
View file
Expand Up @@ -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);
Expand All @@ -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)
Expand Down

0 comments on commit d9f0d6c

Please sign in to comment.