3 Detection and deconvolution

3.1 Introduction

The PSF is not needed with MVM. This is an advantage when the PSF is unknown, or difficult to estimate, which happens relatively often when it is space-variant. However, when the PSF is well-determined, it becomes a drawback because known information is not used for the object reconstruction. This can lead to systematic errors in the photometry, which depends on the PSF and on the source signal to noise ratio. In order to preempt such a bias, a kind of calibration must be performed using simulations (Starck et al. 1999). This section shows how the PSF can be used in the MVM, leading to a deconvolution.

3.2 The deconvolution problem

Consider an image characterized by its intensity distribution (the "data'') I(x,y), corresponding to the observation of a "real image'' O(x,y) through an optical system. If the imaging system is linear and shift-invariant, the relation between the object and the image in the same coordinate frame is a convolution:

 

I(x,y) = (O * P)(x,y) + N(x,y) (3)

P(x,y) is the PSF of the imaging system, and N(x,y)is additive noise. In practice O * P is subject to non-stationary noise. We want to determine O(x,y) knowing I(x,y) and P(x,y). This inverse problem has led to a large amount of work, the main difficulties being the existence of: (i) a cut-off frequency of the PSF, and (ii) the noise. Equations (3) is always an ill-posed problem. This means that there is no unique least-squares solution of minimal norm $\parallel I(x,y) - P(x,y) * O(x,y) \parallel$, and some constraints must be added in order to regularize the problem (Gonzalez 1993; Pratt 1991; Starck 1998). Once the deconvolved image is obtained, it is generally difficult to know what the noise level is in the deconvolved image, and hence it becomes impossible to detect the objects with a confidence interval. For this reason, astronomers generally prefer to apply the detection on non-deconvolved images. Another argument against deconvolution is that is has been shown that some regularization methods like MEM (Narayan & Nityananda 1986) affect the photometry.

3.3 Object reconstruction using the PSF

A reconstructed and deconvolved object can be obtained by searching for a signal O such that the wavelet coefficients of P*O are the same as those of the detected structures. If $\cal T$ describes the wavelet transform operator, and $P_{\rm w}$ the projection operator in the subspace of the detected coefficients, the solution is found by minimization of

 

$$% J(O) = \parallel W - (P_{\rm w} \circ {\cal T}) P * O \parallel$$ (4)

where W represents the detected wavelet coefficients of the data, and P is the PSF. In this approach, each object is deconvolved separately. The flux related to the extent of the PSF will be taken into account. For point sources, the solution will be close to that obtained by PSF fitting. This problem is also different from global deconvolution in the sense that it is well constrained. Except for the positivity of the solution which is always true and must be used, no other constraint is needed. This is due to the fact that the reconstruction is performed from a small set of wavelet coefficients (those above a detection limit). The number of objects are the same as those obtained by the MVM, but the photometry and the morphology is different. The astrometry may also be affected.

3.4 The algorithm

Any minimizing method can be used to obtain the solution O. Since we did not find any problem of convergence, noise amplification, or ringing effect, we chose the van Cittert method on the grounds of its simplicity. For each detected object, we apply the following algorithm:

$$% O^{n+1} = O^{n} + {\cal T}^{-1}(W - (P_{\rm w} \circ {\cal T}) P * O^{n})$$ (5)

where ${\cal T}^{-1}$ is the inverse wavelet transform.

1.

Set n to 0.

2.

Find the initial estimation Oⁿ by applying an inverse wavelet transform to the set W corresponding to the detected wavelet coefficients in the data.

3.

Convolve Oⁿ with the PSF P: Iⁿ = P*Oⁿ.

4.

Determine the wavelet transform W(Iⁿ) of Iⁿ.

5.

Threshold all wavelet coefficients in W(Iⁿ) at position and scales where nothing has been detected (i.e. $P_{\rm w}$ operator). We get $W_{\rm t}(I^n)$.

6.

Determine the residual $W(R) = W - W_{\rm t}(I^n)$.

7.

Reconstruct the residual image Rⁿ by applying an inverse wavelet transform.

8.

Add the residual to the solution: Oⁿ⁺¹ = Oⁿ + Rⁿ.

9.

Threshold negative values in Oⁿ⁺¹.

10.

If $\sigma(R^n) / \sigma(O^0) < \epsilon$ then n = n + 1 and go to step 3.

11.

Oⁿ⁺¹ contains the deconvolved reconstructed object.

In practice, convergence is very fast (less than 20 iterations). The reconstructed image (not deconvolved) can also be obtained just by reconvolving the solution with the PSF.

3.5 Space-variant PSF

Deconvolution methods generally do not take into account the case of space-variant PSF. The standard approach when the PSF varies is to decompose the image into blocks, and to consider the PSF constant inside a given block. Blocks which are too small lead to a problem of computation time (the FFT cannot be used), while blocks which are too large introduce errors due to the use of an incorrect PSF. Blocking artifacts may also appear. Combining source detection and deconvolution opens up an elegant way for deconvolution with a space-variant PSF. Indeed, a straightforward method is derived by just replacing the constant PSF at step 3 of the algorithm with the PSF at the centre of the object. This means that it is not the image which is deconvolved, but its constituent objects.

3.6 Deconvolution and resolution

3.6.1 The intrinsic correlation function

In many cases, there is no sense in trying to deconvolve an image at the resolution of the pixel (especially when the PSF is very large). The idea to limit the resolution is relatively old, because it is already this concept which is used in the CLEAN algorithm (Högbom 1974). Indeed the Clean-Beam fixes the resolution in the final solution. This principle was also developed by Lannes (1987) in a different form. This concept has been re-invented, first by Gull & Skilling (1991) who have called the Clean-Beam the Intrinsic Correlation Function (ICF), and more recently by Magain (1998) and Pijpers (1999).

The ICF is usually a Gaussian, but in some cases it may be useful to take another function. For example, if we want to compare two images I₁ and I₂ which are obtained with two wavelengths or with two different instruments, their PSFs P₁ and P₂ will certainly be different. The classic approach would be to deconvolve I₁ with P₂ and I₂ with P₁, so we are sure that both are at the same resolution. But unfortunately we lose some resolution in doing this. Deconvolving both images is generally not possible because we can never be sure that both solutions O₁ and O₂ will have the same resolution.

A solution would be to deconvolve only the image which has the worse resolution (say I₁), and to limit the deconvolution to the second image resolution (I₂). Then, we just have to take P₂ for the ICF. The deconvolution problem is to find $\tilde O$ (hidden solution) such that:

$$% I_1 = P_1 * P_2 * \tilde O$$ (6)

and our real solution O₁ at the same resolution as I₂ is obtained by convolving $\tilde O$ with P₂. O₁ and I₂ can then be compared.

Introducing an ICF G in the deconvolution equation leads to just considering a new PSF $P^\prime$ which is the convolution of P and G. The deconvolution is carried out using $P^\prime$, and the solution must be reconvolved with G at the end. In this way, the solution has a constrained resolution, but aliasing may occur during the iterative process, and it is not sure that the artifacts will disappear afer the re-convolution with G. Magain (1998) has proposed an original alternative to this problem, by assuming that the PSF can be considered as the convolution product of two terms, the ICF G and an unknown S, P=G*S. Using S instead of P in the deconvolution process, and a sufficiently large FWHM value for G, implies that the Shannon sampling theorem (Shannon 1948) is never violated. But the problem is now to calculate S, knowing P and G, which is again a deconvolution problem. Unfortunately, this delicate point was not discussed in the original paper. Propagation of the error on the S estimation in the final solution has also until now not been investigated, even if this issue seems to be quite important.

3.6.2 ICF calculation

This section describes how to calculate the FWHMfor a given sampling, in order not to violate the Shannon sampling theorem. Gaussian functions are generally chosen for the ICF. The resolution to be achieved is fixed by its standard deviation $\sigma_{\rm G}$, or its FWHM(FWHM=2.34 $\sigma_{\rm G}$). Since the Fourier transform of a Gaussian of standard deviation $\sigma_{\rm G}$ is also a Gaussian of standard deviation $\sigma_\nu = {N \over 2\pi \sigma_{\rm G}}$, (N being the number of pixels), we can estimate the smallest FWHM which does not violate the Shannon sampling theorem. In theory, a Gaussian cannot respect it, but in practice we can consider that values smaller than a given $\epsilon$ have no practical effect, and the Shannon sampling theorem is experimentally respected if

$$% \exp{ -{u^2 \over 2 \sigma_{\nu}^2}} < \epsilon \mbox{ when } u > {N \over 2}.$$ (7)

For $u = {N \over 2}$, we have: $\exp{ -{\pi^2\sigma_{\rm G}^2 \over 2 }} < \epsilon$.
Then the smallest ICF standard deviation $\sigma_{\rm G}$ is given by

$$% \sigma_{\rm G} = \sqrt{-{2 \log{\epsilon} \over \pi^2}}.$$ (8)

Table 1 gives the $\sigma_{\rm G}$ values for different values of $\epsilon$. If the resolution to be achieved is smaller than $\sigma_{\rm G}$, this means that the solution sampling must be fainter than the data sampling.

   

Table 1: ICF standard deviation

$\epsilon$	ICF $\sigma_{\rm G}$	ICF FWHM
10-3	1.18	2.77
10-4	1.37	3.20
10-5	1.53	3.57
10-7	1.81	4.23
10-10	2.16	5.05
10-20	3.05	7.15

3.6.3 Undersampled point spread function

If the PSF is undersampled, it can be used in the same way, but results may not be optimal due to the fact that the sampled PSF varies depending on the position of the source. If an oversampled PSF is available, resulting from theoretical calculation or from a set of observations, it should be used to improve the solution. In this case, each reconstructed object will be oversampled. Equation (4) must be replaced by

 

$$% J(O) = \parallel W - (P_{\rm w} \circ {\cal T} \circ {\cal D}_l) P * O \parallel$$ (9)

where ${\cal D}_l$ is the averaging-decimation operator, consisting of averaging the data in the window of size $l \times l$, and keeping only one average pixel for each $l \times l$ block.