EDP Sciences
Free Access
Issue
A&A
Volume 558, October 2013
Article Number A128
Number of page(s) 9
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201220729
Published online 18 October 2013

© 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 data-sets 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, scale-discretised 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 scale-discretised 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 B-splines 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 S2LET1 code to perform the scale-discretised 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 scale-discretised wavelets (Wiaux et al. 2008a), needlets (Narcowich et al. 2006; Baldi et al. 2009; Marinucci et al. 2008) and B-spline 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 transforms2 (Starck et al. 2006a), invertible filter banks3 (Yeo et al. 2008), needlets (NeedATool4; Pietrobon et al. 2010) and scale-discretised wavelets (S2DW5; Wiaux et al. 2008a). S2LET aims primarily to provide a fast and flexible implementation of the scale-discretised 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 user-friendly 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 scale-discretised axisymmetric wavelets and the corresponding exact scale-discretised 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 scale-discretised 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 band-limited 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 ∈ L2(S2) on the two-dimensional sphere S2 reads (1)where Yℓm are the spherical harmonic functions, which form the canonical orthogonal basis on S2. 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 ω = (θ,φ) ∈ S2 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 band-limited signals in the spherical harmonic basis, with band-limit L if fℓm = 0,  ∀ ≥ L. For band-limited 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 band-limit L, requires the lowest number of samples on the sphere of all sampling theorems, namely (L − 1)(2L − 1) + 1 ~ 2L2 samples (for comparison ~4L2 samples are required by Driscoll & Healy 1994). Fast algorithms to compute the corresponding spherical harmonic transform scale as and are numerically stable to band-limits of at least L = 4096 (McEwen & Wiaux 2011). The GLESP pixelisation scheme (Doroshkevich et al. 2005) also provides a sampling theorem based on the Gauss-Legendre quadrature, and could be used in place of the MW sampling theorem. However, GLESP uses more samples than Gauss-Legendre quadrature requires, which may lead to an overhead when considering large band-limits 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 equal-area pixels.

2.2. Scale-discretised wavelets on the sphere

The scale-discretised wavelet transform allows one to probe spatially localised, scale-dependent content in the signal of interest f ∈ L2(S2). The jth wavelet scale WΨj ∈ L2(S2) is defined as the convolution of f with the wavelet Ψj ∈ L2(S2): (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 orientation6. 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 =  ⟨ f|Yℓ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 low-frequency (low-), approximate information of the signal. The scaling coefficients WΦ ∈ L2(S2) are defined by the convolution of f with the scaling function Φ ∈ L2(S2): (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 J0, 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 band-limited 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 J0 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 band-limit of the signal of interest, yielding J =  ⌈ log λ   (L − 1) ⌉. The choice of the lowest wavelet scale J0 is arbitrary, provided that 0 ≤ J0 < 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.

thumbnail Fig. 1

Wavelets and scaling function constructed with the scale-discretised (SD), needlet and B-spline generating functions (Wiaux et al. 2008a; Marinucci et al. 2008; Starck et al. 2006a) with parameters λ = 3 and J0 = 2 and for band-limit 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 scale-discretised 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 B-spline filters used to construct the isotropic undecimated wavelet transform (Starck et al. 2006a) are also supported, as also shown in Fig. 18. With these three constructions, the wavelets and scaling functions are well-localised 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, scale-dependent features in signals of interest.

thumbnail Fig. 2

Wavelets for scales j ∈  { 2,3,4,5 } and scaling function constructed through a tiling of the harmonic line using scale-discretised functions, with parameters λ = 3 and J0 = 2 and for band-limit 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 band-limit 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 scale-discretised 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 band-limits 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 band-limit k = λj + 1 when using the scale-discretised and needlet kernels, and k = L / λJj−2 when using the B-splines. 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 full-resolution case, as shown in the next section. Figure 3 illustrates the use of the full-resolution and multiresolution transforms on a map of Earth topography data with the scale-discretised 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.

thumbnail Fig. 3

Scale-discretised wavelet transform of a band-limited topography map of the Earth for λ = 3, J0 = 2 and L = 128, i.e. with the scale-discretised wavelets shown in Fig. 2. The wavelet transform decomposes the band-limited signal into wavelet coefficients that extract spatially localised, scale-dependent features. Since the wavelets for different scales j have different band-limits, the wavelet coefficients can be reconstructed at lower resolution on the sphere for lower scales j. Panel a) shows the full-resolution wavelet transform of the topography map. The original Earth topography map is shown in the top-left plot, the scaling coefficients are shown in the top-right plot, while the wavelet coefficients at scales j ∈  { 2,3,4,5 } are shown left-to-right, top-to-bottom 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 full-resolution 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,λ,J0). 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 SSHT9 code. In this case all transforms are theoretically exact and one can analyse and synthesise real and complex signals at floating-point 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 FITS10 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 object-oriented 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 FFTW11 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 scale-discretised wavelet transform of Wiaux et al. (2008a), S2LET also supports the needlet and spline-based 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 scale-discretised and needlet generating functions the jth wavelet scale has compact support in [λj−1,λj+1], whereas the support is much wider with the B-splines, i.e. [0,L / λJj−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 scale-discretised 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 band-limit L and wavelet parameters (λ,J0), recall that the maximum scale is given by J =  ⌈ log λ   (L − 1) ⌉ and hence the wavelet transform (forward or inverse) involves (J − J0 + 3) spherical harmonic transforms (one for the original signal, one for the scaling coefficients and (J − J0 + 1) for the wavelet coefficients). If the scaling coefficients and all wavelet coefficients are reconstructed at full-resolution 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 band-limit 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 band-limit of the signal. The complexity of the overall multiresolution wavelet transform is then dominated by these operations and effectively scales as .

thumbnail Fig. 4

Numerical accuracy and computation time of the scale-discretised wavelet transform computed with S2LET. We consider L = 2i with i ∈  { 2,...,10 }, with parameters λ = 2, J0 = 0. These results are averaged over many realisations of random band-limited signals and were found to be very stable. The scale-discretised transform is either performed at full-resolution (solid lines) or with the multiresolution algorithm (dashed lines). Very good numerical accuracy is achieved by both the full-resolution and multiresolution algorithms (which achieve indistinguishable accuracy), with numerical errors comparable to floating-point 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 full-resolution approach for the band-limits 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 floating-point precision and scales as detailed in the previous section.

We consider band-limits L = 2i 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 band-limited by construction. The complexity is quantified by observing how the computation time tc =  [tsynthesis + tanalysis]/2 scales with band-limit, where the synthesis and analysis computation times are specified by tsynthesis and tanalysis 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 band-limited 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 frec from these coefficients. The final step, which is not included in the computation time either, is to decompose frec into harmonic coefficients in order to compare them with fℓm. The stability of both ϵ and tc is checked by averaging over hundreds of realisations of fℓm for L = 2i 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 band-limit 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 scale-discretised wavelet transform implemented in S2LET are presented in Fig. 4. S2LET achieves very good numerical accuracy, with numerical errors comparable to accumulated floating-point errors only12. Moreover, the full-resolution and multiresolution algorithms are indistinguishable in terms of accuracy. However, the latter is four to five times faster than the former for the band-limits considered since only the wavelet coefficients for j ∈  { J − 1,J } are computed at full-resolution. 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 Nside, since the results were found to be sensitive only to the ratio L / Nside. 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 − J0 + 1) HEALPix spherical harmonic transforms. In the multiresolution case, the results depend on the resolution chosen to reconstruct each wavelet scale.

Table 1

Order of magnitude of the accuracy of the HEALPix spherical harmonic transform, averaged over the parameter Nside.

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 scale-discretised 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 transforms13. All examples were run with the scale-discretised wavelet generating functions.

4.1. Wavelet transform from the command line

S2LET includes ready-to-use high-level 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 full-resolution 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 band-limit 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 band-limited 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 J0) 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 ∈ L2(S2), where the signal of interest s ∈ L2(S2) is contaminated with noise n ∈ L2(S2). We consider zero-mean 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 signal-to-noise 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 ∈ L2(S2), 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. Yj ≡ y  Ψj, Sj ≡ s  Ψj and Nj ≡ n  Ψj. For the zero-mean white Gaussian noise defined by Eq. (19), the noise in wavelet space is also zero-mean and Gaussian, with variance Denoising is performed by hard-thresholding the wavelet coefficients Yj, where the threshold is taken as Tj = 3σj. The denoised wavelet coefficients Dj ≡ d  Ψj are thus given by (23)The denoised signal d ∈ L2(S2) is reconstructed from its wavelet coefficients Dj 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 J0 = 0. For a noisy signal y with S / N(y) = 11.78 dB, the scale-discretised 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 B-spline wavelets while keeping λ and J0 unchanged, the denoised signals have S / N(d) = 14.68 dB and 14.46 dB, respectively.

thumbnail Fig. 5

Wavelet denoising by hard-thresholding, using parameters λ = 2 and J0 = 0 and scale-discretised generating functions. When using needlets and B-spline 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 hard-thresholding of the wavelet coefficients.

5. Summary

In the era of precision astrophysics and cosmology, large and complex data-sets on the sphere must be analysed at high precision in order to confront accurate theoretical predictions. Scale-discretised wavelets are a powerful analysis technique where spatially localised, scale-dependent 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 scale-discretised 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 hard-thresholding 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 scale-discretised wavelet transform on the sphere.


6

As already noted, the extension to directional scale-discretised 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.

7

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).

8

For the B-spline-based construction to probe approximately the same scales as the scale-discretised and needlet ones, we defined the generating functions as

so that the jth filter has (compact) support [0,L / λJj−2] and peaks at the same scales as the jth scale-discretised and needlet filters obtained with the same parameters.

12

The GLESP sampling adopted in MRS also achieves floating point accuracy, although using many more pixels to capture the wavelet scales due to the greater band-limits of the spline-based kernels and the oversampling of the GLESP scheme.

13

Note that the code uses a slightly different notation compared to the equations of this article: B refers to the wavelet scaling parameter (denoted λ herein) and Jmin to the first scale of the transform (denoted J0 herein).

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 Louis-Jeantet foundations.

References

All Tables

Table 1

Order of magnitude of the accuracy of the HEALPix spherical harmonic transform, averaged over the parameter Nside.

All Figures

thumbnail Fig. 1

Wavelets and scaling function constructed with the scale-discretised (SD), needlet and B-spline generating functions (Wiaux et al. 2008a; Marinucci et al. 2008; Starck et al. 2006a) with parameters λ = 3 and J0 = 2 and for band-limit 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
thumbnail Fig. 2

Wavelets for scales j ∈  { 2,3,4,5 } and scaling function constructed through a tiling of the harmonic line using scale-discretised functions, with parameters λ = 3 and J0 = 2 and for band-limit L = 128. This plot was produced with a Matlab demo included in S2LET.

Open with DEXTER
In the text
thumbnail Fig. 3

Scale-discretised wavelet transform of a band-limited topography map of the Earth for λ = 3, J0 = 2 and L = 128, i.e. with the scale-discretised wavelets shown in Fig. 2. The wavelet transform decomposes the band-limited signal into wavelet coefficients that extract spatially localised, scale-dependent features. Since the wavelets for different scales j have different band-limits, the wavelet coefficients can be reconstructed at lower resolution on the sphere for lower scales j. Panel a) shows the full-resolution wavelet transform of the topography map. The original Earth topography map is shown in the top-left plot, the scaling coefficients are shown in the top-right plot, while the wavelet coefficients at scales j ∈  { 2,3,4,5 } are shown left-to-right, top-to-bottom 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
thumbnail Fig. 4

Numerical accuracy and computation time of the scale-discretised wavelet transform computed with S2LET. We consider L = 2i with i ∈  { 2,...,10 }, with parameters λ = 2, J0 = 0. These results are averaged over many realisations of random band-limited signals and were found to be very stable. The scale-discretised transform is either performed at full-resolution (solid lines) or with the multiresolution algorithm (dashed lines). Very good numerical accuracy is achieved by both the full-resolution and multiresolution algorithms (which achieve indistinguishable accuracy), with numerical errors comparable to floating-point 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 full-resolution approach for the band-limits considered.

Open with DEXTER
In the text
thumbnail Fig. 5

Wavelet denoising by hard-thresholding, using parameters λ = 2 and J0 = 0 and scale-discretised generating functions. When using needlets and B-spline 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.