Skip to content

Commit

Permalink
Eliminated load time non-computable initializers in healpix
Browse files Browse the repository at this point in the history
  • Loading branch information
@kbevers
kbevers committed Jul 13, 2016
1 parent b9b9299 commit 196a1bb
Show file tree
Hide file tree
Showing 2 changed files with 250 additions and 236 deletions.
206 changes: 105 additions & 101 deletions src/PJ_healpix.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,8 @@ struct pj_opaque {
};

typedef struct {
int cn; /* An integer 0--3 indicating the position of the polar cap. */
double x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
int cn; /* An integer 0--3 indicating the position of the polar cap. */
double x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
enum Region {north, south, equatorial} region;
} CapMap;

Expand Down Expand Up @@ -160,41 +160,60 @@ static int pnpoly(int nvert, double vert[][2], double testx, double testy) {
int in_image(double x, double y, int proj, int north_square, int south_square) {
if (proj == 0) {
double healpixVertsJit[][2] = {
{-1.0*M_PI- EPS, M_PI/4.0},
{-3.0*M_PI/4.0, M_PI/2.0 + EPS},
{-1.0*M_PI/2.0, M_PI/4.0 + EPS},
{-1.0*M_PI/4.0, M_PI/2.0 + EPS},
{0.0, M_PI/4.0 + EPS},
{M_PI/4.0, M_PI/2.0 + EPS},
{M_PI/2.0, M_PI/4.0 + EPS},
{3.0*M_PI/4.0, M_PI/2.0 + EPS},
{M_PI+ EPS, M_PI/4.0},
{M_PI+ EPS, -1.0*M_PI/4.0},
{3.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
{M_PI/2.0, -1.0*M_PI/4.0 - EPS},
{M_PI/4.0, -1.0*M_PI/2.0 - EPS},
{0.0, -1.0*M_PI/4.0 - EPS},
{-1.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
{-1.0*M_PI/2.0, -1.0*M_PI/4.0 - EPS},
{-3.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
{-1.0*M_PI - EPS, -1.0*M_PI/4.0}
{-1.0*M_PI - EPS, M_FORTPI},
{-3.0*M_FORTPI, M_HALFPI + EPS},
{M_HALFPI, M_FORTPI + EPS},
{M_FORTPI, M_HALFPI + EPS},
{0.0, M_FORTPI + EPS},
{M_FORTPI, M_HALFPI + EPS},
{M_HALFPI, M_FORTPI + EPS},
{3.0*M_FORTPI, M_HALFPI + EPS},
{M_PI + EPS, M_FORTPI},
{M_PI + EPS, -M_FORTPI},
{3.0*M_FORTPI, -1.0*M_HALFPI - EPS},
{M_HALFPI, -1.0*M_FORTPI - EPS},
{M_FORTPI, -1.0*M_HALFPI - EPS},
{0.0, -1.0*M_FORTPI - EPS},
{-1.0*M_FORTPI, -1.0*M_HALFPI - EPS},
{-1.0*M_HALFPI, -1.0*M_FORTPI - EPS},
{-3.0*M_FORTPI, -1.0*M_HALFPI - EPS},
{-1.0*M_PI - EPS, -1.0*M_FORTPI}
};
return pnpoly((int)sizeof(healpixVertsJit)/
sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
} else {
double rhealpixVertsJit[][2] = {
{-1.0*M_PI - EPS, M_PI/4.0 + EPS},
{-1.0*M_PI + north_square*M_PI/2.0- EPS, M_PI/4.0 + EPS},
{-1.0*M_PI + north_square*M_PI/2.0- EPS, 3*M_PI/4.0 + EPS},
{-1.0*M_PI + (north_square + 1.0)*M_PI/2.0 + EPS, 3*M_PI/4.0 + EPS},
{-1.0*M_PI + (north_square + 1.0)*M_PI/2.0 + EPS, M_PI/4.0 + EPS},
{M_PI + EPS, M_PI/4.0 + EPS},
{M_PI + EPS, -1.0*M_PI/4.0 - EPS},
{-1.0*M_PI + (south_square + 1.0)*M_PI/2.0 + EPS, -1.0*M_PI/4.0 - EPS},
{-1.0*M_PI + (south_square + 1.0)*M_PI/2.0 + EPS, -3.0*M_PI/4.0 - EPS},
{-1.0*M_PI + south_square*M_PI/2.0 - EPS, -3.0*M_PI/4.0 - EPS},
{-1.0*M_PI + south_square*M_PI/2.0 - EPS, -1.0*M_PI/4.0 - EPS},
{-1.0*M_PI - EPS, -1.0*M_PI/4.0 - EPS}};
/**
* Assigning each element by index to avoid warnings such as
* 'initializer element is not computable at load time'.
* Before C99 this was not allowed and to keep as portable as
* possible we do it the C89 way here.
**/
double rhealpixVertsJit[12][2];
rhealpixVertsJit[0][0] = -1.0*M_PI - EPS;
rhealpixVertsJit[0][1] = M_FORTPI + EPS;
rhealpixVertsJit[1][0] = -1.0*M_PI + north_square*M_HALFPI- EPS;
rhealpixVertsJit[1][2] = M_FORTPI + EPS;
rhealpixVertsJit[2][0] = -1.0*M_PI + north_square*M_HALFPI- EPS;
rhealpixVertsJit[2][1] = 3*M_FORTPI + EPS;
rhealpixVertsJit[3][0] = -1.0*M_PI + (north_square + 1.0)*M_HALFPI + EPS;
rhealpixVertsJit[3][1] = 3*M_FORTPI + EPS;
rhealpixVertsJit[4][0] = -1.0*M_PI + (north_square + 1.0)*M_HALFPI + EPS;
rhealpixVertsJit[4][1] = M_FORTPI + EPS;
rhealpixVertsJit[5][0] = M_PI + EPS;
rhealpixVertsJit[5][1] = M_FORTPI + EPS;
rhealpixVertsJit[6][0] = M_PI + EPS;
rhealpixVertsJit[6][1] = -1.0*M_FORTPI - EPS;
rhealpixVertsJit[7][0] = -1.0*M_PI + (south_square + 1.0)*M_HALFPI + EPS;
rhealpixVertsJit[7][1] = -1.0*M_FORTPI - EPS;
rhealpixVertsJit[8][0] = -1.0*M_PI + (south_square + 1.0)*M_HALFPI + EPS;
rhealpixVertsJit[8][1] = -3.0*M_FORTPI - EPS;
rhealpixVertsJit[9][0] = -1.0*M_PI + south_square*M_HALFPI - EPS;
rhealpixVertsJit[9][1] = -3.0*M_FORTPI - EPS;
rhealpixVertsJit[10][0] = -1.0*M_PI + south_square*M_HALFPI - EPS;
rhealpixVertsJit[10][1] = -1.0*M_FORTPI - EPS;
rhealpixVertsJit[11][0] = -1.0*M_PI - EPS;
rhealpixVertsJit[11][1] = -1.0*M_FORTPI - EPS;

return pnpoly((int)sizeof(rhealpixVertsJit)/
sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
}
Expand Down Expand Up @@ -247,9 +266,9 @@ XY healpix_sphere(LP lp) {
if (cn >= 4) {
cn = 3;
}
lamc = -3*M_PI/4 + (M_PI/2)*cn;
lamc = -3*M_FORTPI + (M_HALFPI)*cn;
xy.x = lamc + (lam - lamc)*sigma;
xy.y = pj_sign(phi)*M_PI/4*(2 - sigma);
xy.y = pj_sign(phi)*M_FORTPI*(2 - sigma);
}
return xy;
}
Expand All @@ -262,25 +281,25 @@ LP healpix_sphere_inverse(XY xy) {
LP lp;
double x = xy.x;
double y = xy.y;
double y0 = M_PI/4.0;
double y0 = M_FORTPI;

/* Equatorial region. */
if (fabsl(y) <= y0) {
lp.lam = x;
lp.phi = asin(8.0*y/(3.0*M_PI));
} else if (fabsl(y) < M_PI/2.0) {
} else if (fabsl(y) < M_HALFPI) {
double cn = floor(2.0*x/M_PI + 2.0);
double xc, tau;
if (cn >= 4) {
cn = 3;
}
xc = -3.0*M_PI/4.0 + (M_PI/2.0)*cn;
xc = -3.0*M_FORTPI + (M_HALFPI)*cn;
tau = 2.0 - 4.0*fabsl(y)/M_PI;
lp.lam = xc + (x - xc)/tau;
lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
} else {
lp.lam = -1.0*M_PI;
lp.phi = pj_sign(y)*M_PI/2.0;
lp.phi = pj_sign(y)*M_HALFPI;
}
return (lp);
}
Expand Down Expand Up @@ -343,48 +362,48 @@ static CapMap get_cap(double x, double y, int north_square, int south_square,
capmap.x = x;
capmap.y = y;
if (inverse == 0) {
if (y > M_PI/4.0) {
if (y > M_FORTPI) {
capmap.region = north;
c = M_PI/2.0;
} else if (y < -1*M_PI/4.0) {
c = M_HALFPI;
} else if (y < -1*M_FORTPI) {
capmap.region = south;
c = -1*M_PI/2.0;
c = -1*M_HALFPI;
} else {
capmap.region = equatorial;
capmap.cn = 0;
return capmap;
}
/* polar region */
if (x < -1*M_PI/2.0) {
if (x < -1*M_HALFPI) {
capmap.cn = 0;
capmap.x = (-1*3.0*M_PI/4.0);
capmap.x = (-1*3.0*M_FORTPI);
capmap.y = c;
} else if (x >= -1*M_PI/2.0 && x < 0) {
} else if (x >= -1*M_HALFPI && x < 0) {
capmap.cn = 1;
capmap.x = -1*M_PI/4.0;
capmap.x = -1*M_FORTPI;
capmap.y = c;
} else if (x >= 0 && x < M_PI/2.0) {
} else if (x >= 0 && x < M_HALFPI) {
capmap.cn = 2;
capmap.x = M_PI/4.0;
capmap.x = M_FORTPI;
capmap.y = c;
} else {
capmap.cn = 3;
capmap.x = 3.0*M_PI/4.0;
capmap.x = 3.0*M_FORTPI;
capmap.y = c;
}
return capmap;
} else {
double eps;
if (y > M_PI/4.0) {
if (y > M_FORTPI) {
capmap.region = north;
capmap.x = (-3.0*M_PI/4.0 + north_square*M_PI/2.0);
capmap.y = M_PI/2.0;
x = x - north_square*M_PI/2.0;
} else if (y < -1*M_PI/4.0) {
capmap.x = (-3.0*M_FORTPI + north_square*M_HALFPI);
capmap.y = M_HALFPI;
x = x - north_square*M_HALFPI;
} else if (y < -1*M_FORTPI) {
capmap.region = south;
capmap.x = (-3.0*M_PI/4.0 + south_square*M_PI/2);
capmap.y = -1*M_PI/2.0;
x = x - south_square*M_PI/2.0;
capmap.x = (-3.0*M_FORTPI + south_square*M_HALFPI);
capmap.y = -1*M_HALFPI;
x = x - south_square*M_HALFPI;
} else {
capmap.region = equatorial;
capmap.cn = 0;
Expand All @@ -394,21 +413,21 @@ static CapMap get_cap(double x, double y, int north_square, int south_square,
x, y moves to when rHEALPix polar square is disassembled. */
eps = 1e-15; /* Kludge. Fuzz to avoid some rounding errors. */
if (capmap.region == north) {
if (y >= -1*x - M_PI/4.0 - eps && y < x + 5.0*M_PI/4.0 - eps) {
if (y >= -1*x - M_FORTPI - eps && y < x + 5.0*M_FORTPI - eps) {
capmap.cn = (north_square + 1) % 4;
} else if (y > -1*x -1*M_PI/4.0 + eps && y >= x + 5.0*M_PI/4.0 - eps) {
} else if (y > -1*x -1*M_FORTPI + eps && y >= x + 5.0*M_FORTPI - eps) {
capmap.cn = (north_square + 2) % 4;
} else if (y <= -1*x -1*M_PI/4.0 + eps && y > x + 5.0*M_PI/4.0 + eps) {
} else if (y <= -1*x -1*M_FORTPI + eps && y > x + 5.0*M_FORTPI + eps) {
capmap.cn = (north_square + 3) % 4;
} else {
capmap.cn = north_square;
}
} else if (capmap.region == south) {
if (y <= x + M_PI/4.0 + eps && y > -1*x - 5.0*M_PI/4 + eps) {
if (y <= x + M_FORTPI + eps && y > -1*x - 5.0*M_FORTPI + eps) {
capmap.cn = (south_square + 1) % 4;
} else if (y < x + M_PI/4.0 - eps && y <= -1*x - 5.0*M_PI/4.0 + eps) {
} else if (y < x + M_FORTPI - eps && y <= -1*x - 5.0*M_FORTPI + eps) {
capmap.cn = (south_square + 2) % 4;
} else if (y >= x + M_PI/4.0 - eps && y < -1*x - 5.0*M_PI/4.0 - eps) {
} else if (y >= x + M_FORTPI - eps && y < -1*x - 5.0*M_FORTPI - eps) {
capmap.cn = (south_square + 3) % 4;
} else {
capmap.cn = south_square;
Expand All @@ -433,72 +452,57 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
XY xy;
double v[2];
double a[2];
double c[2];
double vector[2];
double v_min_c[2];
double ret_dot[2];
double (*tmpRot)[2];
int pole = 0;

CapMap capmap = get_cap(x, y, north_square, south_square, inverse);
if (capmap.region == equatorial) {
xy.x = capmap.x;
xy.y = capmap.y;
return xy;
}
v[0] = x;
v[1] = y;

v[0] = x; v[1] = y;
c[0] = capmap.x; c[1] = capmap.y;

if (inverse == 0) {
/* Rotate (x, y) about its polar cap tip and then translate it to
north_square or south_square. */
int pole = 0;
double (*tmpRot)[2];
double c[2] = {capmap.x, capmap.y};
a[0] = -3.0*M_FORTPI + pole*M_HALFPI;
a[1] = M_HALFPI;
if (capmap.region == north) {
pole = north_square;
a[0] = (-3.0*M_PI/4.0 + pole*M_PI/2);
a[1] = (M_PI/2.0 + pole*0);
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
vector_sub(v, c, v_min_c);
dot_product(tmpRot, v_min_c, ret_dot);
vector_add(ret_dot, a, vector);
} else {
pole = south_square;
a[0] = (-3.0*M_PI/4.0 + pole*M_PI/2);
a[1] = (M_PI/-2.0 + pole*0);

This comment has been minimized.

Sign in to view

Copy link
@pelson

pelson Dec 21, 2018

Contributor

Think we've lost a sign here.

This comment has been minimized.

Sign in to view

Copy link
@pelson

pelson Dec 21, 2018

Contributor

Raised in #1206
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
vector_sub(v, c, v_min_c);
dot_product(tmpRot, v_min_c, ret_dot);
vector_add(ret_dot, a, vector);
}
xy.x = vector[0];
xy.y = vector[1];
return xy;
} else {
/* Inverse function.
Unrotate (x, y) and then translate it back. */
int pole = 0;
double (*tmpRot)[2];
double c[2] = {capmap.x, capmap.y};
a[0] = -3.0*M_FORTPI + capmap.cn*M_HALFPI;
a[1] = M_HALFPI;
/* disassemble */
if (capmap.region == north) {
pole = north_square;
a[0] = (-3.0*M_PI/4.0 + capmap.cn*M_PI/2);
a[1] = (M_PI/2.0 + capmap.cn*0);
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
vector_sub(v, c, v_min_c);
dot_product(tmpRot, v_min_c, ret_dot);
vector_add(ret_dot, a, vector);
} else {
pole = south_square;
a[0] = (-3.0*M_PI/4.0 + capmap.cn*M_PI/2);
a[1] = (M_PI/-2.0 + capmap.cn*0);
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
vector_sub(v, c, v_min_c);
dot_product(tmpRot, v_min_c, ret_dot);
vector_add(ret_dot, a, vector);
}
xy.x = vector[0];
xy.y = vector[1];
return xy;
}

vector_sub(v, c, v_min_c);
dot_product(tmpRot, v_min_c, ret_dot);
vector_add(ret_dot, a, vector);

xy.x = vector[0];
xy.y = vector[1];
return xy;
}


Expand Down Expand Up @@ -621,9 +625,9 @@ PJ *PROJECTION(healpix) {
P->opaque = Q;

if (P->es) {
Q->apa = pj_authset(P->es); /* For auth_lat(). */
Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
Q->apa = pj_authset(P->es); /* For auth_lat(). */
Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
P->ra = 1.0/P->a;
P->fwd = e_healpix_forward;
P->inv = e_healpix_inverse;
Expand Down
Loading

0 comments on commit 196a1bb

Please sign in to comment.