EDP Sciences
Free Access
Issue
A&A
Volume 589, May 2016
Article Number A2
Number of page(s) 10
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201527923
Published online 04 April 2016

© ESO, 2016

1. Introduction

Astronomical objects are often seen merged or blended on the plane of the sky. This blending can be apparent, because objects at different distances are seen in projection on the plane of the sky or, real, because different objects at the same distance are physically overlapping.

Whatever the reason for the blending, reliable deblending techniques are mandatory for astrophysical projects to meet their scientific objectives. Among the many possible examples, blends of galaxies of different colours can impact performances of photometric redshift algorithms (e.g. Pérez-González et al. 2010; Bellagamba et al. 2012; Parker et al. 2012; Hsu et al. 2014) and conclusions of stellar populations studies (e.g. Yan et al. 2014). Obviously, blending also affects the determination of morphological properties of astronomical objects, for example the shape measurement of faint galaxies in weak lensing cosmological surveys (e.g. Chang et al. 2013; Arneson 2013). In strong gravitational lensing, deblending of the foreground lensing object from the background lensed sources is essential, for example as shown at galaxy scale by Gavazzi et al. (2007) and at cluster-scale by Massey et al. (2015). This is true for at least two reasons. First, one needs to map the visible mass in the lensing object precisely, either to use it as a prior to guide the lens modelling or to infer the mass-to-light ratio in the lens. Second, the image of lensed source must be isolated in the best possible way. Any faint extended arc-like structure, clump, or star-forming region must be seen precisely with minimum light contamination from the lensing object. Our ability to constrain the mass model is completely driven by the amount of details seen in the lensed source, which represent as many observational constraints.

Many of the current techniques to deblend astronomical objects are limited to analytical modelling of their light distribution either in single band (e.g. PSFex, Bertin 2011) or multi bands, sometimes including a simultaneous fit of many overlapping objects (Megamorph; Vika et al. 2013, 2015). Alternatively, some methods make use of high resolution images to flag blended objects and then measure them at different wavelengths using images of lower spatial resolution (e.g. Laidler et al. 2007). Popular softwares like Sextractor (Bertin & Arnouts 1996) use image segmentation to separate blends, which is a technique that was further improved by Zheng et al. (2015). Other techniques include machine learning, recently used in the area of strong gravitational lensing to subtract the light of bright galaxies and to unveil possible lensed background objects without invoking any analytical representation of the galaxies to subtract. This technique is based on a principal component decomposition of the galaxy images using large samples of single-band imaging data (Joseph et al. 2014). A step forwards is to use multi-band images to separate the objects in the lens plane and source plane, also using the colour information. Recent methods have started to make use of multi-band information and combine source position from high resolution images with profile fitting to deblend lower resolution bands (Merlin et al. 2015). Another example is given in Hocking et al. (2015) where neural networks are used to identify objects with different colours in multi-band images.

With the current burst of wide-field sky surveys (DES, KIDS, HSC, Euclid, LSST and WFIRST), data in many optical and near-IR bands will become available. The present paper describes a technique taking advantage of these multi-band data to address the deblending problem using both colour information and a spatial prior, but not involving any analytical modelling of the light distribution of the lensed and lensing objects. Our work is based on a multi-channel extension of the morphological component analysis (MCA) presented in Starck et al. (2004). We illustrate the performances of the algorithm with numerical experiments and with real Hubble Space Telescope (HST) data of strong gravitational lenses.

This paper is organised as follows: in Sect. 2, we introduce the mathematical model we use to understand and separate objects with different colours. In Sect. 3, we describe the mathematical technique used to solve the problem of colour separation in our approach, that is to say morphological component analysis. In Sect. 4, we detail our implementation of the MuSCADeT algorithm. Section 5 shows the performance of our algorithm on simulations that test realistic problems encountered in deblending. We apply our method to real astronomical images from the HST in Sect. 6 with the galaxy clusters MACS J1149+2223 and Abell 2744. We compare our results with current model fitting methods. Section A provides useful information to reproduce our results from the code we made freely available.

2. Deblending and strong lensing

2.1. The source separation problem

We assume the observed data { yi } i = 1,..,Nb in the band i can be represented as (1)where oj are the different observed sources, si,j is the contribution of the jth source in the observation yi, No the number of sources, Nb the number of bands, zi is an additive Gaussian noise, k is the pixel index (k = 1...Np), and Np is the number of pixels.

The parameter s,j corresponds to the spectral energy distribution (SED) of the source oj. The deblending problem consists in finding the different objects oj, which is somewhat complicated since their SED are not known and even the number of objects is not known.

However, several galaxies may have similar colour properties and, therefore, share the same SED, so we can simplify Eq. (1) by considering the data containing only Ns groups of sources, such as, for instance early- and late-type galaxies, and we can restrict the deblending problem to only extract these two groups. We note xj (j = 1..Ns) the image which contains the sum of all objects belonging to the group j, i.e. , where is the number of sources in the group j. We can write (2)Even if this equation looks very similar to Eq. (1), it is in fact simpler since Ns is smaller than No. As a given component xj contains several astrophysical sources, it also gives us more statistics to derive its SED. This linear mixture model can be recast in the following matrix form (3)where Y is a Nb × Np matrix, A is the SED mixing matrix, and X is the Ns × Np matrix, which contains the components xj.

To sum up, we consider that each band is a weighted mix of ns colour components. In the statistical literature, each component is called a source (and it should not be mixed with an astrophysical source). This general problem is called a blind source separation problem (BSS), i.e. estimating both A and X knowing only Y. The weight for a given source is the value of the associated SED at the corresponding wavelength. Figure 1 illustrates the BSS in the case of two sources relative to two populations of galaxies (red and blue galaxies in the figure).

thumbnail Fig. 1

Illustration of the blind source separation in the case of two sources. To make the figure simple, the images in each band are represented as lines in the Y matrix. Sources are lines in the X matrix. On the sketch we figure a red object in the first source and two blue objects in the second source. Matrix A contains the mixing coefficients that allow various combinations of elements of X to produce Y.

Open with DEXTER

thumbnail Fig. 2

Illustration of the PCA colour selection. Left: HST image of the “Refsdal lens” in the galaxy cluster MACS J1149.6+2223. Middle: distribution of the first two PCA coefficients. The red and green dots correspond to the coefficients attributed to the first and second sources by MuSCADeT, respectively. Blue dots are rejected coefficients. Right: corresponding spatial distribution of colours as detected via PCA.

Open with DEXTER

2.2. Determination of the mixing matrix A

The A matrix is central to modelling of multi-band data as it describes the contribution of the different sources to the images taken at different wavelengths. In practice, the elements of A are the SEDs of the objects in the sources X. They can be assumed, for example as template spectra for objects of a given type, or they can be measured because spectroscopic data are available for at least some of the objects in the field of view. In most cases, however, the matrix A needs to be estimated solely from the multi-band imaging data, Y. In order to do this, we use a method based on a principal component analysis (PCA; Jolliffe 1986) of the data.

We consider the multi-band data as an ensemble of vectors { yi = 1..Nb [k] }, where the pixel values in band i at the spatial location k are stored in yi [k], as previously. In other words, { yi = 1..Nb [k] } is the measured SED at location k.

Following our definition of a source, two pixels belonging to a given source x have proportional SEDs. To build the mixing matrix A, one solution is to preselect obvious objects belonging to a same source and to average their SEDs, thus approximating the mixing coefficients corresponding to their source. A more subtle way to do this is to perform a principal component analysis of all SEDs belonging to bright objects and to look for proportional vectors. The details of this procedure can be summarised as follows:

  • We select the brightest objects in all bands and perform the PCA of the SEDs at pixel locations with high signal to noise (S/N). In practice, this is done by applying a wavelet filtering of all bands. Also, to save computation time, we rebin images to 64 × 64 pixels in size.

  • We perform a clustering analysis of the first two PCA components: the linearity and orthogonality of the PCA decomposition implies that proportional vectors see their respective PCA coefficients distributed along the same hyperplane in the PCA space. In other words, vectors with proportional SEDs have their first two components, PC1 and PC2, distributed along lines in the PC1-PC2 space, as illustrated in Fig. 2,

  • We identify the SEDs that have proportional PC1 and PC2 coefficients and average them. Coefficients that are judged too faint or too ambiguous (i.e. they could be a mix of both sources) are rejected,

  • We store the resulting mean SEDs as a column in the A matrix.

This algorithm shows good results in identifying objects with different dominant colours (Fig. 2), but its capabilities in terms of deblending are rather limited when distinct sources spatially overlap. The above PCA analysis is only a spectral analysis. The MCA method proposed in this paper combines the strengths of the morphological analysis and spectral analysis to design a reliable deblending algorithm.

3. Morphological component analysis

The morphological component analysis (MCA) method (Starck et al. 2004, 2010) allows us to separate several components in a given image based on their morphological diversity. Indeed, it was shown that it is possible to disentangle two (or more) signals mixed into one observable, this is only based on the fact that each of those signals can be sparsely represented in their respective data representation domains, called dictionaries, but not in the other’s. For instance, one could separate a periodic signal from a Gaussian profile in an image based on Fourier transform (associated with the periodic signal) and the wavelet transform (for the Gaussian profile). The projection of the mixing over the Fourier dictionary shows enhanced contribution from the periodic component whereas wavelet space shows higher coefficients for the Gaussian component.

3.1. Separation using a sparsity prior

MCA is based on the concept of sparse signal representation. A signal is sparse in a dictionary Φ when it can be well represented, such that x = Φα = ∑ iφiαi and only few coefficients α are different from zero. Some dictionaries, such as Fourier or wavelets, have implicit fast transformation and reconstruction operators that allows us to derive the coefficients α from x (and also to derive x from α) efficiently without having the elements of matrix Φ in memory. In inverse problems, a sparse solution is imposed by adding an 0-norm penalisation term to the data fidelity attachment term. MCA is an iterative algorithm, which separates a single image Y into J components xj, by solving (4)where Φx = α and σ is the noise standard deviation of the noise. Full details can be found in Starck et al. (2010). This method has been used to extract filamentary clouds in Herschel data (André et al. 2010) or, more recently, to improve SNIa detection in the SuperNova Legacy Survey data set (Möller et al. 2015).

3.2. Multi-band dictionaries

As we have multi-band data, we need to use morpho-spectral diversity. The dictionary Φi related to a given component xi is therefore a tensorial product of a spectral dictionary with the spatial dictionary Ψj, i.e. .

3.2.1. Spatial dictionary

In the case of strong lensing, the diversity between the components is mainly related to a different spectral morphology. Therefore we can reasonably use the same spatial dictionary, and in this strong lensing application we use the starlet dictionary (Starck et al. 2007). Starlet transform is an isotropic, undecimated wavelet transform that is computed using consecutive convolutions by a B-spline profile as a scaling function (Starck & Murtagh 2007). The resulting starlet representation is an overcomplete set of coefficients that represent variations in an image at different scales, and is particularly suited to represent astronomical images (Starck & Murtagh 2006).

3.2.2. Morpho-spectral dictionary

We note Ψ the starlet dictionary, we have xj = Ψαj where αj are the starlet coefficients of the jth source. A good choice for the morpho-spectral dictionary is to take Φj = ajΨ, where aj is the jth column of the matrix A. The data attachment term for multichannel data can be written as (5)The multichannel MCA hence consists in changing the data attachment term only, and we need to solve (6)

3.2.3. Lagrangian form and positivity

The sparse recovery problem can be formulated under an augmented Lagrangian form (7)and we can add a positivity constraint to the solution, so we need to solve (8)where λj accounts for the sparsity of each component xj in its own morpho-spectral dictionary. The next section describes how this equation can be solved.

4. MuSCADeT algorithm

The strong lens sparse deblending iterative algorithm (MuSCADeT) is an extension of the MCA algorithm, consisting in applying at each iteration three main steps:

  • 1.

    Perform a gradient step: U = X(n) + μAt(YAX(n)).

  • 2.

    Solve for each j: , and set to zero negative entries in .

  • 3.

    Decrease λj.

In this algorithm μ is the gradient step, derived from the 2-norm of matrix A (Higham 1992) such that μ = 2/ ||A|| 2.

This algorithm is also directly related to the proximal forward-backward algorithm (Combettes & Wajs 2005). A very nice aspect of this algorithm is that the minimisation involved in the second step do not require any iteration, and is obtained by , where ΔΨ,λj is the operator which performs the starlet transform, hard thresholds the starlet coefficients, and reconstructs an image from the thresholded coefficients. Here, λj is the threshold that allows us to select only coefficients that are significant enough to represent the signal. Thresholds are updated at each iteration as described in the following paragraphs. Full details can be found in Starck et al. (2010). Pseudo-algorithm 1 shows the principle of this iterative scheme. MuSCADeT is an iterative process that alternates between a gradient step (line 6) and a filtering of the components in transformed space through iterative hard thresholding (line 9).

4.1. Thresholding strategy

Thresholding aims at selecting the coefficients in transformed space that allow us to reconstruct the desired signal only. In this case, this means that for a given component, we want to select coefficients above noise level that accounts for this component and not for the others. It is therefore crucial to devise an adequate method to adapt thresholds at each iteration of algorithm 1.

Since each iteration moves our solution for components closer to a good separation, the thresholds have to be decreased to capture fainter and fainter structures. A classical way is to operate a linear decrease for instance, where values for λj are linearly sampled between an initial threshold chosen high above noise levels and a sensitivity value K. In general, K is chosen between three and five. The sensitivity value three allows for good completeness of detected coefficients and five ensures a selection that is free from noise-related coefficients. Noise levels are computed using median absolute deviation (Donoho & Johnstone 1994). Although linear or exponential laws are well suited to such problems (Starck et al. 2004), we choose here to rely on a more adaptive strategy based on minimum of maximums (MOM, see Bobin et al. 2007).

At each iteration, we simply estimate the maximum coefficient of each component in its own morpho-spectral dictionary and choose the smallest maxima plus a margin as a threshold. If the result is smaller than the threshold given by a linear decrease, the threshold is updated with the value estimated from MOM as illustrated in line 7 of the MuSCADeT algorithm (1). Full details on this thresholding scheme can be found in Starck et al. (2010).

thumbnail Fig. 3

Source separation with MuSCADeT in the case of a simple colour separation with SEDs estimated from PCA. From left to right: original simulated image of colour sources, first and second components (elements of X in Eq. (3)) as extracted by MuSCADeT, and the residual image after subtracting both components from the original image.

Open with DEXTER

thumbnail Fig. 4

Separation of blended sources with SEDs estimated from PCA (as in Fig. 3).

Open with DEXTER

5. Tests on simulations

In this section, we present several tests conducted on simulated multi-band images that emulate realistic problems. We show that MuSCADeT is able to separate highly blended sources and is robust to the approximation made that several astronomical sources (e.g. galaxies in a cluster) with similar but, not identical, SEDs can be considered a single source.

5.1. Simulations

In all following simulations we generate objects with Gaussian profiles and realistic SEDs. For the first two simulations, all SEDs were extracted from real HST observations of the galaxy cluster MACS J1149+2223 cluster (Kelly et al. 2015) using the PCA method described in Sect. 2.2. Each simulation comprises seven bands. We also add a white Gaussian noise with standard deviation σ = 0.01. An example of the impact of higher noise levels on our image separation is given in Appendix B.

We first apply MuSCADeT to data with four red objects with exactly the same SED and four blue objects with exactly the same SED. All objects have elliptical Gaussian profiles. Three of the simulated bands are used to produce the colour image in Fig. 3. Results of the separation after 100 iterations of MuSCADeT are shown in the other panels of Fig. 3. We observe no contamination between sources since no structure contaminates either component. Since the residuals present no structure, we also conclude that each component has been fully reconstructed.

The second simulation tests the deblending capacities of our algorithm in the ideal case where there is no SED variation between objects of a same component. We generate an extended elliptical Gaussian profile at the centre of the image affected with a red SED. Thirteen blue profiles are evenly spread between the centre and edges of the image such that profiles at the centre are completely blended. Results of the separation are shown in Fig. 4 along with a colour image of our simulated objects. No contamination between sources is visible either and residuals show no structure. We see in particular no difference between blue sources at the centre of the image (the most blended ones) and blue sources on the edges (less blended), meaning that the deblending of each profile is successful whether a profile is highly blended or not. We used 200 iterations of MuSCADeT to produce this result.

thumbnail Fig. 5

Left: part of an HST image of the galaxy cluster MACS J1149+2223, where the four objects used to extract SEDs are indicated with red contours. Middle: extracted SEDs. Each curve corresponds to the SEDs of the galaxies circled in red in the first panel. Right: simulated SEDs. The seven red SEDs are used to produce the upper row of galaxies in Fig. 6. The blue SEDs correspond to the lower row of galaxies (see text).

Open with DEXTER

thumbnail Fig. 6

Results of a separation of sources with poorly known SEDs (as in Fig. 3).

Open with DEXTER

The last simulation tests the robustness of our algorithm to SED variation across objects from a same component. To account for realistic SED variation, we extracted the SEDs from four red galaxies in the MACS J1149+2223 cluster (see Fig. 5), which appear to have the same colour by integrating the flux in each galaxy profile. The resulting SEDs can be seen in the middle panel of Fig. 5. We recorded the slopes of the SEDs and generated a set of eight slopes linearly spread between the maximum and minimum slope estimated. These slopes are then applied to SEDs extracted from cluster MACS J1149+2223 via PCA (see Fig. 5). This way, we have two sets of SEDs that account for red and blue sources and that mimic a range of variations as observed in real images. Sixteen Gaussian profiles (eight red, eight blue) are then generated and each of them is associated with one of the previously generated SEDs. The left panel of Fig. 6 shows three bands of the simulated images as RGB images. Figure 6 shows a colour image of our simulated images and the result of a separation by MuSCADeT. We see again that no structure appears in the residuals and no contamination is found between components. However, the great similarity between SEDs from different components forced us to increase the number of iterations of MuSCADeT to 5000 in order to obtain such results, thus increasing computation time. In each, we used the full algorithm described in this paper including an automated estimation of SEDs through PCA.

thumbnail Fig. 7

Application of MuSCADeT to the MACSJ1149+2223 cluster. Top left: colour image generated using the F475W, F606W, and F814W bands. Our goal is to separate the foreground cluster galaxies (red) from the background lensed galaxy (blue). Middle: the two sources extracted with 2000 iterations of MuSCADeT. Bottom left: original image minus the blue source found by MuSCADeT. Bottom right: same as the bottom left, with the red source subtracted. Top right: residual image obtained after subtraction of the two estimated sources from the original data.

Open with DEXTER

6. Application to real data

6.1. Lens-source separation on MACS J1149+2223

We now apply MuSCADeT to real multi-band data and we compare the performances to traditional model fitting.

We use the deep HST data set of the galaxy cluster MACS J1149+2223 (“Refsdal SN”; Kelly et al. 2015) to carry out this experiment, we show bright cluster members producing strongly lensed images of a distant spiral galaxies with clumpy structures. Our goal is to separate the data into two sources containing the foreground lens galaxies and background lensed object(s). As the cluster contains many member galaxies and as the background galaxy has complex structure, the deblending task is challenging, making this data set a good test for our method.

MACS J1149+2223 has been observed with the ACS in seven bands: F435W, F475W, F555W, F606W, F625W, F775W and F814W (proposal ID: 12068, principal investigator: M. Postman), providing a good spectral coverage for MuSCADeT to work. The data are publicly available from the STScI website1 and drizzled so that the combined frames in each band have the same orientation and pixel scale.

We estimate both SEDs in the mixing matrix using our PCA technique. The result after 2000 iterations of MuSCADeT is shown in Fig. 7, where the two separated sources are shown in the middle panels. In the two lower panels of the figure, we also show the result of the subtraction of source 1 and source 2 from the original data. The colour scale in these images is the same as in the original data. The overall residual image, i.e. with both sources subtracted from the data, is shown with ± 5σ cut levels.

The source separation works very well, with the exception of a few objects with "green" SEDs. One such object is visible in the centre of the image, resulting in a signal in both red and blue sources. This is an intrinsic limitation to our algorithm that separates objects using a limited number of sources, each one with its own SED. Although the SEDs do not need to be known perfectly (as shown in our tests with simulated images), objects with SEDs falling “in between” the SEDs allocated to each source, may lead to inaccurate separation. A possible mitigation strategy is to add extra sources to the decomposition with for example a blue, a green, and a red SED. Deciding whether to do this or not depends on the exact scientific application.

In our example, we do not take the PSF convolution into account. A variation of the PSF with wavelength can introduce artefacts in the separation, especially for objects whose angular size is comparable with the size of the PSF. The central parts of galaxies show such structures and, indeed, in the present data small structures are seen in the residual image at the location of the foreground galaxies. Introducing the PSF convolution can be done in principle, and it might be needed at least for some applications, but at a cost of increased computation time and complexity in the minimisation process.

6.1.1. Comparison with profile fitting

thumbnail Fig. 8

Comparison between the MuSCADeT and that galfit separations. The original colour image of MACS J1149+2223 is shown on the left, followed by the galfit subtraction of the galaxy members and the subtraction using MuSCADeT. The colour cuts are exactly the same in all three images. The red ellipse in the middle panel indicates an area where the galfit fit leads to under-subtraction of a galaxy halo.

Open with DEXTER

A popular way to carry out source separation in galaxy clusters is to fit two-dimensional profiles to the cluster members and to subtract from the data. We apply such a procedure to the HST data of the MACS J1149+2223 cluster and compare with MuSCADeT. In this exercise, we identified nine red objects with elliptical profiles that we fit using the galfit (Peng et al. 2002) software. Each red object is fitted with either one or two Sersic profiles depending on its morphology. Fitting is performed in each of the three bands used to build the colour image (F475W, F606W, and F814W) using the same initial conditions. Since no PSF was used in applying MuSCADeT to this cluster, we used a Dirac function as an input PSF for galfit. The result is presented in the middle panel of Fig. 8 along with the MuSCADeT result. We note that:

  • Profile fitting leaves significant artefacts in the central parts of the fitted galaxies even when a double Sersic model is used.

  • Profile fitting does not model the extended halos of galaxies well, as shown in the circled region in middle panel of Fig. 8. This is both because extended halos are not simple analytical profiles and because the many background lensed structures influence the fit. In principle, these structures can be masked but 1- designing the mask can be difficult and is time consuming, 2- in some parts of the image there is no clean area with no contamination by the background object. The mask would take the majority of the data.

  • The top left lensed image appears to be more blue after galfit’s run than with MuSCADeT. This could mean that MuSCADeT performs badly at extracting extended halo. However, it would imply that too much signal from the red source has been attributed to the blue source. Therefore, when subtracting the blue source from the original images one would see holes in the extended profiles of the red galaxies that are not observed here (see lower left panel of Fig. 7). This implies that galfit is overfitting the extend halo to compensate for blue structures, as pointed out in the previous note, thus removing part of the flux from blue sources.

  • The human time involved in profile fitting can be a limiting factor for large data sets. The user has to decide where to put a galaxy and to find a strategy to estimate the initial guesses for the many parameters involved in the fit. MuSCADeT is fully automated procedure with only one parameter to be chosen, i.e. the sensitivity value K involved in the threshoding scheme.

6.2. Bulge-disk separation on spiral galaxy in Abell 2744

Another possible application of our algorithm is the separation of coloured disk and bulge components in a spiral galaxy. To illustrate this, we use a spiral galaxy in the galaxy cluster Abell 2744, which was imaged with the HST as part of the Hubble Frontier Fields programme.

Abell 2744 was observed with the ACS in three bands: F435W, F606W, and F814W, resulting in the colour image in Fig. 9 (Proposal ID: 11689, principal investigator: R. Dupke). The data are publicly available from the STScI website2 and drizzled so that the combined frames in each band have the same orientation and pixel scale.

Spiral galaxies represent an obvious test bench for our algorithm as they are composed of a red bulge of older stars and of a blue disk dominated by young stars. The ability to separate both colour components in an unbiased way allows us to trace stellar populations in galaxies and/or highlight morphological relations between bulge and disk shapes. Numerous examples can be found in the literature of studies where bulge and disk components are separated morphologically, using profile-fitting techniques featuring Sersic, exponential, or DeVaucouleurs functions in either single or multiple bands (e.g. Vika et al. 2014; Maltby et al. 2012; Pastrav et al. 2013; Kendall et al. 2011).

With MusCADeT, we are able to capture the morphology of the young and old stellar populations of the spiral galaxy in Abell 2744 under the single assumption that they do not display the same colour. Our decomposition in Fig. 9 shows an elliptical red bulge with extended and smooth features along the spiral arms, which elliptical profile-fitting methods would fail to model. The separation of the blended blue and red components in Fig. 9 is overall excellent except for slight cross-contamination of the two colour channels due to features smaller than the PSF size. Adding the PSF in our separation technique is under development. However, even without this refinement, the decomposition presented in our example would be impossible to achieve with profile-fitting techniques given the morphological complexity of the galaxy.

thumbnail Fig. 9

Separation of the red and blue stellar populations in a spiral galaxy. From left to right we show the original HST colour image of a spiral galaxy in the galaxy cluster Abell 2744, the same colour image after subtraction of the blue component estimated with MuSCADeT, the colour image after subtraction of the red component estimated with MuSCADeT, and the residual image after subtraction of both components. The three first images have the same colour cuts. The residuals are shown with cut levels set to five times the noise level in each colour channel.

Open with DEXTER

7. Conclusion

We have developed a new model-independent tool to deblend the images of astronomical objects based on their colour contrast. This is made possible via modern techniques for image analysis based on sparsity of objects in an appropriate dictionary.

More specifically, we created a morpho-spectral dictionary combining Starlets and SEDs to separate images of galaxies with different colours. We show from simulated data that our algorithm is robust against very strong blending as well as against SED variations among objects belonging to the same colour component (source) of our model. SED variations across the field result in an increase of the computing time by one order of magnitude with respect to the ideal case where all objects in a given source have the same SED. This does not hamper a reliable source separation, however.

The method is successfully applied to the deep HST images of the galaxy cluster MACS J1149+2223, in which we separate the foreground red galaxies from background blue clumpy lensed galaxies, and to Abell 2744, where we separate the red and blue stellar populations of a spiral galaxy. This is done in an automated way and with better efficiency that with standard profile fitting. All codes used to produce the results presented here are made freely available for the sake of reproducible research and to make our tool usable for the community.

Future developments of our method include accounting for the PSF in each band and including explicit SED variations across the field of view. These SED variations multiply the complexity of the problem by Np × Ns × Nb, but the effect of the increased complexity can likely be minimised by using sparse priors and physical constrains on the SED profiles. Also, the increased computation time resulting from the extra complexity should be partially compensated by a reduced number of iterations.

The deblending method described here was devised specifically to address the problem of object deblending in the case of strong lensing systems. Deblending is essential in this case to see small and faint details in the lensed structures, which are free of contamination by the bright foreground. The range of applications of the method is nevertheless much broader. Among the numerous possible applications are the identification of star-forming regions in galaxies, model-independent bulge-disk decompositions of galaxies, or even the improvement of photometric redshift in large sky surveys where blending can be a serious issue given the depth of the data.


Acknowledgments

This work is supported by the Swiss National Science Foundation (SNSF), and by the European Community through grant PHySIS (contract no. 640174) and DEDALE (contract No. 665044) within the H2020 Framework Program. The HST images used in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. Support for MAST for non-HST data is provided by the NASA Office of Space Science via grant NNX09AF08G and by other grants and contracts. We would like to thank the members of the CosmoStat team for useful discussions on inverse problem solving using sparse prior.

References

Appendix A: Reproducible research

In the spirit of carrying out reproducible research, we make pubic all codes and resulting products describes in this paper. Table A.1 lists all products that will be made available along with this paper. The MuSCADeT itself is made available as well as input files and routines needed for all benchmark tests and for the application to real data. The routines provide simple examples of how to execute the MuSCADeT algorithm. We encourage potential users to modify them as they wish for their own scientific applications.

Table A.1

Products made available in this paper in the spirit of reproducible research.

Appendix B: Deblending noisy data with known SED

We show here that the MuSCADeT algorithm is capable of separating blended sources even in the case of high noise levels. We

generate simulations of blended objects as in Fig. 4, but with a noise level ten times larger. The high noise levels along with the strong blending make it hard for the PCA estimator to estimate a good SED for the separation.

Although the PCA technique fails in such conditions, the main feature of our algorithm, which is the morphological component analysis-based inversion, still manages to estimate good sources. The price to pay in this case is that the SEDs must be known. Figure B.1 shows the result of a separation performed by MuSCADeT on very noisy data using known SEDs, showing that our algorithm is still able to separate sources. Our PCA SED estimator might replaced in the near future to cope with noisy data. For the present example, we decided to show the residuals after separation of both colour components (middle panels of Fig. B.1) to show that we estimate the sources down to the noise level.

thumbnail Fig. B.1

Separation of blended sources at low S/N. From left to right are shown the original simulated images, the original image after subtraction of the blue component as estimated from MuSCADeT, the original image after subtraction of the red component and the residual image after subtraction of both components.

Open with DEXTER

All Tables

Table A.1

Products made available in this paper in the spirit of reproducible research.

All Figures

thumbnail Fig. 1

Illustration of the blind source separation in the case of two sources. To make the figure simple, the images in each band are represented as lines in the Y matrix. Sources are lines in the X matrix. On the sketch we figure a red object in the first source and two blue objects in the second source. Matrix A contains the mixing coefficients that allow various combinations of elements of X to produce Y.

Open with DEXTER
In the text
thumbnail Fig. 2

Illustration of the PCA colour selection. Left: HST image of the “Refsdal lens” in the galaxy cluster MACS J1149.6+2223. Middle: distribution of the first two PCA coefficients. The red and green dots correspond to the coefficients attributed to the first and second sources by MuSCADeT, respectively. Blue dots are rejected coefficients. Right: corresponding spatial distribution of colours as detected via PCA.

Open with DEXTER
In the text
thumbnail Fig. 3

Source separation with MuSCADeT in the case of a simple colour separation with SEDs estimated from PCA. From left to right: original simulated image of colour sources, first and second components (elements of X in Eq. (3)) as extracted by MuSCADeT, and the residual image after subtracting both components from the original image.

Open with DEXTER
In the text
thumbnail Fig. 4

Separation of blended sources with SEDs estimated from PCA (as in Fig. 3).

Open with DEXTER
In the text
thumbnail Fig. 5

Left: part of an HST image of the galaxy cluster MACS J1149+2223, where the four objects used to extract SEDs are indicated with red contours. Middle: extracted SEDs. Each curve corresponds to the SEDs of the galaxies circled in red in the first panel. Right: simulated SEDs. The seven red SEDs are used to produce the upper row of galaxies in Fig. 6. The blue SEDs correspond to the lower row of galaxies (see text).

Open with DEXTER
In the text
thumbnail Fig. 6

Results of a separation of sources with poorly known SEDs (as in Fig. 3).

Open with DEXTER
In the text
thumbnail Fig. 7

Application of MuSCADeT to the MACSJ1149+2223 cluster. Top left: colour image generated using the F475W, F606W, and F814W bands. Our goal is to separate the foreground cluster galaxies (red) from the background lensed galaxy (blue). Middle: the two sources extracted with 2000 iterations of MuSCADeT. Bottom left: original image minus the blue source found by MuSCADeT. Bottom right: same as the bottom left, with the red source subtracted. Top right: residual image obtained after subtraction of the two estimated sources from the original data.

Open with DEXTER
In the text
thumbnail Fig. 8

Comparison between the MuSCADeT and that galfit separations. The original colour image of MACS J1149+2223 is shown on the left, followed by the galfit subtraction of the galaxy members and the subtraction using MuSCADeT. The colour cuts are exactly the same in all three images. The red ellipse in the middle panel indicates an area where the galfit fit leads to under-subtraction of a galaxy halo.

Open with DEXTER
In the text
thumbnail Fig. 9

Separation of the red and blue stellar populations in a spiral galaxy. From left to right we show the original HST colour image of a spiral galaxy in the galaxy cluster Abell 2744, the same colour image after subtraction of the blue component estimated with MuSCADeT, the colour image after subtraction of the red component estimated with MuSCADeT, and the residual image after subtraction of both components. The three first images have the same colour cuts. The residuals are shown with cut levels set to five times the noise level in each colour channel.

Open with DEXTER
In the text
thumbnail Fig. B.1

Separation of blended sources at low S/N. From left to right are shown the original simulated images, the original image after subtraction of the blue component as estimated from MuSCADeT, the original image after subtraction of the red component and the residual image after subtraction of both components.

Open with DEXTER
In the text

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

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

Initial download of the metrics may take a while.