Discretisation of Functions on the Sphere for High Resolution Applications: a Motivation for HEALPix

Numerical analysis of functions on the sphere involves (1) a class of mathematical operations, whose objects are (2) discretised maps, i.e. quantizations of arbitrary functions according to a chosen tessellation (exhaustive partition of the sphere into finite area elements). Hereafter we mostly specialise our discussion to CMB related applications of HEALPix, but all our statements hold true generally for any relevant deterministic and random functions on the sphere.

Considering point (1): Standard operations of numerical analysis which one might wish to execute on the sphere include convolutions with local and global kernels, Fourier analysis with spherical harmonics and power spectrum estimation, wavelet decomposition, nearest-neighbour searches, topological analysis, including searches for extrema or zero-crossings, computing Minkowski functionals, extraction of patches and finite differencing for solving partial differential equations. Some of these operations become prohibitively slow if the sampling of functions on the sphere, and the related structure of the discrete data set, are not designed carefully.

Regarding point (2): Typically, a whole sky map rendered by a CMB experiment contains (i) signals coming from the sky, which are by design strongly band-width limited (in the sense of spatial Fourier decomposition) by the instrument's angular response function, and (ii) a projection into the elements of a discrete map, or pixels, of the observing instrument's noise; this pixel noise should be random, and white, at least near the discretisation scale, with a band-width significantly exceeding that of all the signals.

With these considerations in mind we propose the following list of desiderata for the mathematical structure of discretised full sky maps:

1. Hierarchical structure of the data base. This is recognised as essential for very large data bases, and was postulated in construction of the Quadrilateralized Spherical Cube (or quad-sphere, see https://lambda.gsfc.nasa.gov/product/cobe/skymap_info_new.cfm), which was used for the COBE data. An argument in favour of this proposition states that the data elements which are nearby in a multi-dimensional configuration space (here, on the surface of a sphere), are also nearby in the tree structure of the data base, hence the near-neighbour searches are conducted optimally in the data storage medium or computer RAM. This property, especially when implemented with a small number of base resolution elements, facilitates various topological methods of analysis, and allows easy construction of wavelet transforms on quadrilateral (and also triangular) grids. Figure 1 shows how a hierarchical partition with quadrilateral structure naturally allows for a binary vector indexing of the data base.

Figure 1: Quadrilateral tree pixel numbering. The coarsely pixelated coordinate patch on the left consists of four pixels. Two bits suffice to label the pixels. To increase the resolution, every pixel splits into 4 daughter pixels shown on the right. These daughters inherit the pixel index of their parent (boxed) and acquire two new bits to give the new pixel index. Several such curvilinearly mapped coordinate patches (12 in the case of HEALPix, and 6 in the case of the COBE quad-sphere) are joined at the boundaries to cover the sphere. All pixels indices carry a prefix (here omitted for clarity) which identifies which base-resolution pixel they belong to.

2. Equal areas of discrete elements of partition. This is advantageous because (i) white noise generated by the signal receiver gets integrated exactly into white noise in the pixel space, and (ii) sky signals are sampled without regional dependence, except for the dependence on pixel shapes, which is unavoidable with tessellations of the sphere. Hence, as much as possible given the experimental details, the pixel size should be made sufficiently small compared to the instrument's resolution to avoid any excessive, and pixel shape dependent, signal smoothing.

3. Iso-Latitude distribution of discrete area elements on a sphere. This property is critical for computational speed of all operations involving evaluation of spherical harmonics. Since the associated Legendre polynomial components of spherical harmonics are evaluated via slow recursions, and can not be simply handled in an analogous way to the trigonometric Fast Fourier Transform, any deviations in the sampling grid from an iso-latitude distribution result in a prohibitive loss of computational performance with the growing number of sampling points, or increasing map resolution. It is precisely this property that the COBE quad-sphere is lacking, and this renders it impractical for applications to high resolution data.

A number of tessellations have been used for discretisation and analysis of functions on the sphere (for example, see Driscoll & Healy (1994), Muciaccia, Natoli & Vittorio (1998), Doroshkevich et al. (2005) — rectangular grids, Baumgardner & Frederickson (1985), Tegmark (1996) — icosahedral grids, Saff & Kuijlaars (1997), Crittenden & Turok (1998) — `igloo' grids, and Szalay & Brunner (1998) — a triangular grid), but none satisfies simultaneously all three stated requirements.

All three requirements formulated above are satisfied by construction with the Hierarchical Equal Area, iso-Latitude Pixelation (HEALPix) of the sphere, which is shown in Figure 2. A more detailed description of HEALPix, its motivations, and applications can be found in Górski et al. (2005).

Figure 2: Orthographic view of HEALPix partition of the sphere. Overplot of equator and meridians illustrates the octahedral symmetry of HEALPix. Light-gray shading shows one of the eight (four north, and four south) identical polar base-resolution pixels. Dark-gray shading shows one of the four identical equatorial base-resolution pixels. Moving clockwise from the upper left panel the grid is hierarchically subdivided with the grid resolution parameter equal to $N_{\mathrm{side}}=$ 1, 2, 4, 8, and the total number of pixels equal to $N_{\mathrm{pix}}= 12 \times N_{\mathrm{side}}^2$ = 12, 48, 192, 768. All pixel centers are located on $N_{\mathrm{ring}} = 4 \times N_{\mathrm{side}}- 1$ rings of constant latitude. Within each panel the areas of all pixels are identical.