An improved method for polarimetric image restoration

An improved method for polarimetric image restoration in interferometry

Abstract

Interferometric radio astronomy data require the effects of limited coverage in the Fourier plane to be accounted for via a deconvolution process. For the last 40 years this process, known as ‘cleaning’, has been performed almost exclusively on all Stokes parameters individually as if they were independent scalar images. However, here we demonstrate for the case of the linear polarisation , this approach fails to properly account for the complex vector nature resulting in a process which is dependant on the axis under which the deconvolution is performed. We present here an improved method, ‘Generalised Complex CLEAN’, which properly accounts for the complex vector nature of polarised emission and is invariant under rotations of the deconvolution axis. We use two Australia Telescope Compact Array datasets to test standard and complex CLEAN versions of the Högbom and SDI CLEAN algorithms. We show that in general the Complex CLEAN version of each algorithm produces more accurate clean components with fewer spurious detections and lower computation cost due to reduced iterations than the current methods. In particular we find that the Complex SDI CLEAN produces the best results for diffuse polarised sources as compared with standard CLEAN algorithms and other Complex CLEAN algorithms. Given the move to widefield, high resolution polarimetric imaging with future telescopes such as the Square Kilometre Array, we suggest that Generalised Complex CLEAN should be adopted as the deconvolution method for all future polarimetric surveys and in particular that the complex version of a SDI CLEAN should be used.

keywords:
techniques: interferometric – techniques: polarimetric – polarisation

1 Introduction

The challenge of deconvolution is wide-spread in radio astronomy, ranging from aperture synthesis (Ryle & Hewish, 1960) to rotation measure synthesis (Brentjens & de Bruyn, 2005). In particular, the technique of aperture synthesis allows observers to use radio interferometers to observe the sky at higher resolution and sensitivity than possible with the largest individual radio telescopes. Such advantages provide the motivation for the next-generation interferometric radio telescopes, such as the LOw Frequency ARray (LOFAR; van Haarlem et al., 2013), the Murchison Widefield Array (MWA; Tingay et al., 2013), the Australian Square Kilometre Array Pathfinder (ASKAP; Hotan et al., 2014) and the Square Kilometre Array (SKA). However, to accurately interpret the observations of a radio interferometer, a method of image reconstruction needs to be applied. Image reconstruction in radio interferometry is often referred to as deconvolution, because it is used to deconvolve the effects of limited uv-coverage on the observation.

A number of possible methods can be used to solve the deconvolution problem including CLEAN, maximum entropy and compressed sensing. Maximum entropy methods are in particularly useful for deconvolution of mosaic images (Cornwell & Evans, 1985; Cornwell, 1988). Additionally modern methods of deconvolution, such as compressed sensing, have begun to be investigated for next-generation radio interferometers (Li et al., 2011; Carrillo et al., 2014; Garsden et al., 2015), however they are still in their infancy and currently limited to very narrow observational bandwidths.

Historically, the CLEAN algorithm and its variations have been successfully applied as a standard method of deconvolution in radio interferometry. The CLEAN algorithm was first introduced by Högbom (1974), and is commonly referred to as Högbom CLEAN, and presents a statistically correct least-squares fit of the measured -data (Schwarz, 1978). Since then, variations of the CLEAN method have been developed to improve the speed and quality of the algorithm (Clark, 1980; Schwab, 1984). Furthermore, variations have been developed to improve deconvolution of resolved and extended structures, such as SDI (Steer-Dwedney-Ito) CLEAN (Steer et al., 1984) and Multi-Scale CLEAN (Cornwell, 2008). With the recent introduction of next-generation radio telescopes with non-coplanar baselines (LOFAR, MWA), wide-field imaging variations of CLEAN have also been developed (Offringa et al., 2014).

One of the strengths of radio astronomy is that polarisation can be measured accurately at radio wavelengths. This enables science that is not accessible at other wavelengths, such as the study of cosmic magnetism which is a key driver of next-generation telescopes like the SKA (Johnston-Hollitt et al., 2015). Additionally, the observed polarisation is vital to calibration, which can effect image quality in total intensity images.

To date, the main motivation for deconvolution and for the CLEAN algorithm has been the deconvolution of total intensity (Stokes ) images. In practice, the same CLEAN algorithms are applied to both the intensity and the polarimetric image, Stokes , , and which respectively represent the two linear and one circular polarisations. The polarised Stokes images themselves are intermediate products produced to extract either the total linear polarisation or the total polarisation which are complex valued vector quantities derived from the component Stokes parameters. In the majority of cases, celestial sources have no detectable circular polarisation, while detectable linear polarisation is common. Thus, of particular interest is the total linear polarisation, , which is given by:

(1)

Because deconvolution has been motivated for Stokes , the Stokes and images are typically deconvolved individually, rather than together as a complex linear polarisation image as suggested by Equation 1. As a result the CLEAN methods used for the past 40 years do not correctly account for the complex vector nature of 1.

All observations with interferometers must be referenced to a standard coordinate system, which can be either fixed with respect to some external coordinate system (basis) or derived from the physical orientation of the crossed dipoles or circular feed horns as they are either fixed on the Earth (e.g. for the MWA), or in the case of a dish array, when the instrument is stowed (e.g. for the Australia Telescope Compact Array and Jansky Very Large Array for linear or circular feeds, respectively). The position angle of the polarisation, , of an electric field is measure by an interferometer as components of the two linear polarisations Stokes and made with respect to the coordinate axis such that:

(2)

For a given polarisation position angle, the magnitude of the measured Stokes and components will change when measured in a rotated coordinate system. However, the final complex polarisation vector reconstructed from the measured linear polarisations is rotationally invariant. The physical linear polarisation vector as a geometric object does not depend on the coordinate system used to measure it. For this reason, rotation of the dipoles or circular feeds with respect to the source during Earth rotation synthesis must be corrected and transformed back to the original coordinate system during calibration via introduction of a phase rotation term 2. Having taken care to correctly account for the rotation of the dipole or circular feed as sources are tracked across the sky, it follows that one should apply a deconvolution process that produces results which are rotationally invariant. The observation of linear polarisation vectors restored from the deconvolution process should also be geometric objects that are independent of coordinate system. Otherwise, the linear polarization images may depend on the choice of coordinate systems. The standard CLEAN methods are not invariant under rotations of when and are CLEANed separately, and the final sky model constructed by standard CLEAN algorithms depends on the chosen orientation of the and axes. In this work, we demonstrate that the CLEANed model will depend on the orientation chosen to perform the deconvolution, then provide an alternative CLEAN deconvolution method that constructs a model independent of the orientation of the and axes.

The alternative method is called ‘Generalised Complex CLEAN’, and is a modification of the CLEAN algorithm. Generalised Complex CLEAN works by locating residual peaks in linear polarisation, rather than locating the peaks independently in Stokes and . Once the peaks are found, an appropriate model is subtracted according to the CLEAN algorithm variant. Here we demonstrate that the deconvolution process is invariant under rotations of , meaning that the residuals do not depend on the chosen axis for and . We then show that complex CLEANing will provide an advantage over traditional CLEANing by detecting more CLEAN components in sources near the noise level, whilst simultaneously detecting less spurious components scattered across the image.

The paper is arranged as follows: first we show using high quality Australia Telescope Compact Array (ATCA) observations that the model constructed by CLEANing and independently is biased to the chosen axes for deconvolution. Then we introduce the Complex Högbom and SDI CLEAN methods, and describe how they can be applied to deconvolve images of complex linear polarisation . Finally, we then demonstrate that complex CLEANing produces a CLEAN model that does not depend on the choice of axes used for deconvolution and show that it produces superior results to current methods.

By showing that the Högbom and SDI CLEAN methods can be extended to not depend on the choice of axes used for deconvolution, we set the stage for extending other methods of deconvolution from real valued images to complex valued images.

2 Observational Data

In this work, we use two ATCA observations to compare the effects of using a standard and Generalised Complex CLEAN on synthesized linear polarisation images. In this section, we present the details of the two observations.

Figure 1: UV coverage of PKS J0334-3900 in units of kilo-. Note that this is the uv-coverage based on the full bandwidth.
Figure 2: 1.4 GHz image of PKS J0334-3900 used in this work. The RMS noise of the Stokes image is 140 Jy b and the beam is 1.2”5.8”, pa = 11.6.

The first data set we use is a 1.4 GHz observation of the radio galaxy PKS J0334-3900 located in the galaxy cluster Abell 3135. PKS J0334-3900 is an intermediate FR I/II radio galaxy with diffuse linearly polarised jets classified as a bent tailed source. The full details of the observations and polarisation analysis can be found in our previously published work (Pratley et al., 2013) and the observations used here are given in Table 1. However, since that work, the data has been reprocessed using a standard self-calibration routine, improving the sensitivity of the 1.4 GHz images, and it is these reprocessed data that are used here. As before the synthesized images were made using uniform weighting which results in a synthesized beam of 1.2 x 5.8 arcsec. A plot of the coverage can be found in Figure 1 and Figure 2 shows the 1.4 GHz image.

Source Config. Time RMS noise Beam
(min) (Jy/beam) (arcsec, deg.)
PKS J0334-3900 6A 59 160 1.155.82, 11.62
1.5A 76
750A 79.7
375 75.4
PKS B1637-771 6A 589.8 35 5.475.05, -36.73
6B 571.2
6C 637.8
6E 602.4
750D 723.6
Table 1: Observation time per baseline configuration, and RMS noise in the final Stokes and images generated from co-added uv-data from all available configurations.
Figure 3: UV coverage of PKS B1637-771 in units of kilo-. Note that this is the uv-coverage based on the full bandwidth.
Figure 4: 1.4 GHz image of PKS B1637-771 used in this work. The RMS noise of the Stokes image is 35 Jy b and the beam is approximately circular and 5” in diameter.

The second dataset is from an observation of the radio galaxy PKS B1637-771 at 1.4 GHz. PKS B1637-771 is a diffuse double lobed radio galaxy, previously stated to be an FRII galaxy (Sadler et al., 2014). A low-surface brightness cocoon of emission surrounds the jets of the radio galaxy, and the complex wispy structure of the cocoon suggests some disruption in the environment or jets. As the jets and cocoon are linearly polarised, this is an excellent example of a highly complex source with diffuse polarised emission and is thus a good case on which to test polarisation cleaning methods.

All ATCA observations of PKS B1637-771 were extracted from the archive (PI: Leahy, project code C888). They were performed between 2000 and 2001 using a phase centre of RA = 16h 44m 20.00s, Dec = -77 15’ 30.00”. The observations were centred at at frequency of 1384 MHz with a bandwidth of 128 MHz over 32 channels. The channels on the edge of the band were removed, leaving an effective bandwidth of 96 MHz. Due to self-generated interference in the 1384 MHz observation, the 8 MHz interval centred at 1408 MHz was also removed.

The observation used five baseline configurations (see Table 1), providing excellent UV-coverage. Taking all configurations into account, the total integration time was 50.78 hours (3046.8 minutes). The flux density scale was set relative to the unresolved calibrator PKS B1934-638, using an assumed total flux density at 1384 MHz of 14.94 0.01 (Reynolds, 1994). The time variations in complex antenna gains and bandpass were calibrated using observations of the unresolved source PKS B1718-649. The observations were calibrated using MIRIAD, while following the standard ATCA pre-CABB reduction method (Sault et al., 1995). Synthesized images were made using uniform weighting of the calibrated data, resulting in a synthesized beam with full width half maximum of 5.5 x 5.0 arcsec. A plot of the coverage can be found in Figure 3 and Figure 4 shows the final Stokes image.

The and root-mean-squared (RMS) noise for each observation can be found in Table 1.

The data from PKS J0334-3900 were used to compare the standard and complex Högbom CLEANs, whilst the observation of PKS B1637-771 was used to compare both the standard and complex SDI CLEANs. In the case of PKS J0334-3900 the standard Högbom CLEAN on the polarised data was found to produce publishable images (e.g. Pratley et al., 2013), however for PKS B1637-771 Högbom CLEAN produces a low quality deconvolution of the observation for this diffuse source, and is plagued by artefacts illustrating the well known poor performance of Högbom CLEAN on low surface brightness, complex emission. For such sources an SDI CLEAN performs better when inspecting the residuals, and has no visible artefacts. The need to accurately CLEAN diffuse sources partly motivates the need for developing the complex SDI CLEAN method.

Symbol Description
, , , Stokes parameters
Linear polarisation,
Synthesized beam, point spread function
Synthesized image of respective intensity
Linear polarisation synthesized image,
Stokes and Gaussian noise
Stokes and CLEAN component images, model images
Stokes and residual images
Stokes and final residual images
Linear polarisation final residual image
Linear polarisation CLEAN component image, model image
Model generated from SDI CLEAN clip.
Synthesized model .
Gain used in CLEAN
Clip value used in SDI CLEAN
Damping factor used in SDI CLEAN
Standard deviation of Gaussian noise in Stokes and
Image of average complex linear polarisation
Constraint used as a cutoff for CLEAN, i.e.
Chosen coordinate axes/frame for Stokes and
, Linear Stokes parameters represented in a rotated frame.
Table 2: List of symbols.

3 Standard CLEAN method applied to linear polarisation data

The standard clean procedure picks the highest peaks in the and synthesized images and , respectively, and then iteratively removes the synthesized beam associated with those peaks via deconvolution. This is the deconvolution used in all major radio astronomy imaging software including AIPS, Newstar (Ikeda et al., 2001), MIRIAD (Sault et al., 1995), CASA (McMullin et al., 2007), WSClean (Offringa et al., 2014) and the pipelines based on these packages for instruments such as LOFAR, MWA and ASKAP. Even in more recent packages such as CASA and WSClean where there is the ability to deconvolve in different ways, the standard method of searching for peaks in and independently in and , respectively, is most often used. 3

However, and are treated as individual scalar images, which doesn’t account for the vector nature of linear polarisation, meaning that the arbitrarily chosen linear polarisation axes directions will be preferentially CLEANed. Here we will lay out the mathematical foundations for the CLEANing process and demonstrate its rotational dependance on the two datasets mentioned above by rotating the coordinate system before CLEANing.

First we define a standard set of symbols in Table LABEL:tab:symbols. In this work, the hypothesis is that deconvolution process, represented by the function , is not invariant under rotations of linear polarisation,

(3)

Saying that is rotationally invariant, is the same as the mathematical statement

(4)

being true for all angles .

If we define linear polarisation in a rotated frame as

(5)

then a rotationally invariant CLEAN method would provide CLEAN components that can be related to the non-rotated frame by a simple rotation

(6)

However, if CLEAN is not rotationally invariant, then

(7)

Rotational invariance of CLEAN is critical for physical results. For instance, the orientation of the dipole could be used as the reference frame of and , if CLEAN is not rotationally invariant then rotating the dipole of the antenna could produce different scientific results.

Next, we present the effects of deconvolving synthesized Stokes and images individually using standard Högbom and SDI CLEANs. The images were CLEANed until the absolute maximum residual in and reached the cutoff of 4 times the root-mean-squared (RMS) noise level. To analyse the CLEAN component images for Stokes and , the polarisation vectors of each component in the model are plotted in a plane. Plots are shown in Figures 5 & 6, showing the linear polarisation distribution of the CLEAN components for the original and rotated frames for a Högbom CLEAN of the PKS J0334-3900 dataset and an SDI CLEAN of the PKS B1637-771 dataset, respectively. For the rotated frames (right panels of Figures 5 & 6) the deconvolution process is performed in the and frame, where and then after deconvolution, the CLEAN components in the and frame are rotated back to the original frame. If the CLEAN processes tested here were rotationally invariant both the original (left) and rotated (right) panels of Figures 5 and 6 would be identical. However, we see from these figures that polarisation distribution of the model is biased to lie along the chosen or axis used for deconvolution and that . Thus, the standard CLEAN methods are not invariant under rotations of linear polarisation and as a result CLEANing using the present methods is not optimal for selecting the most physically important structures.

Figure 5: Plots of the linear polarisation states of the CLEAN components for PKS J0334-3900, generated from Högbom CLEANing Stokes and separately. Left: Polarisation states when the Högbom CLEAN components are generated in the and frame. Right: Polarisation states when the Högbom CLEAN components are generated in the and frame, where is the frame rotated by from the axis and then rotated back to the original coordinate frame. In each plot, the polarisation distribution of the model is biased to lie along the chosen and axes used for deconvolution.
Figure 6: Plots of the linear polarisation states of the CLEAN components for PKS B1637-771, generated from SDI CLEANing Stokes and separately. Left: Polarisation states when the SDI CLEAN components are generated in the and frame. Right: Polarisation states when the SDI CLEAN components are generated in the and frame, where is rotated by from the axis and then rotated back to the original coordinate frame. In each plot, the polarisation distribution of the model is biased by the chosen and used for deconvolution.

Furthermore, we can more clearly demonstrate the lack of rotational invariance by considering the relation and splitting into its real and imaginary components and comparing these for the original and rotated CLEAN processes. Thus in Figures 7 & 8 we show versus and versus , for the Högbom and SDI CLEANs, respectively. If these processes were rotationally invariant all of the components would line along the one-to-one line such that and , shown in blue on the figures. However we see this is far from the case and components are scattered across the plot. In these examples, the standard Högbom CLEAN is slightly more prone to generating biased components than the standard SDI CLEAN, but neither are particularly good.

Figure 7: Plots comparing Högbom CLEAN components generated in the rotated frame and in the frame, for (left) and (right). If standard Högbom CLEAN method was rotationally invariant, components would lie along the relation, shown as the blue line. However, the standard CLEAN method does not follow this relation. The green line represents the least squares straightline fit of the data for and for , showing the difference between the ideal fit.
Figure 8: Plots comparing SDI CLEAN components generated in the rotated frame and in the frame, for (left) and (right). If the complex CLEAN method was rotationally invariant, components would lie along the relation, shown as the blue line. While it follows the line more closely than the standard Högbom CLEAN method in Figure 7, it still does not follow the line. The green line represents the least squares straightline fit of the data for and for , showing the difference between the ideal fit.

4 Generalised Complex Clean

Figure 9: Flow-diagram of the standard Högbom CLEAN algorithm, compared to the Complex Högbom CLEAN algorithm. The standard Högbom CLEAN method deconvolves and individually, but the Complex Högbom CLEAN deconvolves them together. is the synthesized beam, is the residual image, is the CLEAN component image. are the complex linear polarisation synthesized image and final restored image. is the gain parameter, typically .
Figure 10: Flow-diagram of the standard SDI CLEAN algorithm, compared to the complex SDI CLEAN algorithm. The standard SDI CLEAN method deconvolves and individually, but the complex SDI CLEAN deconvolves them together. is the synthesized beam, is the residual image, is the CLEAN component image. are the Stokes synthesized image and final restored image. are the complex linear polarisation synthesized image and final restored image. is the gain parameter, typically . is the clip parameter, typically , and is the model of the clipped residual image, and is the synthesized model. is the scaling of the synthesized model due to the beam volume, where .

In this section we extend the standard CLEAN methods to images of complex linear polarisation. While the methods in this section are simplest to describe using complex numbers, the operations can be translated to operations on the real and imaginary parts in parallel. Where CLEAN assumes that an image is a collection of point sources, the Generalised Complex CLEAN assumes that the image is a collection of point sources with a complex amplitude.

In this work, the synthesized complex polarisation image is defined as

(8)

where is the observed complex linear polarisation, is the instruments point spread function (synthesized beam), is the synthesized complex linear polarisation image, and . and are images of Gaussian noise of the same distribution, with variance . (complex valued) represents any correlation between the noise in and , we expect this term to be close to zero after calibration, because the term is typically due to cross talk between linear or circular feeds. Given that linear polarisation is defined in terms of Stokes and where , it follows that

(9)

and

(10)

The objective of the Generalised Complex CLEAN is to recover the observed linear polarisation given the synthesized complex linear polarisation image and instrumental point spread function (synthesized beam) . Following the traditional CLEAN methods, this can be done by generating a clean component image such that

(11)

where is complex valued, and . In the previous section, the standard method to recover in this context is to use the standard CLEAN methods on and separately to obtain and independently.

The Generalised Complex CLEAN methods start with the residual image equal to the synthesized image, . Each iteration starts by calculating . The brightest sources are located in , which are then modeled and subtracted from . The location and peak of each source is recorded as a CLEAN component in .

When the total number of iterations have been reached, or the residual peak is less than the cutoff , the CLEAN component image is complete and the final residual image has been generated. The restored image is then defined by

(12)

where is a Gaussian beam determined by the resolving power of the observation4. It is also clear that and can be restored separately using the CLEAN components generated from the Generalised Complex CLEAN

(13)

where and . Flow-diagrams for the Complex Högbom and SDI CLEAN algorithms can be found in Figures 9 & 10, where they are compared to their standard counterparts. In the following subsections, we describe the process for a given Complex CLEAN iteration for each of the Complex Högbom and Complex SDI CLEANs.

4.1 Complex Högbom CLEAN

Each iteration starts by calculating , then the position of the peak pixel in is located as . If , or the total number of iterations has been reached, then the Complex CLEAN is complete. Otherwise, the image of the synthesized beam is translated by and scaled by . To improve convergence, the synthesized beam is then further scaled by the gain parameter (typically ). is added to the CLEAN component . The translated and scaled synthesized beam is then subtracted from the residual image .

4.2 Complex SDI CLEAN

First, the value of the maximum pixel in , , is measured. If , or the total number of iterations has been reached, then the Complex CLEAN is complete. A new image, , is generated by performing a clip on the values in , using as the clip parameter, values with a modulus squared below are set to zero. The clipped image is then convolved with the synthesized beam .

Because will have multiple components, needs to be scaled by a damping factor (beam volume factor), . We choose to calculate following the implementation of SDI CLEAN in MIRIAD, however, adapting the calculation for complex numbers. To estimate , the damping factor is calculated by correlating with (a lower limit on may be needed to ensure stability). This correlation is explicitly written as

(14)

It is possible to choose . However, may be close to zero and cause a semi-infinite loop of iterations, so we choose to follow MIRIAD by putting a lower limit of magnitude on , by defining

(15)

This choice of calculation is selected so that is consistent when choosing different and frames. The lower limit of 0.02 on was chosen empirically, and is used in MIRIAD’s standard SDI CLEAN method. Using a lower limit of 0.02 may not always be the optimal choice, in which case one may want to choose a different value for the gain .

Then the gain parameter is used to further scale and . is added to the CLEAN component image , and is subtracted from the residual .

5 Rotational Invariance of Complex CLEAN

In this section, we compare the standard and complex CLEAN methods to show that the Generalised Complex CLEAN is invariant under rotations of linear polarisation.

As with the standard CLEAN procedures discussed in Section 3, we begin by plotting the linear polarisation distribution of the CLEAN components for the original and rotated frames for a Complex Högbom CLEAN of the PKS J0334-3900 dataset and a Complex SDI CLEAN of the PKS B1637-771 dataset in Figures 11 and 12, respectively. Again the for the rotated frames (right panels of Figures 11 and 12) the deconvolution process is performed in the and frame, where and then after deconvolution, the CLEAN components in the and frame are rotated back to the original frame. Compared to the standard Högbom and SDI plots of the same data in Figures 5 & 6, here we see firstly that there is no preferred alignment of components with the deconvolution axes and that both the original and rotated plots show almost identical components. This shows that the complex versions are indeed rotationally invariant and thus correcting picking up physically meaningful components regardless of the axes chosen.

Figure 11: Plots of the linear polarisation states of pixels in the restored source models for PKS J0334-391. The source models were generated from Complex Högbom CLEAN components. Left: Polarisation states of the model when the Complex Högbom CLEAN components are generated in the and frame. Right: Polarisation states of the model when the Complex Högbom CLEAN components are generated in the and frame, where is the frame rotated by from the axis. While there are small variations between the models, there is no preference related to the axes used for deconvolution.
Figure 12: Plots of the linear polarisation states of pixels in the restored source models for PKS B1637-771. The source models were generated from Complex SDI CLEAN components. Left: Polarisation states of the model when the Complex SDI CLEAN components are generated in the and frame. Right: Polarisation states of the model when the Complex SDI CLEAN components are generated in the and frame, where is the frame rotated by from the axis. While there are small variations between the models, there is no preference related to the axes used for deconvolution.

As before we also demonstrate the rotational invariance by considering the relation and splitting into its real and imaginary components and comparing these for the original and rotated Complex CLEAN processes. Figures 13 & 14 show versus and versus , for the Complex Högbom and Complex SDI CLEANs, respectively. Unlike Figures 7 & 8 which use the same data cleaned with the standard methods, here we see the Complex CLEAN components line along the one-to-one line such that and , shown in blue on the Figures. We plot the least squares fit to the data on both plots in green and find it to be in excellent agreement with the one-to-one relation in the case of the Complex Högbom CLEAN and identical for the Complex SDI CLEAN. Again, this demonstrates that the Complex CLEAN methods are rotationally invariant.

Figure 13: Plots comparing the Complex Högbom CLEAN components of PKS J0334-3900 in the and frames, for (left) and (right). It is expected that the Complex CLEAN method should be rotationally invariant and the components should lie along the relation, shown as the blue line. However, deviations from this relation occur as a result of cumulative numerical error in subtractions after many iterations. The green line represents the least squares straightline fit of the data for and for , showing the difference between the ideal fit.
Figure 14: Plots comparing the Complex SDI CLEAN components of PKS B1637-771 generated in the and frames, for (left) and (right). It is expected that the complex CLEAN method should be rotationally invariant and the components should lie along the relation, shown as the blue line. It is known that the SDI CLEAN is less prone to cumulative numerical error when compared to the Högbom CLEAN, because it subtracts many components in a single iteration. This is demonstrated by the components following the relation more closely than in Figure 13. The green line represents the least squares straightline fit of the data for and for , showing the difference between the ideal fit.

Here the scatter from the ideal case starts to become physically meaningful as a method for characterisation of numerical error in the CLEANing process. From the plots it is clear that the Complex SDI CLEAN has less scatter than the Complex Högbom CLEAN and this is consistent with the fact that an SDI CLEAN will produce less numerical error (Steer et al., 1984) after performing a large number of iterations, an advantage over the Högbom CLEAN. To further test this we also performed a Complex SDI clean of PKS J0334-3900 and provide the versus and versus plots in Figure 15. Comparing Figures 13 and 15 shows a stark difference in the scatter which is due to increasing numerical errors associated with underlying single component iterative approach of the Högbom CLEAN over the multiple component cycles used in a SDI clean. Numerically we find that the sample standard deviation (scatter) about the one-to-one line is of order 130Jy for both the real and imaginary parts using the standard Högbom CLEAN of PKS J0334-3900. This improves to around 14Jy when performing a Complex Högbom CLEAN of PKS J0334-3900 and improves further still to around 5Jy for a Complex SDI CLEAN on the same data. We see a similar result for the more complex source, PKS B1637-771 which improves from 40Jy to 4Jy in comparing a standard to a Complex SDI CLEAN. Given that this is measuring scatter in the CLEAN components, this is demonstrating the improvement in the accuracy of the magnitude of the components considered.

Figure 15: Plots comparing Complex SDI CLEAN components of PKS J0334-3900 generated in the and frames, for (left) and (right). The blue line is the ideal relation and the green line represents the least squares straight line fit of the data. It is clear that the Complex SDI CLEAN method has less scatter than the Complex Högbom CLEAN method in Figure 13, since it more closely follows the ideal relation. The green line follows the least squares fit for and for .

In order to illustrate the effect of rotational invariance on the pixel values in the final images of (Re[]) and (Re[]) rather than the CLEAN components, we plot the pixel values in the resultant images in the original and rotated frames against each other. Here we do this for a standard Högbom CLEAN of PKS J0334-3900 (Figure 16), a Complex Högbom CLEAN of PKS J0334-3900 (Figure 17), a standard SDI CLEAN of PKS B1637-771 (Figure 18), and the Complex SDI CLEAN of PKS B1637-771 (Figure 19).

As with the CLEAN component comparison we see a very large scatter about the one-to-one relation for pixel values in both the standard Högbom and SDI CLEANs with the Högbom CLEAN showing the greatest change in pixel values simply due to rotation of the frame in which the CLEAN process is performed. With the Complex Högbom and SDI cleans there is considerable improvement and the SDI pixel values are identical between the rotated and non-rotated frames, while the Högbom CLEAN values show some scatter which is again due to the increased numerical errors of the method.

Figure 16: Plots comparing restored images generated from Högbom CLEAN components of PKS J0334-3900 in the rotated frame and in the frame, for (left) and (right). If standard Högbom CLEAN method was rotationally invariant, all the data points would lie along the relation, shown as the blue line. It is clear that the data points are scattered along this line.
Figure 17: Plots comparing restored images generated from SDI CLEAN components of PKS B1637-771 in the rotated frame and in the frame, for (left) and (right). If the complex CLEAN method was rotationally invariant, all points would lie along the relation, shown as the blue line. It is clear that there is less scatter than the standard Högbom CLEAN method in Figure 16, but the points still do not follow the line closely.
Figure 18: Plots comparing images restored from Complex Högbom CLEAN components of PKS J0334-3900 generated in the and frames, for (left) and (right). It is expected that the complex CLEAN method should be rotationally invariant, and should lie along the relation, shown as the blue line. It is clear that the Complex Högbom CLEAN has less scatter than the non-complex method.
Figure 19: Plots comparing images restored from complex SDI CLEAN components of PKS B1637-771 generated in the and frames, for (left) and (right). It is expected that the complex CLEAN method should be rotationally invariant, and should lie along the relation, shown as the blue line. It is clear that the Complex SDI CLEAN has less scatter than the non-complex method and less scatter than the Complex Högbom CLEAN.

6 Advantages in the Selection of Cut-off values for Components in Complex CLEAN

Having shown that the Generalised Complex CLEAN methods are rotationally invariant, here we discuss the implications of the rotational invariance to application of cut-off values for CLEAN components in the Generalised Complex CLEAN methods and how this differs from the current methods.

The standard CLEAN methods are used to deconvolve images of complex linear polarisation, by CLEANing and restoring the real and imaginary images separately. Unlike Stokes , Stokes and , are expected to have sources of positive and negative values, so the standard CLEAN methods are required to locate peaks of the absolute maximum. In this case, we choose the cutoff to be , where is the standard deviation of the noise in Stokes and , which is typically the same value. Both Stokes and are cleaned until and in each residual image. Furthermore, it is possible for CLEAN to detect a different numbers of components in and , since the component images are generated separately.

As described in Section 4, the Generalised Complex CLEAN methods will generate the component images and simultaneously, and each image will have the same amount of CLEANed components. Since the CLEAN components are located in the peaks of , the Generalised Complex CLEAN will run until .

The advantage of the Generalised Complex CLEAN is that the deconvolution method is rotationally invariant in meaning that in each iteration, the peaks are located in which is also rotationally invariant under rotations ). Thus, the peak amplitude and location for each iteration will also be invariant under rotations and furthermore the cutoff requirement for the Generalised Complex CLEAN is also rotationally invariant, since is the same in the and .

Additionally, if a component satisfies , then one can always rotate to a frame where . Therefore, when is chosen as a cutoff requirement, the Generalised Complex CLEAN subtracts components in and above in all possible frames. By comparison the standard CLEAN methods have a cut-off that does depend on the choice of . Specifically, the standard CLEAN method will not detect components when both

(16)

However, these components are detected by the Generalised Complex CLEAN. Figure 20 shows these constrains against states of possible linear polarisation, where the Generalised Complex CLEAN cut-off is compared to the standard CLEAN cut-off. The effect of this is that when using a standard CLEAN, in order to reach components which are significant (e.g. ) but in regions which are unreachable with a 3 cut off in and , the cut-off is typically relaxed by increasing the number of iterations and hence components. However, this also results in collecting a very large number of spurious CLEAN components which have . In practice one will not actually be able to collect all of the components for which without having an unworkably large number of spurious components and so typically something in the middle is reached were some significant components are missing and a large number of spurious components are included. This scenario both increases the computational requirements and reduces the quality of the final image as compared to a Generalised Complex CLEAN.

6.1 Cut-off in the case of correlations

When there are correlations between Stokes and , which can be caused by instrumental polarisation, the Generalised Complex CLEAN cut-off will change as a function of polariastion angle. While the cut-off is determined by the RMS noise of the signal in the image, it is the variance of the noise that is linear under addition and subtraction of signals. For this reason, the covariance matrix is used to determine how the RMS noise changes under rotations of . We can write the non-diagonal covariance matrix as

(17)

Physically, the eigenvectors of the covariance matrix can be used to determine the direction of the instrumental polarisation, assuming the instrumental polarisation causes the signal’s polarisation to lie along a preferential direction.

We can change the basis of the covariance matrix to see how the coefficients vary as a function of angle

{widetext}
(18)

For any given component at angle , we can rotate our basis so that it lies along the new axis. In this frame, the RMS noise in is

(19)

which is used to determine if the component is above three times the RMS noise level. Therefore, the cut-off as a function of the component’s polarisation angle is

(20)
Figure 20: a) and b) show the phase space of (). a) The black circle of radius represents polarisation states below the detection threshold of . The surface of the circle represents the cutoff for subtracting complex CLEAN components. The square represents polarisation states at the detection level of 3 in or . The polarisation states outside the square are CLEANed using the standard CLEAN method, while the thatched region and the region outside the square are CLEANed using the complex CLEAN method. b) When we rotate the () frame by an angle to view the constraints in the frame it is clear that the standard CLEAN cutoff is not the same for () as for .

7 Complex Clean applied to Linear Polarisation Images

Having demonstrated the rotational invariance of Generalised Complex CLEAN we now present a case study to highlight the advantages of the process over the traditional methods by examining the images which result from a standard versus Complex SDI CLEAN of the PKS B1637-771 data. The Stokes and CLEAN component images for the standard and Complex SDI CLEAN are shown in Figure 21 while the respective residual images are shown in Figure 22. As shown in the figures the Complex SDI clean finds more components associated with the diffuse structure in the source and almost no components outside the source. By comparison the standard method misses many components associated with the diffuse emission and inserts spurious components associated with noise. This is consistent with the differences in accessibility of components in the plane discussed in Section 6. Additionally, for the Generalised Complex CLEAN we see the residual values change smoothly from positive to negative through a range of values correctly representing the source structure, as opposed to effectively being single valued with a hard transition between the positive and negative domains. These hard transitions between domains result in the so-called canals (similar to ‘depolarisation’ canals’) seen in the residual image of the polarised intensity for the standard CLEAN shown in the left panel of Figure 23. Furthermore, as these features are a result of the rotational invariance the canals and domains change with rotations of linear polarisation before deconvolution. By comparison the residual image for the total linear polarisation derived from the Generalised Complex CLEAN shown in the right panel of Figure 23 is both smooth and shows less remaining flux. The final restored total linear polarisation images for both methods are shown in Figure 24, though it is difficult to tell in the greyscale image more diffuse polarised flux is recovered by the Generalised Complex CLEAN and less iterations are required.

Figure 21: Stokes (top row) and Stokes (bottom row) CLEAN component images using the standard (left column) and Complex (right column) SDI CLEAN methods. The Complex SDI CLEAN creates more components related to the smooth structure of the source, while the standard SDI CLEAN creates more CLEAN components associated with the noise.
Figure 22: Stokes (top row) and Stokes (bottom row) residual images made using the standard (left column) and Complex (right column) SDI CLEAN methods. The residuals made using the Complex SDI CLEAN are lower than the residuals made using the standard SDI CLEAN. Furthermore, the residuals made using the standard SDI CLEAN are smooth over the positive and negative domains, but discontinuous at the domain boundaries (seen as canals in the polarisation residual intensity). The residuals for the Complex SDI CLEAN are smooth across these boundaries and correctly preserve the overall structure of the source.
Figure 23: Linear polarisation residual intensity image of PKS B1637-771 using made the standard (left) and Complex (right) SDI CLEAN methods. The complex method shows a smooth residual, and the standard method shows canals within the residuals. These canals are nulls in the residuals of and , due to a change in sign.
Figure 24: Total linear polarisation intensity images of PKS B1637-771 made using the standard (left) and Complex (right) SDI CLEAN methods.

8 Distribution of Code

Python code for Complex Högbom and SDI CLEANs are available from the authors. Additionally, we have written a patch to allow a Complex SDI CLEAN task to run within MIRIAD. The patch is also available from the authors and we hope that CSIRO will consider including it in the standard MIRIAD release in the future. At present the method is designed to work for single pointing data, however, we have formulated an extension of the Complex SDI CLEAN for mosaiced data (see Appendix) which may be incorporated in MIRIAD in the future.

9 Conclusions

In this work, we have presented a new method of CLEANing polarimetric data. We show that our method, ‘Generalised Complex CLEAN’, provides a more physically consistent CLEAN model than the currently used methodologies, particularly for complex, diffuse sources. Generalised Complex CLEAN overcomes many of the problems with current CLEAN by properly accounting for the complex vector nature of polarisation. By looking for peaks in instead of separately in and , Generalised Complex CLEAN provides a method that is rotationally invariant and thus detects more true components at low sign-to-noise and fewer spurious ones and does not suffer from canals in the residual images, thereby recovering more flux. Furthermore, we show that by using Generalised Complex CLEAN the magnitude of the CLEAN components is more accurate than for standard CLEAN algorithms. By examining the relationship between the CLEAN components generated in a rotated and non-rotated frame we are able to measure the effects of numerical precision for different CLEAN methods. We find that the Complex SDI method is considerably better than a Complex Högbom CLEAN for preserving properties of components under rotations. This is believed to be the result of less numerical error in the SDI process. If this is the case, then other deconvolution methods which use fewer iterative subtractions may also be more suitable for large images. Both the Python code and a MIRIAD patch to run Generalised Complex CLEAN algorithms will be made available to the community.

Data from a number of current centimetre-wavelength instruments should benefit from the application of Generalised Complex CLEAN including the ATCA, SMA, ALMA, JVLA, GMRT and WRST. Furthermore, as we move to higher resolutions with instruments such as the the Square Kilometre Array which will generate more and more CLEAN components, it will be vital to consider the effects of cumulative numerical errors. In particular SKA1-MID, which will operate in the centimetre wavelength regime, is expected to generate images with of order 35,000 35,000 pixels requiring vast numbers of CLEAN components. Considering the results presented here it would seem prudent for future instruments to adopt a Complex SDI CLEAN in order to both have more physically meaningful results and reduce cumulative numerical errors. At meter-wavelengths instruments like SKA1-Low, which will produce wider and deeper images than ever before with of order 80,000 80,000 pixels (Dehghan et al., 2016), will certainly suffer from cumulative errors if deconvolution processes with large numbers of iterative steps are used. Such low frequency instruments may also benefit from deconvolution via a Complex SDI clean, but the issues of direction dependent calibration, and primary and synthesized beam variations across the large fields of view make this a more complex and challenging and case. Thus, more work will be needed in these areas to assess the advantages of a Generalised Complex CLEAN on widefield, low frequency data.

Acknowledgements

MJ-H and LP acknowledge support from the Marsden Fund administered by the Royal Society of New Zealand. We thank Dr Michiel Brentjens and Dr Andre Offringa for clarifying the CLEANING approach and possibilities in CASA and WSClean. We also acknowledge Dr Roman Klapaukh for input to preliminary discussions on the computation impact of the complex deconvolution. Finally we thank the referee and the editors for their constructive suggestions which helped expand the applicability of this work.

Appendix A Complex SDI CLEAN adapted for mosaics

Sault et al. (1996) explained that only minimal changes are needed to extend SDI CLEAN to operate on mosaics. Similarly, Complex SDI CLEAN can be extended to deconvolve polarisation mosaics. In this section, we describe the minimal changes needed to extend complex SDI CLEAN to mosaics.

In particular, for mosaicing the residual is weighted by the RMS image of the image when locating the peak. We adopt the same approach, and define the linear polarisation RMS image as , where is the RMS of the complex linear polarisation at , calculated using the variance

(21)

and is the mean in complex linear polarisation and is the total number of samples used to calculate the variance.

First the maximum of is found, then is constructed by clipping for values .

The final difference is the calculation of , the factor taking into account the scaling needed for the synthesized model . This is now calculated with weighting each image by the RMS

(22)

Then, one can scale the synthesized model by .

All other steps and calculations of the complex SDI method remain the same. 5

Footnotes

  1. Several older papers within the VLBI literature mention the task CXCLN (Cotton et al., 1984; Cotton, 1992), in classic AIPS, which uses a Högbom algorithm described as a ‘complex CLEAN’. This task is functionally different, since it was developed to account for asymmetric uv-coverage of Stokes and in VLBI observations when interferometers using circular feeds have only measured one of the circular correlations LR or RL (Cotton, 1992). The understanding of the VLBI community is that this approach is only needed when only half the correlations are measured, and that it is identical to the traditional approach in cases where all 4 correlations (RR, RL, LR and LL) are measured (Aaron, 1997), but as we will show here, this is incorrect. Note we have adopted the name ‘Generalised Complex CLEAN’ for our method to avoid confusion with CXCLN.
  2. This will be the case for all antennas unless the instrument has triaxial feeds such as ASKAP (Hotan et al., 2014) which keep the dipoles at a constant position relative to the source.
  3. CASA allows three options beyond the standard cleaning approach when using a Clark CLEAN only: i) searching for peaks in , ii) searching for peaks in and simultaneously, or iii) and total polarisation simultaneously. In all cases it then produces individual clean component images for each Stokes parameter. This latter approach is of most relevance here but this is designed to constrain peaks so as to select the most highly polarised components associated with a continuum (Stokes ) source first. Constraining components to have both a strong peak in and total polarisation is in fact worse for some science applications than just searching for peaks in e.g. when looking at diffuse polarisation in the Galactic Plane. In such cases the spatial scales of emission in Stokes and polarisation is different and the interferometer can resolve out the Stokes whilst still detecting polarised structures resulting in strong polarised signals with no matching continuum emission (see Gaensler et al., 2001; de Bruyn & Brentjens, 2005 for examples). For this reason LOFAR which runs AWImager only cleans in the Stokes parameters independently. The most flexible current CLEANing package is WSClean (Offringa et al., 2014) which will allow searches for peaks in the sum of squares of any sensible combination of polarisations, without requiring a Stokes counterpart. While WSClean can provide this functionality through its “jointpolarizations” option, there are some issues when also implementing multi-scale clean and cleaning jointly in and is not usually how the software is run.
  4. Here we have assumed that the and beams are the same. In some situations, such as for significant flagging, or for time switched polarization measurement e.g. quarter wave plates, this may not be true. In these cases there will be an additional bias to spatial scales where the coverage is mismatched. For point sources, the effect is unimportant and more generally speaking, as there is typically not much difference in the and beam, the overall effect will be small.
  5. The RMS image cannot be generated from and because . However, MIRIAD does not have the function for calculating , so may be used as an approximation.

References

  1. Aaron S., 1997, EVN Memo, 15
  2. Brentjens M. A., de Bruyn A. G., 2005, A&A, 441, 1217
  3. Carrillo R. E., McEwen J. D., Wiaux Y., 2014, MNRAS, 439, 3591
  4. Clark B. G., 1980, A&A, 89, 377
  5. Cornwell T. J., 1988, A&A, 202, 316
  6. Cornwell T. J., 2008, IEEE Journal of Selected Topics in Signal Processing, 2, 793
  7. Cornwell T. J., Evans K. F., 1985, A&A, 143, 77
  8. Cotton W. D., 1992, AIPS Memo, 79
  9. Cotton W. D., Geldzahler B. J., Marcaide J. M., Shapiro I. I., Sanroma M., Rius A., 1984, ApJ, 286, 503
  10. de Bruyn A. G., Brentjens M. A., 2005, A&A, 441, 931
  11. Dehghan S., Johnston-Hollitt M., Hollitt C., 2016, ArXiv e-prints: 1601.00067
  12. Gaensler B. M., Dickey J. M., McClure-Griffiths N. M., Green A. J., Wieringa M. H., Haynes R. F., 2001, ApJ, 549, 959
  13. Garsden H., et al., 2015, A&A, 575, A90
  14. Högbom J. A., 1974, A&AS, 15, 417
  15. Hotan A. W., et al., 2014, PASA, 31, e041
  16. Ikeda M., Nishiyama K., Ohishi M., Tatematsu K., 2001, in Harnden Jr. F. R., Primini F. A., Payne H. E., eds, Astronomical Data Analysis Software and Systems X Vol. 238 of Astronomical Society of the Pacific Conference Series, Development of Radio Astronomical Data Reduction Software NEWSTAR. p. 522
  17. Johnston-Hollitt M., et al., 2015, Advancing Astrophysics with the Square Kilometre Array (AASKA14), p. 92
  18. Li F., Cornwell T. J., de Hoog F., 2011, A&A, 528, A31
  19. McMullin J. P., Waters B., Schiebel D., Young W., Golap K., 2007, in Shaw R. A., Hill F., Bell D. J., eds, Astronomical Data Analysis Software and Systems XVI Vol. 376 of Astronomical Society of the Pacific Conference Series, CASA Architecture and Applications. p. 127
  20. Offringa A. R., et al., 2014, MNRAS, 444, 606
  21. Pratley L., Johnston-Hollitt M., Dehghan S., Sun M., 2013, MNRAS, 432, 243
  22. Reynolds J., 1994, ATNF Technical Memo Series, p. 39.3/040
  23. Ryle M., Hewish A., 1960, MNRAS, 120, 220
  24. Sadler E. M., Ekers R. D., Mahony E. K., Mauch T., Murphy T., 2014, MNRAS, 438, 796
  25. Sault R. J., Staveley-Smith L., Brouw W. N., 1996, A&AS, 120, 375
  26. Sault R. J., Teuben P. J., Wright M. C. H., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds, Astronomical Data Analysis Software and Systems IV Vol. 77 of Astronomical Society of the Pacific Conference Series, A Retrospective View of MIRIAD. p. 433
  27. Schwab F. R., 1984, AJ, 89, 1076
  28. Schwarz U. J., 1978, A&A, 65, 345
  29. Steer D. G., Dewdney P. E., Ito M. R., 1984, A&A, 137, 159
  30. Tingay S. J., et al., 2013, PASA, 30, 7
  31. van Haarlem M. P., et al., 2013, A&A, 556, A2
103230
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
Edit
-  
Unpublish
""
The feedback cannot be empty
Submit
Cancel
Comments 0
""
The feedback cannot be empty
   
Add comment
Cancel

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.