5 Morphological component analysis

5.1 Introduction

The content of an image is often complex, and there is not a single transform which is optimal to represent all the contained features. For example, the Fourier transform better represents some textures, while the wavelet transform better represents singularities. Even if we limit our class of transforms to the wavelet one, decision have to be taken between an isotropic wavelet transform which produce good results for isotropic objects (such stars and galaxies in astronomical images, cells in biological images, etc.), or an orthogonal wavelet transform, which is better for images with edges. This has motivated the development of different methods (Chen et al. 1998; Meyer et al. 1998; Huo 1999), and the two most frequently discussed approaches are the Matching Pursuit (MP) (Mallat & Zhang 1993) and the Basis pursuit (BP) (Chen et al. 1998). A dictionary ${\cal D}$ being defined as a collection of waveforms $(\varphi_{\gamma})_{\gamma \in \Gamma}$, the general principe consists in representing a signal s as a "sparse'' linear combination of a small number of basis such that:

$$% s = \sum_{\gamma} a_{\gamma} \varphi_{\gamma}$$ (11)

or an approximate decomposition

$$% s = \sum_{i=1}^m a_{\gamma_i} \varphi_{\gamma_i} + R^{(m)} .$$ (12)

Matching pursuit (Mallat & Zhang 1993; Mallat 1998) method (MP) uses a greedy algorithm which adaptively refines the signal approximation with an iterative procedure:

• Set s⁰ = 0 and R⁰ = 0.
• Find the element $\alpha_k \varphi_{\gamma_k}$ which best correlates with the residual.
• Update s and R:
  

  $$s^{k+1} = s^k + \alpha_k \varphi_{\gamma_k}$$

  R^k+1 = s - s^k . (13)

  
  

In case of non orthogonal dictionaries, it has been shown (Chen et al. 1998) that MP may spend most of the time correcting mistakes made in the first few terms, and therefore is suboptimal in term of sparsity.

Basis pursuit method (Chen et al. 1998) (BP) is a global procedure which synthesizes an approximation $\tilde{s}$ to s by minimizing a functional of the type

$$% \Vert s - \tilde{s}\Vert _{\ell_2} ^2 + \lambda \cdot \Vert\alpha\Vert _{\ell_1}, \quad \tilde{s} = \Phi \alpha.$$ (14)

Between all possible solutions, the chosen one has the minimum l¹ norm. This choice of l₁ norm is very important. A l₂ norm, as used in the method of frames (Daubechies 1988), does not preserve the sparsity (Chen et al. 1998).

In many cases, BP or MP synthesis algorithms are computationally very expensive. We present in the following an alternative approach, that we call Combined Transforms Method (CTM), which combines the different available transforms in order to benefit of the advantages of each of them.

5.2 The combined transformation

Depending on the content of the data, several transforms can be combined in order to get an optimal representation of all features contained in our data set. In addition to the ridgelet and the curvelet transform, we may want to use the à trous algorithm which is very well suited to astronomical data, or the undecimated wavelet transform which is commonly used in the signal processing domain.

Other transform such wavelet packets, the Fourier transform, the Pyramidal median transform (Starck et al. 1998), or other multiscale morphological transforms, could also be considered. However, we found that in practice, these four transforms (i.e. curvelet, ridgelet, à trous algorithm, and undecimated wavelet transform) furnishes a very large panel of waveforms which is generally large enough to well represents all features contained in the data.

In general, suppose that we are given K linear transforms $T_1, \ldots, T_K$ and let $\alpha_k$ be the coefficient sequence of an object x after applying the transform T_k, i.e. $\alpha_k = T_k x$. We will suppose that for each transform T_k we have available a reconstruction rule that we will denote by T^-1_k although this is clearly an abuse of notations.

Therefore, we search a vector ${\bf\alpha} = {\alpha_1, \dots, \alpha_{K}}$ such that

$$% s = \Phi {\bf\alpha}$$ (15)

where $\Phi {\bf\alpha} = \sum_{k=1}^K T_k^{-1} \alpha_k$. As our dictionary is overcomplete, there is an infinity of vectors verifing this condition, and we need to solve the following optimization problem:

$$% {\rm min} \parallel s - \phi {\bf\alpha} \parallel^2 + {\cal C}(\bf\alpha)$$ (16)

where ${\cal C}$ is a penalty term. We easily see that chosing ${\cal C}({\bf\alpha})$ = $\parallel {\bf\alpha} \parallel_{l_1}$ leads to the BP method, where the dictionary ${\cal D}$ is only composed of the basis elements of the chosen transforms.

Two iterative methods, soft-CTM and hard-CTM, allowing us to realize such a combined transform, are described in this section.

Figure 10: Top, original image containing lines and gaussians. Botton left, reconstructed image for the à trous wavelet coefficient, bottom right, reconstructed image from the ridgelet coefficients.

   
5.3 Soft-CTM

Noting T₁, ..., T_K the K transform operators, a solution ${\bf\alpha}$ is obtained by minimizing a functional of the form:

 

$$% J({\bf\alpha}) = \parallel s - \sum_{k=1}^{K} T_k^{-1} \alpha_k \parallel_2^2 + \lambda \sum_k \parallel \alpha_k \parallel_1$$ (17)

where s is the original signal, and $\alpha_k$ are the coefficients obtained with the transform T_k.

Figure 11: Upper left, galaxy SBS 0335-052 (10 $\mu$m), upper middle, upper right, and bottom left, reconstruction respectively from the ridgelet, the curvelet and wavelet coefficients. Bottom middle, residual image. Bottom right, artifact free image.

Figure 12: Upper left, galaxy SBS 0335-052 (20 $\mu$m), upper right, addition of the reconstructed images from both the ridgelet and the curvelet coefficients, bottom left, reconstruction from the wavelet coefficients, and bottom right, residual image.

An simple algorithm to achieve such an solution is:

1.

Initialize $L_{\max}$, the number of iterations N_i, $\lambda = L_{\max}$, and  $\delta_{\lambda} = \frac{L_{\max}}{N_i}$.

2.

While $\lambda >= 0$ do

3.

For k = 1, ..., K do

• Calculate the residual $R = s - \sum_k T_k^{-1} \alpha_k$.
• Calculate the transform T_k of the residual: r_k = T_k R.
• Add the residual to $\alpha_k$: $\alpha_k = \alpha_k + r_k$.
• Soft threshold the coefficient $\alpha_k$ with the $\lambda$ threshold.

4.

$\lambda = \lambda - \delta$, and goto 2.

Figure 10 illustrates the result in the case where the input image contains only lines and Gaussians. In this experiment, we have initialized $L_{\max}$ to 20, and $\delta$ to 2 (10 iterations). Two transform operators were used, the à trous wavelet transform and the ridgelet transform. The first is well adapted to the detection of Gaussian due to the isotropy of the wavelet function (Starck et al. 1998), while the second is optimal to represent lines (Candès & Donoho 1999b). Figure 10 top, bottom left, and bottom right represents respectively the original image, the reconstructed image from the à trous wavelet coefficient, and the reconstructed image from the ridgelet coefficient. The addition of both reconstructed images reproduces the original one.

In some specific cases where the data are sparse in all bases, it has been shown (Huo 1999; Donoho & Huo 2001) that the solution is identical to the solution when using a $\parallel . \parallel_0$ penalty term. This is however generally not the case. The problem we met in image restoration applications, when minimizing Eq. (17), is that both the signal and noise are split into the bases. The way the noise is distributed in the coefficients $\alpha_k$ is not known, and leads to the problem that we do not know at which level we should threshold the coefficients. Using the threshold we would have used with a single transform makes a strong over-filtering of the data. Using the l¹ optimization for data restoration implies to first study how the noise is distributed in the coefficients. The hard-CTM method does not present this drawback.

Figure 13: Top left, HST image of A370, top right coadded image from the reconstructions from the ridgelet and the curvelet coefficients, bottom left reconstruction from the à trous wavelet coefficients, and bottom right addition of the three reconstructed images.

   
5.4 Hard-CTM

The following algorithm consists in hard thresholding the residual successively on the different bases.

1.

For noise filtering, estimate the noise standard deviation $\sigma$, and set $L_{\min} = k_{\sigma}$. Otherwise, set $\sigma=1$ and $L_{\min} = 0$.

2.

Initialize $L_{\max}$, the number of iterations N_i, $\lambda = L_{\max}$ and  $\delta_{\lambda} = \frac{L_{\max} - L_{\min}}{N_i}$.

3.

Set all coefficients $\alpha_k$ to 0.

4.

While $\lambda >= L_{\min}$ do

5.

for k = 1, .., K do

• Calculate the residual $R = s - \sum_k T_k^{-1} \alpha_k$.
• Calculate the transform T_k of the residual: r_k = T_k R.
• For all coefficients $\alpha_{k,i}$ do
  • Update the coefficients: if $\alpha_{k,i} \ne 0$ or $\mid r_{k,i} \mid > \lambda \sigma$ then $\alpha_{k,i} = \alpha_{k,i} + r_{k,i}$.

6.

$\lambda = \lambda - \delta_{\lambda}$, and goto 5.

For an exact representation of the data, $k_{\sigma}$ must be set to 0. Choosing $k_{\sigma} > 0$ introduces a filtering. If a single transform is used, it corresponds to the standard k-sigma hard thresholding.

It seems that starting with a high enough $L_{\max}$ and a high number of iterations would lead to the l⁰ optimization solution, but this remains to be proved.

Figure 14: Top left, galaxy NGC 2997, top right reconstructed image from the à trous wavelet coefficients, bottom left, reconstruction from the ridgelet coefficients, and bottom right addition of both reconstructed images.

5.5 Experiments

5.5.1 Experiment 1: Infrared Gemini data

Figure 11 upper left shows a compact blue galaxy located at 53 Mpc. The data have been obtained on ground with the GEMINI-OSCIR instrument at 10 $\mu$m. The pixel field of view is $0.089^{\prime\prime}$/pix, and the source was observed during 1500s. The data are contaminated by a noise and a stripping artifact due to the instrument electronic. The same kind of artifact pattern were observed with the ISOCAM instrument (Starck et al. 1999).

This image, noted D₁₀, has been decomposed using wavelets, ridgelets, and curvelets. Figure 11 upper middle, upper right, and bottom left show the three images R₁₀, C₁₀, W₁₀reconstructed respectively from the ridgelets, the curvelets, and the wavelets. Image in Fig. 11 bottom middle shows the residual, i.e. e₁₀ = D₁₀ - (R₁₀ + C₁₀ + W₁₀). Another interesting image is the artifact free one, obtained by subtracting R₁₀ and C₁₀ from the input data (see Fig. 11 bottom right). The galaxy has well been detected in the wavelet space, while all stripping artifact have been capted by the ridgelets and curvelets.

Figure 12 upper left shows the same galaxy, but at 20 $\mu$m. We have applied the same decomposition on D₂₀. Figure 12 upper right shows the coadded image R₂₀ + C₂₀, and we can see bottom left and right the wavelet reconstruction W₂₀ and the residudal e₂₀ = D₂₀ - (R₂₀ + C₂₀ + W₂₀).

5.5.2 Experiment 2: A370

Figure 13 upper left shows the HST A370 image. It contains many anisotropic features such the gravitationnal arc, and the arclets. The image has been decomposed using three transforms: the ridgelet transform, the curvelet transform, and the à trous wavelet transform. Three images have then been reconstructed from the coefficients of the three basis. Figure 13 upper right shows the coaddition of the ridgelet and curvelet reconstructed images. The à trous reconstructed image is displayed in Fig. 13 lower left, and the coaddition of the three images can be seen in Fig. 13 lower right. The gravitational arc and the arclets are all represented in the ridgelet and the curvelet basis, while all isotropic features are better represented in the wavelet basis.

5.5.3 Experiment 3: Elongated-point like object separation in astronomical images.

Figure 14 shows the result of a decomposition of a spiral galaxy (NGC 2997). This image (Fig. 14 top left) contains many compact structures (stars and HII region), more or less isotropic, and large scale elongated features (NGC 2997 spiral part). Compact objects are well represented by isotropic wavelets, and the elongated features are better represented by a ridgelet basis. In order to benefit of the optimal data representation of both transforms, the image has been decomposed on both the "à trous'' wavelet transform and on the ridgelet transform by using the same method as described in Sect. 5.4. When the functional is minimized, we get two images, and their coaddition is the filtered version of the original image. The reconstructions from the à trous coefficients, and from the ridgelet coefficients can be seen in Fig. 14 top right and bottom left. The addition of both images is presented in Fig. 14 bottom right.

We can see that this Morphological Component Analysis (MGA) allows us to separate automatically features in an image which have different morphological aspects. It is very different from other techniques such as Principal Component Analysis or Independent Component Analysis (Cardoso 1998) where the separation is performed via statistical properties.

Acknowledgements

We are grateful to the referee for helpful comments on an earlier version.