next previous
Up: CLEAN and WIPE

1. Introduction

WIPE| is a Fourier synthesis method recently developed in radio imaging and optical interferometry (Lannes et al. 1994, 1996). The name of WIPE is associated with that of ÇLEAN|, a deconvolution method intensively used in astronomy (Högbom 1974; Schwarz 1978). (WIPE| can also be used as a deconvolution technique.) As a matter of fact, WIPE| was devised, quite independently, on the grounds of well known properties in harmonic analysis and band-limited interpolation (Lannes et al. 1987a). There also exists a probabilistic interpretation of this method (Maréchal & Lannes 1997). Nevertheless, to some extent, WIPE can equally well be regarded as an updated version of ÇLEAN|. This paper is aimed precisely at clarifying the relationship between these two methods.

In Sect. 2, we first show that CLEAN lies in the family of "matching pursuit" techniques (Mallat & Zhang 1993). In other words, the principle of the "iterative harmonic analysis'' of ÇLEAN|, as it is exhibited in Schwarz (1978), can be regarded as a special case of a more general approach. We then present the reconstruction methods that can be implemented for solving the related optimization problems, first, without any control of the propagation of errors (Sect. 3), and afterwards, with a complete control of the stability parameters (Sect. 4). The first weak point of CLEAN can thus be exhibited: in situations of astrophysical interest (for example, when observing extended complex objects with dilute apertures), the matching process of CLEAN is tex2html_wrap2636. As shown in Sect. 5, the regularization principle of WIPE|, which can be applied to CLEAN as it is, remedies this lack of robustness. Unfortunately, CLEAN has another weak point: the basis functions used in the matching pursuit process, the `clean beams,' are not well suited for reconstructing the boundaries of the structuring entities of the images at the selected resolution level. As a result, CLEAN must be interrupted before the best possible fit is reached. As indicated in Sect. 5, the basic version of WIPE overcomes this difficulty in a simple and efficient manner. The multiresolution extension of WIPE|, underlying the main aspects of this paper, reinforces the arguments justifying this methodological choice.

The guiding idea of our analysis is based on the theoretical framework presented in Sect. 1.1. (The reader who is not familiar with the related basic notions is invited to consult Appendices 1 to 4 of Lannes et al. 1987a). As CLEAN was first used as a Fourier synthesis technique, it was natural to compare the principles of CLEAN and WIPE| in this context. Such a comparative study can only be started from a common statement of the problem. The corresponding formulation, which makes the theoretical framework more attractive, is presented in Sect. 1.2. We then also specify the general conditions of the numerical simulations in support of our analysis.

1.1. Theoretical framework

In many inverse problems, the "reconstructed object" is described in an "object representation space"  E generated by  m vectors tex2html_wrap_inline2650 selected among a family tex2html_wrap_inline2652 of M  vectors tex2html_wrap_inline2656. The latter may be regarded as the "atoms" of the object representation under consideration. The linear space generated by all the vectors of tex2html_wrap_inline2658, tex2html_wrap_inline2660, is a subspace of some tex2html_wrap2870 space: the "object space" tex2html_wrap_inline2662. The "data vector" tex2html_wrap_inline2664 lies in another tex2html_wrap2872 space: the "data space" tex2html_wrap_inline2666. (When the data are complex quantities, it is always possible to work in the real linear space underlying the complex linear space directly involved in the analysis.) To a first approximation, tex2html_wrap_inline2668 is related to the object (or image) to be reconstructed via a linear operator  A from tex2html_wrap_inline2672 into tex2html_wrap_inline2674. The basic problem is to choose the vectors tex2html_wrap_inline2676 in tex2html_wrap_inline2678, and thereby to construct the object representation space  E . The solutions are then obtained by minimizing on  E the quadratic functional:
equation272

Clearly tex2html_wrap_inline2684 is the norm on tex2html_wrap_inline2686; the scalar products on tex2html_wrap_inline2688 and tex2html_wrap_inline2690, as well as the norms, will be distinguished by the subscripts o and d, respectively.

Let F be the "image" of E by  A (the space of the tex2html_wrap_inline2698's, tex2html_wrap_inline2700 spanning  E ), tex2html_wrap_inline2704 be the operator from  E onto  F induced by  A , and tex2html_wrap_inline2712 be the projection of tex2html_wrap_inline2714 onto  F (see Fig. 1 (click here)). The vectors tex2html_wrap_inline2718 minimizing  q on  E , the solutions of the problem, are such that tex2html_wrap2874. They are identical up to a vector lying in the kernel of tex2html_wrap_inline2726. Note that tex2html_wrap_inline2728. (By definition, the kernel of A , tex2html_wrap_inline2732, also referred to as the null space of  A , is the set of vectors tex2html_wrap_inline2736 such that tex2html_wrap2876). The condition tex2html_wrap2878 is a necessary condition for tex2html_wrap_inline2742 to be reduced to tex2html_wrap_inline2744.

As tex2html_wrap_inline2746 is orthogonal to  F , the solutions tex2html_wrap_inline2750 of the problem are characterized by the property:
displaymath2638

On denoting by tex2html_wrap_inline2752 the adjoint of A , this property can also be written in the form


displaymath2639
where r is regarded as a residue. The solutions tex2html_wrap_inline2758 are therefore characterized by the fact that the corresponding residue is orthogonal to  E , i.e.,
equation294

This condition is of course equivalent to tex2html_wrap_inline2762, where tex2html_wrap_inline2764 is the projection (operator) of tex2html_wrap_inline2766 onto  E . The solutions of the problem are therefore the solutions of the equation on  E :


equation298

For any tex2html_wrap_inline2772 and any tex2html_wrap_inline2774, we have
displaymath2640
hence tex2html_wrap_inline2776. This explicitly shows that Eq. (3) is the "normal equation":
equation305

When the problem is well posed, tex2html_wrap_inline2778 is a one-to-one map (tex2html_wrap_inline2780); the solution is then unique: there exists only one vector tex2html_wrap_inline2782 such that tex2html_wrap_inline2784. This vector, tex2html_wrap_inline2786, is said to be the least-squares solution of the equation tex2html_wrap2880.

In this case, let tex2html_wrap_inline2790 be a variation of tex2html_wrap_inline2792 in  F , and tex2html_wrap_inline2796 be the corresponding variation of tex2html_wrap_inline2798 (see Fig. 1 (click here)). As shown in Appendix 1, the robustness of the reconstruction process is then governed by the inequality


equation311
where tex2html_wrap_inline2800 is the condition number of tex2html_wrap_inline2802:
equation317
tex2html_wrap_inline2804 and tex2html_wrap_inline2806 respectively denote the greatest lower bound and the least upper bound of tex2html_wrap_inline2808 for the tex2html_wrap_inline2810's with norm unity in E:
equation322

The closer to 1 is the condition number, the easier and the more robust is the reconstruction process. As
displaymath2641
tex2html_wrap_inline2814 and tex2html_wrap_inline2816 are the smallest and largest eigenvalues of tex2html_wrap_inline2818, respectively.

   Figure 1: Uniqueness of the solution and robustness of the reconstruction process. Operator  A  is an operator from the object space tex2html_wrap_inline2822 into the data space tex2html_wrap_inline2824. The object representation space  E is a particular subspace of tex2html_wrap_inline2828 (see text). The image of  E by  A , the range of tex2html_wrap_inline2834, is denoted by  F . In this Euclidean representation, tex2html_wrap_inline2838 is the projection of the data vector tex2html_wrap_inline2840 onto  F . The inverse problem must be stated so that tex2html_wrap_inline2844 is a one-to-one map from  E onto  F , the condition number tex2html_wrap_inline2850 having a reasonable value (see Eq. (5))

In the analysis of the different methods that can be devised and implemented for solving the problem, three main aspects must be examined:

1) the precise definition of tex2html_wrap_inline2852, A and E ;

2) the selected minimization technique;

3) the robustness of the reconstruction process.

The first point is essential for the interpretation of the results. It is related to the statement of the problem and, if need be, to its "regularization;" tex2html_wrap_inline2858 is not necessarily the vector of the experimental data. In such a case, the definition of A must of course take into account the corresponding manipulations. As regards  E , which may also be involved in the regularization of the problem, it is important to note that this space may be constructed, in a global manner or step by step, interactively or automatically. As suggested by the statement of point 2, many different techniques can be used for minimizing  q on  E ; some of these are certainly more efficient than others, but in this case, the choice is not crucial. The last point concerns the study of the propagation of errors. The part played by inequality (5) in the development of the corresponding error analysis shows that a good reconstruction procedure must also provide, in particular, the condition number tex2html_wrap_inline2868.

1.2. Formulation of the problems of Fourier synthesis in astronomy

In the problems of Fourier synthesis encountered in astronomy, the "object function" of interest, tex2html_wrap_inline2900, is a real-valued function of an angular position variable tex2html_wrap_inline2902. The "object model variable" tex2html_wrap_inline2904 lies in some object space tex2html_wrap_inline2906 whose vectors, the functions tex2html_wrap_inline2908, are defined at a very high level of resolution. In this sense, one may say that tex2html_wrap_inline2910 emulates the Hilbert space of square integrable real-valued functions tex2html_wrap_inline2912. More precisely, tex2html_wrap_inline2914 is characterized by two key parameters: the extension tex2html_wrap_inline2916 of its field, and its resolution scale tex2html_wrap_inline2918; tex2html_wrap_inline2920 is therefore much smaller than the resolution limit that can be reasonably expected at the end of the reconstruction process: tex2html_wrap_inline2922, where N  is some power of 2. (The larger is  N , the more oversampled is the object field.) To define the object space more explicitly, we introduce the finite grid (see Fig. 2 (click here)):


displaymath2886

Let tex2html_wrap_inline2928 be an integer vector lying in tex2html_wrap_inline2930: tex2html_wrap_inline2932. Clearly, the functions


equation350
where
equation354
form an orthogonal set with
equation361

In our comparative analysis of CLEAN and WIPE|, tex2html_wrap_inline2934 is the tex2html_wrap3136 space generated by the basis vectors tex2html_wrap_inline2936, tex2html_wrap_inline2938 spanning  G (see Fig. 2 (click here)). The functions tex2html_wrap_inline2942 lying in this space can therefore be expanded in the form:
equation370

In the general context of this paper, the functions tex2html_wrap_inline2944, which play the role of interpolation or scaling functions, can be regarded as the "elementary particles" of the object representation. Evidently, other orthogonal sets of scaling functions can be used (cf. Sect. 4.3 of Lannes et al. 1994). Let us finally note that tex2html_wrap_inline2946 is a subspace of tex2html_wrap_inline2948; its inner product can be explicitly expressed in the form (cf. Eqs. (8)-(11)):
displaymath2887
Thus,
displaymath2888

   Figure 2: Object grid tex2html_wrap_inline2950 for N = 8

The Fourier transform of tex2html_wrap_inline2954 is defined by the relationship
displaymath2889
where tex2html_wrap_inline2956 is a two-dimensional angular spatial frequency: tex2html_wrap3144. According to Eqs. (11), (8) and (9), we therefore have:
displaymath2890
where
displaymath2891
with tex2html_wrap_inline2960 and
displaymath2892

The "experimental data" tex2html_wrap_inline2962 are blurred values of tex2html_wrap_inline2964 on a finite list of frequencies in the Fourier domain: tex2html_wrap_inline2966.

As tex2html_wrap_inline2968 is a tex2html_wrap3152 function, it is natural to define tex2html_wrap_inline2970 as the complex conjugate of tex2html_wrap_inline2972. The `experimental frequency list' tex2html_wrap_inline2974 is defined consequently: if tex2html_wrap_inline2976, then tex2html_wrap_inline2978 (except for tex2html_wrap3154: in the convention adopted in the sequel, either the null frequency does not lie in tex2html_wrap_inline2982, or there exists only one occurrence of this point). The experimental frequency coverage generated by tex2html_wrap_inline2984 is therefore centrosymmetric (see Fig. 3 (click here)a).

The experimental data vector tex2html_wrap_inline2986 lies in the "experimental data space" tex2html_wrap_inline2988. By definition, tex2html_wrap_inline2990 is the real Euclidean space underlying the space of tex2html_wrap3156 functions on tex2html_wrap_inline2992 such that tex2html_wrap3158. This space, whose dimension is equal to tex2html_wrap_inline2996, is endowed with the scalar product
equation434
where tex2html_wrap_inline2998, and W is a given weighting function (see for instance Eqs. (17) and (18) further on). Note that this scalar product is real: tex2html_wrap_inline3002. In ÇLEAN|, the data space tex2html_wrap_inline3004 reduces to tex2html_wrap_inline3006, whereas in WIPE|, tex2html_wrap_inline3008 is a natural extension of tex2html_wrap_inline3010 (cf. Sect. 5.1).

Let tex2html_wrap_inline3012 be the Fourier domain: tex2html_wrap_inline3014. In Fourier synthesis, the frequency coverage to be synthesized is a centrosymmetric region tex2html_wrap_inline3016. In Fig. 3 (click here)a, this region is the disc centred on the origin. CLEAN and WIPE share a common objective, that of the image to be reconstructed. This image, tex2html_wrap_inline3018, is defined so that its Fourier transform is quadratically negligible outside tex2html_wrap_inline3020. More explicitly, tex2html_wrap_inline3022 is defined by a convolution relation of the form:
equation447

The "synthetic beam" s is a function resulting from the choice of tex2html_wrap_inline3026: the "clean beam" in ÇLEAN|, the "neat beam" in WIPE (see Fig. 3 (click here)b). The Fourier data corresponding to tex2html_wrap_inline3028 are then defined by the relationship:
equation452

Clearly, tex2html_wrap_inline3032 lies in tex2html_wrap_inline3034.

The neat beam can be regarded as a sort of optimal clean beam: the optimal apodized point spread function that can be designed within the limits of the tex2html_wrap3162 principle. (Apodization theory is concerned with the control of the distribution of light over the exit pupil of an optical system in order to achieve a suppression of side lobes of the diffraction pattern -- Slepian 1965). More precisely, the neat beam  s is a centrosymmetric function lying in the object space tex2html_wrap_inline3038, and satisfying the following properties:

In the terminology adopted in this paper, an atom such as  s is a finite linear combination of elementary particles tex2html_wrap_inline3096. The integers tex2html_wrap_inline3098 involved in this linear combination lie in a subset tex2html_wrap_inline3100 of  G (see Figs. 2 (click here) and 3 (click here)b). As explicitly shown in Sect. 2 of Lannes et al. (1996), the computation of the corresponding coefficients does not raise any particular difficulty.

In a first approach to Fourier synthesis, Eq. (13) suggests that the basis functions of the object representation space  E should be translated versions of  s : a finite number of atoms centred on the nodes of grid  G . This very natural idea, which is exploited as it is, in the matching pursuit principle of CLEAN (cf. Sect. 2.2), may be completely relaxed in WIPE|. More precisely, the matching pursuit process may be performed at the level of the elementary particles. The resolution limit of the reconstruction process is then kept under control thanks to an appropriate regularization principle (cf. Sect. 5.1).

The simulations presented in this paper correspond to the general conditions defined in Fig. 3 (click here). The "object function" tex2html_wrap_inline3110 was defined as a set of "tex2html_wrap_inline3112 functions" centred on the nodes of grid tex2html_wrap_inline3114 tex2html_wrap3184. The image to be reconstructed (tex2html_wrap_inline3118) then lies in tex2html_wrap_inline3120. Figure 4 (click here) gives an idea of the complexity of this image.

  figure509
Figure 4: Image to be reconstructed, tex2html_wrap_inline3122, at the resolution level defined in Fig. 3 (click here). Note that the portion of the field represented here is twice as large as that defined in Fig. 3 (click here)b

The Fourier data were blurred by adding a Gaussian noise. More precisely, for all the tex2html_wrap3188 (see Fig. 3 (click here)a), the standard deviation of tex2html_wrap_inline3126 was set equal to 5% of the total flux of the object: tex2html_wrap3190. The weighting function  W introduced in Eq. (12) was explicitly defined by the formula
equation521
where
equation528


next previous
Up: CLEAN and WIPE

Copyright by the European Southern Observatory (ESO)
web@ed-phys.fr