---

• Outline
• Glossary
• 6. Deconvolution in Imaging
  • Previous: 6.1 Sky Models
  • Next: 6.3 CLEAN Implementations

---

Import standard modules:

Import section specific modules:

---

6.2 Interative Deconvolution with Point Sources (CLEAN)

If we carry on with the logic from the previous section ➞ that we can use point sources as a generic unit for building sky models then we are ready to introduce the standard deconvolution method in radio interfereometry, CLEAN. Deconvolution is a broad term, that is defined by different, discipline-specific techniques to separate out two functions which have previously been (approximately) convolved together. Such as in the case of the array point spread function and the sky.

This problem of deconvolution seems fairly simple at first glance. We know how to convolve two functions $f$ and $g$ to produce an new function $h$ where

$$h = f \circ g \qquad.$$

Given, $h$ and one of the other functions, say $g$ then $f$ can be deconvolved by using the convolution theorem

$$\mathscr{F}\{f\} = \frac{\mathscr{F}\{h\}}{\mathscr{F}\{g\}}$$

Easy, we have solved deconvolution, unfortunately this is only true for a noise-less system. And, as we will see in the next chapter ➞, our interferometric observation definitely contains noise in the visibility measurements. When we are measuring the sky visibility function we are convolving the sky with the array response and some instrumental noise $\epsilon$.

$$I^D = \text{PSF} \circ I_{true} + \epsilon$$

If this $\epsilon$ is small compared to the sky signal then we can approximate $\epsilon$ as 0 and do this simple deconvolution method, which is called inverse filtering. Unfortunately, in radio astronomy, and most cases where deconvolution is required, $\epsilon$ is not small. Applying an inverse filter to a noisy signal will result in a poor estimate of the original signal. As an example of why inverse filtering fails in synthesis imaging let us apply the method to the simple KAT-7 observation we have been using in previous sections. As this is simulated data we know the exact sky model and produce a dirty image which is free of any calibration issues, i.e. the observation is full known.

Left: Dirty image of a point source field observed by KAT-7. Right: PSF response of KAT-7 for a 6 hour observation at a declination of $-30^{\circ}$.

Left: a noise sky model derived by applying an inverse filter. Right: The true point source sky model.

We compute the inverse filtered deconvolved sky model using the array PSF and dirty image:

$$I_{true} = \mathscr{F}^{-1} \left \{ \frac{\mathscr{F}\{I^D\}}{\mathscr{F}\{\textrm{PSF}\}} \right \} $$

The resulting sky model in the left figure above is just noise with some bright structure along the edges due to the edge effects. As we can see, we do not recover any sources in the true sky model (right figure above). We need a different method to deal with the noise.

Perhaps we can modify this inverse filter method by only selecting a part of the image which is clearly made up of sources and not noise. We can do this by setting a minimum threshold on the dirty image and applying the inverse filter to the thresholded image.

Left: Dirty image thresholded to a minimum of 2. Right: The sky model recovered by applying inverse filtering to the thresholded dirty image.

By applying a threshold to the dirty image (figure above on the left) we are effectively selecting out only signal so then we should be able to recover a sky model of at least the brightest sources. Unfortunately, we do not, as can be seen in the figure above on the right. Instead we see rings around the bright sources. This is because we are removing the PSF structure around the sources in the dirty image. When we apply the inverse PSF filter those missing sidelobes result in source edge effects, i.e. rings. Now, this is not terrible as we could go through the image detecting circles with a Hough transform and come up with a sky model of positions (but not necessarily flux). But, surely we could do something even better.

First, we know the PSF response of the array. Second, we know that every source in the sky is convolved with (approximately) that PSF. We just need to find sources, which are not due to noise, and subtract the PSF multiplied by a factor (to account for different flux values) from that position in the image, and then add that position and flux multiplier to a sky model. Further, since each pixel contains some amount of noise or possible flux due to the PSF from some other source, we do not want to subtract all the flux in a single pixel in one go. Instead, we could iteratively subtract a portion of the flux in a pixel until some noise threshold is reached. This is called iterative deconvolution, a type of greedy algorithm, which is the basis what we call the class of CLEAN deconvoloution algorithms in synthesis imaging.

To understand the details of CLEAN lets do a swift review of the visibility sampling function and the array PSF.

6.2.1 Sampling, the PSF and the dirty image

In $\S$ 5.2 ➞ we discussed the effects of sampling in imaging. When we observe with a radio interferometric array we collect visibility data on a range of $uv$ coordinates, which we then Fourier transform into an image. The $uv$ coordinates of the visibility data don't fully cover the $uv$-plane. We saw examples of this in plots of KAT-7 $uv$ converage. For example:

Figure: uv coverage from a declination $-30^{\circ}$, 6 hour, 10 MHz observation with KAT-7.

Conceptually, we can think of this partial $uv$ coverage as the product of the complete visibility function with full $uv$ coverage filling the UV-plane $V(u,v)$, with a sampling function $S(u,v)$ that is unity (for uniform weighting) where we have $uv$ samples and zero where we don't:

$$ S(u,v) \cdot V_{true}(u,v) $$

From the van Cittert-Zernike Theorem we know our image on the sky is the Fourier transform of the conceptual complete visibility function:

$$ I_{true}(l,m) = \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} V_{true}(u,v) \; e^{-2\pi i(ul+vm)} \; du \, dv$$

But what we get from the interferometer isn't $V_{true}(u,v)$, it is the incomplete $uv$ coverage product $S(u,v) \cdot V_{true}(u,v)$:

$$ \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} S(u,v) \cdot V_{true}(u,v) \; e^{-2\pi i(ul+vm)} \; du \, dv \;=\; ?? $$

How do we deal with this? The Convolution Theorem states that multiplication in one domain is equivalent to convolution in the Fourier domain:

$$ h \cdot g \;\; \rightleftharpoons \;\; H \; \circ \; G $$

In our incomplete $uv$-coverage situation, we are multiplying the conceptual complete visibility function $V_{true}(u,v)$ by a sampling function $S(u,v)$, which equates to a convolution in the image domain:

$$ S(u,v) \cdot V_{true}(u,v) \;\; \rightleftharpoons \;\; \textrm{PSF}(l,m) \; \circ \; I_{true}(l,m) $$

where $I_{true}(l,m)$ is the true sky image, the Fourier Transform of the complete visibility $V(u,v)$, and $\textrm{PSF}(l,m)$ is the Fourier transform of the sampling function $S(u,v)$, the point spread functions. In older literature the term dirty beam is often used instead of the PSF, they are equivalent. We term the product $\textrm{PSF}(l,m) \; \circ \; I_{true}(l,m)$ the "Dirty Image" $I^D(l,m)$.

Summary

• The Fourier transform of our sampled visibility $V(u,v)$ data gives us a dirty image $I^D(l,m)$.
• The dirty image $I^D(l,m)$ is the convolution of the true image $I_{true}(l,m)$ with the PSF.
• The PSF (or dirty beam $\textrm{PSF}(l,m)$) is the Fourier transform of the sampling function $S(u,v)$

6.2.2 Point Source Assumption

One convenient feature of radio astronomy images is that they consist for the most part of point sources on an empty sky. This is just an approximation of course, but it is an approximation that makes decomvolution much easier, because we can then consider our PSF as a point source response (or impulse response) - every point source on the true image is replaced by a PSF in the dirty image.

Let's set up a little example of this.

First, set up a toy sky, with three point sources and some noise:

Now set up a fake PSF. We will just make up a star shape to be the PSF:

Now we convolve our true sky image $I_{true}(l,m)$ with the PSF, to get the dirty image $I^D(l,m)$

Lets plot all three together to see what we have just done:

6.2.3 CLEAN

Now we are ready to develop a simple version of the iterative deconvolution CLEAN algorithm called Högbom's method ⤴. We will describe Högbom's simple CLEAN method here and discuss advances in CLEAN algorithms in the next section ➞.

In Högbom's CLEAN, the deconvolution is performed by iteratively finding peaks in the image and removing scaled versions of the PSF at the peak pixel positions. It is an image-domain CLEAN, meaning that the deconvolution is done entirely in the image $lm$-plane (covered in the next section is the development of visibility-domain version of CLEAN). Högbom's method is further discussed in the next section ➞. The version of Högbom's method below is intended to be a simlpe example, for a generic implementation of the method see the Högbom's Method notebook ➞.

Högbom's Algorithm (Image-domain CLEAN):

1. Make a copy the dirty image $I^D(l,m)$ called the residual image $I^R(l,m)$.
2. Find the maximum pixel value and position of the maximum in the residual image $I^R(l,m)$.
3. Subtract the PSF multiplied by the peak pixel value $f_{max}$ and a gain factor $g$ from the residual image $I^R(l,m)$ at the position of the peak.
4. Record the position and magnitude of the point source subtracted in a model, i.e. $g \cdot f_{max}$.
5. Go to (Step 2.), unless all remaining pixel values are below some user-specified threshold or the number of iterations have reached some user-specified limit.
6. Convolve the accumulated point source sky model with a restoring beam, termed the CLEAN beam (usually a Gaussian fitted to the main lobe of the dirty beam)
7. Add the remainder of the residual image $I^R(l,m)$ to the CLEAN image formed in (6.) to form the final restored image.

• Input: Dirty image, PSF
• Parameters: gain, iteration limit OR flux threshold
• Output: Sky model, residual image, restored image

Let us use the algorithm steps above to deconvolve the sky model from the dirty image:

Let us plot the original dirty image and the final residual image scaled to see the residuals:

After a few iterations the residual image is beginning to look noise like, with some PSF structure remaining. As the remaining pixel value approach the noise level we need to halt deconvolution. This is a tricky decision and many issues occur when we 'over' or 'under' deconvolve an image. The resulting artefacts are further exasperated by poor calibration, but that will come later. For now though, the process of halting an iterative deconvolution requires setting either a mimimum flux threshold which halts when all pixels in the residual image are at or below the threshold, or we keep track of the number of iterations and halt after a specified number of cycles.

As for the final sky model, we can print the components.

As we can see, many of the components are at or near the same pixel values. We should expect this as we are not removing all the flux from a pixel during an iteration, only a portion of the flux which is determined by the gain $g$ scale factor. These are not separate sources, but in the sky model they are presented as such. An additional step of source finding ($\S$ 6.5 ➞) needs to be applied to the model in order to combine these sources.

Plotting the deconvolved sky model image next to the original true sky model image (below) we can see that the three sources are in the new sky model, though the recovery is not perfect. Perhaps we did under-deconvolved, perhaps we are including some noise into the sky model, there are a number of issues that cause deconvolution to be imperfect. But, what we can see is that iterative deconvolution produces more realistic results compared to the inverse filtering methods we started the section discussing.

The last steps of CLEAN are somewhat optional, they are performed to produce a 'restored' image, this can be thought of an idealized image produced by an interferometric array with the same resolution as the observing array but with all spatial modes (up to the maximum resolution) fully sampled. That is, a restoring beam, which is usually taken to be a 2-D Gaussian of the same width as the main lobe of the PSF, is convoled with the point source sky model. This is done to reintroduce the the resolution limits of the array without including the PSF sidelobe structure. The restored image is the 'pretty' image we like to show off.

In our toy example here it is hardly worth doing a proper Gaussian fit, but let's do it anyway:

And now we can convolve the sky model with the restoring beam and add back in the residuals to produce the final restored image:

Now let us plot the results of the CLEAN:

Notice how the CLEAN is imperfect!

Even in a simple field of distant point sources, with no noise, PSF side-lobes can overlap and add up which leads to wrongly identifying sources with CLEAN. And as soon as we add noise the situation degrades further - when we start CLEANing close to the noise level then noise starts being added to the model as sources. Also individual point sources can be subtracted in multiple components in adjacent pixels.

But having said that - CLEAN usually performs adequately. As long as we are careful to look at the results of our CLEAN, and are aware of image artifacts and noise levels!

There are also lots of parameters to fiddle with, even in our simple CLEAN. Try changing the gain and cutoff threshold and re-run this notebook, to see how they effect the CLEAN. Also try changing the noise level of the original map by changing noise_rms - you will see how the CLEAN performs better when the signal-to-noise ratio is higher.

In the next section we will dive deeper into the different implementations of CLEAN, limitations of point source-based deconovlution, and current developments in advancing deconvolution methods.

---

Next: 6.3 CLEAN Implementations