diff --git a/include/GHealpix.hpp b/include/GHealpix.hpp index d5474ac..d1cfe3a 100644 --- a/include/GHealpix.hpp +++ b/include/GHealpix.hpp @@ -49,9 +49,11 @@ public: // Methods void read(const GFitsHDU* hdu); int nside(void) const; + int order(void) const; int num_pixels(void) const; double omega(void) const; GSkyDir pix2ang(const int& ipix); + int ang2pix(GSkyDir dir) const; private: // Private methods @@ -59,15 +61,20 @@ private: void copy_members(const GHealpix& pixels); void free_members(void); GHealpix* clone(void) const; + int nside2order(int nside); void pix2xy(const int& ipix, int* x, int* y); + int xy2pix(int x, int y) const; void pix2ang_ring(int ipix, double* theta, double* phi); void pix2ang_nest(int ipix, double* theta, double* phi); + int ang2pix_z_phi_ring(double z, double phi) const; + int ang2pix_z_phi_nest(double z, double phi) const; // Private data area int m_nside; //!< Number of divisions of each base pixel (1-8192) int m_npface; //!< int m_ncap; //!< - int m_order; //!< Ordering (0=ring, 1=nested, -1=?) + int m_order; //!< Order + int m_scheme; //!< Ordering scheme (0=ring, 1=nested, -1=?) int m_coordsys; //!< Coordinate system (0=equatorial, 1=galactic, -1=?) int m_num_pixels; //!< Number of pixels int m_size_pixels; //!< Vector size of each pixel diff --git a/pyext/GHealpix.i b/pyext/GHealpix.i index 5be778d..14cb27f 100644 --- a/pyext/GHealpix.i +++ b/pyext/GHealpix.i @@ -40,6 +40,7 @@ public: int num_pixels(void) const; double omega(void) const; GSkyDir pix2ang(const int& ipix); + int ang2pix(GSkyDir dir) const; }; diff --git a/src/GHealpix.cpp b/src/GHealpix.cpp index 411420d..69c83dc 100644 --- a/src/GHealpix.cpp +++ b/src/GHealpix.cpp @@ -19,7 +19,8 @@ /* __ Includes ___________________________________________________________ */ #include -#include +//#include +#include #include "GException.hpp" #include "GTools.hpp" #include "GHealpix.hpp" @@ -37,15 +38,19 @@ /* __ Debug definitions __________________________________________________ */ +/* __ Local prototypes ___________________________________________________ */ +unsigned int isqrt(unsigned int arg); + /* __ Constants __________________________________________________________ */ -const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4}; -const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7}; +const int jrll[12] = {2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4}; +const int jpll[12] = {1, 3, 5, 7, 0, 2, 4, 6, 1, 3, 5, 7}; +const int order_max = 13; +const int ns_max = 1 << order_max; /* __ Static conversion arrays ___________________________________________ */ static short ctab[0x100]; static short utab[0x100]; - /*========================================================================== = = = GHealpix constructors/destructors = @@ -178,6 +183,8 @@ void GHealpix::read(const GFitsHDU* hdu) m_fact2 = 4.0 / m_num_pixels; m_fact1 = 2 * m_nside * m_fact2; m_omega = fourpi / m_num_pixels; + m_order = nside2order(m_nside); + // Get ordering scheme from ORDERING keyword std::string ordering; @@ -185,14 +192,14 @@ void GHealpix::read(const GFitsHDU* hdu) ordering = hdu->card("ORDERING")->string(); } catch (GException::fits_key_not_found &e) { - m_order = -1; // Flag unknown ordering + m_scheme = -1; // Flag unknown ordering } // Decode ordering scheme string if (ordering == "RING") - m_order = 0; + m_scheme = 0; else if (ordering == "NESTED") - m_order = 1; + m_scheme = 1; else throw GException::healpix(G_LOAD, "Invalid ordering"); @@ -268,6 +275,16 @@ int GHealpix::nside(void) const /***********************************************************************//** + * @brief Returns ordering parameter. + ***************************************************************************/ +int GHealpix::order(void) const +{ + // Return order + return m_order; +} + + +/***********************************************************************//** * @brief Returns number of pixels. ***************************************************************************/ int GHealpix::num_pixels(void) const @@ -300,7 +317,7 @@ GSkyDir GHealpix::pix2ang(const int& ipix) double phi = 0.0; // Perform ordering dependent conversion - switch (m_order) { + switch (m_scheme) { case 0: pix2ang_ring(ipix, &theta, &phi); break; @@ -328,6 +345,49 @@ GSkyDir GHealpix::pix2ang(const int& ipix) } +/***********************************************************************//** + * @brief Returns pixel for a given sky direction + * + * @param[in] ipix Pixel number (0,1,...,num_pixels) + ***************************************************************************/ +int GHealpix::ang2pix(GSkyDir dir) const +{ + // Declare result + int ipix; + double z; + double phi; + + // Compute coordinate system dependent (z,phi) + switch (m_coordsys) { + case 0: + z = cos(pihalf-dir.dec()); + phi = dir.ra(); + break; + case 1: + z = cos(pihalf-dir.b()); + phi = dir.l(); + break; + default: + break; + } + + // Perform ordering dependent conversion + switch (m_scheme) { + case 0: + ipix = ang2pix_z_phi_ring(z, phi); + break; + case 1: + ipix = ang2pix_z_phi_nest(z, phi); + break; + default: + break; + } + + // Return pixel index + return ipix; +} + + /*========================================================================== = = = GHealpix private methods = @@ -344,6 +404,7 @@ void GHealpix::init_members(void) m_npface = 0; m_ncap = 0; m_order = 0; + m_scheme = 0; m_coordsys = 0; m_num_pixels = 0; m_size_pixels = 0; @@ -380,6 +441,7 @@ void GHealpix::copy_members(const GHealpix& pixels) m_npface = pixels.m_npface; m_ncap = pixels.m_ncap; m_order = pixels.m_order; + m_scheme = pixels.m_scheme; m_coordsys = pixels.m_coordsys; m_num_pixels = pixels.m_num_pixels; m_size_pixels = pixels.m_size_pixels; @@ -434,6 +496,32 @@ GHealpix* GHealpix::clone(void) const /***********************************************************************//** + * @brief Convert nside to order + * + * @param[in] nside Number of sides. + ***************************************************************************/ +int GHealpix::nside2order(int nside) +{ + // Initialise order + int order = -1; + + // Determine order + for (int m = 0; m <= order_max; ++m) { + int nstest = 1 << m; + if (nside == nstest) { + order = m; + break; + } + if (nside < nstest) + break; + } + + // Return order + return order; +} + + +/***********************************************************************//** * @brief Convert pixel index to (x,y) coordinate * * @param[in] ipix Pixel index for which (x,y) are to be computed. @@ -456,6 +544,19 @@ void GHealpix::pix2xy(const int& ipix, int* x, int* y) /***********************************************************************//** + * @brief Convert (x,y) coordinate to pixel + * + * @param[in] x x coordinate. + * @param[in] y y coordinate. + ***************************************************************************/ +int GHealpix::xy2pix(int x, int y) const +{ + // Return pixel + return utab[x&0xff] | (utab[x>>8]<<16) | (utab[y&0xff]<<1) | (utab[y>>8]<<17); +} + + +/***********************************************************************//** * @brief Convert pixel index to (theta,phi) angles for ring ordering * * @param[in] ipix Pixel index for which (theta,phi) are to be computed. @@ -470,48 +571,34 @@ void GHealpix::pix2ang_ring(int ipix, double* theta, double* phi) if (ipix < 0 || ipix >= m_num_pixels) throw GException::out_of_range(G_PIX2ANG_RING, ipix, 0, m_num_pixels-1); - // Setup - int ipix1 = ipix + 1; // in {1, m_num_pixels} - int nl2 = 2 * m_nside; - int nl4 = 4 * m_nside; - - // Determine points in each polar cap (0 for nside=1) - int ncap = 2 * m_nside * (m_nside-1); - - double fact1 = 1.5 * m_nside; - double fact2 = 3.0 * m_nside * m_nside; - // Handle North Polar cap - if (ipix1 <= ncap) { - double hip = ipix1/2.; - double fihip = floor(hip); - int iring = (int)floor(sqrt(hip - sqrt(fihip)))+1; - int iphi = ipix1 - 2*iring*(iring - 1); - *theta = acos(1. - iring*iring / fact2); - *phi = (1.*iphi - 0.5) * pi/(2.*iring); + if (ipix < m_ncap) { + int iring = int(0.5*(1+isqrt(1+2*ipix))); // counted from North pole + int iphi = (ipix+1) - 2*iring*(iring-1); + *theta = acos(1.0 - iring*iring / m_fact2); + *phi = (iphi - 0.5) * pi/(2.0*iring); } // Handle Equatorial region - else if(ipix1 <= nl2*(5*m_nside+1)) { - int ip = ipix1 - ncap - 1; - int iring = (int)floor(ip/nl4) + m_nside; - int iphi = (int)fmod(ip,nl4) + 1; - double fodd = 0.5 * (1 + fmod((double)(iring+m_nside),2)); - *theta = acos((nl2 - iring) / fact1); - *phi = (1.*iphi - fodd) * pi/(2.*m_nside); + else if (ipix < (m_num_pixels - m_ncap)) { + int ip = ipix - m_ncap; + int iring = ip/(4*m_nside) + m_nside; // counted from North pole + int iphi = ip%(4*m_nside) + 1; + double fodd = ((iring+m_nside)&1) ? 1 : 0.5; + int nl2 = 2*m_nside; + *theta = acos((nl2 - iring) / m_fact1); + *phi = (iphi - fodd) * pi/nl2; } - // Handle South Polar cap + // Handle South Polar cap else { - int ip = m_num_pixels - ipix1 + 1; - double hip = ip/2.; - double fihip = floor(hip); - int iring = (int)floor(sqrt( hip - sqrt(fihip))) + 1; - int iphi = (int)(4.*iring + 1 - (ip - 2.*iring*(iring-1))); - *theta = acos( -1. + iring*iring / fact2 ); - *phi = (1.*iphi - 0.5) * pi/(2.*iring); + int ip = m_num_pixels - ipix; + int iring = int(0.5*(1+isqrt(2*ip-1))); // Counted from South pole + int iphi = 4*iring + 1 - (ip - 2*iring*(iring-1)); + *theta = acos(iring*iring / m_fact2 - 1.0); + *phi = (iphi - 0.5) * pi/(2.*iring); } - + // Return return; } @@ -534,7 +621,7 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) // Get face number and index in face int nl4 = 4 * m_nside; - int face_num = ipix/m_npface; // Face number in {0,11} + int face_num = ipix >> (2*m_order); // Face number in {0,11} int ipf = ipix & (m_npface - 1); // Get pixel coordinates @@ -543,7 +630,7 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) pix2xy(ipf, &ix, &iy); // Computes the z coordinate on the sphere - int jr = jrll[face_num]*m_nside - ix - iy - 1; + int jr = (jrll[face_num] << m_order) - ix - iy - 1; // Declare result variables int nr; @@ -585,6 +672,118 @@ void GHealpix::pix2ang_nest(int ipix, double* theta, double* phi) } +/***********************************************************************//** + * @brief Returns pixels which contains angular coordinates (z,phi) + * + * @param[in] z Cosine of zenith angle - cos(theta). + * @param[in] phi Azimuth angle in radians. + ***************************************************************************/ +int GHealpix::ang2pix_z_phi_ring(double z, double phi) const +{ + // Initialise pixel + int ipix = 0; + + // Setup + double za = abs(z); + double tt = modulo(phi,twopi) * inv_pihalf; // in [0,4) + + // Equatorial region + if (za <= twothird) { + double temp1 = m_nside*(0.5+tt); + double temp2 = m_nside*z*0.75; + int jp = int(temp1-temp2); // index of ascending edge line + int jm = int(temp1+temp2); // index of descending edge line + int ir = m_nside + 1 + jp - jm; // in {1,2n+1} + int kshift = 1 - (ir & 1); // kshift=1 if ir even, 0 otherwise + int ip = (jp+jm-m_nside+kshift+1)/2; // in {0,4n-1} + ip = modulo(ip,4*m_nside); + ipix = m_ncap + (ir-1)*4*m_nside + ip; + } + + // North & South polar caps + else { + double tp = tt - int(tt); + double tmp = m_nside * sqrt(3*(1-za)); + int jp = int(tp*tmp); // increasing edge line index + int jm = int((1.0-tp)*tmp); // decreasing edge line index + int ir = jp + jm + 1; // ring number counted from the closest pole + int ip = int(tt*ir); // in {0,4*ir-1} + ip = modulo(ip,4*ir); + if (z>0) + ipix = 2*ir*(ir-1) + ip; + else + ipix = m_num_pixels - 2*ir*(ir+1) + ip; + } + + // Return pixel + return ipix; +} + + +/***********************************************************************//** + * @brief Returns pixels which contains angular coordinates (z,phi) + * + * @param[in] z Cosine of zenith angle - cos(theta). + * @param[in] phi Azimuth angle in radians. + ***************************************************************************/ +int GHealpix::ang2pix_z_phi_nest(double z, double phi) const +{ + // Initialise face and pixel numbers + int face_num; + int ix; + int iy; + + // Setup + double za = abs(z); + double tt = modulo(phi,twopi) * inv_pihalf; // in [0,4) + + // Equatorial region + if (za <= twothird) { + double temp1 = ns_max*(0.5+tt); + double temp2 = ns_max*z*0.75; + int jp = int(temp1-temp2); // index of ascending edge line + int jm = int(temp1+temp2); // index of descending edge line + int ifp = jp >> order_max; // in {0,4} + int ifm = jm >> order_max; + if (ifp == ifm) // faces 4 to 7 + face_num = (ifp==4) ? 4: ifp+4; + else if (ifp < ifm) // (half-)faces 0 to 3 + face_num = ifp; + else // (half-)faces 8 to 11 + face_num = ifm + 8; + ix = jm & (ns_max-1); + iy = ns_max - (jp & (ns_max-1)) - 1; + } + + // Polar region, za > 2/3 + else { + int ntt = int(tt); + double tp = tt-ntt; + double tmp = ns_max*sqrt(3*(1-za)); + int jp = int(tp*tmp); // increasing edge line index + int jm = int((1.0-tp)*tmp); // decreasing edge line index + if (jp >= ns_max) jp = ns_max-1; // for points too close to the boundary + if (jm >= ns_max) jm = ns_max-1; + if (z >= 0) { + face_num = ntt; // in {0,3} + ix = ns_max - jm - 1; + iy = ns_max - jp - 1; + } + else { + face_num = ntt + 8; // in {8,11} + ix = jp; + iy = jm; + } + } + + // Get pixel + int ipf = xy2pix(ix, iy); + ipf >>= (2*(order_max - m_order)); // in {0, nside**2 - 1} + return ipf + (face_num<<(2*m_order)); // in {0, 12*nside**2 - 1} +} + + + /*========================================================================== = = = GHealpix friends = @@ -601,15 +800,18 @@ ostream& operator<< (ostream& os, const GHealpix& pixels) { // Put header in stream os << "=== GHealpix ===" << endl; - os << " Number of base pixel div. .: " << pixels.m_nside << endl; - os << " Number of pixels ..........: " << pixels.m_num_pixels << endl; + os << " Nside (number of divisions): " << pixels.m_nside << endl; + os << " Npface (pixels per face) ..: " << pixels.m_npface << endl; + os << " Ncap (number of cap pixels): " << pixels.m_ncap << endl; + os << " Npix (number of pixels) ...: " << pixels.m_num_pixels << endl; + os << " Order .....................: " << pixels.m_order << endl; os << " Pixel vector size .........: " << pixels.m_size_pixels << endl; os << " Solid angle ...............: " << scientific << pixels.m_omega - << fixed << endl; + << fixed << " sr" << endl; // Put ordering in stream os << " Ordering ..................: "; - switch (pixels.m_order) { + switch (pixels.m_scheme) { case 0: os << "Ring" << endl; break; @@ -628,22 +830,19 @@ ostream& operator<< (ostream& os, const GHealpix& pixels) os << " Coordinate system .........: "; switch (pixels.m_coordsys) { case 0: - os << "Equatorial (RA,Dec)" << endl; + os << "Equatorial (RA,Dec)"; break; case 1: - os << "Galactic (l,b)" << endl; + os << "Galactic (l,b)"; break; case -1: - os << "*** Unknown ***" << endl; + os << "*** Unknown ***"; break; default: - os << "*** Invalid ***" << endl; + os << "*** Invalid ***"; break; } - for (int i = 0; i < pixels.m_num_pixels; ++i) - os << pixels.m_dir[i] << " "; - // Return output stream return os; } @@ -654,3 +853,16 @@ ostream& operator<< (ostream& os, const GHealpix& pixels) = Other functions used by GHealpix = = = ==========================================================================*/ + +/***********************************************************************//** + * @brief Integer n that fulfills n*n <= arg < (n+1)*(n+1) + * + * @param[in] arg Argument. + * + * Returns the integer \a n, which fulfills \a n*n <= arg < (n+1)*(n+1). + ***************************************************************************/ +unsigned int isqrt(unsigned int arg) +{ + // Return + return unsigned(sqrt(arg+0.5)); +}