A&A 455, 741755 (2006)
DOI: 10.1051/00046361:20053820
SunyaevZel'dovich cluster reconstruction in multiband bolometer camera surveys
S. Pires^{1}  J. B. Juin^{1}  D. Yvon^{1}  Y. Moudden^{1}  S. Anthoine^{2}  E. Pierpaoli^{3}
1 
CEA  CE Saclay, DSM/DAPNIA, 91191 GifsurYvette Cedex, France
2 
PACM, Princeton University, Princeton, NJ 08544, USA
3 
California Institute of Technology, Mail Code 13033, Pasadena, CA 91125, USA
Received 12 July 2005 / Accepted 30 March 2006
Abstract
We present a new method for the reconstruction of SunyaevZel'dovich (SZ)
galaxy clusters in future SZsurvey experiments using multiband bolometer cameras such as Olimpo, APEX, or Planck. Our goal is to optimise SZCluster extraction from our observed noisy maps. None of the algorithms used in the detection chain is tuned using prior knowledge of the SZCluster signal, or other astrophysical sources (Optical Spectrum, Noise Covariance Matrix, or covariance of SZ Cluster wavelet coefficients). First, a blind separation of
the different astrophysical components that contribute to the observations is conducted
using an Independent Component Analysis (ICA) method.
This is a new application of ICA to multichannel astrophysical data analysis.
Then, a recent non linear filtering technique in the wavelet domain, based on multiscale entropy and the False Discovery Rate (FDR) method, is used
to detect and reconstruct the galaxy clusters. We use the Source Extractor software to identify the detected clusters. The proposed method was applied on realistic simulations of observations that we produced as mixtures of synthetic maps of the four brightest light sources in the range 143 GHz to 600 GHz namely the SunyaevZel'dovich effect, the Cosmic Microwave Background (CMB) anisotropies, the extragalactic InfraRed point sources and the Galactic Dust Emission. We also implemented a simple model of optics and noise to account for instrumental effects.
Assuming nominal performance for the near future SZsurvey Olimpo, our detection chain recovers 25% of the cluster of mass larger than
,
with 90% purity.
Our results are compared with those obtained with published algorithms.
This new method has a high global detection efficiency in the highpurity/low completeness region, being however a blind algorithm (i.e. without using any prior assumptions on the data to be extracted).
Key words: cosmology: largescale structure of Universe  cosmology: observations  galaxies: clusters: general  methods: data analysis  methods: numerical
Spectacular advances in cosmology have taken place through observation of the
sky at millimeter and
submillimeter wavelengths. The Boomerang (Netterfield et al. 2002), MAXIMA (Hanany et al. 2000), Archeops (Tristram et al. 2005), WMAP (Hinshaw et al. 2003), CBI (Contaldi et al. 2002) and DASI (Halverson et al. 2002) experiments measured the anisotropies of the Cosmic Microwave Background with high precision. Combined with other observations, such as distant supernovae of Riess et al. (1999); Knop et al. (2003) and/or Clusters (White et al. 1999; Pierpaoli et al. 2001) or the LargeScale Structure Observation (Eisenstein et al. 2005), these experiments allowed one to place tight constraints on parameters of generic cosmological models.
These results gave rise to new questions. Consensus cosmological models assume that the
largescale evolution of the Universe is driven by Dark Matter and Dark Energy. Dark Matter has been sought for 30 yr,
and most candidates have been rejected in the light of experimental results. Dark Energy will be one of the toughest challenges of modern cosmology. One way to learn about the nature of dark energy is to study the evolution of the Universe at late times, during which dark energy is thought to have impacted on the formation of largescale structures. Cluster surveys are an ideal tool for this purpose.
Largescale structures can be studied using several probes: optical galaxy catalogs, Xray cluster surveys such as XMMLSS (Pierre et al. 2004),
SZ cluster surveys such as Olimpo (Masi et al. 2003), APEX, AMI (Jones 2001), AMIBA (AMIBA 20052007) and weak shear such as CFHTLS "wide''.
The mass and redshift distribution of clusters of galaxies,
(where M, z and
are the cluster mass and redshift, and the solid angle), is predicted with fair precision using Halo models and Nbody simulations. SZ cluster, Xray and
optical surveys measure quantities related to cluster mass through complex phenomena, such as interstellar gas heating processes, clustering and processes of stellar formation and evolution that are still not fully understood. The SZ signal from galaxy clusters is believed to be the simplest direct observable of cluster mass, though some uncertainties on extragalactic gas heating mechanisms remain. Measuring the cluster mass function will require these systematics to be understood. Pioneering SZ cluster surveys will obtain their first data this year. The Planck satellite survey will follow (Lamarre et al. 2004). In this paper we prepare for their analysis using simulated data.
Melin (2004) and Schafer (2004), following Herranz et al. (2002a,b), developed a cluster extraction tool based on an optimised filter. The filter is optimised to damp the noise as well as the other astrophysical components. The SZ signal is small compared to the other astrophysical sources, rendering the task of filtering difficult. In addition a filter damps some of the signal at the same time as it rejects parasites. We developed the idea of a 3 step method using as a first step a source separation algorithm. The source separation algorithm sorts the noisy SZ signal from the other astrophysical components and the filtering step is left with the easier task of improving the SZ cluster signal versus noise. A detection algorithm extracts the SZ cluster candidates.
In the 100 to 600 GHz range, the brightest components of the sky are the Cosmic Microwave Background (CMB), the Infrared Point Sources, the Galactic dust emission and, swamped by the previous ones, the SZ clusters. It follows that the true sky map
,
in a given optic band centered on ,
can be modeled as a sum of distinct astrophysical radiations as in



(1) 
where
denote spatial or angular indexes on 2D or spherical maps.
The observed sky map
is the result of the convolution of the true sky with the optical beam of the telescope plus the unavoidable contribution of instrumental noise
.
Multiband signal processing in astrophysical applications deals with exploiting correlations between data maps at different ,
to enhance the estimation of specific objects of interest such as the power spectrum of the CMB spatial fluctuations, or the mass function of galaxy clusters. Considering the case where the radiative properties of the sources are completely isotropic, in the sense that they do not depend on the direction of observation, the above model can be rewritten in the following factored form:

(2) 
where S_{i} is the spatial template (i.e. the convolution of the true component map with the instrumental beam at frequency )
and
the emission law of the
ith astrophysical component. Although this is a coarse approximation in the case of Infrared Point Sources as described in Sect. 2.2, it is mostly valid for the other three components. With observations available in m
channels, Eq. (2) can be written in matrix
form:

(3) 
where
is a vector in
,
A is an
matrix, n is the number of contributing astrophysical components,
is now a vector in
and
in
.
An important simplifying assumption made here is that the instrumental beam does not vary with the optical frequency .
This is acceptable with Olimpo but definitely not for other multichannel experiments such as Planck. Accounting for a varying instrumental beam requires considering a more complex convolutive mixture model for the observed multichannel data either with known or parameterized beams. Some methods have been developed to deal with such mixtures but most of these move to the Fourier domain and consider an instantaneous mixture model in this representation. This is not acceptable here as the component of interest (i.e. SZ clusters) would be below the noise level in all channels in this representation: this is due to the nonlocality of the Fourier basis functions.
Instead, we chose to move the data to the wavelet domain. Wavelets are known to provide sparse representations of images with point singularities in an otherwise smooth background. Thus enhancement of the SZ signal concentration is expected. Then, once computed, in order to reconstruct the astrophysical component maps, the estimated separating filter is applied onto the original data.
Equation (3) shows that the observations consist of linear mixtures of astrophysical components with different weights and additive noise. Assuming no prior knowledge
of A, the problem of recovering the particular component of interest i.e. the SZ effect map of galaxy clusters, can be seen as a blind source separation problem which we approach in this paper using Independent Component Analysis (ICA). Such methods have been previously used quite successfully in the analysis of astrophysical data from present or future multichannel experiments such as WMAP or Planck (Moudden et al. 2005; Kuruoglu et al. 2003; Bobin et al. 2005; Maino et al. 2002) where the focus is on estimating the map of CMB anisotropies. We concentrate here on separating the SZ component, which has very different statistical properties compared to CMB and thus requires the use of other ICA methods. Then, due to the additive instrumental noise, the source separation process has to be followed by a filtering technique.
We use an iterative filtering based on Bayesian methods to filter the noise (Starck et al. 2006).
This method uses a Multiscale Entropy prior which is only defined for nonsignificant wavelet coefficients selected by the False Discovery Rate (FDR) Method (Benjamini & Hochberg 1995). Then, we identify SZ clusters in filtered maps using the Source Extractor software SExtractor (Bertin 2003) which conveniently provides source extraction along with adequate photometry and deblending.
In the following, we first describe our sky model. We use the Olimpo SZcluster survey (Masi et al. 2003) as an example of a future ambitious SZ survey. We model the Olimpo instrument and produce simulated observed maps. With minor parameter changes, these algorithms will also allow us to simulate observed local maps typical of other bolometer camera
SZ surveys (i.e. Planck or APEX observations). Then we explain the data processing methods used to recover the map of SZcluster signal and the "observed'' cluster catalog. We quantify the efficiency of our detection chain and compare it to a recently published method.
This paper is at the interface of astrophysics and applied mathematics.

Figure 1:
Graphical summary of the simulation and detection chain presented in this paper. 
Open with DEXTER 
The input of our simulation is a list of cosmological model parameter values. We use a
cosmological model with
,
,
,
,
,
h=0.7,
,
n_{s}= 1,
and a seed for the random number generator. This simulation was fully written in C++, so that we can mass produce our simulations on a PC farm, when needed.
Big Bang cosmology models assume that the Universe has been evolving from a hot and dense plasma. While it expands, its contents cools. At some time, decoupling occured, when the temperature of the Universe dropped below
eV and the mean free path of photons suddenly became infinite so that the large majority of them have been propagating freely through the Universe ever since. These photons are nowadays observed to have a blackbody spectrum at temperature 2.728 K (Hinshaw et al. 2003): the Cosmic Microwave Background (CMB). Small anisotropies in the CMB are observed and interpreted as resulting from temperature and density fluctuations in the primordial plasma at the decoupling time.
Bolometer cameras measure incoming power variations, while scanning the sky. Thus the main blackbody spectrum (monopole) will not be observed. In the following simulations we also assume that the CMB dipole (Doppler effect due to the Earth's movement) has been subtracted from the data.
In order to simulate a map of CMB anisotropies, we first enter the cosmological parameter list in the cmbeasy code (Doran & Mueller 2004) and compute a Power Spectrum in Spherical Harmonics up to multipole 5000 of the CMB fluctuations. Then using a small field approximation (White et al. 1999), we generate a random spatially correlated Gaussian field of primordial anisotropies with the previous power spectrum. Figure 2 shows a typical simulated CMB anisotropy map in units of .

Figure 2:
The 4 physical components of the sky included in our simulation:
upper right is a map of the CMB's anisotropies in unit of K,
upper left is the SZ Cluster map, in units of the Compton parameter on line of sight y,
lower right is the IR point source map, convolved with a beam of 2 arcmin, in Jy at 350 GHz, lower left is the Galactic dust map in units of MJy/st at 100 m. 
Open with DEXTER 
2.2 IR point sources
In our frequency range, infrared point sources are expected to be a significant contamination. Borys et al. (2003) from the
SCUBA experiment published a list of infrared point sources at 350 GHz and a
law from this data.
We extend this law to lower brightness, and generate a catalog of IR point sources in our field of view.
We assume that their optical spectra follow a grey body law parameterised as:

(4) 
Where Planck stands for Planck's black body spectrum,
stands for frequency, T_{0} is the temperature of the
black body and
is by definition the spectral index. For each IR point source, the spectral index was randomly picked between 1.5 and 2, and T_{0} was set to 30 K. Thus each IR point source has a different optical spectrum. Then the IR point source is placed on our map as a single pixel. IR sources are positioned randomly on the map but a user chosen fraction of them can be positioned within clusters of galaxies, and would make SZCluster photometry more difficult. We obtain an IR point source map in units of Jy at 350 GHz as shown in Fig. 2.

Figure 3:
Simulated maps in Olimpo's four frequency bands:
upper right is the 147 GHz band, upper left is the 217 GHz band, lower right is the 385 GHz band: CMB anisotropies, IR point sources and Galactic Dust blend in this band.
Lower left 500 GHz: IR point sources and Galactic Dust are the dominant features at high frequencies. The SZ cluster signal is dominated by other astrophysical sources at all frequencies. 
Open with DEXTER 
We know from FIRAS and IRAS that the galactic interstellar dust emits in the far infrared. For
these simulations, we wanted to generate truly random galactic dust maps. So, instead of using the maps published by Schlegel et al., we followed
Bouchet & Gispert (1999) and assumed that Galactic Dust was described by a power spectrum of spatial correlations in
.
We then randomly generated a map with the same algorithms as used for the CMB map simulations. We assume a homogeneous optical spectrum over our map given by Eq. (4) with ,
and T_{0}= 20 K. This is a reasonnable assumption since SZCluster observations will take place at high galactic latitudes, and we will be working on small size maps. Finally, we normalise the map rms at its observed values at high galactic latitude:
at
m, the rms flux is observed at 1 MJ/st. We obtain a galactic Dust Emission Map, shown in Fig. 2 in natural units
of MJ/st at 100 m.
2.4 SunyaevZel'dovich clusters
Hot intracluster gas is observed as a plasma at a temperature of 2 to 20 keV. A small fraction of CMB photons traveling through clusters Comptonscatter on electron gas. Since the electron kinetic energy is much higher than the CMB photon energy, CMB photons statistically gain energy when diffused. This effect, known as the thermal SunyaevZel'dovich effect (1980), results in a small distortion in the CMB BlackBody spectrum, in the direction of clusters of galaxies.
The probability that a photon scatters while crossing a cluster is given by:

(5) 
The optical spectral dependance is given by:

(6) 
In these expressions,
refer to the electron temperature, density and mass, respectively;
is the Thomson crosssection and
is the dimensionless frequency of observation in terms of unperturbed CMB temperature T= 2.728 K .
Simulation of SZ Cluster maps requires modelling largescale structure formation.
For this we translated in C++ the ICosmo semianalytic model (Refregier & R. Massey 2004) which, starting from a cosmological parameter set, allows us to compute the cluster density versus mass and redshift:
.
Several fits to data and physical assumptions are involved and implemented in this computation. We run CICosmo using the options of a concordance CDM model and the mass function from
Sheth & Thormen (Sheth et al. 2001). The virial radius of a cluster is computed assuming a spherical collapse (Peebles 1993; Peacok 1999). This is the radius where the average overdensity is
times the critical density (chosen independent of redshift).
Gas heating is modeled according to Pierpaoli et al. (2003), with the additional parameter T^{*}. Then we randomly generate a cluster catalog in our field of view, according to density profile keeping those clusters which have an integrated Compton parameter Y larger than
st (corresponding to a mass cut of 10^{14} Solar Mass).
An ellipticity is attributed to each cluster according to experimental observations (Cooray 2000). We assume that the cluster gas density follows an elliptical beta model. For the time being, clusters are positioned on the map randomly and the orientation of the main axis of each ellipse is also chosen randomly. We project to account for two point correlations in cluster positions in
the near future. A typical SZ cluster map in natural units of yCompton is shown in Fig. 2.
A detailed model of observations can only be developed after data acquisition. In the following, we choose to use a very simple observation model, nevertheless representative of an ongoing project, the Olimpo balloon experiment. The Olimpo camera will observe the sky at 4 frequencies, 143, 217, 385 and 600 GHz. The millimeterwave filter bandwiths will be close to 30 GHz and are assumed to have a tophat shape. The beams at the 4 frequencies are assumed to be symmetric, Gaussian and of 3, 2, 2 and 2 arcmin FHWM respectively. The optical efficiency (mirrors, filters), the photon noise and the bolometer sensitivities and noise contributions are summarized by an equivalent observed noise temperature
on the Sky in unit of
.
Based on the BOOMERanG experiment, typical values in the four channels of Olimpo are expected at 150, 200, 500, and 5000
respectively.
The simplest and most optimistic model is to assume that the observed noise along bolometer timelines will be white, in which
case, after projection (or coaddition, Yvon & Mayet 2005) on 2D maps, pixel noise can be very easly computed as:

(7) 
where
stands for the number of bolometers working in a particular frequency band,
and
is the observation length on this pixel.
High resolution SZ bolometer surveys will cover their large survey adding up small patches of the sky. Noise map patterns are likely to be complex, and will only be known after data has been acquired.
We chose in this first paper to begin with a simple case of homogeneous white noise on the full field of view.
Frequency [GHz] 
143 
217 
385 
600 
FWHM [arcmin] 
3 
2 
2 
2 
White noise level [
] 
150 
200 
500 
5000 
We now have all the tools to compute simulated observed sky maps. For each frequency band, we compute a conversion
factor from the maps in their natural units into the observed unit at the bolometer level namely
.
The exception
is the IR point source map where, because of the random spectral indexes of the point sources, a conversion factor is computed for each
point source. We sum the four physical components into a "true'' sky map which is then convolved with the experimental beam. Then we add the noise. We end up with one map per frequency band (Fig. 3).
These maps would be what the analysis team would recover from the data, after pointing reconstruction, parasite
removal, decorrelation of instrumental systematics in the data, and mapmaking: much analysis work. In the following,
we explain a set of algorithms optimised for recovering SZ cluster signals in the observed maps.
The simulated observations generated according to the method described in the previous section consist of linear mixtures of the four main sources of diffuse radiation as expected in the four channels of the Olimpo instrument. Focusing on separating the SZ map, we want to make the most of the data in the four frequency bands (centered on 143, 217, 385 and 600 GHz) and exploit the fact that each map is actually a different point of view on the same scene: multichannel data requires processing these maps in a coherent manner. We now provide an overview of the general principles of ICA for component separation from multichannel data. Then, we give a more detailed description of JADE, the specific ICA algorithm we used.
Blind Source Separation (BSS) is a problem that occurs in multidimensional data processing.
The overall goal is to recover unobserved signals, images or sources S from
mixtures of these sources X observed typically at the output of an array of m sensors.
The simplest mixture model takes the form of Eq. (3):
where X and S are random vectors of respective sizes
,
and A is an
matrix.
The entries of S are assumed to be independent random variables. Multiplying S by A linearly mixes the n sources into m observed processes.
Independent Component Analysis methods were developed to solve the BSS problem, i.e. given a batch of T observed samples of X, estimate the mixing matrix A or its inverse and reconstruct the corresponding T samples of the source vector S, relying mostly on the statistical independence of the source processes. Inverting the mixing process requires that the number of mixtures m be greater than or equal to the number of sources n. Solving the underdetermined case where m < n requires some more a priori information on the sources. This is an currently active area of research with close ties to the ideas of sparse representation (Lee et al. 1999; Georgiev et al. 2005; Lewicki & Sejnowski 2000).
Note that with the above model, the independent sources can only be recovered up to a multiplication by a nonmixing matrix i.e. up to a permutation and a scaling of the entries of S. Although independence is a strong assumption, it is in many cases physically plausible. This exceeds the simple second order decorrelation obtained for instance using Principal Component Analysis (PCA): decorrelation is not enough to recover the source processes since any rotation of a white random vector remains a white random vector.
Algorithms for blind component separation and mixing
matrix estimation depend on the a priori model used for the probability
distributions of the sources (Cardoso 2001) although rather coarse assumptions can be made (Cardoso 1998; Hyvärinen et al. 2001).
In a first set of techniques,
source separation is achieved in a noiseless setting, based on the
nonGaussian nature of all but possibly one of the components.
Most mainstream ICA techniques belong to this category: JADE (Cardoso 1999),
FastICA, Infomax (Hyvärinen et al. 2001). In a second set of blind techniques, the components
are modeled as Gaussian processes and,
in a given representation (Fourier, wavelet, etc.), separation requires that the sources
have diverse, i.e. non proportional, variance profiles.
For instance, the Spectral Matching ICA method (SMICA) (Moudden et al. 2005; Delabrouille et al. 2003) considers the case of mixed stationary Gaussian components in a noisy context:
moving to a Fourier representation, colored
components can be separated based on the diversity of their power spectra.
In the case where the main component of interest is well modeled by a
peaked distribution with long tails (e.g. a Laplace distribution) as is the case with SZ maps, methods from the first set are expected to yield better results. Next we provide a description of JADE, the nonGaussian ICA method we used.
The Joint Approximate Diagonalization of Eigenmatrices method (JADE) assumes a linear mixture model as in (8) where the independent sources S are nonGaussian
iid^{} random processes with the additional assumption of a high signal to noise ratio (i.e.
).The mixing matrix is assumed to be square and invertible so that (de)mixing is just a change of basis. Although the noiseless assumption may not be true in the problem at hand, the algorithm may still be applied and a proper change of representation can move us closer to such a setting, as discussed in the next Sect. 4.3.
As mentioned above, second order statistics do not retain enough information for source separation in this context: finding a change of basis in which the data covariance matrix is diagonal will not in general enable us to identify the independent sources properly. Nevertheless, decorrelation is half the job (Cardoso 1998) and one may seek the basis in which the data is represented by maximally independent processes among those bases in which the data is decorrelated. This leads to socalled orthogonal algorithms: after a proper whitening of the data by multiplication with the inverse of a square root of the covariance matrix of the data W, one seeks a rotation R (which leaves things white) so that
defined by

(9) 
and
are estimations of the sources and of the inverse of the mixing matrix.

Figure 4:
SZ component map extracted by JADE from the four observed noisy maps. The SZ cluster signal, subdominant at all observed frequencies, now appears clearly. No obvious leftovers from other astrophysical sources are seen. Remaining noise is small, because we prefiltered data before JADE processing, and we simulated the nominal noise levels of an ambitious project: Olimpo. 
Open with DEXTER 
JADE is such an orthogonal ICA method and, like most mainstream ICA techniques, it exploits higher order statistics so as to achieve a non linear decorrelation. In the case of JADE, statistical independence is assessed using fourth order cross cumulants defined by:
where
stands for statistical expectation and the y_{i} are the entries of vector Y modeled as random variables. Then, the correct change of basis (i.e. rotation) is found by diagonalizing the fourth order cumulant tensor. If the y_{i} were independent, all the cumulants with at least two different indices would be zero. As a consequence of the independence assumption of the source processes S and of the whiteness of Y for all rotations R, the fourth order tensor F is well structured: JADE was precisely devised to take advantage of the algebraic properties of F. JADE objective function is given by



(11) 
which can be interpreted as a joint diagonalization criterion. Fast and robust algorithms are available for the minimization of
with respect to R based on Jacobi's method for matrix diagonalization (Pham 2001). More details on JADE can be found in Cardoso (1998,1999); Hyvärinen et al. (2001).
4.3 JADE in wavelet space
We chose to use JADE after a Wavelet transform. Wavelets come into play as a
sparsifying^{}
transform. Moving the data to a wavelet representation does not affect its information content and applying a wavelet transform on both sides of (8) does not affect the mixing matrix and the model structure is preserved. However, the statistical distribution of the data coefficients in the new representation is different: wavelets are known to lead to sparse approximately iid representations of structured data. Further, the local (coefficient wise) signal to noise ratio depends on the choice of a representation. A wavelet transform tends to grab the informative coherence between pixels while averaging the noise contributions, thus enhancing structures in the data.
Although the standard ICA model is for a noiseless setting, the derived methods can be applied to real data. Performance will depend on the detectability of significant coefficients i.e. on the sparsity of the statistical distribution of the coefficients. Moving to a wavelet representation will often lead to more robustness to noise.
We used the biorthogonal wavelet transform based on the Antonini 7/9 filters (Antonini et al. 1992). Although this choice may not be optimal for the source separation objective, it does enhance separation. Finding the optimal representation of the data for component separation is not a trivial task given the diversity of structures encountered in the different component maps. Some very recent work in BSS considers exploiting the morphological diversity of the source processes by resorting to highly overcomplete dictionaries (i.e. unions of several bases) to enhance source separation in some specific cases (Starck et al. 2004,2005; Zibulevsky & Pearlmutter 2001; Bobin et al. 2005). How this could be useful for astrophysical component separation in experiments such as Olimpo or Planck is currently being investigated.
Once the data has been transformed to a proper representation (e.g. wavelets but also ridgelets and curvelets (Starck et al. 2003), if the 2D or 3D data is strongly anisotropic), we apply the standard JADE method to the new multichannel coefficients. Once the mixing matrix is estimated, the initial source maps are obtained using the adequate inverse transform (Fig. 4).
JADE has been used to separate the signals from four mixtures and four sources in the
presence of Gaussian noise. We have added the expected experimental level of Gaussian noise
to each observed mixture map. In order to prevent bias induced by pixelisation in our simulation and detection algorithms, we overpixelised our observed maps, and thus pixel noise is enhanced. But all algorithms for BSS that require whitening are sensitive to additive noise. In order to minimise the impact of this noise at the ICA step, we chose to convolve our observed map before JADE with a Gaussian of optical beams width.
To optimise signal recovery, we perform a filtering to remove the remaining noise. We tested different methods to remove noise. Thanks to the original simulated SZ map (without noise), we can easily compare the results of those filterings in the next part.
A rather common linear filtering technique uses a Gaussian filter,
generally isotropic. The standard method consists of convolving the
observed map
with a Gaussian window G with standard deviation
:



(12) 
The filtering depends strongly on the value of
.
This value is adjusted arbitrarly or based on a priori information.
An alternative to Gaussian filtering is Wiener filtering which is a
method that attempts to minimize the mean squared error between the
original and the restored signal. Wiener filtering convolves the
observed map
with a weight function i.e. by
assigning the following weight to each mode in Fourier Space:

(13) 
where u, v are the spatial frequencies,
is a model of the map power
spectrum and is in practice derived from the data.
The weight function makes it possible to attenuate or to remove part
of the frequencies if the signaltonoise ratio is low.
The filtering depends on the model of the noise.
The Wiener filter is the optimal filter if both the signal and the
noise are well modeled as Gaussian Random Fields.
This condition is not verified for SZ maps
which display highly nonGaussian features. Nevertheless, Wiener filtering
generally outperforms the simple Gaussian filtering.
In the following paragraphs, we have tested a state of the art nonlinear
filtering method in order to improve signal recovery.
The Maximum Entropy Method (MEM) is commonly used in
astronomy for image processing (see Starck & Murtagh 2002; Marshall et al. 2002; Starck et al. 2001,
for a full description). It is based on entropy and Bayesian methods.
The Bayesian approach provides the means to incorporate prior knowledge
in data analysis. Choosing the prior is one of the most critical aspects
of Bayesian analysis. Here an entropic prior is used.
In data filtering, an entropy functional assesses the information content
of a data set. Among several possible definitions of entropy,
the most commonly used in image processing is the Gull & Skilling (1991)
definition:



(14) 
where k and l are respectively the position and the scale index.
m is a model, chosen typically to be the sky background.
has a global maximum at S=m.
However, MEM
does not allow negative values in the solution, which is a problem in
experimental SZ data analysis, where we measure
fluctuations about zero. To overcome this problem, it has been
proposed (Maisinger et al. 2004) to replace
by:



(15) 
where
.
Here m does not play the same role as in (14).
It is a constant fixed to the expected signal rms.
To overcome the difficulties encountered by the MEM to restore images
containing both high and low frequencies, Pantin & Starck (1996) have
suggested a definition of entropy in a multiscale framework which we
describe in the next section. It has been shown that the main drawbacks
of the MEM (i.e. model dependent solution, oversmoothing of compact objects, ...)
disappear.
 1.
 Multiscale Entropy definition: the Multiscale Entropy method is based on the standard MEM prior derived from the wavelet decomposition of a signal. In this section, we use the
"a trous'' algorithm, with the B^{3}Spline filter that matches well our cluster shapes and is fast to compute. The idea is to consider the entropy of a signal as the sum of the information at each scale of its wavelet transform. The information of a wavelet coefficient is related to the probability of its being due to noise.
The Undecimated Isotropic Wavelet Transform (UIWT) decomposes an
image S as in:
where k, l are the position index and j is the scale index.
C_{J} is a coarse or smooth version of the original image S and w_{j} represents the details
in S at scale 2^{j}
(see Starck & Murtagh for details). Thus, the algorithm outputs J+1 subband arrays of size
.
We will use an indexing convention such that j = 1 corresponds to the finest scale (high frequencies).
Denoting as H(S) the information relative to the signal and
h(w_{j}(k,l)) the
information relative to a single wavelet coefficient,
the entropy is now defined as:



(16) 
where l is the number of scales and n is the size of map.
 2.
 Entropy definition:
the function h in (16) assesses the
amount of information carried by a specific wavelet coefficient.
Several functions have been proposed for h.
A discussion and comparison between different entropy definitions
can be found in Starck et al. (2006).
We choose the NOISEMSE entropy Starck et al. (2001) for the SZ
reconstruction problem in which the entropy is derived
using a model of the noise contained in the data:



(17) 
where
P_{n}(w_{j}(k,l)) is the probability that the coefficient
w_{j}(k,l) can be due to the noise:
.
For Gaussian noise, we have:
P_{n}(w_{j}(k,l)) 
= 



= 

(18) 
and Eqs. (17) becomes (19)



(19) 
The NOISEMSE is very close to the l_{1} norm
(i.e. absolute value of the wavelet coefficient) when the coefficient
value is large, which is known to produce good results for the
analysis of piecewise smooth images (Donoho & Elad 2003).
In the case of nonGaussian or spatially inhomogeneous noise, our implementation will take a noise map as input. Then we can easily derive its entropy using the above definition.
 3.
 Signal and noise information: the SZ component obtained by blind separation is swamped by noise.
The following algorithm assumes that the observed map can be decomposed as:



(20) 
Then, we can decompose the information contained
in our image in two components, the first one (H_{s}) corresponding to the noncorrupted part,
and the other one (H_{n}) describing a component which contains no
information for us:



(21) 
For each wavelet coefficient
w_{j}(k,l), we have to estimate the fractions h_{n}and h_{s} of h:



(22) 
5.2.3 Multiscale Entropy Filtering
 1.
 Filtering:
the problem of filtering the observed map
can be expressed as follows.
We look for a filtered map
such that the difference between
and
minimizes the information due to the signal (to recover all
the signal) and such that
minimizes the information due to the noise
(we want no noise). These two requirements are competing. A tradeoff is necessary, because, on the one hand, we want to remove all the noise (heavy filtering) and on the other hand, we want to recover the signal with fidelity. In practice, we minimize for each wavelet coefficient
w_{j}(k,l):
where
w_{j}(k,l) are the wavelet coefficients of the observed map
,
the wavelet coefficients of the filtered map
and
is the socalled regularization (tradeoff) parameter.
 2.
 Selecting significant Wavelet coefficients:
whatever the filtering, the signal is always substantially modified. We want to fully reconstruct significant structures, without imposing strong regularization while eliminating the noise efficiently The introduction of the multiresolution support (Murtagh et al. 1995), helps to do so. It applies the previous regularization (i.e. filtering) only on the nonsignificant (noisy) wavelet coefficients (Pantin & Starck 1996). Thus the choice of the regularization parameter is not a critical point. The other components of the maps are left untouched.
The new Multiscale Entropy becomes:



(23) 
where
,
and M is the "multiresolution support'' defined as:



(24) 
M describes, in a Boolean way, whether the data contains information at a given scale j and at a given position (k,l).
w_{j}(k,l) is said to be significant if the probability that the wavelet coefficient is due to noise is small. In the case of Gaussian noise, a coefficient
w_{j}(k,l) is significant if
,
where
is the noise standard deviation at scale j, and k is a constant. Without an objective method for selecting the threshold, it is adjusted arbitrarly, generally taken between 3 and 5 (Murtagh et al. 1995).
 3.
 Selecting significant Wavelet coefficients using the FDR:
the False Discovery Rate (FDR) is a new statistical procedure by
Benjamini & Hochberg (1995) which offers an effective way to select an adaptative
threshold to compute the multiresolution support . This technique has been described by
Hopkins et al. (2002); Starck et al. (2006); Miller et al. (2001) with several examples of
astrophysical applications. The FDR procedure provides the means to adaptively
control the fraction of false discoveries over total discoveries.
The FDR is given by the ratio (25), that is, the proportion of declared
active pixels which are false positives:



(25) 
where
is the number of pixels truly inactive declared active,
and
is the number of pixels declared active.
The FDR formalism ensures that, on average, the False Discovery Rate is no larger than
which lies between 0 and 1. This procedure guarantees control over the FDR
in the sense that:
The unknown factor
is the proportion of truly inactive pixels. A complete description of the FDR method can be found in Miller et al. (2001). In Hopkins et al. (2002) and Starck et al. (2006), it has been shown that the FDR outperforms standard method for source detection.
In this application, we use the FDR method in a multiresolution framework
(see Starck et al. 2006).
We select a detection threshold
for each scale.
The different values of
at each scale were set as in Starck et al. (2006) and left untouched during our work:
where
.
A wavelet coefficient
w_{j}(k,l) is considered significant if its absolute value
is larger than
as seen below.
 4.
 Multiscale Entropy Filtering algorithm:
assuming Gaussian noise, the Multiscale Entropy restoration method
reduces to finding the image
that minimizes
,
given the map
output of source separation with:
where
is the noise standard deviation in
,
J is the number
of Wavelet scales,
is the Wavelet Transform operator and
is the multiscale entropy only defined for nonsignificant coefficients (outside the support selected by the FDR thresholding). Full details of the minimization algorithm can be found in Starck et al. (2001).
The results presented in the next section are obtained on a SZ map with a uniform Gaussian
white noise but the method still holds for a nonuniform Gaussian noise over the map.
6.1 Extraction algorithm  SExtractor
We use a public source software for extraction of SZCluster candidates from the filtered maps: SExtractor (Bertin 2003). SExtractor turned out to be fast, convenient, and easy to configure. It handles conveniently Cluster Extraction, with observed sizes ranging from a few pixels up to angular diameters of a degree.
We used the Noiseless configuration, assuming that our filters are efficient enough for this assumption to be valid. We mainly use the source identification capability and the deblending algorithms that be very useful when extracting largemass clusters. In order to recover as many candidates as possible, we set
and
(the smallest cluster size after optical beam dilution is over 9 pixels). Photometry of Cluster candidates is done using SExtractor's FLUX_AUTO Mode.
SExtractor photometry performs very well on our data: the recovery error on integrated Compton
flux Y is smaller than 4 percent, for Y larger than
sr (id
arcmin^{2}).
Thus SExtractor outputs a catalog of Cluster candidates with their position, flux and size. This catalog is likely to be contaminated by false detections due to residual noise, or spurious point sources that have survived the Source Separation step. While processing observed data, we will have to live with this contamination. Using simulations, we quantify in the following the SZCluster selection function (Completeness) and this contamination.
6.2 Completeness and Purity of the recovered catalog
In order to identify true SZ Clusters from contamination in our catalog, we used an association criterion: for each candidate cluster, we scan the generated cluster catalog (see Sect. 2.4) for the closest cluster and store the angular distance between the two. Doing so, we observed a distribution with a strong peak at a distance smaller than 4 pixels (3.5 arcmin), corresponding to the true detected clusters, and a long flat tail corresponding to random associations. We then decided to tag candidates as true when their distance to the closest cluster in the simulated catalog is less than 4 pixels, neglecting the small fraction of random associations passing this test. This criterion can and will be enhanced using an improved distance involving recovered cluster flux and spatial distance.
Purity is defined as the ratio of true detections to the total number of detections:

(28) 
Completeness is defined as the ratio of true detections to the total number of clusters in the simulation:

(29) 
Completeness and purity are the main performance criteria of a detection chain and should be maximised. We will use them to compare detection chains in the results section.

Figure 5:
Input SZ cluster map ( upper left) map, and the 3 maps as recovered by Completeness and ( upper right) Gaussian filtering
,
bottom left, Wiener filtering and bottom right MEFDR method. We chose to plot 25 maps, to point out differences that do not appear on 400 maps. 
Open with DEXTER 
In the following, we will quantify the performance of our algorithm chain in two steps. We will focus on the SZ map reconstruction: first the source separation step with JADE and then the additional FDRMultiScale Entropy Method (MEFDR) filtering step. The MEFDR filter will be compared first to the classical Gaussian and Wiener filters. Then we will discuss assumptions and performances of our method relative to the one published in Pierpaoli et al. (2005).
The comparison between two methods for astrophysical image reconstruction can be performed in
several ways.
One may initially investigate the characteristics of residuals. This technique is typically used by mathematicians.
An astrophysicist may be more inclined to assess how well the objects of interest (SZ cluster of galaxies, in this case) are recovered: how many of them are found, and with what precision.
Finally, a cosmologist would want to know to what extent the reconstruction procedure limits the ability to estimate cosmological parameters. This implies consideration of the purity and completeness as functions of cluster mass threshold, as well as precision in the reconstruction of the (central or integrated) Compton parameter for such masses.
In what follows, we will compare the map's reconstruction error. Then we will focus on the
SZ Cluster detection performance of the full detection chain, by comparing the recovered source catalogs to the Simulated SZ cluster catalog, in terms of purity and Completeness as defined in Sect. 6.2. This allows a parameterfree quantification of performances and tradeoffs involved in cluster detection, being mostly independent of photometry details.
Note also that the cluster detection efficiency does not convey information on the spread of
the input/output y parameter. A more detailed description of the selection function, contamination, a discussion on photometry issues and methods and of cosmological implications will be given in a following paper (Juin et al. 2005).
7.1 JADE
Figure 4 shows the recovered SZ map, after source separation. JADE performs very well on our maps: the SZ signal that subdominant at all scales now appears clearly. No obvious residuals from other astrophysical sources are seen. JADE assumes and requires that input data have zero mean. This point is easy to meet in bolometer experiments since our bolometer cameras measure flux variations while the telescope scans the sky.
Second, JADE loses calibration information while processing data. Sky components are separated, but calibration of each output map has to be restored. It is of common observational procedure to calibrate a survey on the brightest sources in the field (which can be observed in follow up experiments). We chose in the following to calibrate our SZ maps using the 100 brightest clusters in the field. We average their recovered integrated flux (SExtractor FLUX_AUTO mode), and scale the map to match the average integrated flux of the 100 brightest simulated clusters. In the following, before computing statistical tests, all the maps have been normalised in this way. We will suppress this feature by taking into account prior knowledge (the optical spectral dependance of the SZ component) in future work.
Third, it is very important to minimise noise in the observed maps before source separation (JADE). JADE was designed to run on noiseless data. And it is quite sensitive to noise. We chose to apply to the four observed maps, before JADE, a simple Gaussian filter with the widths of the optical beams. If the noise level is not reduced before JADE, then the recovered mixing matrix will be inaccurate, and the SZ Cluster map will be polluted by the remains of the other astrophysical components, inducing efficiency loses. Once the mixing matrix has been carefully estimated, one can apply it to the unfiltered observed maps to extract the observed SZ cluster map, and then apply an optimised filter and cluster detection algorithms to extract the cluster catalog.
We characterize in our simulations the residual error on the estimated SZ map after component separation. The computed error map has a variance of 0.888, a skewness of 0.0019, a kurtosis
of 0.0156 and is essentially white noise with no obvious visible structures. The foreground contribution to the noise map has a variance of 0.00273 (skewness of 0.0019, kurtosis of 0.0156) and shows that the component separation is very efficient on our simulated data, and strengthens the assumption of Gaussian noise made in the subsequent nonlinear filtering and data processing steps. Should the experimental noise be significantly nonGaussian, the presented detection chain could still be used. The price to pay will be a lack of statistical efficiency. The overall performance of the detection chain has to be assessed by means of simulations as shown in our paper.
Additional filtering is then applied to recover smaller clusters. We now quantify performances of the MEFDR filter relative to the simpler filtering methods (Gaussian and Wiener). We chose to show results obtained with maps processed by JADE running in wavelet space, our best method. We will quantify the results by computing an error map, and its properties. Then we will compare the output of the extraction procedures using the three filters.
Figure 5 shows the maps after two classical filtering methods and our FDR multiscale entropy method. All maps look similar in nominal noise conditions of such ambitious experiments as future SZ Cluster survey, only statistical tests can discriminate between the 3 filtering methods. We computed the error map (filtered map, input map subtracted) for each filtered option. Computing the rms of the error maps divided by the rms of the simulated SZ map (dimensionless rms) leads to
,
and
.
The CPU time on a
pixels map was
s,
s, and
s. Thus our filtering algorithm is reasonably fast.
A more accurate test is to plot (see Fig. 6) the dimensionless standard deviation of error maps (filtered map, input map subtracted) at each scale of the wavelet decomposition:

(30) 
where
and
,
are the standard deviation of the error map and the input SZ map, selecting the scale L of our wavelet transform. MEFDR method performs better than simpler filters at all scales. Running our algorithms while changing the noise level showed that the higher the noise level, the larger the advantage in using our MEFDR filtering compared to simple Gaussian, or Wiener filters.
Then considering our goal of detecting clusters, the relevant test is to compare the recovered
catalog to the simulated cluster catalog, the input of the simulated map. Observed cluster catalogs are extracted and processed as explained in Sect. 6.2.
In a future paper, we will present the selection function, contamination, a discussion on photometry issues and constraints on cosmology.
In this paper where we discuss detection procedures the relevant information is the curve of purity versus completeness while raising the threshold on the cluster observed flux
(Fig. 7). As expected, the higher the threshold, the higher the purity, and the lower the completeness.
We see once more that MEFDR out performs the other filter options, allowing greater completeness at any required level of purity.
In what follows, we will compare the method described above with the one presented in
Pierpaoli et al. (2005).
We will compare the general performance of the two methods in reconstructing the SZ images and
then assess the global purity and completeness, as defined in Sect. 6.2.
As we are not assessing here the goodness of the Y parameter reconstruction, a full comparison with the Pierpaoli et al. (2005) results (i.e. purity
and completeness as a function of cluster mass threshold, spread of
the input/output Y parameter relation ) is not possible at this time.
We will first describe the Pierpaoli et al. (2005) method, and then preform the comparison
according to the mentioned criteria.

Figure 6:
Dimensionless standard deviation of error maps (filtered map, input SZ map subtracted) at each scale of the wavelet decomposition, for the three studied filters options. MEFDR filtering method (red squares) is a significant improvement compared to simpler filtering methods (Gaussian, black dots and Wiener, blue triangles). 
Open with DEXTER 
The method presented in Pierpaoli et al. (2005) is formally a nonblind component separation that has been optimized to recover galaxy clusters.
In this method component separation and the deconvolution of the beam effect are done in one step by computing the Bayes Least Square estimate under a Gaussian scale mixture model of
"neighborhoods'' of wavelet coefficients. These "neighborhoods'' are sets of wavelet coefficients associated with the same location and that behave in a coherent manner. The optical frequency dependency of each component, the beam size and the noise level for each observation are assumed to be known. The procedure relies on the possibility to discriminate between the different components (CMB, clusters, infrared point sources and galaxy dust) by modeling the joint probability of these neighborhoods by Gaussian scale mixtures (see Pierpaoli et al. 2005 for details).
A Gaussian scale mixture is a random vector
,
where u is a Gaussian vector independent of the scalar random variable z  allowed to be
nonGaussian.
The model is determined by the probability distribution p_{z} and the covariance matrices of u. The covariance matrices have to be adjusted at each scale of the wavelet decomposition as a function of the resolution of the observed maps. They are computed prior to running the algorithm on simulated maps of
each component at the correct resolution. The probability distribution p_{z} is also componentdependent. Pierpaoli et al. (2005) investigate several different possibilities for p_{z}, showing
that different distributions lead to very different results in the mass reconstruction.
For example, the Gaussian assumption (
in the following) is most able to recover the maximum number of clusters but is inefficient in reconstructing the appropriate central intensity of these clusters; a nonGaussian distribution modeled on the simulated SZ map (
)
is most able to recover the input Comptonisation parameter for bright clusters while missing the reconstruction of very low intensity ones.

Figure 7:
Purity versus completeness for the three filtering options presented in this paper. Completeness + MEFDR (red squares) out performs the other filter options Gauss (black dots) and Wiener (blue triangles), allowing more completeness at any required level of purity. In the example presented here we simulated 9942 clusters according to
distribution, with a threshold value on
of
st: note that the completeness is normalized to this total number of simulated clusters. The insert is a zoom on the high purity region of the graph. 
Open with DEXTER 

Figure 8:
Maps computed by two of the algorithms (25 deg^{2} size). Left is the input SZ Cluster map. Centered is the map output of completeness after MEFDR filtering. Right is the map output of the Pierpaoli et al. method, using the Gaussian assumption for the probability distribution. 
Open with DEXTER 
In Pierpaoli et al. (2005) the focus is on reconstructing the signal of the most massive clusters,
since those are the ones which are likely to lead to the most precise constraints on cosmological parameters. To this aim, the distribution
is preferable to
.
As this paper is more focused on cluster detection efficiency (also including lowintensity clusters) and disregards the precision of the reconstruction, we will consider here the Gaussian prior for p_{z} (
in what follows).
A Gaussian prior for p_{z} reduces the estimator to a local Wiener filtering in
Wavelet space, as all information about the signal nonGaussianity is lost.
The result of Pierpaoli et al. (2005) is a set of beamdeconvolved maps (one for each physical component considered) which can be directly
compared with the input ones, or with other method's.
The cluster's y map is used here for comparison. As for cluster detection, we will present results obtained with the code used in Pierpaoli et al. (2005), as well as the SExtractor code adopted here.

Figure 9:
Standard deviation versus scale for SZ cluster maps recovered by the completeness MEFDR method (red squares) and Pierpaoli et al.
method (blue triangles). The top axis gives the wavelet's typical scale (in arcmin) corresponding to the index on the bottom axis. 
Open with DEXTER 
In Fig. 8 we present the maps that will be used for the comparison.
While both maps recover well the intense clusters, the
processed map shows more lowintensity structure than the MEFDR map.
This could be due to a better reconstruction of lowintensity clusters, as well as undesirable noise.
Computing the
error map rms, we find
,
which should be compared MEFDR.
An analysis of the residuals' normalised standard deviation for different wavelet scales (shown in Fig. 9) shows that the two methods mainly differ at angular scales equal and above 14 arcmin. The major contribution to the difference in the total rms is therefore not associated with the
typical (low mass) cluster scale, but to much larger ones.

Figure 10:
Purity versus completeness curves for the four software detection chains studied in the following. We use red squares for CompletenessFDR+SExtractor, black dots for
and blue triangles for
with peak detection, black diamonds for
using SExtractor for source detection. The inserted plot is a zoom on the high purity region of the graph. 
Open with DEXTER 
We compare here the detection efficiency defined in Sect. 6.2.
For consistency with our detection chain (see Sect. 7.1), we recalibrate
the Pierpaoli maps so that the average integrated Compton flux of the 100 brightest clusters matches the value observed in the input SZ cluster map (small correction).
Figure 10 presents the curves of completeness versus purity for completeness
FDRME + SExtractor,
,
using the peak finding algorithm presented in Pierpaoli et al. (2005) to detect clusters and the
method using SExtractor (PierNew).
As expected,
is not as good as the other methods at all purity, since the
probability p_{z}used here is optimized for accurate recovery of the most massive cluster's central intensity
and not for cluster detection.
(the Gaussian distribution) improves
results, especially when SExtractor is used to detect the
clusters (PierNew).
At low threshold (low purity) PierNew recovers more clusters than MEFDR.
The lowintensity structures visible in Fig. 8 contain
a sizable number of clusters. As completeness FDRME is designed to filter out noise better, it also filters out some lowintensity clusters which then no longer contribute to the detection rate.
At high purities, PierNew and completeness FDRME provide equivalent performance, given the statistical uncertainly of the small number of clusters involved in the highthreshold cut.
We simulated observed sky maps at frequencies of a typical large multiband bolometric SZCluster survey. This simulation includes the main features of the "background'': point sources with varying spectral index and brightest sky components. We implemented a complete software detection chain, working in 3 steps. First we use a source separation algorithm that is based on a Wavelet transform and the JADEICA algorithm. Then, we filter the SZCluster map using an FDRMultiscale Entropy method. Finally, we detect the clusters on the filtered maps using the SExtractor software in a noiseless configuration. Despite several simplifying assumptions detailed in the paper, our detection chain proves to be very efficient, yielding 25% of clusters of mass larger than
detected with 90% purity.
We compare our detection algorithm to a previously published method based on wavelets (Pierpaoli et al. 2005).
We restrict our comparison criteria to the global detection efficiency, as defined in Sect. 6.2. We find a detection efficiency in the
highpurity region comparable to the Gaussian probability case (
method), which is the model that provides the best performance for this comparison among those presented in Pierpaoli et al. (2005).
The MEFDR detection efficiency slightly degrades at lower cluster intensity than the
method, as MEFDR filters out lowintensity clusters during the denoising procedure.
These results, however, are impressive as MEFDR, unlike Pierpaoli et al. (2005) methods, is a blind algorithm that makes no assumption on the physical properties of the signal to be recovered.
There remain few open issues concerning the use of the proposed detection chain to process upcoming real observations. Once configured, the detection chain was used without
"close supervision''. Also, the wavelet transform used jointly with JADE was chosen a priori. JADE as well as the MEFDR filter were run as black box algorithms, while the SExtractor parameters were set to fairly standard values, since the input data maps were prefiltered. We also left out the smallest scales in the MEFDR maps to reduce high frequency spatial noise. If the experimental data are not as clean as in our simulations, how will these simplifications affect the result?
If the noise level is increased or if the noise displays strong "striping'' (as a result of low frequency noise in the bolometer timelines), JADE is no longer appropriate for the component separation step. Hence the experimental data processing will have to include a good "destriper'' in the mapmaking calculations. The proposed MEFDR filter proved to be more robust than the Gaussian or Wiener filter, when the noise level is increased.
Our last concern is that the Galactic Dust Maps in our simulations are generated as spatially Gaussian, which is not true of observed "Galactic dust clouds''. But the optical spectrum of Galactic Dust being very different from SZ clusters, we foresee that ICA will separate these features efficiently.
We will focus on improving the source separation step. We will include optical spectrum knowledge in the source separation algorithms. This will allow us to overcome the intrinsic scaling indeterminacy of the blind linear mixture model and so prevent loosing track of map calibration as in
Sect. 7.1. Additionally, JADE must run on noisefree data. We will design our source separation algorithm to handle noisy data. We studied the performance of our method in terms of detection efficiency. A further assessment in terms of reconstruction accuracy and photometry is necessary since these are valuable for cosmology. Of course detection efficiency doesn't provide all the information we want to know in order to do cosmology, as the accuracy of the reconstruction and photometry issues are also important. In a following paper (Juin et al. 2005), we will use these algorithms to compute selection functions, contaminations and constraints on cosmology foreseen for the upcoming bolometric SZcluster surveys.
Acknowledgements
The authors wish to thank A. Réfrégier, R. Teyssier, J. Rich, C. Magneville, J. P. Pansart, J. L. Starck (CEA/DAPNIA), J. B. Melin and J. Bartlett (Univ.P7, APC), E. Bertin (IAP), for many private discussions and the reviewer for helpful comments on an earlier version. E.P. is an ADVANCE fellow (NSF grant AST0340648), also supported by NASA grant NAG511489.
 AMIBA 20052007,
AMIBA experiment, see: http://amiba.asiaa.sinica.edu.tw/

Antonini, M., Barlaud, M., Mathieu, P., & Daubechies, I. 1992,
Vol. I
(In the text)

Benjamini, Y., & Hochberg, Y. 1995, J. R. Stat. Soc. B, 57,
289 (In the text)
 Bertin, E.
2003, http://terapix.iap.fr/
(In the text)
 Bobin, J.,
Moudden, Y., Starck, J.L., & Elad, M. 2005, in SPARS'05,
Signal Processing with Adaptative Sparse Structured
Representations, Vol. 5914, Rennes, France
 Borys, C.,
Chapman, S. C., Halpern, M., & Scott, D. 2003, MNRAS, 344,
385 [NASA ADS] [CrossRef] (In the text)
 Bouchet,
F. R., & Gispert, R. 1999, New Astron., 4, 443 [NASA ADS]
[arXiv:astroph/9903176]
(In the text)
 Cardoso,
J.F. 1998, Proc. IEEE, Special issue on blind identification and
estimation, 9, 2009
 Cardoso, J. F.
1999, Neural Computation, 11, 157 [CrossRef] (In the text)
 Cardoso,
J. F. 2001, in Proc. ICA 2001, San Diego
(In the text)
 Contaldi, C. R.,
Bond, J. R., Pogosyan, D., et al. 2002, Proc. XVIII IAP
Coll. On the nature of dark energy, Paris, 15 July 2002
[arXiv:astroph/0210303]
(In the text)
 Cooray, A. R. 2000,
MNRAS
(In the text)
 Delabrouille,
J., Cardoso, J.F., & Patanchon, G. 2003, MNRAS, 346, 1089 [NASA ADS] [CrossRef]
 Donoho,
D. L., & Elad, M. 2003, Proc. Nat. Aca. Sci., 100,
2197 [NASA ADS] (In the text)
 Doran, M.,
& Mueller, C. M. 2004, Journal of Cosmology and
Astroparticle Physics, 0409, 003 [CrossRef]
[arXiv:astroph/0311311]
(In the text)

Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al.
2005, ApJ, 633, 560 [NASA ADS] [CrossRef] (In the text)
 Georgiev,
P. G., Theis, F., & Cichocki, A. 2005, IEEE Transactions
on Neural Networks, 16, 992 [CrossRef]
 Gull,
S., & Skilling, J. 1991, MEMSYS5 Quantified Maximum Entropy
User's Manual, Royston, England
(In the text)
 Halverson,
N. W., Leitch, E. M., Pryke, C., et al. 2002, ApJ,
568, 38 [NASA ADS] [CrossRef] (In the text)
 Hanany, S.,
Ade, P., Balbi, A., et al. 2000, ApJ, 545, l5
(In the text)

Herranz, D., Sanz, J., Hobson, M. P., et al. 2002a, MNRAS,
336, 1057 [NASA ADS] [CrossRef]

Herranz, D., Sanz, J. L., Barreiro, R. B., &
MartinezGonzález, E. 2002b, ApJ, 580, 610 [NASA ADS] [CrossRef]
 Hinshaw, G.,
Spergel, D. N., Benett, C. L., et al. 2003, ApJS, 148,
97 [NASA ADS] [CrossRef] (In the text)

Hopkins, A. M., Miller, C. J., Connolly, A. J.,
et al. 2002, AJ, 123, 1086 [NASA ADS] [CrossRef]

Hyvärinen, A., Karhunen, J., & Oja, E. 2001, Independent
Component Analysis (New York: John Wiley), 481
 Jones, M. E.
2001, MNRAS, in 8th Taipei Astrophysics Workshop AMiBA 2001
(In the text)
 Juin, J. B., Yvon, D., &
Refregier, A. 2005, A&A, submitted
[arXiv:astroph/0512378]
(In the text)
 Knop, R. A.,
Aldering, G., Amanullah, R., et al. 2003, ApJ, 598, 102 [NASA ADS] [CrossRef] (In the text)
 Kuruoglu, E.,
Bedini, L., Paratore, M., Salerno, E., & Tonazzini, A. 2003,
Neural Networks, 16, 479 [CrossRef]
 Lamarre, J.,
Puget, J., Bouchet, F., et al. 2004, New Astronomy Reviews,
Proc. of the workshop on The Cosmic Microwave Background and its
Polarization, ed. S. Hanany, & R. A. Olive
(In the text)
 Lee,
T.W., Lewicki, M., Girolami, M., & Sejnowski, T. 1999, IEEE
Signal Processing Letters, 6, 87 [CrossRef]
 Lewicki,
M. S., & Sejnowski, T. J. 2000, Neural Computation,
12, 337 [CrossRef]
 Maino,
D., Farusi, A., Baccigalupi, C., Perrotta, F., & Banday, A.
2002, MNRAS, 334, 53 [NASA ADS] [CrossRef]

Maisinger, K., Hobson, M. P., & Lasenby, A. N. 2004,
MNRAS, 347, 339 [NASA ADS] [CrossRef] (In the text)

Marshall, P. J., Hobson, M. P., Gull, S. F., &
Bridle, S. L. 2002, MNRAS, 335, 1037 [NASA ADS] [CrossRef]
 Masi, S., Ade,
P., de Bernardis, P., et al. 2003, Mem. Soc. Astron. It.,
74, 96 [NASA ADS] (In the text)

Melin, J.B. 2004, Amas de galaxies et effet SunyaevZel'dovich:
observations et étude des effets de sélection des
sondages, Univ. Denis Diderot, Paris VII
(In the text)
 Miller,
C. J., Genovese, C., Nichol, R. C., et al. 2001, AJ,
112, 3492 [NASA ADS] [CrossRef]
 Moudden, Y.,
Cardoso, J.F., Starck, J.L., & Delabrouille, J. 2005, Eurasip
Journal on Applied Signal Processing, 15, 2437 [CrossRef]

Murtagh, F., Starck, J.L., & Bijaoui, A. 1995, A&AS, 112,
179 [NASA ADS] (In the text)

Netterfield, C., Ade, P., Bock, J., et al. 2002, ApJ, 571,
604 [NASA ADS] [CrossRef] (In the text)
 Pantin, E.,
& Starck, J.L. 1996, A&AS, 315, 575 [NASA ADS] (In the text)
 Peacok,
J. 1999, MNRAS, 284, 885 [NASA ADS]

Peebles, P. 1993 (Princeton: Princeton Univ. Press)
 Pham, D.T.
2001, SIAM Journal on Matrix Analysis and Applications, 22,
1136 [CrossRef] (In the text)

Pierpaoli, E., Anthoine, S., Huufenberger, K., & Daubechies, I.
2005, MNRAS, 359, 261 [NASA ADS] [CrossRef] (In the text)
 Pierpaoli, E.,
Scott, D., & White, M. 2001, MNRAS, 325, 77 [NASA ADS] [CrossRef]

Pierpaoli, E., Borgani, S., Scott, D., & White, M. 2003,
MNRAS
(In the text)
 Pierre, M.,
Valtchanov, I., et al. 2004, Journal of Cosmology and
Astroparticle Physics, 09, 11 [NASA ADS] [CrossRef] (In the text)
 Refregier,
A., & Massey, R. J. R. 2004, AJ, 127, 3102 [NASA ADS] [CrossRef] (In the text)
 Riess,
A. G., Filippenko, A. V., Li, W., & Schmidt,
B. P. 1999, AJ, 118, 2668 [NASA ADS] [CrossRef] (In the text)

Schafer, B. 2004, MNRAS, submitted
(In the text)
 Sheth, R. K., Mo,
H. J., & Tormen, G. 2001, MNRAS
(In the text)
 Starck,
J.L. & Murtagh, F. 2002, Astronomical Image and Data Analysis
(SpringerVerlag)
 Starck,
J.L., Murtagh, F., Querre, P., & Bonnarel, F. 2001, A&A,
368, 730 [EDP Sciences] [NASA ADS] [CrossRef]
 Starck,
J.L., Candes, E., & Donoho, D. 2003, A&A, 398, 785 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
 Starck,
J.L., Elad, M., & Donoho, D. 2004, Advances in Imaging and
Electron Physics, 132
 Starck, J.L.,
Moudden, Y., Bobin, J., Elad, M., & Donoho, D. 2005, in SPIE
Conference on Wavelets, 5914, San Diego, CA, USA
 Starck,
J.L., Pires, S., & Refregier, A. 2006, A&A, 451, 1139 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
 Sunyaev,
R., & Zeldovich, I. B. 1980, ARA&A, 18, 537 [NASA ADS] (In the text)
 Tristram,
M., Patanchon, G., MaciasPerez, Y. F., et al. 2005, A&A,
436, 785 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)

White, M., Carlstrom, J. E., Dragovan, M., & Holzapfel,
W. L. 1999, ApJ, 514, 12 [NASA ADS] [CrossRef]
 Yvon, D.,
& Mayet, F. 2005, A&A, 436, 729 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)

Zibulevsky, M., & Pearlmutter, B. 2001, NeuralComputation, 13,
863 [CrossRef]
Copyright ESO 2006