S2LET: A code to perform fast wavelet analysis on the sphere
^{1} Department of Physics and
Astronomy, University College London, London WC1E
6BT, UK
email: boris.leistedt.11@ucl.ac.uk
^{2} Mullard Space Science Laboratory
(MSSL), University College London, Surrey
RH5 6NT,
UK
email: jason.mcewen@ucl.ac.uk
^{3} Institute of Electrical Engineering,
Ecole Polytechnique Fédérale de Lausanne (EPFL), 1015
Lausanne,
Switzerland
email: pierre.vandergheynst@epfl.ch
^{4} Department of Radiology and Medical
Informatics, University of Geneva (UniGE), 1211
Geneva,
Switzerland
^{5} Department of Medical Radiology,
Lausanne University Hospital (CHUV), 1011
Lausanne,
Switzerland
email: yves.wiaux@epfl.ch
^{6} Institute of Sensors, Signals
& Systems, Heriot Watt University, Edinburgh
EH14 4AS,
UK
Received:
13
November
2012
Accepted:
27
August
2013
We describe S2LET, a fast and robust implementation of the scalediscretised wavelet transform on the sphere. Wavelets are constructed through a tiling of the harmonic line and can be used to probe spatially localised, scaledependent features of signals on the sphere. The reconstruction of a signal from its wavelets coefficients is made exact here through the use of a sampling theorem on the sphere. Moreover, a multiresolution algorithm is presented to capture all information of each wavelet scale in the minimal number of samples on the sphere. In addition S2LET supports the HEALPix pixelisation scheme, in which case the transform is not exact but nevertheless achieves good numerical accuracy. The core routines of S2LET are written in C and have interfaces in Matlab, IDL and Java. Real signals can be written to and read from FITS files and plotted as Mollweide projections. The S2LET code is made publicly available, is extensively documented, and ships with several examples in the four languages supported. At present the code is restricted to axisymmetric wavelets but will be extended to directional, steerable wavelets in a future release.
Key words: techniques: high angular resolution / methods: data analysis / techniques: image processing
© ESO, 2013
1. Introduction
Signals defined or measured on the sphere arise in numerous disciplines, where analysis techniques defined explicitly on the sphere are now in common use. In particular, wavelets on the sphere (Antoine & Vandergheynst 1998, 1999; Baldi et al. 2009; Marinucci et al. 2008; McEwen et al. 2006a; Narcowich et al. 2006; Starck et al. 2006a; Wiaux et al. 2006a, 2005, 2007, 2008a; Yeo et al. 2008) have been applied very successfully to problems in astrophysics and cosmology, where datasets are increasingly large and need to be analysed at high resolution in order to confront accurate theoretical predictions (e.g. Barreiro et al. 2000; Basak & Delabrouille 2012; Cayón et al. 2001; Deriaz et al. 2012; Labatie et al. 2012; Lan & Marinucci 2008; McEwen et al. 2006b, 2007a,b, 2008; Pietrobon et al. 2008; Starck et al. 2006b; Schmitt et al. 2010; Vielva et al. 2004, 2006a,b, 2007; Wiaux et al. 2006b, 2008b).
While wavelet theory is well established in Euclidean space (see e.g. Daubechies 1992), multiple wavelet frameworks have been developed on the sphere, only a fraction of which lead to exact transforms in both the continuous and discrete settings. In fact, discrete methodologies (Schröder & Sweldens 1995; Sweldens 1996, 1997) achieve exactness in practice but may not lead to a stable basis on the sphere (Sweldens 1997). In the continuous setting several constructions are theoretically exact, and have been combined with sampling theorems on the sphere to enable exact reconstruction in the discrete setting also. In particular, scalediscretised wavelets (Wiaux et al. 2008a) lean on a tiling of the harmonic line to yield an exact wavelet transform in both the continuous and discrete settings. In the axisymmetric case, the scalediscretised wavelets reduce to needlets (Narcowich et al. 2006; Baldi et al. 2009; Marinucci et al. 2008), which were developed independently using an analogous tiling of the harmonic line. Similarly, the isotropic undecimated wavelet transform (UWT) developed by (Starck et al. 2006a) exploits Bsplines of order 3 to cover the harmonic line with filters with greater overlap but nevertheless compact support.
In this paper we describe the new publicly available S2LET^{1} code to perform the scalediscretised wavelet transform of complex signals on the sphere. At present S2LET is restricted to axisymmetric wavelets (i.e. azimuthally symmetric when centred on the poles) and includes generating functions for axisymmetric scalediscretised wavelets (Wiaux et al. 2008a), needlets (Narcowich et al. 2006; Baldi et al. 2009; Marinucci et al. 2008) and Bspline wavelets (Starck et al. 2006a). We intend to extend the code to directional, steerable wavelets and spin functions in a future release. The core routines of S2LET are written in C, exploit fast algorithms on the sphere, and have interfaces in Matlab, IDL and Java.
We note that many very useful public codes are already available to compute wavelet transforms on the sphere, including isotropic undecimated wavelet, ridgelet and curvelet transforms^{2} (Starck et al. 2006a), invertible filter banks^{3} (Yeo et al. 2008), needlets (NeedATool^{4}; Pietrobon et al. 2010) and scalediscretised wavelets (S2DW^{5}; Wiaux et al. 2008a). S2LET aims primarily to provide a fast and flexible implementation of the scalediscretised transform with exact reconstruction on the sphere using the sampling theorem of McEwen & Wiaux (2011), although it has also been extended to support some of the features of these other codes. Furthermore, particular attention has been paid in the development of S2LET to prove a userfriendly code, supporting multiple programming languages, and which is extensively documented.
The remainder of this article is organised as follows. In Sect. 2 we detail the construction of scalediscretised axisymmetric wavelets and the corresponding exact scalediscretised wavelet transform on the sphere. In Sect. 3 we describe the S2LET code, including implementation details, computational complexity and numerical performance. We present a number of simple examples using S2LET in Sect. 4, along with the code to execute them. We conclude in Sect. 5.
2. Wavelets on the sphere
We review the construction of scalediscretised wavelets on the sphere through tiling of the harmonic line (Wiaux et al. 2008a). Directional, steerable wavelets were also considered by Wiaux et al. (2008a), however we restrict our attention to axisymmetric wavelets here. Furthermore, the use of a sampling theorem on the sphere guarantees that spherical harmonic coefficients capture all the information content of bandlimited signals, resulting in theoretically exact harmonic and wavelet transforms. One may alternatively adopt samplings of the sphere for which exact quadrature rules do not exist, such as HEALPix (Górski et al. 2005), but which nevertheless exhibit other useful properties, leading to numerically accurate but not theoretically exact transforms.
2.1. Harmonic analysis on the sphere
The spherical harmonic decomposition of a square integrable signal f ∈ L^{2}(S^{2}) on the twodimensional sphere S^{2} reads (1)where Y_{ℓm} are the spherical harmonic functions, which form the canonical orthogonal basis on S^{2}. The spherical harmonic coefficients f_{ℓm}, with ℓ ∈N and m ∈Z such that  m ≤ ℓ, form a dual representation of the signal f in the harmonic basis on the sphere. The angular position ω = (θ,φ) ∈ S^{2} is specified by colatitude θ ∈ [0,π] and longitude φ ∈ [0,2π). The spherical harmonic coefficients are given by (2)with the surface element dΩ(ω) = sinθdθdφ. We consider bandlimited signals in the spherical harmonic basis, with bandlimit L if f_{ℓm} = 0, ∀ℓ ≥ L. For bandlimited signals sampling theorems can be invoked so that both forward and inverse transforms can be reduced to finite summations that are theoretically exact. Sampling theorems effectively encode a quadrature rule for the exact evaluation of integrals on the sphere from a finite set of sampling nodes. Various sampling theorems exist in the literature (e.g. Driscoll & Healy 1994; Healy et al. 1996; McEwen & Wiaux 2011). In this work we adopt the McEwen & Wiaux (2011) sampling theorem (hereafter MW), which is based on an equiangular sampling scheme and, for a given bandlimit L, requires the lowest number of samples on the sphere of all sampling theorems, namely (L − 1)(2L − 1) + 1 ~ 2L^{2} samples (for comparison ~4L^{2} samples are required by Driscoll & Healy 1994). Fast algorithms to compute the corresponding spherical harmonic transform scale as and are numerically stable to bandlimits of at least L = 4096 (McEwen & Wiaux 2011). The GLESP pixelisation scheme (Doroshkevich et al. 2005) also provides a sampling theorem based on the GaussLegendre quadrature, and could be used in place of the MW sampling theorem. However, GLESP uses more samples than GaussLegendre quadrature requires, which may lead to an overhead when considering large bandlimits and numerous wavelet scales. We focus on the MW sampling scheme to obtain a theoretically exact transform. Alternative sampling schemes that are not based on sampling theorems also exist such as HEALPix (Górski et al. 2005), which is supported by S2LET, MRS and Needatool. HEALPix does not lead to exact transforms on the sphere but the resulting approximate transforms nevertheless achieve good accuracy and benefit from other practical advantages, such as equalarea pixels.
2.2. Scalediscretised wavelets on the sphere
The scalediscretised wavelet transform allows one to probe spatially localised, scaledependent content in the signal of interest f ∈ L^{2}(S^{2}). The jth wavelet scale W^{Ψj} ∈ L^{2}(S^{2}) is defined as the convolution of f with the wavelet Ψ^{j} ∈ L^{2}(S^{2}): (3)where * denotes complex conjugation. Convolution on the sphere is defined by the inner product of f with the rotated wavelet ℛ_{ω}Ψ^{j}. We restrict our attention to axisymmetric wavelets, i.e. wavelets that are azimuthally symmetric when centred on the poles. Consequently, the rotation operator ℛ_{ω} is parameterised by angular position ω = (θ,φ) only and not also orientation^{6}. For the axisymmetric case the spherical harmonic decomposition of W^{Ψj} is then simply given by a weighted product in harmonic space: (4)where , f_{ℓm} = ⟨ fY_{ℓm} ⟩ and , and where δ_{m0} is the Kronecker delta symbol.
The wavelet coefficients extract the detail information of the signal only; a scaling function and corresponding scaling coefficients must be introduced to represent the lowfrequency (lowℓ), approximate information of the signal. The scaling coefficients W^{Φ} ∈ L^{2}(S^{2}) are defined by the convolution of f with the scaling function Φ ∈ L^{2}(S^{2}): (5)or in harmonic space, (6)where and Φ_{ℓ0}δ_{m0} = ⟨ ΦY_{ℓm} ⟩.
Provided the wavelets and scaling function satisfy an admissibility property (defined below), the function f may be reconstructed exactly from its wavelet and scaling coefficients by (7)or equivalently in harmonic space by (8)The parameters J_{0}, J define the lowest and highest scales j of the wavelet decomposition and must be defined consistently to extract and reconstruct all the information content of f. These parameters depend on the construction of the wavelets and scaling function and are defined explicitly in the next paragraphs. The admissibility condition under which a bandlimited function f can be decomposed and reconstructed exactly is given by the following resolution of the identity: (9)We are now in a position to define wavelets and a scaling function that satisfy the admissibility property. In this paper, we use the smooth generating functions defined by Wiaux et al. (2008a) in order to tile the harmonic line. Alternative definitions are also supported by S2LET, as presented at the end of this section. Consider the C^{∞} Schwartz function with compact support on [−1,1]: (10)for t ∈R. We introduce the positive real parameter to map s(t) to (11)which has compact support in [1 / λ,1]. We then define the smoothly decreasing function k_{λ} by (12)which is unity for t < 1/λ, zero for t > 1, and is smoothly decreasing from unity to zero for t ∈ [1 / λ,1]. We finally define the wavelet generating function by (13)and the scaling function generating function by (14)The wavelets and scaling function are constructed from their generating functions to satisfy the admissibility condition given by Eq. (9). A natural approach is to define from the generating functions κ_{λ} to have support on [λ^{j −1},λ^{j + 1}], yielding (15)For these wavelets Eq. (9) is satisfied for ℓ ≥ λ^{J0}, where J_{0} is the lowest wavelet scale used in the decomposition. The scaling function Φ is constructed to extract the modes that cannot be probed by the wavelets (i.e. modes with ℓ < λ^{J0}): (16)To satisfy exact reconstruction, J is set to ensure the wavelets reach the bandlimit of the signal of interest, yielding J = ⌈ log _{λ} (L − 1) ⌉. The choice of the lowest wavelet scale J_{0} is arbitrary, provided that 0 ≤ J_{0} < J. The wavelets and scaling function may then be reconstructed on the sphere through an inverse spherical harmonic transform. The harmonic tiling and real space representation of these wavelets are shown in Figs. 1 and 2 respectively.
Fig. 1 Wavelets and scaling function constructed with the scalediscretised (SD), needlet and Bspline generating functions (Wiaux et al. 2008a; Marinucci et al. 2008; Starck et al. 2006a) with parameters λ = 3 and J_{0} = 2 and for bandlimit L = 128. The tiling is shown in the top panel, and the profiles of the reconstructed wavelets in the bottom panel. 

Open with DEXTER 
In addition to the scalediscretised generating functions (Wiaux et al. 2008a), S2LET also supports the needlet functions (Marinucci et al. 2008)^{7}, which yield a similar tiling of the harmonic line, as shown in Fig. 1. The Bspline filters used to construct the isotropic undecimated wavelet transform (Starck et al. 2006a) are also supported, as also shown in Fig. 1^{8}. With these three constructions, the wavelets and scaling functions are welllocalised both spatially on the sphere and also in harmonic space. Consequently, the associated wavelet transforms on the sphere can be used to extract spatially localised, scaledependent features in signals of interest.
Fig. 2 Wavelets for scales j ∈ { 2,3,4,5 } and scaling function constructed through a tiling of the harmonic line using scalediscretised functions, with parameters λ = 3 and J_{0} = 2 and for bandlimit L = 128. This plot was produced with a Matlab demo included in S2LET. 

Open with DEXTER 
3. The S2LET code
In this section we describe the S2LET code. We first introduce a multiresolution algorithm to capture each wavelet scale in the minimum number of samples on the sphere, which follows by taking advantage of the reduced bandlimit of the wavelets for scales j < J − 1. This multiresolution algorithm reduces the computation cost of the transform considerably. We then provide details of the implementation, the computational complexity and the numerical accuracy of the scalediscretised wavelet transform supported in S2LET. We finally outline planned future extensions of the code.
3.1. Multiresolution algorithm
In harmonic space, the wavelet coefficients are simply given by the weighted product of the spherical harmonic coefficients of f and the wavelets, as expressed in Eq. (4). Although the wavelet coefficients can be analysed at the same resolution as the signal f (i.e., at full resolution), by construction they have different bandlimits for different scales j, as shown in Fig. 1. The reconstruction can thus be performed at lower resolution, without any loss of information if a sampling theorem is used. This approach yields a multiresolution algorithm where the wavelet coefficients are reconstructed with the minimal number of samples on the sphere: the jth wavelet coefficients have bandlimit k = λ^{j + 1} when using the scalediscretised and needlet kernels, and k = L / λ^{J−j−2} when using the Bsplines. When the MW sampling theorem is used, the wavelets are recovered on (k − 1)(2k − 1) + 1 samples on the sphere. This approach leads to significant improvements in terms of speed and memory use compared to the fullresolution case, as shown in the next section. Figure 3 illustrates the use of the fullresolution and multiresolution transforms on a map of Earth topography data with the scalediscretised filters and the MW scheme. When adopting the HEALPix sampling of the sphere, multiresolution can also be used. However HEALPix does not rely on a sampling theorem and therefore the resolution for the reconstruction of each wavelet scale must be chosen heuristically and adapted to the desired accuracy. For example, in the MRS code (Starck et al. 2006a) it is chosen such that . More detail on the accuracy of the wavelet transform with HEALPix are provided below.
Fig. 3 Scalediscretised wavelet transform of a bandlimited topography map of the Earth for λ = 3, J_{0} = 2 and L = 128, i.e. with the scalediscretised wavelets shown in Fig. 2. The wavelet transform decomposes the bandlimited signal into wavelet coefficients that extract spatially localised, scaledependent features. Since the wavelets for different scales j have different bandlimits, the wavelet coefficients can be reconstructed at lower resolution on the sphere for lower scales j. Panel a) shows the fullresolution wavelet transform of the topography map. The original Earth topography map is shown in the topleft plot, the scaling coefficients are shown in the topright plot, while the wavelet coefficients at scales j ∈ { 2,3,4,5 } are shown lefttoright, toptobottom respectively in the remaining plots. Panel b) shows the same decomposition but using the multiresolution algorithm. The signals shown in panel b) contain the same information as in panel a) but represented in the minimal number of samples on the sphere. These plots were produced by one of the many Matlab demos provided with S2LET. 

Open with DEXTER 
3.2. Implementation
The core numerical routines of S2LET are implemented in C. By adopting a low level programming language such as C for the implementation of the core algorithms, computational efficiency is optimised. The C library includes the fullresolution and multiresolution wavelet transforms, with specific optimisations for real signals in order to take advantage of all symmetries of the spherical harmonic transform. The wavelet transform is computed in harmonic space through Eqs. (4) and (6), for the input parameters (L,λ,J_{0}). To reconstruct signals on the sphere, by default S2LET uses the exact spherical harmonic transform of the MW sampling theorem (McEwen & Wiaux 2011) implemented in the SSHT^{9} code. In this case all transforms are theoretically exact and one can analyse and synthesise real and complex signals at floatingpoint precision. S2LET has been extended to also support the HEALPix sampling scheme, in which case the transform is not theoretically exact but nevertheless achieves good numerical accuracy.
We provide interfaces for the C library in three languages: Matlab, IDL and Java. The Matlab and IDL codes also include routines to read/write signals on the sphere stored in either HEALPix FITS^{10} files or the FITS file format used to stored MW sampled signals. In addition, functionality to plot the Mollweide projection of real signals for both MW or HEALPix samplings is included. The Java interface includes an objectoriented representation of sampled maps, spherical harmonics and wavelet transforms. All routines and interfaces are well documented and illustrated with several examples for both the MW and HEALPix samplings. These examples cover multiple combinations of parameters and types of signals. S2LET requires SSHT, which implements fast and exact algorithms to perform the forward and inverse spherical harmonic transforms corresponding to the MW sampling theorem (McEwen & Wiaux 2011). SSHT in turn requires the FFTW^{11} package for the computation of fast Fourier transforms. The fast spherical harmonic transforms implemented in SSHT compute Wigner functions, and thus the spherical harmonic functions, through efficient recursion using either the method of Trapani & Navaza (2006) or Risbo (1996). Here we present results using the recursion of Risbo (1996). The fast spherical harmonic transform algorithms implemented in SSHT scale as (McEwen & Wiaux 2011).
Although primarily intended to perform the scalediscretised wavelet transform of Wiaux et al. (2008a), S2LET also supports the needlet and splinebased wavelet transforms developed by Marinucci et al. (2008) and Starck et al. (2006a). As shown in Fig. 1, these generating functions yield the same number of wavelet scales (for the parameter choices described previously). However, with the scalediscretised and needlet generating functions the jth wavelet scale has compact support in [λ^{j−1},λ^{j+1}], whereas the support is much wider with the Bsplines, i.e. [0,L / λ^{J−j−2}] in the S2LET implementation. As a consequence, when using the multiresolution algorithm the wavelet coefficients must be captured on a greater number of pixels than with the scalediscretised or needlet kernels, while probing approximately the same scales, as shown in Fig. 1.
The complexity of the axisymmetric wavelet transform is dominated by spherical harmonic transforms since the wavelet transforms are computed efficiently in harmonic space, through Eqs. (4) and (6) for the forward transform and through Eq. (8) for the inverse transform. Given a bandlimit L and wavelet parameters (λ,J_{0}), recall that the maximum scale is given by J = ⌈ log _{λ} (L − 1) ⌉ and hence the wavelet transform (forward or inverse) involves (J − J_{0} + 3) spherical harmonic transforms (one for the original signal, one for the scaling coefficients and (J − J_{0} + 1) for the wavelet coefficients). If the scaling coefficients and all wavelet coefficients are reconstructed at fullresolution in real space, the axisymmetric wavelet transform scales as . However, in the previous section we established a multiresolution algorithm that takes advantage of the reduced bandlimit of the wavelets for scales j < J − 1. With the multiresolution algorithm with a sampling theorem, only the finest wavelet scales j ∈ { J − 1,J } are computed at maximal resolution corresponding to the bandlimit of the signal. The complexity of the overall multiresolution wavelet transform is then dominated by these operations and effectively scales as .
Fig. 4 Numerical accuracy and computation time of the scalediscretised wavelet transform computed with S2LET. We consider L = 2^{i} with i ∈ { 2,...,10 }, with parameters λ = 2, J_{0} = 0. These results are averaged over many realisations of random bandlimited signals and were found to be very stable. The scalediscretised transform is either performed at fullresolution (solid lines) or with the multiresolution algorithm (dashed lines). Very good numerical accuracy is achieved by both the fullresolution and multiresolution algorithms (which achieve indistinguishable accuracy), with numerical errors comparable to floatingpoint precision, found empirically to scale as as shown by the red line in panel a). Computation time scales as for both algorithms as shown by the red line in panel b), in agreement with theory. The multiresolution algorithm is four to five times faster than the fullresolution approach for the bandlimits considered. 

Open with DEXTER 
3.3. Numerical validation
We first evaluate the performance of S2LET in terms of accuracy and complexity using the MW sampling theorem, for which all transforms are theoretically exact. We show that S2LET achieves floatingpoint precision and scales as detailed in the previous section.
We consider bandlimits L = 2^{i} with i ∈ { 2,...,10 } and generate sets of spherical harmonic coefficients f_{ℓm} following independent Gaussian distributions . We then perform the wavelet decomposition and reconstruct the harmonic coefficients, denoted by . We evaluate the accuracy of the transform using the error metric , which is theoretically zero since all signals are bandlimited by construction. The complexity is quantified by observing how the computation time t_{c} = [t_{synthesis} + t_{analysis}]/2 scales with bandlimit, where the synthesis and analysis computation times are specified by t_{synthesis} and t_{analysis} respectively. Since we evaluate the wavelet transform in real space, a preliminary step is required to reconstruct the signal f from the randomly generated f_{ℓm}. This step is not included in the computation time since its only purpose is to generate a valid bandlimited test signal on the sphere. The analysis then denotes the decomposition of f into wavelet coefficients W^{Ψj} and scaling coefficients W^{Φ} on the sphere. The synthesis refers to recovering the signal f^{rec} from these coefficients. The final step, which is not included in the computation time either, is to decompose f^{rec} into harmonic coefficients in order to compare them with f_{ℓm}. The stability of both ϵ and t_{c} is checked by averaging over hundreds of realisations of f_{ℓm} for L = 2^{i} with i ∈ { 2,...,8 } and a few realisations with i ∈ { 9,10 }. The results proved to be very stable, i.e. the variances of the error and timing metrics are lower than 5%. Recall that for given bandlimit L the number of samples on the sphere required by the exact quadrature is (2L − 1)(L − 1) + 1. All tests were run on an Intel 2.0 GHz Core i7 processor with 8GB of RAM. On this machine, precision of floating point numbers is of the order of ~10^{16}, and errors are expected to add up and accumulate when considering linear operations such as the spherical harmonic and wavelet transforms. The accuracy and timing performance of the scalediscretised wavelet transform implemented in S2LET are presented in Fig. 4. S2LET achieves very good numerical accuracy, with numerical errors comparable to accumulated floatingpoint errors only^{12}. Moreover, the fullresolution and multiresolution algorithms are indistinguishable in terms of accuracy. However, the latter is four to five times faster than the former for the bandlimits considered since only the wavelet coefficients for j ∈ { J − 1,J } are computed at fullresolution. As shown in Fig. 4, computation time scales as for both algorithms, in agreement with theory.
S2LET can also be used with HEALPix, in which case the accuracy of the spherical harmonic transform is critical to the accuracy of the wavelet transform (since HEALpix does not rely on a sampling theorem it does not exhibit theoretically exact harmonic transforms, unlike SSHT or GLESP). The performances of the spherical harmonic transforms in HEALPix and GLESP have been widely studied in the past (see, e.g., Reinecke 2011; Doroshkevich et al. 2011; Reinecke & Seljebotn 2013), and that of the MW sampling were presented in McEwen & Wiaux (2011). We do not compile the entirety of these results here, but we have reproduced the essential results on our machine; Table 1 summarises the orders of accuracy of the HEALPix iterative spherical harmonic transform. Using the same setup as previously, we calculated the maximum error on the spherical harmonic coefficients when performing the transform back and forth, averaged over the values of N_{side}, since the results were found to be sensitive only to the ratio L / N_{side}. Even with several iterations, which multiplies the number of transforms and thus computation time, the spherical harmonic transform in HEALPix remains at least an order of magnitude less accurate than the MW and GLESP counterparts (which, being both theoretically exact, achieve comparable performances, see Reinecke 2011; Reinecke & Seljebotn 2013; Doroshkevich et al. 2011; McEwen & Wiaux 2011). Since the wavelet transforms implemented in MRS and Needatool are also computed in harmonic space, their complexity and accuracy are dominated by that of the underlying spherical harmonic transforms. As a consequence, when adopting the HEALPix scheme, S2LET, MRS and Needatool achieve similar performances, resulting from the computation time and accumulated errors of (J − J_{0} + 1) HEALPix spherical harmonic transforms. In the multiresolution case, the results depend on the resolution chosen to reconstruct each wavelet scale.
Order of magnitude of the accuracy of the HEALPix spherical harmonic transform, averaged over the parameter N_{side}.
3.4. Future extensions
In future work we plan to extend S2LET to support directional, steerable wavelets on the sphere (Wiaux et al. 2008a). We also plan to exploit recent ideas leading to fast (spin) spherical harmonic transforms (McEwen & Wiaux 2011) to yield faster algorithms than those developed by Wiaux et al. (2008a) and McEwen et al. (2013) to compute directional wavelet transforms on the sphere. Finally, we intend to add support to analyse spin signals on the sphere (cf. Geller et al. 2008; Geller & Marinucci 2011; Starck et al. 2009). In a future release, the code will also be parallelised, which will lead to further speed improvements. The S2LET code will thus be under active development with future releases forthcoming. In any case, we hope this first version of the S2LET code will prove useful for axisymmetric scalediscretised wavelet analysis on the sphere. Indeed, the code has already been used as an integral part of the new exact flaglet wavelet transform on the ball (Leistedt & McEwen 2012), the spherical space constructed by augmenting the sphere with the radial line.
4. Examples
The S2LET code is extensively documented and ships with several examples in the four languages supported. In this section we present a subset of short examples, along with the code to execute them in order to demonstrate the ease of using S2LET to perform wavelet transforms^{13}. All examples were run with the scalediscretised wavelet generating functions.
4.1. Wavelet transform from the command line
S2LET includes readytouse highlevel programs to directly decompose a real signal into wavelet coefficients. The inputs are a FITS file containing the signal of interest and the parameters for the transform. The program writes the output coefficients in FITS files in the same directory as the input file and with a consistent naming scheme. These commands are available for both HEALPix and MW sampling schemes. For the MW sampling case illustrated in Example 1, the wavelet transform is theoretically exact and the band limit corresponds to the resolution of the input map, which will be read automatically from the file. The transform may be performed in fullresolution or multiresolution by adjusting the multiresolution flag specified by the last parameter (respectively 0 and 1), and the output wavelet coefficients are computed at full and minimal resolution accordingly. For the case of a HEALPix map, as illustrated in Example 2, the bandlimit must be supplied as the last parameter in the command. The output scaling and wavelet coefficients of a HEALPix map are reconstructed and stored in FITS files at the same resolution as the input map. For both MW and HEALPix samplings the output coefficients may be read and plotted using the Matlab or IDL routines.
Example 1. Performing the forward (analysis) and inverse (synthesis) wavelet transform of a real signal (MW sampling) from the command line.
Example 2. Performing the forward (analysis) and inverse (synthesis) wavelet transform of a real signal (HEALPix sampling) from the command line.
4.2. Wavelet transform in Matlab and IDL
Examples 3 and 4 read real signals on the sphere from FITS files, calculate the wavelet coefficients and plot them using a Mollweide projection. The first case is a Matlab example where the input map is a simulation of the cosmic microwave background in the HEALPix sampling. The second case is a IDL example where the input map is a topography map of the Earth in MW sampling. S2LET ships with versions of these two examples in C, Matlab and IDL.
Example 3. Performing the wavelet transform of a real signal (HEALPix sampling) using the Matlab interface.
Example 4. Performing the wavelet transform of a real signal (MW sampling) using the IDL interface.
4.3. Wavelet denoising in C
Example 5 illustrates the use of the wavelet transform to denoise a signal on the sphere. The input noisy map is a bandlimited topography map of the Earth in MW sampling at resolution L = 128. It is read from a FITS file, decomposed into wavelet coefficients (for given parameters λ and J_{0}) which are then denoised by thresholding. The denoised signal is reconstructed from the denoised wavelet coefficients and written to a FITS file.
In this example we consider a noisy signal y = s + n ∈ L^{2}(S^{2}), where the signal of interest s ∈ L^{2}(S^{2}) is contaminated with noise n ∈ L^{2}(S^{2}). We consider zeromean white Gaussian noise on the sphere, where the variance of the harmonic coefficients of the noise is specified by (19)A simple way to evaluate the fidelity of the observed signal y is through the signaltonoise ratio (S/N), define on the sphere by (20)where the signal energy is defined by (21)We seek a denoised version of y, denoted by d ∈ L^{2}(S^{2}), with large S / N(d) so that d isolates the informative signal s. When taking the wavelet transform of the noisy signal y, one expects the energy of the informative part to be concentrated in a small number of wavelet coefficients, whereas the noise energy should be spread over various wavelet scales. In this particular toy example, the signal has significant power on large scales, as shown in Fig. 3, which are well described in the wavelet basis and less affected by the random white noise. Since the transform is linear, the wavelet coefficients of the jth scale are simply given by the sum of the individual contributions: (22)where capital letters denote the wavelet coefficients, i.e. Y^{j} ≡ y ⋆ Ψ^{j}, S^{j} ≡ s ⋆ Ψ^{j} and N^{j} ≡ n ⋆ Ψ^{j}. For the zeromean white Gaussian noise defined by Eq. (19), the noise in wavelet space is also zeromean and Gaussian, with variance Denoising is performed by hardthresholding the wavelet coefficients Y^{j}, where the threshold is taken as T^{j} = 3σ^{j}. The denoised wavelet coefficients D^{j} ≡ d ⋆ Ψ^{j} are thus given by (23)The denoised signal d ∈ L^{2}(S^{2}) is reconstructed from its wavelet coefficients D^{j} and the scaling coefficients of y, which are not thresholded. The denoising procedure outlined above is particularly simple and more sophisticated denoising strategies can be developed; we adopt this simple denoising strategy merely to illustrate the use of the S2LET code. In this example we perform the wavelet transform with parameters λ = 2 and J_{0} = 0. For a noisy signal y with S / N(y) = 11.78 dB, the scalediscretised wavelet denoising recovers a denoised signal d with S / N(d) = 14.66 dB. The initial, noisy and denoised maps are shown in Fig. 5. When switching to needlets and Bspline wavelets while keeping λ and J_{0} unchanged, the denoised signals have S / N(d) = 14.68 dB and 14.46 dB, respectively.
Fig. 5 Wavelet denoising by hardthresholding, using parameters λ = 2 and J_{0} = 0 and scalediscretised generating functions. When using needlets and Bspline wavelets, the denoised signals have S / N(d) = 14.68 dB and 14.46 dB respectively. This example is included in S2LET as a documented demo program. 

Open with DEXTER 
Example 5. Denoising a real signal (MW sampling) in C through hardthresholding of the wavelet coefficients.
5. Summary
In the era of precision astrophysics and cosmology, large and complex datasets on the sphere must be analysed at high precision in order to confront accurate theoretical predictions. Scalediscretised wavelets are a powerful analysis technique where spatially localised, scaledependent signal features of interest can be extracted and analysed. Combined with a sampling theorem, this framework leads to an exact multiresolution wavelet analysis, where signals on the sphere can be reconstructed from their scaling and wavelet coefficients exactly.
We have described S2LET, a fast and robust implementation of the scalediscretised wavelet transform. Although the first public release of S2LET is restricted to axisymmetric wavelets, the generalisation to directional, steerable wavelets will be made available in a future release. The core numerical routines of S2LET are written in C and have interfaces in Matlab, IDL and Java. Both MW and HEALPix pixelisation schemes are supported. In this article we have presented a number of examples to illustrate the ease of use of S2LET for performing wavelet transform of real signals stored as FITS files and to plot scaling and wavelet coefficients on Mollweide projections of the sphere. We have also detailed a denoising example where denoising is performed through simple hardthresholding in wavelet space. Although only a simple denoising strategy was performed to illustrate the use of the S2LET code, it nevertheless performed very well, highlighting the effectiveness of the scalediscretised wavelet transform on the sphere.
As already noted, the extension to directional scalediscretised wavelets has been derived by Wiaux et al. (2008a). At present the S2LET code supports axisymmetric wavelets only; directional wavelets will be added in a later release.
In our implementation of needlets we use a scaling function to represent the approximate information in the signal, which is not always included (e.g., NeedAtool; Pietrobon et al. 2010).
For the Bsplinebased construction to probe approximately the same scales as the scalediscretised and needlet ones, we defined the generating functions as
so that the jth filter has (compact) support [0,L / λ^{J−j−2}] and peaks at the same scales as the jth scalediscretised and needlet filters obtained with the same parameters.
Acknowledgments
B.L. is supported by the Perren Fund and the Impact Fund. J.D.M. is supported in part by a Newton International Fellowship from the Royal Society and the British Academy. Y.W. is supported by the Center for Biomedical Imaging (CIBM) of the Geneva and Lausanne Universities, EPFL and the Leenaards and LouisJeantet foundations.
References
 Antoine, J.P., & Vandergheynst, P. 1998, J. Math. Phys., 39, 3987 [NASA ADS] [CrossRef] [Google Scholar]
 Antoine, J.P., & Vandergheynst, P. 1999, Appl. Comput. Harm. Anal., 7, 1 [CrossRef] [Google Scholar]
 Baldi, P., Kerkyacharian, G., Marinucci, D., & Picard, D. 2009, Ann. Stat., 37, 1150 [NASA ADS] [CrossRef] [Google Scholar]
 Barreiro, R. B., Hobson, M. P., Lasenby, A. N., et al. 2000, MNRAS, 318, 475 [NASA ADS] [CrossRef] [Google Scholar]
 Basak, S., & Delabrouille, J. 2012, MNRAS, 419, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Cayón, L., Sanz, J. L., MartínezGonzález, E., et al. 2001, MNRAS, 326, 1243 [NASA ADS] [CrossRef] [Google Scholar]
 Daubechies, I. 1992, Ten Lectures on Wavelets, Soc. Ind. Appl. Math. [Google Scholar]
 Deriaz, E., Starck, J.L., & Pires, S. 2012, A&A, 540, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doroshkevich, A. G., Naselsky, P. D., Verkhodanov, O. V., et al. 2005, Int. J. Mod. Phys. D, 14, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Doroshkevich, A. G., Verkhodanov, O. V., Naselsky, P. D., et al. 2011, Int. J. Mod. Phys. D, 20, 1053 [NASA ADS] [CrossRef] [Google Scholar]
 Driscoll, J. R., & Healy, D. M. J. 1994, Adv. Appl. Math., 15, 202 [CrossRef] [Google Scholar]
 Geller, D., & Marinucci, D. 2011, J. Math. Anal. Appl., 375, 610 [CrossRef] [Google Scholar]
 Geller, D., Hansen, F. K., Marinucci, D., Kerkyacharian, G., & Picard, D. 2008, Phys. Rev. D., 78, 123533 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Healy, D., Jr., Rockmore, D., Kostelec, P. J., & Moore, S. S. B. 1996, The J. Fourier Anal. Appl., 9, 341 [CrossRef] [Google Scholar]
 Labatie, A., Starck, J. L., & LachièzeRey, M. 2012, AJ, 746, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Lan, X., & Marinucci, D. 2008, Electron. J. Stat., 2, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Leistedt, B., & McEwen, J. D. 2012, IEEE Trans. Sig. Proc., 60 [Google Scholar]
 Marinucci, D., Pietrobon, D., Balbi, A., et al. 2008, MNRAS, 383, 539 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., & Wiaux, Y. 2011, IEEE Trans. Sig. Proc., 59, 5876 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Hobson, M. P., & Lasenby, A. N. 2006a, unpublished [arXiv:astroph/0609159] [Google Scholar]
 McEwen, J. D., Hobson, M. P., Lasenby, A. N., & Mortlock, D. J. 2006b, MNRAS, 371, L50 [NASA ADS] [Google Scholar]
 McEwen, J. D., Vielva, P., Hobson, M. P., MartínezGonzález, E., & Lasenby, A. N. 2007a, MNRAS, 376, 1211 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Vielva, P., Wiaux, Y., et al. 2007b, J. Fourier Anal. Appl., 13, 495 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Wiaux, Y., Hobson, M. P., Vandergheynst, P., & Lasenby, A. N. 2008, MNRAS, 384, 1289 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Vandergheynst, P., & Wiaux, Y. 2013, in Wavelets and Sparsity XIV, SPIE international symposium on optics and photonics (invited contribution) [Google Scholar]
 Narcowich, F. J., Petrushev, P., & Ward, J. D. 2006, SIAM J. Math. Anal., 38, 574 [CrossRef] [Google Scholar]
 Pietrobon, D., Amblard, A., Balbi, A., et al. 2008, Phys. Rev. D., 78, 103504 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Balbi, A., Cabella, P., & Gorski, K. M. 2010, ApJ, 723, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Reinecke, M. 2011, A&A, 526, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Reinecke, M., & Seljebotn, D. S. 2013, A&A, 554, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Risbo, T. 1996, J. Geodesy, 70, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Schmitt, J., Starck, J. L., Casandjian, J. M., Fadili, J., & Grenier, I. 2010, A&A, 517, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Schröder, P., & Sweldens, W. 1995, Computer Graphics Proceedings (SIGGRAPH 95), 161 [Google Scholar]
 Starck, J.L., Moudden, Y., Abrial, P., & Nguyen, M. 2006a, A&A, 446, 1191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J.L., Pires, S., & Réfrégier, A. 2006b, A&A, 451, 1139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J.L., Moudden, Y., & Bobin, J. 2009, A&A, 497, 931 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sweldens, W. 1996, Appl. Comput. Harmon. Anal., 3, 186 [CrossRef] [MathSciNet] [Google Scholar]
 Sweldens, W. 1997, SIAM J. Math. Anal., 29, 511 [CrossRef] [Google Scholar]
 Trapani, S., & Navaza, J. 2006, Acta Crystallog. Sect. A, 62, 262 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., MartínezGonzález, E., Barreiro, R. B., Sanz, J. L., & Cayón, L. 2004, ApJ, 609, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., MartínezGonzález, E., & Tucci, M. 2006a, MNRAS, 365, 891 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., Wiaux, Y., MartínezGonzález, E., & Vandergheynst, P. 2006b, New Astron. Rev., 50, 880 [NASA ADS] [CrossRef] [Google Scholar]
 Vielva, P., Wiaux, Y., MartínezGonzález, E., & Vandergheynst, P. 2007, MNRAS, 381, 932 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., Jacques, L., & Vandergheynst, P. 2005, ApJ, 632, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., Jacques, L., Vielva, P., & Vandergheynst, P. 2006a, ApJ, 652, 820 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., Vielva, P., MartínezGonzález, E., & Vandergheynst, P. 2006b, Phys. Rev. Lett., 96, 151303 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Wiaux, Y., McEwen, J. D., & Vielva, P. 2007, J. Fourier Anal. Appl., 13, 477 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., McEwen, J. D., Vandergheynst, P., & Blanc, O. 2008a, MNRAS, 388, 770 [NASA ADS] [CrossRef] [Google Scholar]
 Wiaux, Y., Vielva, P., Barreiro, R. B., MartínezGonzález, E., & Vandergheynst, P. 2008b, MNRAS, 385, 939 [NASA ADS] [CrossRef] [Google Scholar]
 Yeo, B., Ou, W., & Golland, P. 2008, IEEE Transactions on Image Processing, 17, 283 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Order of magnitude of the accuracy of the HEALPix spherical harmonic transform, averaged over the parameter N_{side}.
All Figures
Fig. 1 Wavelets and scaling function constructed with the scalediscretised (SD), needlet and Bspline generating functions (Wiaux et al. 2008a; Marinucci et al. 2008; Starck et al. 2006a) with parameters λ = 3 and J_{0} = 2 and for bandlimit L = 128. The tiling is shown in the top panel, and the profiles of the reconstructed wavelets in the bottom panel. 

Open with DEXTER  
In the text 
Fig. 2 Wavelets for scales j ∈ { 2,3,4,5 } and scaling function constructed through a tiling of the harmonic line using scalediscretised functions, with parameters λ = 3 and J_{0} = 2 and for bandlimit L = 128. This plot was produced with a Matlab demo included in S2LET. 

Open with DEXTER  
In the text 
Fig. 3 Scalediscretised wavelet transform of a bandlimited topography map of the Earth for λ = 3, J_{0} = 2 and L = 128, i.e. with the scalediscretised wavelets shown in Fig. 2. The wavelet transform decomposes the bandlimited signal into wavelet coefficients that extract spatially localised, scaledependent features. Since the wavelets for different scales j have different bandlimits, the wavelet coefficients can be reconstructed at lower resolution on the sphere for lower scales j. Panel a) shows the fullresolution wavelet transform of the topography map. The original Earth topography map is shown in the topleft plot, the scaling coefficients are shown in the topright plot, while the wavelet coefficients at scales j ∈ { 2,3,4,5 } are shown lefttoright, toptobottom respectively in the remaining plots. Panel b) shows the same decomposition but using the multiresolution algorithm. The signals shown in panel b) contain the same information as in panel a) but represented in the minimal number of samples on the sphere. These plots were produced by one of the many Matlab demos provided with S2LET. 

Open with DEXTER  
In the text 
Fig. 4 Numerical accuracy and computation time of the scalediscretised wavelet transform computed with S2LET. We consider L = 2^{i} with i ∈ { 2,...,10 }, with parameters λ = 2, J_{0} = 0. These results are averaged over many realisations of random bandlimited signals and were found to be very stable. The scalediscretised transform is either performed at fullresolution (solid lines) or with the multiresolution algorithm (dashed lines). Very good numerical accuracy is achieved by both the fullresolution and multiresolution algorithms (which achieve indistinguishable accuracy), with numerical errors comparable to floatingpoint precision, found empirically to scale as as shown by the red line in panel a). Computation time scales as for both algorithms as shown by the red line in panel b), in agreement with theory. The multiresolution algorithm is four to five times faster than the fullresolution approach for the bandlimits considered. 

Open with DEXTER  
In the text 
Fig. 5 Wavelet denoising by hardthresholding, using parameters λ = 2 and J_{0} = 0 and scalediscretised generating functions. When using needlets and Bspline wavelets, the denoised signals have S / N(d) = 14.68 dB and 14.46 dB respectively. This example is included in S2LET as a documented demo program. 

Open with DEXTER  
In the text 