Polarized wavelets and curvelets on the sphere
J.L. Starck  Y. Moudden  J. Bobin
Laboratoire AIM (UMR 7158), CEA/DSMCNRSUniversité Paris Diderot, IRFU, SEDISAP, Service d'Astrophysique, Centre de Saclay, 91191 GifSurYvette 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 fullsky 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 QU or EB wavelet transforms and QU or EB 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 baryonphoton 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 socalled acoustic oscillations in the baryonphoton 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 fullsky 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 nonGaussianity 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 socalled ``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 nonGaussianity 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 QU or EB wavelet transforms, and QU or EB 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
Figure 1: QU orthogonal wavelet transform. 

Open with DEXTER 
Fullsky 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:
(1)  
which can also be written which by definition expresses the fact that the quantities are spin2 fields on the sphere. The suitable generalization of the Fourier representation for such fields is the spin2 spherical harmonics basis denoted , in which we can expand:
Figure 2: Top: examples of backprojections of Qwavelet coefficients. Bottom: examples of backprojections of Uwavelet coefficients. 

Open with DEXTER 
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 biorthogonal wavelet transform on each face of the Healpix map, separately for Q and U.
Figure 1 shows the flowgraph of this QU orthogonal wavelet transform (QUOWT). 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^{2} 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:
(3)  
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 nonzero coefficient allows us to view the different atoms in the dictionary related to the QUOWT 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 Modulephase non linear multiscale transform
3.1 Introduction
Given a polarized map in the standard QU representation, consider a different point of view and define the modulus M and phase P maps as follows:(5)  
(6) 
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 nonnegative and ii) the phase field takes its values on the unit circle S^{1}. 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 MPmultiscale 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 logmap of p_{1} onto the tangent space of at p_{0}. The backprojection is obtained using the inverse of the logmap Exp_{p0}^{}.
For instance, if we choose
then
and
.
The
and
maps are then defined as follows:
(7)  
(8) 
The multiscale transform for manifold valued data introduced in Rahman et al. (2005) is equivalent to a twostep interpolationrefinement 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 nonlinearity 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:w_{j+1}^{P}  =  (11)  
c_{j+1,k}^{P}  =  (12) 
This multiscale transform is invertible and its inverse is computed as follows:
(13)  
(14) 
The picture in Fig. 3 features some examples of backprojections of MPmultiscale coefficients.
Figure 3: Examples of MPmultiscale coefficients backprojection. 

Open with DEXTER 
3.3 Undecimated MPmultiscale 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:(17)  
(18) 
where:
(19) 
This multiscale transform is invertible and its inverse is computed as follows:
(20) 
Figure 4 top shows a simulated polarized field of the synchrotron emission and its noisy version. We have applied the MPmultiscale 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).
Figure 4: Polarized field smoothing  top left: simulated synchroton emission. Top right: same field corrupted by additive noise. Bottom left: MPmultiscale reconstruction after setting to zero all coefficients from the three first scales. Bottom right: MPmultiscale reconstruction after setting to zero all coefficients from the five first scales. 

Open with DEXTER 
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:
(21) 
A possible scaling function (Starck et al. 1998), defined in the spherical harmonics representation, is where B_{3} is the cubic Bspline compactly supported over [2, 2]. Denoting a rescaled version of with cutoff frequency , a multiresolution decomposition of f on a dyadic scale is obtained recursively:
c_{0}  =  
c_{j}  =  (22) 
where the zonal low pass filters h_{j} are defined by
=  
=  (23) 
Figure 5: Qisotropic wavelet transform backprojection ( left) and Uisotropic wavelet backprojection ( right). 

Open with DEXTER 
The cutoff 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_{0}. 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 Qwavelet coefficient (left) and a Uwavelet 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
Figure 6: Simulated observations on the sphere of the polarized galactic dust emission. 

Open with DEXTER 
Figure 7: QUundecimated wavelet transform of the simulated polarized map of galactic dust emission shown in Fig. 6. 

Open with DEXTER 
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 socalled first generation discrete curvelet described in Donoho & Duncan (2000) and Starck et al. (2002) consists in applying the ridgelet transform to subimages of a wavelet decomposition of the original image. By construction, the subimages 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 PlanckSurveyor 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 Qcurvelet coefficient and Ucurvelet 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 spin2 spherical harmonics basis denoted
:
(28) 
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 pseudoscalar 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.
Figure 8: Top, Qcurvelet backprojection ( left) and zoom ( right). Bottom, Ucurvelet backprojection ( left) and zoom. 

Open with DEXTER 
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.
Figure 9: Top: Q and U CMB polarization data maps from channel ka of the WMAP experiment. Left: low pass and wavelet coefficients in three scales of the formal E mode. Right: low pass and wavelet coefficients in three scales of the formal B mode. 

Open with DEXTER 
Reconstruction
Figure 10: Eisotropic wavelet transform backprojection ( left) and Bisotropic wavelet backprojection ( right). 

Open with DEXTER 
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):
(30) 
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:
Q  =  
=  (31)  
U  =  
=  (32) 
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 Ewavelet coefficients, and, on the right, backprojections of Bwavelet coefficients on the right hand side at different scales.
E  B curvelet
Figure 11: Ecurvelet coefficient backprojection. 

Open with DEXTER 
Figure 12: Bcurvelet coefficient backprojection. 

Open with DEXTER 
Similarly to the EBwavelet constructions, we can easily construct an EBcurvelet 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 Ecurvelet coefficient and Fig. 12 shows the backprojectionof a Bcurvelet coefficient.
Figure 13: Top, simulated input polarized image ( left) and noisy polarized field ( right). Bottom, denoising of the polarized field using the EBisotropic undecimated wavelets ( left) and the EBcurvelets ( right). 

Open with DEXTER 
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.
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 QUfield was filtered by thresholding either the EBisotropic wavelet coefficients of the polarized dust map (Fig. 13 bottom left) or the EBisotropic 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 noisefree image.
On each noisy image, we applied six different transformations, thresholded the coefficients and reconstructed the filtered images. The transforms we used are the QUwavelets (QUUWT),
the EBwavelets (EBUWT), the QUcurvelet (QUCUR), the EBcurvelet (EBCUR), the biorthogonal wavelet transform (OWT) and the modulusphase undecimated multiscale transform (MPUWT).
For each filtered image Q (resp. U), we computed the SignaltoNoise Ratio (SNR) in dB:
(35) 
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 noisefree 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, QUwavelets 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 biorthogonal 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 ModulusPhase 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 Modphase multiscale transform can be used for restoration purposes. However, for other applications such as nonGaussianity 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 modulusphase 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 pseudoscalar 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 OliveiraCosta, 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: SaintMalo 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: SaintMalo 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ínezGonzá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 eprints
 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 eprints (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 eprints
 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ínezGonzález, E., Barreiro, R. B., Sanz, J. L., & Cayón, L. 2004, ApJ, 609, 22 [NASA ADS] [CrossRef]
 Vielva, P., Wiaux, Y., MartínezGonzá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_{0}, Exp_{p0}(s) is the map which takes a vector s of the tangent space of at p_{0} and provides the point p_{1} by travelling along the geodesic starting at p_{0} 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
Figure 1: QU orthogonal wavelet transform. 

Open with DEXTER  
In the text 
Figure 2: Top: examples of backprojections of Qwavelet coefficients. Bottom: examples of backprojections of Uwavelet coefficients. 

Open with DEXTER  
In the text 
Figure 3: Examples of MPmultiscale coefficients backprojection. 

Open with DEXTER  
In the text 
Figure 4: Polarized field smoothing  top left: simulated synchroton emission. Top right: same field corrupted by additive noise. Bottom left: MPmultiscale reconstruction after setting to zero all coefficients from the three first scales. Bottom right: MPmultiscale reconstruction after setting to zero all coefficients from the five first scales. 

Open with DEXTER  
In the text 
Figure 5: Qisotropic wavelet transform backprojection ( left) and Uisotropic wavelet backprojection ( right). 

Open with DEXTER  
In the text 
Figure 6: Simulated observations on the sphere of the polarized galactic dust emission. 

Open with DEXTER  
In the text 
Figure 7: QUundecimated wavelet transform of the simulated polarized map of galactic dust emission shown in Fig. 6. 

Open with DEXTER  
In the text 
Figure 8: Top, Qcurvelet backprojection ( left) and zoom ( right). Bottom, Ucurvelet backprojection ( left) and zoom. 

Open with DEXTER  
In the text 
Figure 9: Top: Q and U CMB polarization data maps from channel ka of the WMAP experiment. Left: low pass and wavelet coefficients in three scales of the formal E mode. Right: low pass and wavelet coefficients in three scales of the formal B mode. 

Open with DEXTER  
In the text 
Figure 10: Eisotropic wavelet transform backprojection ( left) and Bisotropic wavelet backprojection ( right). 

Open with DEXTER  
In the text 
Figure 11: Ecurvelet coefficient backprojection. 

Open with DEXTER  
In the text 
Figure 12: Bcurvelet coefficient backprojection. 

Open with DEXTER  
In the text 
Figure 13: Top, simulated input polarized image ( left) and noisy polarized field ( right). Bottom, denoising of the polarized field using the EBisotropic undecimated wavelets ( left) and the EBcurvelets ( right). 

Open with DEXTER  
In the text 
Copyright ESO 2009