HEALPix is an interesting framework that can be used to store information about spherical surfaces in a digital format as pixels that have nice properties (hierarchical, equal-area and isolatitude), and it is used by NASA and SETI@home (among others). The detail of HEALPix is presented in a paper by Górski et al. (2005), which was published in The Astrophysical Journal and is titled “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere” (DOI:10.1086/427976). Despite being an excellent paper that is filled with mathematical info and descriptions, it seems like the authors ran out of steam as they got to Equation 18, since this equation apparently does not accurately reflect the general case for calculating the longitude index in the nested scheme of pixel indexing. The following describes a revised version of Equation 18, along with some additional comments. (Be sure to have a copy of the paper at hand — if you cannot download the original paper through the DOI link above, see the Harvard Astronomy Abstract Service for alternative links or visit arXiv:astro-ph/0409513.)
Having wrecked my brain on this, it seems that Equation 18 (on page 765) actually only describes the case for pixels in the equatorial zone. The general equation that describes the longitude index also for polar pixels needs to look slightly different, namely:
j = (F2 . nr + h + s) / 2 ,
where nr is introduced shortly and replaces the Nside constant in the published Equation 18. The variable nr represents the number of pixels in one quadrant of a specific constant-latitude ring and is therefore a function of ring index i: for the north polar cap (where i < Nside) we have nr = i, for the entire equatorial belt (where Nside <= i <= 3 . Nside) we use nr = Nside (which again produces the original Equation 18 for this region), and for the south polar cap (where i > 3 . Nside) we need to define nr = 4 . Nside - i (thus, using mirror symmetry with respect to the equator).
For verification purposes, nested pixel index 29 makes for a good test candidate along with Nside = 4, which is also conveniently illustrated in the bottom panel of Figure 4 (on page 764). For this case the pixel is thoroughly in the northern polar region and we have: pn = 29, pn' = 13, f = 1, frow = 0, F1 = 2, F2 = 3, x = 3, y = 2, v = 5, h = 1, i = 2, s = 1, and z = cos(theta) = 11 / 12 = 0.9167. Using the original Equation 18, it incorrectly produces j = 7, resulting in phi / pi = 1.6250 when substituted into Equation 5. Visual inspection of Figure 4 shows that the phi / pi value should be somewhere between 0.75 and 1.0 for this pixel (considering the pixel’s boundaries), and, since this is the 4th pixel in the second ring, the value of j should be somewhere in the region of 3 or 4, but not 7. Even if one were to attempt substituting j = 7 into Equation 9 (actually meant for the equatorial belt), this would still produce an incorrect result of phi / pi = 0.8125 (although it would be more in line with the expected range). But if one were to use the revised equation stated above for the north polar cap, we would get j = 4 that would result in phi / pi = 0.8750 from Equation 5 (which is also exactly in the center of our expected range). This last result is also the value returned by the reference implementation of the HEALPix software, version 2.12a (available from SourceForge).
A few minor points before I get to an interesting observation in the last section:
- This is really just a technicality, but on page 761 (roughly in the center of the right column) where it states “for odd and even values of Ntheta, respectively”, I suspect that the words “odd and even” should have been reversed (if I interpret this correctly).
- If you are wondering what the variables p and ph mean, which are used for the first time on page 763 above Equation 2 without introduction, this is how I have come to think of them: variable p is simply an integer pixel index that runs from 0 up to and including 12 . Nside2 - 1. The variable ph can then be interpreted as the half-step-size pixel-index counterpart of p, which runs from 0.5 up to and including 6 . Nside2 in steps of half units. Furthermore, the p' variable is simply the equatorial pixel-index, starting at zero for the first pixel in the equatorial band. Also, although p is used to index pixels in a linear fashion, similar to the ring scheme, it is probably best not to interpret p as the ring-scheme index per se.
- Although the paper does not state this explicitly (and I was wondering about it as I got to the start of page 765), one can calculate the base-resolution pixel index number, f, from the nested index, pn, using f = I( pn / Nside2 ), where I(x) represents the integer floor function as it does in the rest of the paper. This may be obvious, but it is nice to know for sure that you are using the right equation.
Finally, and I found this quite refreshing, if you want to calculate the longitude, phi, of a pixel center from its nested pixel index, pn, you actually need not calculate either s or j to arrive at phi. Substitution of the revised Equation 18 into Equations 5 and 9 (for the different regions) produces a unified general equation that can be applied anywhere on the sphere:
phi = ( pi / 4 ) . ( F2 . nr + h ) / nr .
This last equation could likely be used to optimize some of the functions in the reference implementation of the HEALPix software; for example, one may eliminate calculations and usages of the internal variables named
kshift (the boolean complement of variable s) in the
pix2ang_nest functions of the C, Java, Fortran and IDL code, as well as the
pix2ang_z_phi implementation in C++, and possibly elsewhere.