--- libwcs/proj.c | 1365 +++++++++++++++++++++++++++++++++++++++++++++++++++++- libwcs/wcs.c | 17 +- libwcs/wcs.h | 19 +- libwcs/wcsinit.c | 8 + libwcs/wcslib.h | 11 +- 5 files changed, 1401 insertions(+), 19 deletions(-) diff --git a/libwcs/proj.c b/libwcs/proj.c index 8460329..87d7c53 100644 --- a/libwcs/proj.c +++ b/libwcs/proj.c @@ -1,3 +1,5 @@ +#include + /*============================================================================ * * WCSLIB - an implementation of the FITS WCS proposal. @@ -67,6 +69,9 @@ * tscset tscfwd tscrev TSC: tangential spherical cube * cscset cscfwd cscrev CSC: COBE quadrilateralized spherical cube * qscset qscfwd qscrev QSC: quadrilateralized spherical cube +* hpxset hpxfwd hpxrev HPX: HEALPix projection +* xpfset xpffwd xpfrev XPH: HEALPix polar, aka "butterfly" +* toaset toafwd toarev TOA: TOAST projection * * * Driver routines; prjset(), prjfwd() & prjrev() @@ -237,11 +242,11 @@ #include #include "wcslib.h" -int npcode = 26; -char pcodes[26][4] = +int npcode = 29; +char pcodes[29][4] = {"AZP", "SZP", "TAN", "STG", "SIN", "ARC", "ZPN", "ZEA", "AIR", "CYP", "CEA", "CAR", "MER", "COP", "COE", "COD", "COO", "SFL", "PAR", "MOL", - "AIT", "BON", "PCO", "TSC", "CSC", "QSC"}; + "AIT", "BON", "PCO", "TSC", "CSC", "QSC", "HPX", "XPH", "TOA"}; const int AZP = 101; const int SZP = 102; @@ -269,6 +274,9 @@ const int PCO = 602; const int TSC = 701; const int CSC = 702; const int QSC = 703; +const int HPX = 801; +const int XPH = 802; +const int TOA = 803; /* Map error number to error message for each function. */ const char *prjset_errmsg[] = { @@ -288,6 +296,31 @@ const char *prjrev_errmsg[] = { #define copysgn(X, Y) ((Y) < 0.0 ? -fabs(X) : fabs(X)) #define copysgni(X, Y) ((Y) < 0 ? -abs(X) : abs(X)) + +/* Vector functions for TOAST partitioning */ + +const double deg2rad = PI / 180.0; + +typedef struct vec +{ + double lon, lat; + double x, y, z; +} +Vec; + +void vCopy (Vec *v, Vec *c); +void vCalcRADec(Vec *v); +void vCalcXYZ (Vec *v); +void vMidpoint (Vec *a, Vec *b, Vec *c); +void vPixCenter(Vec *a, Vec *b, Vec *c, Vec *d, Vec *v); +int vCross (Vec *a, Vec *b, Vec *c); +double vDot (Vec *a, Vec *b); +double vNormalize(Vec *a); +double vPrint (Vec *v, char *lbl); + +void splitIndex(unsigned long index, int level, int *x, int *y); + + /*==========================================================================*/ int prjset(pcode, prj) @@ -349,6 +382,12 @@ struct prjprm *prj; cscset(prj); } else if (strcmp(pcode, "QSC") == 0) { qscset(prj); + } else if (strcmp(pcode, "HPX") == 0) { + hpxset(prj); + } else if (strcmp(pcode, "XPH") == 0) { + xphset(prj); + } else if (strcmp(pcode, "TOA") == 0) { + toaset(prj); } else { /* Unrecognized projection code. */ return 1; @@ -4364,6 +4403,1100 @@ double *phi, *theta; return 0; } +/*============================================================================ +* HPX: HEALPix projection. +* +* Given: +* prj->p[1] H - the number of facets in longitude. +* prj->p[2] K - the number of facets in latitude +* +* Given and/or returned: +* prj->r0 Reset to 180/pi if 0. +* prj->phi0 Reset to 0.0 +* prj->theta0 Reset to 0.0 +* +* Returned: +* prj->flag HPX +* prj->code "HPX" +* prj->n True if K is odd. +* prj->w[0] r0*(pi/180) +* prj->w[1] (180/pi)/r0 +* prj->w[2] (K-1)/K +* prj->w[3] 90*K/H +* prj->w[4] (K+1)/2 +* prj->w[5] 90*(K-1)/H +* prj->w[6] 180/H +* prj->w[7] H/360 +* prj->w[8] (90*K/H)*r0*(pi/180) +* prj->w[9] (180/H)*r0*(pi/180) +* prj->prjfwd Pointer to hpxfwd(). +* prj->prjrev Pointer to hpxrev(). + + +*===========================================================================*/ + +int hpxset(prj) + +struct prjprm *prj; + +{ + strcpy(prj->code, "HPX"); + prj->flag = HPX; + prj->phi0 = 0.0; + prj->theta0 = 0.0; + + prj->n = ((int)prj->p[2])%2; + + if (prj->r0 == 0.0) { + prj->r0 = R2D; + prj->w[0] = 1.0; + prj->w[1] = 1.0; + } else { + prj->w[0] = prj->r0*D2R; + prj->w[1] = R2D/prj->r0; + } + + prj->w[2] = (prj->p[2] - 1.0) / prj->p[2]; + prj->w[3] = 90.0 * prj->p[2] / prj->p[1]; + prj->w[4] = (prj->p[2] + 1.0) / 2.0; + prj->w[5] = 90.0 * (prj->p[2] - 1.0) / prj->p[1]; + prj->w[6] = 180.0 / prj->p[1]; + prj->w[7] = prj->p[1] / 360.0; + prj->w[8] = prj->w[3] * prj->w[0]; + prj->w[9] = prj->w[6] * prj->w[0]; + + prj->prjfwd = hpxfwd; + prj->prjrev = hpxrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int hpxfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct prjprm *prj; +double *x, *y; + +{ + double abssin, sigma, sinthe, phic; + int hodd; + + if( prj->flag != HPX ) { + if( hpxset( prj ) ) return 1; + } + + sinthe = sindeg( theta ); + abssin = fabs( sinthe ); + +/* Equatorial zone */ + if( abssin <= prj->w[2] ) { + *x = prj->w[0] * phi; + *y = prj->w[8] * sinthe; + +/* Polar zone */ + } else { + +/* DSB - The expression for phic is conditioned differently to the + WCSLIB code in order to improve accuracy of the floor function for + arguments very slightly below an integer value. */ + hodd = ((int)prj->p[1]) % 2; + if( !prj->n && theta <= 0.0 ) hodd = 1 - hodd; + if( hodd ) { + phic = -180.0 + (2.0*floor( prj->w[7] * phi + 1/2 ) + prj->p[1] ) * prj->w[6]; + } else { + phic = -180.0 + (2.0*floor( prj->w[7] * phi ) + prj->p[1] + 1 ) * prj->w[6]; + } + + sigma = sqrt( prj->p[2]*( 1.0 - abssin )); + + *x = prj->w[0] *( phic + ( phi - phic )*sigma ); + + *y = prj->w[9] * ( prj->w[4] - sigma ); + if( theta < 0 ) *y = -*y; + + } + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int hpxrev(x, y, prj, phi, theta) + +const double x, y; +struct prjprm *prj; +double *phi, *theta; + +{ + double absy, sigma, t, yr, xc; + int hodd; + + if (prj->flag != HPX) { + if (hpxset(prj)) return 1; + } + + yr = prj->w[1]*y; + absy = fabs( yr ); + +/* Equatorial zone */ + if( absy <= prj->w[5] ) { + *phi = prj->w[1] * x; + t = yr/prj->w[3]; + if( t < -1.0 || t > 1.0 ) { + return 2; + } else { + *theta = asindeg( t ); + } + +/* Polar zone */ + } else if( absy <= 90 ){ + + hodd = ((int)prj->p[1]) % 2; + if( !prj->n && yr <= 0.0 ) hodd = 1 - hodd; + if( hodd ) { + xc = -180.0 + (2.0*floor( prj->w[7] * x + 1/2 ) + prj->p[1] ) * prj->w[6]; + } else { + xc = -180.0 + (2.0*floor( prj->w[7] * x ) + prj->p[1] + 1 ) * prj->w[6]; + } + + sigma = prj->w[4] - absy / prj->w[6]; + + if( sigma == 0.0 ) { + return 2; + } else { + + t = ( x - xc )/sigma; + if( fabs( t ) <= prj->w[6] ) { + *phi = prj->w[1] *( xc + t ); + } else { + return 2; + } + } + + t = 1.0 - sigma*sigma/prj->p[2]; + if( t < -1.0 || t > 1.0 ) { + return 2; + } else { + *theta = asindeg ( t ); + if( y < 0 ) *theta = -*theta; + } + + } else { + return 2; + } + + return 0; +} + +/*============================================================================ +* XPH: HEALPix polar, aka "butterfly" projection. +* +* Given and/or returned: +* prj->r0 Reset to 180/pi if 0. +* prj->phi0 Reset to 0.0 if undefined. +* prj->theta0 Reset to 0.0 if undefined. +* +* Returned: +* prj->flag XPH +* prj->code "XPH" +* prj->w[0] r0*(pi/180)/sqrt(2) +* prj->w[1] (180/pi)/r0/sqrt(2) +* prj->w[2] 2/3 +* prj->w[3] tol (= 1e-4) +* prj->w[4] sqrt(2/3)*(180/pi) +* prj->w[5] 90 - tol*sqrt(2/3)*(180/pi) +* prj->w[6] sqrt(3/2)*(pi/180) +* prj->prjfwd Pointer to xphfwd(). +* prj->prjrev Pointer to xphrev(). +*===========================================================================*/ + +int xphset(prj) + +struct prjprm *prj; + +{ + strcpy(prj->code, "XPH"); + prj->flag = XPH; + + if (prj->r0 == 0.0) { + prj->r0 = R2D; + prj->w[0] = 1.0; + prj->w[1] = 1.0; + } else { + prj->w[0] = prj->r0*D2R; + prj->w[1] = R2D/prj->r0; + } + + prj->w[0] /= sqrt(2.0); + prj->w[1] /= sqrt(2.0); + prj->w[2] = 2.0/3.0; + prj->w[3] = 1e-4; + prj->w[4] = sqrt(prj->w[2])*R2D; + prj->w[5] = 90.0 - prj->w[3]*prj->w[4]; + prj->w[6] = sqrt(1.5)*D2R; + + prj->prjfwd = xphfwd; + prj->prjrev = xphrev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int xphfwd(phi, theta, prj, x, y) + +const double phi, theta; +struct prjprm *prj; +double *x, *y; + +{ + double abssin, chi, eta, psi, sigma, sinthe, xi; + + if (prj->flag != XPH) { + if (xphset(prj)) return 1; + } + + /* Do phi dependence. */ + chi = phi; + if (180.0 <= fabs(chi)) { + chi = fmod(chi, 360.0); + if (chi < -180.0) { + chi += 360.0; + } else if (180.0 <= chi) { + chi -= 360.0; + } + } + + /* phi is also recomputed from chi to avoid rounding problems. */ + chi += 180.0; + psi = fmod(chi, 90.0); + + /* y is used to hold phi (rounded). */ + *x = psi; + *y = chi - 180.0; + + /* Do theta dependence. */ + sinthe = sindeg(theta); + abssin = fabs(sinthe); + + if (abssin <= prj->w[2]) { + /* Equatorial regime. */ + xi = *x; + eta = 67.5 * sinthe; + + } else { + /* Polar regime. */ + if (theta < prj->w[5]) { + sigma = sqrt(3.0*(1.0 - abssin)); + } else { + sigma = (90.0 - theta)*prj->w[6]; + } + + xi = 45.0 + (*x - 45.0)*sigma; + eta = 45.0 * (2.0 - sigma); + if (theta < 0.0) eta = -eta; + } + + xi -= 45.0; + eta -= 90.0; + + /* Recall that y holds phi. */ + if (*y < -90.0) { + *x = prj->w[0]*(-xi + eta); + *y = prj->w[0]*(-xi - eta); + + } else if (*y < 0.0) { + *x = prj->w[0]*(+xi + eta); + *y = prj->w[0]*(-xi + eta); + + } else if (*y < 90.0) { + *x = prj->w[0]*( xi - eta); + *y = prj->w[0]*( xi + eta); + + } else { + *x = prj->w[0]*(-xi - eta); + *y = prj->w[0]*( xi - eta); + } + + return 0; + +} + +/*--------------------------------------------------------------------------*/ + +int xphrev(x, y, prj, phi, theta) + +const double x, y; +struct prjprm *prj; +double *phi, *theta; + +{ + double abseta, eta, eta1, sigma, xi, xi1, xr, yr; + const double tol = 1.0e-12; + + if (prj->flag != XPH) { + if (xphset(prj)) return 1; + } + + + xr = x*prj->w[1]; + yr = y*prj->w[1]; + if (xr <= 0.0 && 0.0 < yr) { + xi1 = -xr - yr; + eta1 = xr - yr; + *phi = -180.0; + } else if (xr < 0.0 && yr <= 0.0) { + xi1 = xr - yr; + eta1 = xr + yr; + *phi = -90.0; + } else if (0.0 <= xr && yr < 0.0) { + xi1 = xr + yr; + eta1 = -xr + yr; + *phi = 0.0; + } else { + xi1 = -xr + yr; + eta1 = -xr - yr; + *phi = 90.0; + } + + xi = xi1 + 45.0; + eta = eta1 + 90.0; + abseta = fabs(eta); + + if (abseta <= 90.0) { + if (abseta <= 45.0) { + /* Equatorial regime. */ + *phi += xi; + *theta = asindeg (eta/67.5); + + /* Bounds checking. */ + if (45.0+tol < fabs(xi1)) return 2; + + } else { + /* Polar regime. */ + sigma = (90.0 - abseta) / 45.0; + + /* Ensure an exact result for points on the boundary. */ + if (xr == 0.0) { + if (yr <= 0.0) { + *phi = 0.0; + } else { + *phi = 180.0; + } + } else if (yr == 0.0) { + if (xr < 0.0) { + *phi = -90.0; + } else { + *phi = 90.0; + } + } else { + *phi += 45.0 + xi1/sigma; + } + + if (sigma < prj->w[3]) { + *theta = 90.0 - sigma*prj->w[4]; + } else { + *theta = asindeg (1.0 - sigma*sigma/3.0); + } + if (eta < 0.0) *theta = -(*theta); + + /* Bounds checking. */ + if (eta < -45.0 && eta+90.0+tol < fabs(xi1)) return 2; + } + + } else { + /* Beyond latitude range. */ + *phi = 0.0; + *theta = 0.0; + return 2; + } + + return 0; +} + + +/*============================================================================ +* TOA: TOAST projection. +* +* Given: +* prj->p[1] The HTM level (number of iterative subdivisions of space) +* +* Returned: +* prj->flag TOA +* prj->code "TOA" +* prj->prjfwd Pointer to toafwd(). +* prj->prjrev Pointer to toarev(). + + +*===========================================================================*/ + +int toaset(prj) + +struct prjprm *prj; + +{ + strcpy(prj->code, "TOA"); + prj->flag = TOA; + + prj->prjfwd = toafwd; + prj->prjrev = toarev; + + return 0; +} + +/*--------------------------------------------------------------------------*/ + +int toafwd(phi, theta, prj, x, y) + +const double phi, theta; +struct prjprm *prj; +double *x, *y; +{ + Vec ref; + Vec center; + Vec corner [4]; + Vec midpoint[4]; + Vec normal [4]; + + unsigned long index, maxindex; + + int debug = 0; + + int i, level, maxlevel, npix; + int xindex, yindex; + int prime, opposite; + + double direction[4]; + double size; + + + if(debug > 1) + printf("\nTOAFWD> LON = %10.6f LAT = %10.6f\n\n", phi, theta); + + + if( prj->flag != TOA ) + if( toaset( prj ) ) + return 1; + + + // We calculate down to level 27 to allow + // for fraction pixel coordinates + + maxlevel = 27; + + maxindex = 0b1 << maxlevel; + + + // Get location of interest from command line + + ref.lon = phi; + ref.lat = theta; + + while(ref.lon < 0.) ref.lon += 360.; + while(ref.lon >= 360.) ref.lon -= 360.; + + if(ref.lat > 90.) ref.lat = 90.; + if(ref.lat < -90.) ref.lat = -90.; + + vCalcXYZ(&ref); + + + // Get the level of the map + // This determines how many pixels there are in the map + + level = (int)prj->p[1]; + + npix = 0b1 << (level+8); + + + // The first level has to be done by hand + + if(ref.lon < 90.) // quadrant 11 [3] + { + if(debug > 1) + printf("TOAFWD> Quadrant 11\n\n"); + + index = 0b11; + + corner[0].lon = 0.; + corner[0].lat = 90.; + + corner[1].lon = 0.; + corner[1].lat = 0.; + + corner[2].lon = 0.; + corner[2].lat = -90.; + + corner[3].lon = 90.; + corner[3].lat = 0.; + + prime = 1; + } + + else if(ref.lon < 180.) // quadrant 10 [2] + { + if(debug > 1) + printf("TOAFWD> Quadrant 10\n\n"); + + index = 0b10; + + corner[0].lon = 180.; + corner[0].lat = 0.; + + corner[1].lon = 0.; + corner[1].lat = 90.; + + corner[2].lon = 90.; + corner[2].lat = 0.; + + corner[3].lon = 0.; + corner[3].lat = -90.; + + prime = 0; + } + + else if(ref.lon < 270.) // quadrant 00 [0] + { + if(debug >1) + printf("TOAFWD> Quadrant 00\n\n"); + + index = 0b00; + + corner[0].lon = 0.; + corner[0].lat = -90.; + + corner[1].lon = 270.; + corner[1].lat = 0.; + + corner[2].lon = 0.; + corner[2].lat = 90.; + + corner[3].lon = 180.; + corner[3].lat = 0.; + + prime = 1; + } + + else // quadrant 01 [1] + { + if(debug > 1) + printf("TOAFWD> Quadrant 01\n\n"); + + index = 0b01; + + corner[0].lon = 270.; + corner[0].lat = 0.; + + corner[1].lon = 0.; + corner[1].lat = -90.; + + corner[2].lon = 0.; + corner[2].lat = 0.; + + corner[3].lon = 0.; + corner[3].lat = 90.; + + prime = 0; + } + + level = 1; + + vCalcXYZ(&corner[0]); + vCalcXYZ(&corner[1]); + vCalcXYZ(&corner[2]); + vCalcXYZ(&corner[3]); + + if(debug > 1) + { + vPrint(&corner[0], "corner[0]"); + vPrint(&corner[1], "corner[1]"); + vPrint(&corner[2], "corner[2]"); + vPrint(&corner[3], "corner[3]"); + } + + + // Drill down level by level + + while(level < maxlevel) + { + // Find cell edge midpoints + + vMidpoint(&corner[1], &corner[0], &midpoint[0]); + vMidpoint(&corner[2], &corner[1], &midpoint[1]); + vMidpoint(&corner[3], &corner[2], &midpoint[2]); + vMidpoint(&corner[0], &corner[3], &midpoint[3]); + + if(debug > 1) + { + vPrint(&midpoint[0], "midpoint[0]"); + vPrint(&midpoint[1], "midpoint[1]"); + vPrint(&midpoint[2], "midpoint[2]"); + vPrint(&midpoint[3], "midpoint[3]"); + } + + + // We also need the center point (midpoint of the HTM diagonal) + + opposite = (prime + 2) % 4; + + vMidpoint(&corner[prime], &corner[opposite], ¢er); + + if(debug > 1) + { + vPrint(&corner[prime], "Prime"); + vPrint(&corner[opposite], "Opposite"); + vPrint(¢er, "center"); + } + + + // Find the lines connecting the center and the edge midpoints + // (actually the normal vectors to those planes) + + vCross(&midpoint[0], ¢er, &normal[0]); + vNormalize(&normal[0]); + + vCross(¢er, &midpoint[2], &normal[2]); + vNormalize(&normal[0]); + + vCross(&midpoint[1], ¢er, &normal[1]); + vNormalize(&normal[0]); + + vCross(¢er, &midpoint[3], &normal[3]); + vNormalize(&normal[0]); + + + // The dot product of these with the point of interest vector + // tells us which half of the cell that point is in + + if(debug > 1) + { + vPrint(&ref, "ref"); + vPrint(&normal[0], "normal[0]"); + vPrint(&normal[1], "normal[1]"); + vPrint(&normal[2], "normal[2]"); + vPrint(&normal[3], "normal[3]"); + } + + direction[0] = vDot(&normal[0], &ref); + direction[1] = vDot(&normal[1], &ref); + direction[2] = vDot(&normal[2], &ref); + direction[3] = vDot(&normal[3], &ref); + + if(debug > 1) + printf("DIRECTIONS> %.4f %.4f %.4d %.4f\n", direction[0], direction[1], direction[2], direction[3]); + + + // Use this to define the next level down + + if(direction[0] >= 0. && direction[3] >= 0.) // Upper left (00) + { + if(debug > 1) + printf("TOAFWD> Level %2d, Subquadrant 00\n", level); + + index = (index << 2) + 0b00; + + if(prime == 0) prime = 0; + else if(opposite == 0) prime = 0; + else prime = 1; + + //-- + vCopy(&midpoint[0], &corner[1]); + vCopy( ¢er, &corner[2]); + vCopy(&midpoint[3], &corner[3]); + } + + if(direction[0] < 0. && direction[1] >= 0.) // Upper right (01) + { + if(debug > 1) + printf("TOAFWD> Level %2d, Subquadrant 01\n", level); + + index = (index << 2) + 0b01; + + if(prime == 1) prime = 1; + else if(opposite == 1) prime = 1; + else prime = 0; + + vCopy(&midpoint[0], &corner[0]); + //-- + vCopy(&midpoint[1], &corner[2]); + vCopy( ¢er, &corner[3]); + } + + if(direction[2] >= 0. && direction[3] < 0.) // Lower left (10) + { + if(debug > 1) + printf("TOAFWD> Level %2d, Subquadrant 10\n", level); + + index = (index << 2) + 0b10; + + if(prime == 3) prime = 3; + else if(opposite == 3) prime = 3; + else prime = 0; + + vCopy(&midpoint[3], &corner[0]); + vCopy( ¢er, &corner[1]); + vCopy(&midpoint[2], &corner[2]); + //-- + } + + if(direction[2] < 0. && direction[1] < 0.) // Lower right (11) + { + if(debug > 1) + printf("TOAFWD> Level %2d: Subquadrant 11\n", level); + + index = (index << 2) + 0b11; + + if(prime == 2) prime = 2; + else if(opposite == 2) prime = 2; + else prime = 1; + + vCopy( ¢er, &corner[0]); + vCopy(&midpoint[1], &corner[1]); + //-- + vCopy(&midpoint[2], &corner[3]); + } + + if(debug > 1) + { + vPrint(&corner[0], "corner[0]"); + vPrint(&corner[1], "corner[1]"); + vPrint(&corner[2], "corner[2]"); + vPrint(&corner[3], "corner[3]"); + + printf("prime = %d\n", prime); + } + + ++level; + } + + + splitIndex(index, level, &xindex, &yindex); + + *x = (double)xindex * (double)npix / (double)maxindex; + *y = (double)yindex * (double)npix / (double)maxindex; + + if(debug) + { + printf("TOAFWD> (lon, lat):(%10.6f, %10.6f) -> (xpix, ypix):(%9.4f, %9.4f) fraction of image:[%.3f,%.3f]\n", + ref.lon, ref.lat, *x+0.5, *y+0.5, *x/npix, *y/npix); + fflush(stdout); + } + + + // FITS pixel coordinates start at 0.5 + // (so the center of the first pixel is at 1.0) + + *x += 1.0; + *y += 1.0; + + return(0); +} + +/*--------------------------------------------------------------------------*/ + +int toarev(x, y, prj, phi, theta) + +const double x, y; +struct prjprm *prj; +double *phi, *theta; +{ + Vec ref; + Vec center; + Vec corner [4]; + Vec midpoint[4]; + + unsigned long index, maxindex, currindex; + + int debug = 0; + + int i, npix, level, maxlevel; + int xindex, yindex; + int xsplit, ysplit; + int prime, opposite; + + double size; + + if(debug > 1) + printf("\nTOAREV> X = %9.4f Y = %9.4f\n\n", x, y); + + if( prj->flag != TOA ) + if( toaset( prj ) ) + return 1; + + + // We calculate down to level 27 to allow + // for fraction pixel coordinates + + maxlevel = 27; + + maxindex = 0b1 << maxlevel; + + + // Get the level of the map + // This determines how many pixels there are in the map + + level = (int)prj->p[1]; + + npix = 0b1 << (level + 8); + + if(x < 0) return 2; + if(x > npix+1) return 2; + if(y < 0) return 2; + if(y > npix+1) return 2; + + xindex = (x-1.0) / npix * maxindex + 0.5; + yindex = (y-1.0) / npix * maxindex + 0.5; + + if(debug > 1) + { + printf("TOAREV> xindex = %o\n", xindex); + printf("TOAREV> yindex = %o\n\n", yindex); + } + + + currindex = maxlevel; + + xsplit = xindex >> (currindex-1) & 0b1; + ysplit = yindex >> (currindex-1) & 0b1; + + + // The first level has to be done by hand + + if(xsplit == 0 && ysplit == 0) // quadrant 00 [0] + { + if(debug > 1) + printf("TOAREV> Quadrant 00\n\n"); + + corner[0].lon = 0.; + corner[0].lat = -90.; + + corner[1].lon = 270.; + corner[1].lat = 0.; + + corner[2].lon = 0.; + corner[2].lat = 90.; + + corner[3].lon = 180.; + corner[3].lat = 0.; + + prime = 1; + } + + else if(xsplit == 1 && ysplit == 0) // quadrant 01 [1] + { + if(debug > 1) + printf("TOAREV> Quadrant 01\n\n"); + + corner[0].lon = 270.; + corner[0].lat = 0.; + + corner[1].lon = 0.; + corner[1].lat = -90.; + + corner[2].lon = 0.; + corner[2].lat = 0.; + + corner[3].lon = 0.; + corner[3].lat = 90.; + + prime = 0; + } + + else if(xsplit == 0 && ysplit == 1) // quadrant 10 [2] + { + if(debug > 1) + printf("TOAREV> Quadrant 10\n\n"); + + corner[0].lon = 180.; + corner[0].lat = 0.; + + corner[1].lon = 0.; + corner[1].lat = 90.; + + corner[2].lon = 90.; + corner[2].lat = 0.; + + corner[3].lon = 0.; + corner[3].lat = -90.; + + prime = 0; + } + + else // quadrant 11 [3] + { + if(debug > 1) + printf("TOAREV> Quadrant 11\n\n"); + + corner[0].lon = 0.; + corner[0].lat = 90.; + + corner[1].lon = 0.; + corner[1].lat = 0.; + + corner[2].lon = 0.; + corner[2].lat = -90.; + + corner[3].lon = 90.; + corner[3].lat = 0.; + + prime = 1; + } + + level = 1; + + vCalcXYZ(&corner[0]); + vCalcXYZ(&corner[1]); + vCalcXYZ(&corner[2]); + vCalcXYZ(&corner[3]); + + if(debug > 1) + { + vPrint(&corner[0], "corner[0]"); + vPrint(&corner[1], "corner[1]"); + vPrint(&corner[2], "corner[2]"); + vPrint(&corner[3], "corner[3]"); + } + + // Drill down level by level + + while(level < maxlevel) + { + // Find cell edge midpoints + + vMidpoint(&corner[1], &corner[0], &midpoint[0]); + vMidpoint(&corner[2], &corner[1], &midpoint[1]); + vMidpoint(&corner[3], &corner[2], &midpoint[2]); + vMidpoint(&corner[0], &corner[3], &midpoint[3]); + + if(debug > 1) + { + vPrint(&midpoint[0], "midpoint[0]"); + vPrint(&midpoint[1], "midpoint[1]"); + vPrint(&midpoint[2], "midpoint[2]"); + vPrint(&midpoint[3], "midpoint[3]"); + } + + + // We also need the center point (midpoint of the HTM diagonal) + + opposite = (prime + 2) % 4; + + vMidpoint(&corner[prime], &corner[opposite], ¢er); + + if(debug > 1) + { + vPrint(&corner[prime], "Prime"); + vPrint(&corner[opposite], "Opposite"); + vPrint(¢er, "center"); + } + + + // The next level up in xindex and yindex tells + // us which half of the cell that point is in + + currindex -= 1; + + xsplit = xindex >> (currindex-1) & 0b1; + ysplit = yindex >> (currindex-1) & 0b1; + + + // Use this to define the next level down + + if(ysplit == 0 && xsplit == 0) // Upper left (00) + { + if(debug > 1) + printf("TOAREV> Level %2d, Subquadrant 00\n", level); + + if(prime == 0) prime = 0; + else if(opposite == 0) prime = 0; + else prime = 1; + + //-- + vCopy(&midpoint[0], &corner[1]); + vCopy( ¢er, &corner[2]); + vCopy(&midpoint[3], &corner[3]); + } + + else if(ysplit == 0 && xsplit == 1) // Upper right (01) + { + if(debug > 1) + printf("TOAREV> Level %2d, Subquadrant 01\n", level); + + if(prime == 1) prime = 1; + else if(opposite == 1) prime = 1; + else prime = 0; + + vCopy(&midpoint[0], &corner[0]); + //-- + vCopy(&midpoint[1], &corner[2]); + vCopy( ¢er, &corner[3]); + } + + else if(ysplit == 1 && xsplit == 0) // Lower left (10) + { + if(debug > 1) + printf("TOAREV> Level %2d, Subquadrant 10\n", level); + + if(prime == 3) prime = 3; + else if(opposite == 3) prime = 3; + else prime = 0; + + vCopy(&midpoint[3], &corner[0]); + vCopy( ¢er, &corner[1]); + vCopy(&midpoint[2], &corner[2]); + //-- + } + + else // Lower right (11) + { + if(debug > 1) + printf("TOAREV> Level %2d, Subquadrant 11\n", level); + + if(prime == 2) prime = 2; + else if(opposite == 2) prime = 2; + else prime = 1; + + vCopy( ¢er, &corner[0]); + vCopy(&midpoint[1], &corner[1]); + //-- + vCopy(&midpoint[2], &corner[3]); + } + + if(debug > 1) + { + vPrint(&corner[0], "corner[0]"); + vPrint(&corner[1], "corner[1]"); + vPrint(&corner[2], "corner[2]"); + vPrint(&corner[3], "corner[3]"); + + printf("prime = %d\n", prime); + } + + ++level; + } + + + vPixCenter(&corner[0], &corner[1], &corner[2], &corner[3], &ref); + + vCalcRADec(&ref); + + *phi = ref.lon; + *theta = ref.lat; + + if(debug) + { + printf("TOAREV> (lon, lat):(%10.6f, %10.6f) <- (xpix, ypix):(%9.4f, %9.4f)\n", + *phi, *theta, x, y); + fflush(stdout); + } + + return(0); +} + + + /* This routine comes from E. Bertin sextractor-2.8.6 */ int @@ -4506,6 +5639,230 @@ poly_end: return 0; } +/***************************************************/ +/* */ +/* vCopy() */ +/* */ +/* Copy the contents of one vector to another */ +/* */ +/***************************************************/ + +void vCopy(Vec *v, Vec *c) +{ + c->lon = v->lon; + c->lat = v->lat; + + c->x = v->x; + c->y = v->y; + c->z = v->z; + + return; +} + + + +/***************************************************/ +/* */ +/* vCalcRADec() */ +/* */ +/* Update vector with RA and Dec based on x,y,z */ +/* */ +/***************************************************/ + +void vCalcRADec(Vec *v) +{ + v->lon = atan2(v->y, v->x)/deg2rad; + v->lat = asin(v->z)/deg2rad; + + while(v->lon >= 360.) v->lon -= 360.; + while(v->lon < 0.) v->lon += 360.; + + return; +} + + + +/***************************************************/ +/* */ +/* vCalcXYZ() */ +/* */ +/* Update vector with x,y,z based on RA,Dec */ +/* */ +/***************************************************/ + +void vCalcXYZ(Vec *v) +{ + v->x = cos(v->lat * deg2rad) * cos(v->lon * deg2rad); + v->y = cos(v->lat * deg2rad) * sin(v->lon * deg2rad); + v->z = sin(v->lat * deg2rad); + + return; +} + + + +/***************************************************/ +/* */ +/* vMidpoint() */ +/* */ +/* Finds the midpoint between two points on the */ +/* sky. */ +/* */ +/***************************************************/ + +void vMidpoint(Vec *a, Vec *b, Vec *c) +{ + c->x = a->x + b->x; + c->y = a->y + b->y; + c->z = a->z + b->z; + + vNormalize(c); + + return; +} + + + +/***************************************************/ +/* */ +/* vPixCenter() */ +/* */ +/* Finds the center of a a pixel (four corners) */ +/* on the sky. */ +/* */ +/***************************************************/ + +void vPixCenter(Vec *a, Vec *b, Vec *c, Vec *d, Vec *v) +{ + v->x = a->x + b->x + c->x + d->x; + v->y = a->y + b->y + c->y + d->y; + v->z = a->z + b->z + c->z + d->z; + + vNormalize(v); + + return; +} + + + +/***************************************************/ +/* */ +/* vCross() */ +/* */ +/* Vector cross product. */ +/* */ +/***************************************************/ + +int vCross(Vec *v1, Vec *v2, Vec *v3) +{ + v3->x = v1->y*v2->z - v2->y*v1->z; + v3->y = -v1->x*v2->z + v2->x*v1->z; + v3->z = v1->x*v2->y - v2->x*v1->y; + + if(v3->x == 0. + && v3->y == 0. + && v3->z == 0.) + return 0; + + return 1; +} + + +/***************************************************/ +/* */ +/* vDot() */ +/* */ +/* Vector dot product. */ +/* */ +/***************************************************/ + +double vDot(Vec *a, Vec *b) +{ + double sum; + + sum = a->x * b->x + + a->y * b->y + + a->z * b->z; + + return sum; +} + + +/***************************************************/ +/* */ +/* vNormalize() */ +/* */ +/* Normalize the vector */ +/* */ +/***************************************************/ + +double vNormalize(Vec *v) +{ + double len; + + len = 0.; + + len = sqrt(v->x * v->x + v->y * v->y + v->z * v->z); + + if(len == 0.) + len = 1.; + + v->x = v->x / len; + v->y = v->y / len; + v->z = v->z / len; + + return len; +} + + +/***************************************************/ +/* */ +/* vPrint() */ +/* */ +/* Print out vector (for debugging) */ +/* */ +/***************************************************/ + +double vPrint(Vec *v, char *label) +{ + vCalcRADec(v); + + printf("VECTOR> %9.6f %9.6f %9.6f -> %10.6f %10.6f (%s)\n", + v->x, v->y, v->z, v->lon, v->lat, label); + fflush(stdout); +} + + +/***************************************************/ +/* */ +/* splitIndex() */ +/* */ +/* Cell indices are Z-order binary constructs. */ +/* The x and y pixel offsets are constructed by */ +/* extracting the pattern made by every other bit */ +/* to make new binary numbers. */ +/* */ +/***************************************************/ + +void splitIndex(unsigned long index, int level, int *x, int *y) +{ + int i; + unsigned long val; + + val = index; + + *x = 0; + *y = 0; + + for(i=0; i> (2*i)) & 0b1) << i); + *y = *y + (((val >> (2*i+1)) & 0b1) << i); + } + + return; +} + /* Dec 20 1999 Doug Mink - Change cosd() and sind() to cosdeg() and sindeg() * Dec 20 1999 Doug Mink - Include wcslib.h, which includes proj.h, wcsmath.h * Dec 20 1999 Doug Mink - Define copysign only if it is not defined @@ -4524,4 +5881,6 @@ poly_end: * * Mar 14 2011 Doug Mink - If no coefficients in ZPN, make ARC * Mar 14 2011 Doug Mink - Add Emmanuel Bertin's TAN polynomial from Ed Los + * Dec 22 2016 John Good - Support for HEALPix and TOAST */ + diff --git a/libwcs/wcs.c b/libwcs/wcs.c index 84cfb16..2d43d69 100644 --- a/libwcs/wcs.c +++ b/libwcs/wcs.c @@ -375,13 +375,16 @@ char *ctype2; /* FITS WCS projection for axis 2 */ strcpy (ctypes[24], "CSC"); strcpy (ctypes[25], "QSC"); strcpy (ctypes[26], "TSC"); - strcpy (ctypes[27], "NCP"); - strcpy (ctypes[28], "GLS"); - strcpy (ctypes[29], "DSS"); - strcpy (ctypes[30], "PLT"); - strcpy (ctypes[31], "TNX"); - strcpy (ctypes[32], "ZPX"); - strcpy (ctypes[33], "TPV"); + strcpy (ctypes[27], "HPX"); + strcpy (ctypes[28], "XPH"); + strcpy (ctypes[29], "NCP"); + strcpy (ctypes[30], "GLS"); + strcpy (ctypes[31], "DSS"); + strcpy (ctypes[32], "PLT"); + strcpy (ctypes[33], "TNX"); + strcpy (ctypes[34], "ZPX"); + strcpy (ctypes[35], "TPV"); + strcpy (ctypes[36], "TOA"); /* Initialize distortion types */ strcpy (dtypes[1], "SIP"); diff --git a/libwcs/wcs.h b/libwcs/wcs.h index 0d69049..3ef1695 100644 --- a/libwcs/wcs.h +++ b/libwcs/wcs.h @@ -192,14 +192,17 @@ struct WorldCoor { #define WCS_CSC 24 /* COBE quadrilateralized Spherical Cube */ #define WCS_QSC 25 /* Quadrilateralized Spherical Cube */ #define WCS_TSC 26 /* Tangential Spherical Cube */ -#define WCS_NCP 27 /* Special case of SIN from AIPS*/ -#define WCS_GLS 28 /* Same as SFL from AIPS*/ -#define WCS_DSS 29 /* Digitized Sky Survey plate solution */ -#define WCS_PLT 30 /* Plate fit polynomials (SAO) */ -#define WCS_TNX 31 /* Tangent Plane (NOAO corrections) */ -#define WCS_ZPX 32 /* Zenithal Azimuthal Polynomial (NOAO corrections) */ -#define WCS_TPV 33 /* Tangent Plane (SCAMP corrections) */ -#define NWCSTYPE 34 /* Number of WCS types (-1 really means no WCS) */ +#define WCS_HPX 27 /* Tangential Spherical Cube */ +#define WCS_XPH 28 /* Tangential Spherical Cube */ +#define WCS_NCP 29 /* Special case of SIN from AIPS*/ +#define WCS_GLS 30 /* Same as SFL from AIPS*/ +#define WCS_DSS 31 /* Digitized Sky Survey plate solution */ +#define WCS_PLT 32 /* Plate fit polynomials (SAO) */ +#define WCS_TNX 33 /* Tangent Plane (NOAO corrections) */ +#define WCS_ZPX 34 /* Zenithal Azimuthal Polynomial (NOAO corrections) */ +#define WCS_TPV 35 /* Tangent Plane (SCAMP corrections) */ +#define WCS_TOA 36 /* TOAST */ +#define NWCSTYPE 37 /* Number of WCS types (-1 really means no WCS) */ /* Coordinate systems */ #define WCS_J2000 1 /* J2000(FK5) right ascension and declination */ diff --git a/libwcs/wcsinit.c b/libwcs/wcsinit.c index c184ca0..e8d40db 100644 --- a/libwcs/wcsinit.c +++ b/libwcs/wcsinit.c @@ -536,6 +536,14 @@ char *wchar; /* Suffix character for one of multiple WCS */ hgetr8c (hstring, keyword, &mchar, &wcs->prj.p[i]); } } + else if (wcs->prjcode == WCS_HPX) { + hgetr8c (hstring, "PV2_1", &mchar, &wcs->prj.p[1]); + hgetr8c (hstring, "PV2_2", &mchar, &wcs->prj.p[2]); + } + + else if (wcs->prjcode == WCS_TOA) { + hgetr8c (hstring, "PV2_1", &mchar, &wcs->prj.p[1]); + } /* Initialize TNX, defaulting to TAN if there is a problem */ if (wcs->prjcode == WCS_TNX) { diff --git a/libwcs/wcslib.h b/libwcs/wcslib.h index 7b10f9a..0603da8 100644 --- a/libwcs/wcslib.h +++ b/libwcs/wcslib.h @@ -109,7 +109,7 @@ extern void poly_addcste(polystruct *poly, double *cste), #endif extern int npcode; -extern char pcodes[26][4]; +extern char pcodes[29][4]; struct prjprm { char code[4]; @@ -219,6 +219,15 @@ struct prjprm { int qscset(struct prjprm *); int qscfwd(const double, const double, struct prjprm *, double *, double *); int qscrev(const double, const double, struct prjprm *, double *, double *); + int hpxset(struct prjprm *); + int hpxfwd(const double, const double, struct prjprm *, double *, double *); + int hpxrev(const double, const double, struct prjprm *, double *, double *); + int xphset(struct prjprm *); + int xphfwd(const double, const double, struct prjprm *, double *, double *); + int xphrev(const double, const double, struct prjprm *, double *, double *); + int toaset(struct prjprm *); + int toafwd(const double, const double, struct prjprm *, double *, double *); + int toarev(const double, const double, struct prjprm *, double *, double *); int raw_to_pv(struct prjprm *prj, double x, double y, double *xo, double *yo); #else int prjset(), prjfwd(), prjrev(); -- 2.11.0