15.7 Deconvolution

The first imaging step presented before leads to a convolution equation whose solution is the convolution product of the sky brightness distribution (apodized by the interferometer primary beam) by the dirty beam.

To derive the astronomically meaningful result, i.e. ideally the sky brightness, a deconvolution is required. Deconvolution is always a non linear process, and requires (in one way or another) to impose some constraints on the solution, or in other words to add some information, to better select plausible solutions. Such additional constraints can be explicit (e.g. positivity, or user specified finite support) or qualitative.

15.7.1 The CLEAN method

The standard deconvolution technique, CLEAN relies on such a qualitative constraint: it assumes that the sky brightness is essentially a ensemble of point sources (the sky is dark, but full of stars). The algorithm which derives from such an assumption is straightforward. It is a simple ``matching pursuit''

1. Initialize a Residual map to the Dirty map
2. Initialize a Clean component list to zero.
3. Assume strongest feature in Residual map originates from a point source
4. Add a fraction $\gamma$ (the Loop Gain) of this point source to the Clean component list, remove the same fraction, convolved with the dirty beam, from the Residual map.
5. If the strongest feature in the Residual map is larger than some threshold, go back to point 3 (each such step is called an iteration).
6. If the strongest feature is below threshold, of if the number of iterations $N_{\mathrm{iter}}$ is too large, go to point 7.
7. Convolve the Clean component list by a properly chosen Clean Beam (this is called the restoration step).
8. add to the result the Residual map to obtain the Clean Map.

The CLEAN algorithm as a number of free parameters. The loop gain controls the convergence of the method. In theory, $0 < \gamma < 2$, but in practice one should use $\gamma \simeq 0.1 - 0.2$, depending on sidelobe levels, source structure and dynamic range. While high values of $\gamma$ would in principle give faster convergence, since the remaining flux is $\propto (1-\gamma)^{N_{\mathrm{iter}}}$ if the object is made of a single point source, deviations from an ideal convolution equation force to use significantly lower values in order to avoid non linear amplifications of errors. Such deviations from the ideal convolution equation are unavoidable because of thermal noise, and also of phase and amplitude errors which distort the dirty beam.

The threshold for convergence and number of iterations define to which accuracy the deconvolution proceeds. It is common practice to CLEAN down to about the noise level or slightly below. However, in case of strong sources, the residuals may be dominated by dynamic range limitations rather than by noise.

The clean beam used in the restoration step plays an important role. It is usually selected as a 2-D Gaussian, which allows the convolution to be computed by a simple Fourier transform, although other choices could be possible. The size of the clean beam is a key parameter. It should be selected to match the (inner part of) the dirty beam, otherwise the flux density estimates may be incorrect. To understand this problem, let us note first that the units of the dirty image are undefined. Simply, a 1 Jy isolated point source appears with a peak value of 1 in the dirty map. This is no longer true (because of sidelobes) if there is more than one point source, or a fortiori, an extended source. The unit of the clean image is well defined: it is Jy per beam, which can easily be converted to brightness temperature from the effective clean beam solid angle and the observing wavelength. Now, assume the source being observed is just composed of 2 separate point sources of equal flux, and that the dirty beam is essentially a Gaussian. Let us clean the dirty image in such a way that only 1 of the 2 point sources is actually included in the clean component list. If we restore the clean image with a clean beam which is, e.g. twice smaller than the original dirty beam, the final result will undoubtedly be odd. The second source would appear extended and have a larger flux than the first one. No such problem appears if the clean beam matches the dirty beam. Admittedly, the above example shows a problem which results from a combination of two effects: an inappropriate choice for the clean beam, and an insufficient deconvolution. However, the second problem always exists to some extent, because of noise in the original data set. Hence, to minimize errors, it is important to match the clean and dirty beams.

Note that in some circumstances, there may be no proper choice. An example is a dirty beam with narrow central peak on top of a broad ``shoulder''. Small scale structures will be properly reconstructed, but larger ones not.

The last step in the CLEAN method plays a double role. On one hand, it protects against insufficient deconvolution. Furthermore, since the residual image should be essentially noise if the deconvolution has converged, it allows noise estimate on the cleaned image.

15.7.2 Interpretation of CLEAN

If CLEAN converged, the Clean component list is a plausible solution of the measurement equation (within the noise), but it is not unique... Hence, because of convolution by the clean beam, the clean image is not a solution. However, besides allowing a reasonable definition of the image unit in case of incomplete convergence, there are two reasons to convolve by a clean beam. First, convolution by the clean beam smears out artifacts due to extrapolation beyond the measured area of the $uv$ plane. This is an a posteriori regularization. Second, the clean components are forced to reside on the grid defined by the image. This discrete representation has a number of limitations (e.g. necessity of negative clean components, limited accuracy due to the finite size of the component list), which are reduced by convolution by the clean beam, because the clean image then has finite resolution and can be properly represented on a discrete grid provided the Nyquist sampling is preserved.

An important property of CLEAN is that (to first order) only the inner quarter of the dirty image can be properly cleaned. This is easily understood when dirty beam and dirty images are computed on the same grid size, since a source at one edge of the inner quarter requires knowledge of the dirty beam sidelobes beyond the map size to be deconvolved from the opposite edge. However, this also remains true if one computes the dirty beam on a twice larger grid than the dirty image: more than the inner quarter can be deconvolved, but because of aliasing, the map edges can never be.

Finally, CLEAN offers a very simple way to impose further constraints on the class of solution which is acceptable, by allowing definition of a support. This can be the standard (simple or multiple) Clean Box available in many non interactive implementations, or a user defined mask in interactive implementations. The search region can even be modified from iteration to iteration to help clean convergence. Such a flexible support is available inside the MAPPING program. Note however that the Clean Box or support should not be too limited: cleaning the noise is necessary too (as well as incorporating negative Clean component).

15.7.3 The CLEAN variants

The original CLEAN method is due to [Hogbom 1974]. Several variants exist.

One of the most popular (CLARK) is due to [Clark 1980], and involves minor and major cycles. In Minor cycles, an Hogbom CLEAN is performed, but with a truncated dirty beam, and only on the list of brightest pixels. This search is fast, because of the dirty beam truncation and because of the limited support. The Clean components identified during the minor cycles are removed at once by a FFT during a Major cycle. Because removal is done by FFT, slightly more than the inner map quarter can be cleaned.

A second variant, called MX, due to [Cotton & Schwab 1984], is similar to the CLARK method, except that the Clean components are removed from the $uv$ table at the Major cycle stage (and thus the imaging process is repeated at each major cycle). This avoid aliasing of sidelobes, allows to clean more than the inner quarter, but is relatively slow because of the re-imaging at major cycles. Unless disk storage is a real problem, a faster result of equal quality is obtained by standard Clean with a twice larger map.

The next variant, called SDI (from [Steer et al. 1984]), is again like the CLARK method, but in Minor cycles, no deconvolution is performed, but only a selection of the strongest components down to some threshold. Major cycles are identical to those of the CLARK method. Although the principle is simple, the implementation is not easy because of normalization subtleties in the minor cycle stage. This method is reasonably well suited for more extended structures, but could become unstable if the threshold is inappropriate.

The Multi Resolution Clean (MRC, [Wakker & Schwartz 1988]) separates the problem in a smooth map and a difference map. Since the measurement equation is linear, both maps can be Cleaned (with Hogbom or Clark method) independently. This is faster than the standard CLEAN because the smooth map can be compressed by pixel averaging, and only fine structure left in difference map, so fewer Clean components are required.

15.7.4 The GILDAS implementation

All the above variants are implemented in the GILDAS software. All of them, except MX, are implemented both as tasks and as interactive commands in the MAPPING program. The later implementation allows definition of a flexible support constraint. The default method is CLARK. SDI & MRC are usually not necessary for Plateau de Bure, because of the small ratio between the field of view (primary beam) and the resolution ($< 30$).

MX is implemented only as a task, and not recommended because of its relatively slow speed. Since Plateau de Bure images are relatively small ( $128 \times 128$), it is easier to use a standard clean on larger images.

The GILDAS software does not include any implementation of the Maximum Entropy Method, MEM. The main reason is that MEM is not suited for limited $uv$ coverage. But MEM also has some undesirable properties, among which its attempt to give a unique solution, with no physical justification, the noise dependent resolution, and the definition of a global criterium for adjustment to data. Furthermore, no noise estimate is possible on MEM deconvolved images.