---

Polarized wavelets and curvelets on the sphere

J.-L. Starck - Y. Moudden - J. Bobin

Laboratoire AIM (UMR 7158), CEA/DSM-CNRS-Université Paris Diderot, IRFU, SEDI-SAP, Service d'Astrophysique, Centre de Saclay, 91191 Gif-Sur-Yvette Cedex, France

Received 13 November 2008 / Accepted 16 January 2009

Abstract
The statistics of the temperature anisotropies in the primordial cosmic microwave background radiation field provide a wealth of information for cosmology and for estimating cosmological parameters. An even more acute inference should stem from the study of maps of the polarization state of the CMB radiation. Measuring the extremely weak CMB polarization signal requires very sensitive instruments. The full-sky maps of both temperature and polarization anisotropies of the CMB to be delivered by the upcoming Planck Surveyor satellite experiment are hence being awaited with excitement. Multiscale methods, such as isotropic wavelets, steerable wavelets, or curvelets, have been proposed in the past to analyze the CMB temperature map. In this paper, we contribute to enlarging the set of available transforms for polarized data on the sphere. We describe a set of new multiscale decompositions for polarized data on the sphere, including decimated and undecimated Q-U or E-B wavelet transforms and Q-U or E-B curvelets. The proposed transforms are invertible and so allow for applications in data restoration and denoising.

Key words: cosmology: cosmic microwave background - methods: data analysis - techniques: image processing

1 Introduction

The statistical analysis of the slight intensity fluctuations in the primordial cosmic microwave background radiation field, for which evidence was found for the first time in the early 1990's in the observations made by COBE (Smoot et al. 1992), is a major issue in modern cosmology as these are strongly related to the cosmological scenarios describing the properties and evolution of our Universe. In the Big Bang model, the observed CMB anisotropies are an imprint of primordial fluctuations in baryon-photon density from a time when the temperature of the Universe was high enough above 3000 K for matter and radiation to be tightly coupled. At that time, the attraction of gravity and the repulsive radiation pressure were opposed, thus generating so-called acoustic oscillations in the baryon-photon fluid, causing peaks and troughs to appear in the power spectrum of the spatial anisotropies of the CMB. With the Universe cooling down as it expanded, matter and radiation finally decoupled. Photons were set free in a nearly transparent Universe, while the density fluctuations collapsed under the effect of gravity into large scale structures such as galaxies or clusters of galaxies. Due to the expansion of the Universe, the CMB photons are now observed in the microwave range but are still distributed according to an almost perfect black body emission law. Another major result was the measurement of the polarization state and anisotropies of the CMB radiation field by DASI (Kovac et al. 2002). Only a fraction of the total CMB radiation is polarized so that extremely sensitive instruments are needed. Polarization of the CMB radiation is a consequence of the Thomson scattering of photons on electrons. But for the outgoing population of photons to be polarized, the radiation incident on the scatterer needs to be anisotropic and have a quadrupole moment. The statistics of the CMB polarization anisotropies are also a source of information for cosmology. Inference of cosmological parameters from the joint statistics of the CMB anisotropies should benefit from both the complementarity and the redundancy of the information carried by the additional measurement of CMB polarization. Hence the full-sky maps with unprecedented sensitivity and angular resolution of both temperature and polarization anisotropies of the CMB to be delivered by the upcoming Planck Surveyor satellite experiment are awaited with excitement.

Multiscale transforms on the sphere

Wavelets have been used in many studies of CMB data analysis, especially for non-Gaussianity detection (Vielva et al. 2004,2006; Jin et al. 2005; McEwen et al. 2007; Hobson et al. 1999; Starck et al. 2004a; Aghanim & Forni 1999). Continuous wavelet transforms on the sphere (Tenorio et al. 1999; Antoine 1999; Holschneider 1996; Cayón et al. 2001) are the most used transforms. Recently, Starck et al. (2006) proposed an invertible isotropic undecimated wavelet transform (UWT) on the sphere based on spherical harmonics. A similar wavelet construction has been published in Marinucci et al. (2008), Faÿ & Guilloux (2008) and Faÿ et al. (2008) using so-called ``needlet filters''. Wiaux et al. (2008) have also proposed an algorithm which allows us to reconstruct an image from its steerable wavelet transform. Since reconstruction algorithms are available, these new tools can be used for other applications than non-Gaussianity detection, such as denoising, deconvolution, component separation (Delabrouille et al. 2008; Moudden et al. 2005; Bobin et al. 2008) or inpainting (Abrial et al. 2008,2007). Other multiscale transforms on the sphere such as ridgelets and curvelets have been developed (Starck et al. 2006), and are well adapted to detect anisotropic features. It has been shown that such constructions are very useful for the detection of cosmic strings (Jin et al. 2005; Hammond et al. 2008; Starck et al. 2004a). More generally, each transform can be characterized by a given dictionary , and the transformation consists in representing the input data y as a linear combination of the dictionary elements: , where are the coefficients. In the case of a wavelet transform, the dictionary is composed of wavelet functions, and are the wavelet coefficients. Depending on the content of the data, a given dictionary may be more adapted to detect the signal of interest. A dictionary is considered as well designed for a class of signals if the transformation of any of these signals leads to a sparse representation, i.e. if most of the coefficients are equal or close to zero, while only few coefficients have a high amplitude. If the morphological signature of the features to be detected are not known, we should not restrict ourselves to analyze the data with wavelets or curvelets only, but rather use all available tools (Starck et al. 2004a).

In this paper, we describe a set of new multiscale decompositions for polarized data on the sphere, including decimated and undecimated Q-U or E-B wavelet transforms, and Q-U or E-B curvelets. Choosing the right transform will depend on the application and on the typical structures occurring in the data to be analyzed. Indeed, as illustrated next, the different transforms are associated with different waveform dictionaries. Section 2 presents an orthogonal decomposition for polarized data on the sphere. Section 3 introduces a new decomposition where the modulus and the phase are processed independently. Sections 4 and 5 extend the isotropic wavelet transform and the curvelet transform on the sphere to the case of polarized data, and Sect. 6 introduces two other wavelet and curvelet decompositions which are based on the E/B mode separation. The proposed transforms are invertible and so allow for applications in data restoration and denoising. Section 7 reports on denoising experiments using these new polarized transforms, thus demonstrating their usefulness in practice.

2 Orthogonal representation of polarized data

2.1 Introduction

Full-sky CMB polarization data, as expected from the upcoming Planck experiment, consists of measurements of the Stokes parameters so that in addition to the temperature T map, Q and U maps are given as well. The fourth Stokes parameter commonly denoted V is a measure of circular polarization. In the case of CMB which is not expected to have circularly polarized anisotropies, V vanishes. The former three quantities, T, Q and U then fully describe the linear polarization state of the CMB radiation incident along some radial line of sight: T is the total incoming intensity, Q is the difference between the intensities transmitted by two perfect orthogonal polarizers the directions of which define a reference frame in the tangent plane, and U is the same as Q but with polarizers rotated 45 degrees in that tangent plane. Clearly, Q and U are not invariant through a rotation of angle of the local reference frame around the line of sight. In fact, it is easily shown that:

which can also be written which by definition expresses the fact that the quantities are spin-2 fields on the sphere. The suitable generalization of the Fourier representation for such fields is the spin-2 spherical harmonics basis denoted , in which we can expand:

2.2 Multiscale representation

The easiest way to build a multiscale transform for polarized data is to use the Healpix representation (Górski et al. 2005), and to apply a bi-orthogonal wavelet transform on each face of the Healpix map, separately for Q and U. Figure 1 shows the flow-graph of this Q-U orthogonal wavelet transform (QU-OWT). Recall that the base resolution of the Healpix representation divides the sphere into twelve curvilinear quadrilateral faces of equal area placed on three rings around the poles and equator. Each face is subsequently divided into nside² pixels of exactly equal surface but with varying shape. It follows that Q and U are reconstructed at position k from their wavelet coefficients w_j,p^Q, w_j,p^U, c^Q_J,p and c^U_J,p according to:

which can also be written as:

This wavelet transform is not redundant i.e. the decomposition has the same number of coefficients as the input data, and it is invertible so that the Q and U maps can be reconstructed exactly.

When we apply such a decomposition, we implicitly use a dictionary on which we project the data. As discussed previously, the shape of the dictionary elements, also called atoms, is very important to have an efficient analysis of the data. In the case of polarized data, it is not straightforward to imagine these shapes from Eq. (4). In order to visualize them, we can perform a backprojection i.e. we apply the inverse wavelet transform to sets of wavelet coefficients where only one coefficient is different from zero. Repeating the same experiment, changing only the scale and position of the non-zero coefficient allows us to view the different atoms in the dictionary related to the QU-OWT transform that we use. Figure 2 shows examples of backprojections of Q wavelet coefficients (top) and backprojections of U wavelet coefficients (bottom). The shapes of the individual atoms do not look close to the astronomical patterns we would expect in our data. Therefore, this decomposition may not be optimal to analyze polarized astronomical data, although this would need to be confirmed in practice. The following sections describe other polarized wavelet transforms with different morphological properties.

3 Module-phase non linear multiscale transform

3.1 Introduction

Given a polarized map in the standard Q-U representation, consider a different point of view and define the modulus M and phase P maps as follows:

Because the smoothness of the Q and U maps should result in some smoothness of the modulus map M and the phase map P, one may consider devising a multiscale modulus/phase decomposition of the spin 2 field .

The specificity of the modulus/phase decomposition of is twofold: i) the modulus field is non-negative and ii) the phase field takes its values on the unit circle S¹. Recently, Rahman et al. (2005) introduced a multiscale analysis technique for manifold valued data that will be described in the following paragraph. We then define the modulus/phase (MP) multiscale transform as follows:

1.

apply a classical multiscale transform (i.e. wavelets) to the modulus map M;

2.

apply the multiscale analysis technique for manifold valued data described in Rahman et al. (2005) to the phase map P.

3.2 Decimated MP-multiscale transform

Let us provide some essential notation: we assume that the entries of the phase map P lie in a manifold (e.g. ). According to Rahman et al. (2005), take and define as the log-map of p₁ onto the tangent space of at p₀. The back-projection is obtained using the inverse of the log-map Exp_p0.

For instance, if we choose then and . The and maps are then defined as follows:

The multiscale transform for manifold valued data introduced in Rahman et al. (2005) is equivalent to a two-step interpolation-refinement scheme similar to the lifting scheme described in Daubechies & Sweldens (1998). The wavelet coefficients and low pass approximation pixels are then computed as follows at each scale j and pixel k:

The wavelet coefficient w_j+1,k^P at pixel k and scale j is the projection of its prediction/interpolation onto the tangent space of at c_j,2k+1^P. The low pass approximation c_j+1,k^P at scale j+1 is computed by updating c_j,2k^P from the wavelet coefficient  w_j+1,k^P.

The main advantage of this scheme is its ability to capture local regularities while guaranteeing the low pass approximation to belong to the manifold . Indeed, the wavelet coefficient  w_j+1,k^P at pixel k and scale j+1 is computed as the Exp map at  c_j,2k+1^P of an approximation of c_j,2k^P.

Note also that even if the definitions of the and maps involve the absolute phase (i.e. ), at least they only require the computation of differences of phases values thus avoiding the explicit manipulation of an absolute phase.

However the non-linearity of the proposed transform is a major drawback when considering denoising and restoration applications.

Illustration:

In the case of polarized data, the entries of the phase map P lie in . In the following experiments, and are chosen such that:

This multiscale transform is invertible and its inverse is computed as follows:

The picture in Fig. 3 features some examples of backprojections of MP-multiscale coefficients.

3.3 Undecimated MP-multiscale transform

For image restoration purposes, the use of undecimated multiscale transforms has been shown to provide better results than decimated transforms (Starck & Murtagh 2006; Starck et al. 1998). The aforementioned modulus/phase multiscale analysis can be extended to an undecimated scheme consisting in: i) applying an undecimated wavelet transform to the modulus map; ii) analyzing the phase map P using an extension to the undecimated case of the multiscale transform described in Rahman et al. (2005). In that case, Eqs. (9) and (10) are replaced with the following equations:

where with and h_l > 0. The low pass approximation c_j+1,k^P is then computed from a linear combination (linear filter) of a neighborhood of c_j,k weighted by the positive scalars . Note that from scale j to scale j+1, the spatial size of the neighborhood increases by a factor 2 which would be equivalent to downsize by a factor 2 the band pass filter of the classical wavelet decomposition scheme.

3.4 Example

In the case of polarized data, the entries of the phase map P lie in . In the following experiments, is chosen such that:

where:

This multiscale transform is invertible and its inverse is computed as follows:

Figure 4 top shows a simulated polarized field of the synchrotron emission and its noisy version. We have applied the MP-multiscale transform and we remove the first three scales (i.e. we put all coefficients to zero) before reconstructing. The resulting image is shown on the bottom left of Fig. 4. The bottom right of Fig. 4 corresponds to the same experiment, but by removing the five first scales. We can see that the field is smoother and smoother, but respecting the large scale structure of the field. This transform will be very well suited to CMB studies where the phase is analyzed independently of the modulus, such as in Dineen et al. (2005) and Naselsky et al. (2005).

4 Polarized wavelet transform using spherical harmonics

4.1 Isotropic undecimated wavelet transform on the sphere (UWTS)

The undecimated isotropic transform on the sphere described in Starck et al. (2006) is similar in many respects to the usual à trous isotropic wavelet transform. It is obtained using a zonal scaling function which depends only on colatitude and is invariant with respect to a change in longitude . It follows that the spherical harmonic coefficients of vanish when which makes it simple to compute spherical harmonic coefficients of where * stands for convolution:

A possible scaling function (Starck et al. 1998), defined in the spherical harmonics representation, is where B₃ is the cubic B-spline compactly supported over [-2, 2]. Denoting a rescaled version of with cut-off frequency , a multi-resolution decomposition of f on a dyadic scale is obtained recursively:

where the zonal low pass filters h_j are defined by

The cut-off frequency is reduced by a factor of 2 at each step so that in applications where this is useful such as compression, the number of samples could be reduced adequately. Using a pixelization scheme such as Healpix (Górski et al. 2005), this can easily be done by dividing by 2 the Healpix nside parameter when computing the inverse spherical harmonics transform. As in the à trous algorithm, the wavelet coefficients can be defined as the difference between two consecutive resolutions, . This defines a zonal wavelet function as

With this particular choice of wavelet function, the decomposition is readily inverted by summing the coefficient maps on all wavelet scales

where we have made the simplifying assumption that f is equal to c₀. Obviously, other wavelet functions could be used just as well, such as the needlet function (Marinucci et al. 2008).

4.2 Extension to polarized data

By applying the above scalar isotropic wavelet transform to each component T, Q, U of a polarized map on the sphere, we have:

where c_J^X stands for the low resolution approximation to component X and w_j^X is the map of wavelet coefficients of that component on scale j. This leads to the following decomposition:

Figure 5 shows the backprojection of a Q-wavelet coefficient (left) and a U-wavelet coefficient (right). Figure 7 shows the undecimated isotropic polarized wavelet transform of the dust image shown in Fig. 6 using six scales, i.e. five wavelet scales and the coarse approximation.

5 Polarized curvelet transform

The 2D ridgelet transform (Candès & Donoho 1999a) was developed in an attempt to overcome some limitations inherent in former multiscale methods e.g. the 2D wavelet, when handling smooth images with edges i.e. singularities along smooth curves. Ridgelets are translation invariant ridge functions with a wavelet profile in the normal direction. Although ridgelets provide sparse representations of smooth images with straight edges, they fail to efficiently handle edges along curved lines. This is the framework for curvelets which were given a first mathematical description in Candès & Donoho (1999b). Basically, the curvelet dictionary is a multiscale pyramid of localized directional functions with anisotropic support obeying a specific parabolic scaling such that at scale 2^-j, its length is 2^-j/2 and its width is 2^-j. This is motivated by the parabolic scaling property of smooth curves. Other properties of the curvelet transform as well as decisive optimality results in approximation theory are reported in Candès & Donoho (1999c,b). Notably, curvelets provide optimally sparse representations of manifolds which are smooth away from edge singularities along smooth curves. Several digital curvelet transforms (Starck et al. 2002; Candès et al. 2006; Donoho & Duncan 2000) have been proposed which attempt to preserve the essential properties of the continuous curvelet transform and many papers (Herrmann et al. 2008; Starck et al. 2004b) report on their successful application in image processing experiments. The so-called first generation discrete curvelet described in Donoho & Duncan (2000) and Starck et al. (2002) consists in applying the ridgelet transform to sub-images of a wavelet decomposition of the original image. By construction, the sub-images are well localized in spaceand frequency and the subsequent ridgelet transform provides the necessary directional sensitivity. This latter implementation in combination with the good geometric properties of the Healpix pixelization scheme, inspired the digital curvelet transform on the sphere (Starck et al. 2006). The digital curvelet transform on the sphere is clearly invertible in the sense that each step of the overall transform is itself invertible. The curvelet transform on the sphere has a redundancy factor of 16J +1 when J scales are used, which may be a problem for handling huge data sets such as from the future Planck-Surveyor experiment. This can be reduced by substituting the pyramidal wavelet transform to the undecimated wavelet transform in the above algorithm. More details on the wavelet, ridgelet, curvelet algorithms on the sphere can be found in Starck et al. (2006) and software related to these new transforms is available from the web page http://jstarck.free.fr/mrs.html. As for the isotropic wavelet on the sphere, a straightforward extension to polarized data will consist in applying successively the curvelet transform on the sphere to the three components T, Q and U. Figure 8 shows the backprojection of a Q-curvelet coefficient and U-curvelet coefficient. Clearly, the shapes of these polarized curvelet functions are very different from the polarized wavelet functions. In the next section, we will build very different dictionaries using the E/B mode decomposition.

6 Polarized E/B wavelet and E/B curvelet

6.1 Introduction

We have seen that the generalization of the Fourier representation for polarized data on the sphere is the spin-2 spherical harmonics basis denoted :

At this point, it is convenient (Zaldarriaga 1998) to introduce the two quantities denoted E and B which are defined on the sphere by

where stands for the usual spin 0 spherical harmonics basis functions. The quantities E and B are derived by applying the spin lowering operator twice to and the spin raising operator twice to so that E and B are real scalar fieldson the sphere, invariant through rotations of the local reference frame. The normalization of and chosen in the latter definition is purely conventional but it appears to be rather popular Bunn et al. (2003); Zaldarriaga & Seljak (1997). Still, we could multiply and by some and we would have just as good a representation of the initial polarization maps. Through a change of parity E will remain invariant whereas the sign of pseudo-scalar B will change. The E and B modes defined here are not so different from the gradient (i.e. curl free) and curl (i.e. divergence free) components encountered in the analysis of vector fields. Finally, the spatial anisotropies of the Gaussian CMB temperature and polarization fields are completely characterized in this new linear representation by the power spectra and cross spectra of T, E and B. Thanks to the different parities of T and E on one side and B on the other, the sufficient statistics reduce to only four spectra namely . For a given cosmological model, it is possible to give a theoretical prediction of these spectra. Aiming at inverting the model and inferring the cosmological parameters, an important goal of CMB temperature and polarization data analysis is then to estimate the latter power spectra, based on sampled, noisy sometimes incomplete T,Q,U spherical maps.

6.2 E/B isotropic wavelet

Following the above idea of representing CMB polarization maps by means of E and B modes, we propose a formal extension of the previous undecimated isotropic wavelet transform that will allow us to handle linear polarization data maps T,Q,U on the sphere. Practically, the maps we consider are pixelized using for instance the Healpix pixelization scheme. In fact, we are not concerned at this point with the recovery of E and B modes from pixelized or incomplete data maps which itself is not a trivial task. The extension of the wavelet transform on the sphere we describe here makes use of the E and B representation of polarized maps described above in a formal way. Given polarization data maps T, Q and U, the proposed wavelet transform algorithm consists of the following steps:

1.

apply the spin 2 spherical harmonics transform to and . Practically, the Healpix software package provides an implementation of this transform for maps that use this pixelization scheme. Otherwise, a fast implementation was recently proposed by Wiaux et al. (2007);

2.

combine the decomposition coefficients and from the first step into and and build formal E and B maps associated with Q and U by applying the usual inverse spherical harmonics transform, as in Eq. (29). For numerical and algorithmic purposes, it may be efficient to stay with the spherical harmonics representation of E and B;

3.

apply the undecimated isotropic transform on the sphere described above to map T and to the E, B representation of the polarization maps.

The wavelet coefficient maps w_j^T, w_j^E, w_j^B and the low resolution approximation maps c_J^T, c_J^E, c_J^B are obtained by applying the isotropic undecimated wavelet transform described in Sect. 4.1 to the T, E, B representation of the polarized data. Figure 9 shows the result of applying the proposed transform to the polarized CMB data map ka from the WMAP experiment. The top two images show the initial Q and U maps while the subsequent maps are the low pass and wavelet coefficients' maps in a four scale decomposition. The scaling function we used is a cubic box spline as proposed in Sect. 4.1. The wavelet coefficients were obtained as the difference between two successive low pass approximations of the multiresolution decomposition of the E and B maps. The proper choice for the scaling and wavelet functions will depend on the application and the existence of constraints to be enforced.

Reconstruction

Obviously, the transform described above is invertible and the inverse transform is readily obtained by applying the inverse of each of the three steps in reverse order. If, as in the example decomposition above, we take the wavelet function to be the difference between two successive low pass approximations, the third step is inverted by simply summing the last low pass approximation with the maps of wavelet coefficients from all scales as in Eq. (25):

where c_J^X stands for the low resolution approximation to component X and w_j^X is the map of wavelet coefficients of that component on scale j. Finally, noting that:

the initial representation of the polarized data in terms of T, Q and U maps is reconstructed from its wavelet coefficients using the following equations:

where

with denoting the transpose conjugate of W so that is the scalar dot product of and W while is an operator (or matrix) acting on its left hand side as a projection along W and reconstruction along . In practice, the Healpix software package provides us with an implementation of the forward and inverse spin 0 and spin 2 spherical harmonics transforms which we need to implement the proposed inverse transform given by Eqs. (33) and (34). Clearly, as mentioned earlier in Sect. 4.1, we could have chosen some other wavelet function than merely the difference between two consecutive scaling functions, and the transformation would still be nearly as simple to invert. Figure 10 shows, on the left, backprojections of E-wavelet coefficients, and, on the right, backprojections of B-wavelet coefficients on the right hand side at different scales.

E - B curvelet

Similarly to the EB-wavelet constructions, we can easily construct an EB-curvelet transform by first computing the E and B components using the spin 2 spherical harmonics transform, and then applying a curvelet transform on the sphere separately on each of these two components. Figure 11 shows the backprojection of an E-curvelet coefficient and Fig. 12 shows the backprojectionof a B-curvelet coefficient.

7 Experiments: application to denoising

Thanks to the invertibility of the different proposed transforms for polarization maps on the sphere, it is possible to use these transformations for restoration applications. The denoising algorithm we use here consists in three consecutive steps:

1.

apply the chosen polarized multiscale transform;

2.

set to zero those coefficients whose absolute values are below a given threshold. We have used a threshold equal to five times the noise standard deviation;

3.

reconstruct the denoised field using the inverse transform.

More sophisticated thresholding strategies exist (Starck & Murtagh 2006) which can be used just as well on coefficients of polarized wavelets and curvelets.

Figure 13 illustrates the results of a simple denoising experiment. Noise was added to a simulated dust image (see Fig. 13 top left and right), and the noisy QU-field was filtered by thresholding either the EB-isotropic wavelet coefficients of the polarized dust map (Fig. 13 bottom left) or the EB-isotropic curvelet coefficients (Fig. 13 bottom right). Both decompositions produce nice visual results.

In order to be more quantitative, we used two test images (synchrotron and dust) with different noise levels. The noise levels were taken equal to 0.1,0.2,0.5,1. and 2 times the standard deviation of the noise-free image. On each noisy image, we applied six different transformations, thresholded the coefficients and reconstructed the filtered images. The transforms we used are the QU-wavelets (QU-UWT), the EB-wavelets (EB-UWT), the QU-curvelet (QU-CUR), the EB-curvelet (EB-CUR), the biorthogonal wavelet transform (OWT) and the modulus-phase undecimated multiscale transform (MP-UWT). For each filtered image Q (resp. U), we computed the Signal-to-Noise Ratio (SNR) in dB:

where is the standard deviation of the noise free original image component Q (resp. U) and is the standard deviation of the error image i.e. the difference between the noise-free component Q (resp. U) and the filtered component Q (resp. U). These errors are given in Table 1 for the dust image and in Table 2 for the synchrotron image.

It is clear from the above results that the different transforms described here do not perform as well on this specific numerical experiment. For the dust, QU-wavelets are better when the noise level is not so high, while curvelets do better when the noise and signal levels are of the same order. This can be explained by the fact that structures on small scales are more or less isotropic, and therefore better represented by wavelets, while large scale structures are anisotropic and therefore better analyzed using curvelets. When increasing the noise level, structures on the smallest scales are no longer detected by either of the two dictionaries. Only large scale features are detectable, and curvelets do this job better. Dealing with the polarized synchrotron map, curvelets do better than wavelets in all cases experimented with here. Although the bi-orthogonal wavelet transform is clearly not as good as the others in these experiments, it could nevertheless be very useful in situations where computation time is an important issue. Indeed, since it doesn't make use of the spherical harmonics transform and also because it is not redundant, it is a very fast transform. Finally the worse results were obtained using the Modulus-Phase multiscale transform. This can be explained by the fact that the Gaussian approximation we made for the noise in the wavelet transform of the modulus is not accurate enough. Furthermore, it is not clear what is the best way to correct the phase wavelet coefficients from the noise. A better understanding of the noise behavior after transformation is clearly required before the Mod-phase multiscale transform can be used for restoration purposes. However, for other applications such as non-Gaussianity studies, the latter transform could prove an interesting tool to use as well.

Table 1:   Error (in dB) for both the Q and U components between the original dust image and respectively its noisy version and the results obtained by hard thresholding representations of the noisy data in different dictionaries.

Table 2:   Error (in dB) for both the Q and U components between the original synchrotron image and respectively its noisy version and the results obtained by the hard thresholding using different dictionaries.

8 Conclusion

The present contribution has enlarged the set of available representations of polarized data on the sphere. We described the construction of new multiscale decompositions, namely the modulus-phase transform, the QU and EB isotropic undecimated wavelet and curvelet transforms for polarized data on the sphere. The latter are derived as extensions of the undecimated wavelet and curvelet transforms for scalar maps on the sphere. The proposed extensions use a formal representation of T, Q and U polarization data maps in terms of the scalar T, E and pseudo-scalar B maps. The proposed transforms are invertible allowing for applications in image restoration and denoising as was shown in a preliminary experiment. Ongoing research is concerned with assessing and understanding the merits of the different transforms we described, in such problems as restoring, denoising or inpainting of sparse linearly polarized data, as well as in the blind separation of mixed linearly polarized components.

Software

The software related to this paper, MRS/POL, and its full documentation will be included in the next version of the the MRS software which is available from the following web page: http://jstarck.free.fr/mrs.html

Acknowledgements

Some of the results in this paper have been derived using the Healpix (Górski et al. 2005).

References

• Abrial, P., Moudden, Y., Starck, J., et al. 2007, J. Fourier Analysis and Applications, 13, 729 [CrossRef]
• Abrial, P., Moudden, Y., Starck, J., et al. 2008, Statistical Methodology, 5, 289 [NASA ADS] [CrossRef]
• Aghanim, N., & Forni, O. 1999, A&A, 347, 409 [NASA ADS]
• Antoine, J.-P. 1999, in Wavelets in Physics, 23
• Bobin, J., Moudden, Y., Starck, J. L., Fadili, J., & Aghanim, N. 2008, Statistical Methodology, 5, 307 [NASA ADS] [CrossRef]
• Bunn, E. F., Zaldarriaga, M., Tegmark, M., & de Oliveira-Costa, A. 2003, Phys. Rev. D, 67, 023501 [NASA ADS] [CrossRef]
• Candès, E., & Donoho, D. 1999a, Philosophical Transactions of the Royal Society of London A, 357, 2495 [NASA ADS] [CrossRef] (In the text)
• Candès, E. J., & Donoho, D. L. 1999b, in Curve and Surface Fitting: Saint-Malo 1999, ed. A. Cohen, C. Rabut, & L. Schumaker (Nashville, TN: Vanderbilt University Press) (In the text)
• Candès, E. J., & Donoho, D. L. 1999c, in Curve and Surface Fitting: Saint-Malo 1999, ed. A. Cohen, C. Rabut, & L. Schumaker (Nashville, TN: Vanderbilt University Press)
• Candès, E., Demanet, L., Donoho, D., & Ying, L. 2006, SIAM Multiscale Model. Simul., 5/3, 861
• Cayón, L., Sanz, J. L., Martínez-González, E., et al. 2001, MNRAS, 326, 1243 [NASA ADS] [CrossRef]
• Daubechies, I., & Sweldens, W. 1998, J. Fourier Analysis and Applications, 4, 247 [CrossRef] (In the text)
• Delabrouille, J., Cardoso, J., Le Jeune, M., et al. 2008, ArXiv e-prints
• Dineen, P., Rocha, G., & Coles, P. 2005, MNRAS, 358, 1285 [NASA ADS] [CrossRef] (In the text)
• Donoho, D., & Duncan, M. 2000, in Proc. Aerosense 2000, Wavelet Applications VII, ed. H. Szu, M. Vetterli, W. Campbell, & J. Buss, SPIE, 4056, 12
• Faÿ, G., & Guilloux, F. 2008, ArXiv e-prints (In the text)
• Faÿ, G., Guilloux, F., Betoule, M., et al. 2008, Phys. Rev. D, 78, 083013 [NASA ADS] [CrossRef] (In the text)
• Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] (In the text)
• Hammond, D. K., Wiaux, Y., & Vandergheynst, P. 2008, ArXiv e-prints
• Herrmann, F., Moghaddam, P., & Stolk, C. 2008, ACHA, 24, 150
• Hobson, M. P., Jones, A. W., & Lasenby, A. N. 1999, MNRAS, 309, 125 [NASA ADS] [CrossRef]
• Holschneider, M. 1996, J. Math. Phys., 37, 4156 [NASA ADS] [CrossRef]
• Jin, J., Starck, J.-L., Donoho, D., Aghanim, N., & Forni, O. 2005, Eurasip J. Appl. Signal Proc., 15, 2470 [CrossRef]
• Kovac, J. M., Leitch, E. M., Pryke, C., et al. 2002, Nature, 420, 772 [NASA ADS] [CrossRef] (In the text)
• Marinucci, D., Pietrobon, D., Balbi, A., et al. 2008, MNRAS, 383, 539 [NASA ADS] (In the text)
• McEwen, J. D., Vielva, P., Wiaux, Y., et al. 2007, J. Fourier Analysis and Applications, 13, 495 [NASA ADS] [CrossRef]
• Moudden, Y., Cardoso, J.-F., Starck, J.-L., & Delabrouille, J. 2005, Eurasip J. Appl. Signal Proc., 15, 2437 [CrossRef]
• Naselsky, P., Chiang, L.-Y., Olesen, P., & Novikov, I. 2005, Phys. Rev. D, 72, 063512 [NASA ADS] [CrossRef] (In the text)
• Rahman, I. U., Drori, I., Stodden, V. C., Donoho, D. L., & Schröder, P. 2005, Multiscale Modeling & Simulation, 4, 1201 [CrossRef] (In the text)
• Smoot, G. F., Bennett, C. L., Kogut, A., et al. 1992, ApJ, 396, L1 [NASA ADS] [CrossRef] (In the text)
• Starck, J.-L., & Murtagh, F. 2006, Astronomical Image and Data Analysis, 2nd Ed. (Berlin: Springer), A&AL
• Starck, J.-L., Murtagh, F., & Bijaoui, A. 1998, Image Processing and Data Analysis: The Multiscale Approach (Cambridge University Press)
• Starck, J.-L., Candès, E., & Donoho, D. 2002, IEEE Transactions on Image Processing, 11, 131
• Starck, J.-L., Aghanim, N., & Forni, O. 2004a, A&A, 416, 9 [NASA ADS] [CrossRef] [EDP Sciences]
• Starck, J.-L., Elad, M., & Donoho, D. 2004b, Advances in Imaging and Electron Physics, 132
• Starck, J.-L., Moudden, Y., Abrial, P., & Nguyen, M. 2006, A&A, 446, 1191 [NASA ADS] [CrossRef] [EDP Sciences] (In the text)
• Tenorio, L., Jaffe, A. H., Hanany, S., & Lineweaver, C. H. 1999, MNRAS, 310, 823 [NASA ADS] [CrossRef]
• Vielva, P., Martínez-González, E., Barreiro, R. B., Sanz, J. L., & Cayón, L. 2004, ApJ, 609, 22 [NASA ADS] [CrossRef]
• Vielva, P., Wiaux, Y., Martínez-González, E., & Vandergheynst, P. 2006, New Astron. Rev., 50, 880 [NASA ADS] [CrossRef]
• Wiaux, Y., Jacques, L., & Vandergheynst, P. 2007, J. Comput. Phys., 226, 2359 [NASA ADS] [CrossRef] (In the text)
• Wiaux, Y., McEwen, J. D., Vandergheynst, P., & Blanc, O. 2008, MNRAS, 388, 770 [NASA ADS] [CrossRef] (In the text)
• Zaldarriaga, M. 1998, ApJ, 503, 1 [NASA ADS] [CrossRef] (In the text)
• Zaldarriaga, M., & Seljak, U. 1997, Phys. Rev. D, 55, 1830 [NASA ADS] [CrossRef]

Footnotes

... Healpix

http://healpix.jpl.nasa.gov

...p

In differential geometry, the Exp map and Log map are generalizations of the usual exponential and logarithm function. In this paper, the manifold is a Riemannian manifold. In that case, the Exp map at point p₀, Exp_p0(s) is the map which takes a vector s of the tangent space of at p₀ and provides the point p₁ by travelling along the geodesic starting at p₀ in the direction s.

...ka

Available at http://lambda.gsfc.nasa.gov/product/map/current/.

All Tables

Table 1:   Error (in dB) for both the Q and U components between the original dust image and respectively its noisy version and the results obtained by hard thresholding representations of the noisy data in different dictionaries.

Table 2:   Error (in dB) for both the Q and U components between the original synchrotron image and respectively its noisy version and the results obtained by the hard thresholding using different dictionaries.

All Figures

---

Copyright ESO 2009