Richardson–Lucy Deconvolution with a Spatially Variant Point-spread Function of Chandra: Supernova Remnant Cassiopeia A as an Example

The RL algorithm iteratively estimates a true image from an observed image using Bayesian inference. It generally assumes that the PSF does not change with position in the image. The RL algorithm is expressed by

$$\begin{eqnarray}&&{W}_{i}^{(r+1)}={W}_{i}^{(r)}\displaystyle \sum _{k}\displaystyle \frac{{P}_{{ik}}{H}_{k}}{{\sum }_{j}{P}_{{jk}}{W}_{j}^{(r)}},\end{eqnarray} \tag{ 1 }$$

where i and j are mapping the image in the sky, and k is mapping the image on the detector. The indices of the summation run through all the pixels. W^(r) is the restored image after r iterations, and H is the observed image on the ACIS detector. P_jk is the probability that a photon emitted in sky W bin j is measured in data space H bin k, or P(H_k ∣W_j ).

Previous Chandra image deconvolution approaches (e.g., Thimmappa et al. 2020; Sobolenko et al. 2022) used a simplified approximation for the P_jk values, i.e., they used the same PSF for each j bin. Here we assume that the PSF changes as a function of the off-axis angle and the roll angle. As a consequence, the Chandra RL algorithm is extended. The formula for RL_sv is obtained by rewriting Equation (1) as

$$\begin{eqnarray}&&{W}_{i}^{(r+1)}={W}_{i}^{(r)}\displaystyle \sum _{k}\displaystyle \frac{{P}_{{iik}}{H}_{k}}{{\sum }_{j}{P}_{{jjk}}{W}_{j}^{(r)}}.\end{eqnarray} \tag{ 2 }$$

P_jjk refers to a PSF at a position j (first index) which returns a probability that an event emitted at W_j (second index) is observed at H_k (third index), or P_j (H_k ∣W_j ). Computational cost and memory requirements need to be minimized for calculating the third-order tensor of the PSF, which is a distinctive feature of the RL_sv algorithm. When H is corrected for slight differences among pixels in effective area and exposures, its normalization can be chosen arbitrarily. Here we use H_k = N_k /A_k , where N_k is the detector count image, and A_k is the Hadamard product of effective area and exposure time.

There are regularization techniques to enhance the RL method, which are also readily available for the RL_sv method. Among these, total variation (TV) regularization (Rudin et al. 1992) is effective in handling statistical errors, which is used in the RL method (Dey et al. 2006). The formula for RL_sv with the regularization is obtained by rewriting Equation (2) as

$$\begin{eqnarray}&&{W}_{i}^{(r+1)}=\displaystyle \frac{{W}_{i}^{(r)}}{1-{\lambda }_{\mathrm{TV}}{\rm{div}}\left(\tfrac{{\rm{\nabla }}{W}_{i}^{(r)}}{| {\rm{\nabla }}{W}_{i}^{(r)}| }\right)}\displaystyle \sum _{k}\displaystyle \frac{{P}_{{iik}}{H}_{k}}{{\sum }_{j}{P}_{{jjk}}{W}_{j}^{(r)}}.\end{eqnarray} \tag{ 3 }$$

The only difference from the RL_sv algorithm of Equation (2) is the regularization term of $1-{\lambda }_{\mathrm{TV}}{\rm{div}}({\rm{\nabla }}{W}_{i}^{(r)}/| {\rm{\nabla }}{W}_{i}^{(r)}| )$, where λ_TV is the regularization parameter, ${\rm{div}}(\cdot )$ is the divergence, and ${\rm{\nabla }}{W}_{i}^{(r)}$ is the gradient of ${W}_{i}^{(r)}$. In this paper, we utilize the parameter λ_TV = 0.002, as proposed by Dey et al. (2006).

Deconvolution methods require an understanding of their applicability to a practical condition, as well as optimization of computation cost and accuracy (for features of various methods see Naik & Sahu 2013). The RL method is well studied and has been used to incorporate regularization (e.g., van Kempen & van Vliet 2000; Dey et al. 2006; Yuan et al. 2008; Yongpan et al. 2010) and the recent trend of deep learning (e.g., Agarwal et al. 2020). For Chandra users, RL deconvolution with a single PSF is frequently used because the method is already implemented as arestore in the Chandra Interactive Analysis of Observations (CIAO; Fruscione et al. 2006), Chandra's standard data processing package.

Compared to other methods, the RL method forces the deconvolved image of each iteration to be non-negative, and its integral value is conserved. Additionally, the method converges to the maximum likelihood solution for a Poisson noise distribution (Shepp & Vardi 1982), which is suitable for Chandra images with noise from counting statistics. Depending on the application, it is less prone to ringing artifacts than inverse PSF-based methods (e.g., Sekko et al. 1999; Neelamani et al. 2004); see the results of the comparison by Dalitz et al. (2015). According to White (1994), it is robust against small errors in the PSF.

Because Cas A is a bright and diffuse X-ray source with a moderately large apparent diameter, it is an ideal target to demonstrate the RL_sv method. It has been observed by Chandra almost every year since 1999. The Chandra data of Cas A used in this paper are listed in Table 1: ACIS-S observation of 2004 in Obs. ID = 4636, 4637, 4639, and 5319. The image size is 1489 × 1488 pixels, or 743'' × 742'' given a unit pixel of 0492. Data processing and analysis were performed using CIAO version 4.13. The data were reprocessed from the level 1 event files by chandra_repro. Since the roll angle and optical axis of the four observations are almost the same (the maximum difference of the optical axis location is about 4 unit pixels), all the events were merged into one by merge_obs. The total exposure time was 428.29 ks.

The Chandra telescope system consists of four pairs of nested reflecting surfaces, configured in the Wolter type I geometry. The high-energy response is achieved by coating the mirrors with iridium. It has attained the highest angular resolution of 0492 among existing X-ray telescopes. Its mirror has been extensively calibrated on the ground and in orbit (Jerius et al. 2000). The Chandra PSF is positional-dependent, mainly due to aberrations. The RL_sv method includes the position dependence of the PSF, which is useful for highly extended X-ray sources.

Because creating a PSF for each position is computationally expensive, it is decimated at some intervals. For reference, creating a PSF takes several seconds, depending on the computational environment and desired accuracy. The sampling interval of the PSFs was chosen to be 35 × 35 pixels (total of 43 × 43 = 1849). The interval was determined empirically by trying several different ranges. In general, the PSFs simulated from each observation should be merged for a precise calculation. Here, the PSF of the Obs. ID of 4636 was used as a representative since its sampling is decimated. The PSFs at the lattice points were generated by CIAO's simulate_psf using the Model of AXAF Response to X-rays (Wise et al. 1997; Davis et al. 2012) at a monochromatic energy of 2.3 keV. They were applied to the observed image with energies from 0.5 to 7.0 keV.

Figure 1 shows all the PSFs sampled every 35 × 35 pixels. The optical axis is located at the northeast in the image, where the spread of the PSF is minimum. As the position is away from the optical axis, the tail of the PSF increases with its gradual shift of the elliptical axis. Although it is a trade-off with photon statistics, it is effective to run the RL_sv method with the optimal monoenergetic PSF for each of the multiple energy decompositions (see Figure 2). This is because, at shorter wavelengths, the effect of diffuse reflection due to the roughness of the mirror surface is not negligible (Jerius et al. 2000).

Zoom In Zoom Out Reset image size

Figure 3(a) is an observed ∼400 ks image using the energy range from 0.5 to 7.0 keV, as explained in Section 3.1. We applied the RL_sv method to the image with a sampling interval of each PSF of 35 × 35 pixels. The number of iterations is 200. Note that the choice of iteration number is discussed in Section 4.1. The result of the RL_sv method is presented in Figure 3(b). The unit of flux and its range in Figure 3(b) is the same as in Figure 3(a). The overall structures in the RL_sv image become more vivid than the original ones. The images around the off-axis are significantly improved compared to those around the optical axis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To make the differences more precise, we present magnified images of the original image in Figures 3(a-1, -2, -3) and the RL_sv-image in Figures 3(b-1, -2, -3). The three regions represent a sharp filament in the northeast, complicated filaments in the north, and a slightly diffuse area in the south. The filamentary structures in the northeast and north become sharper in the RL_sv-image. We will quantify the filament width in detail and discuss the systematic uncertainties associated with the method in Section 4.2.

We considered a way of assessing an appropriate number of iterations, which is one of the issues with the RL method. This is because the method has the property of excessive amplification of noise as the number of iterations increases. We propose a method to suppress convergence using statistical errors during the iteration. The formula of the method is written as

$$\begin{eqnarray}&&{W}_{i}^{(r+1)}={W}_{i}^{(r)}\displaystyle \sum _{k}\displaystyle \frac{{P}_{{iik}}G({N}_{k})/{A}_{k}}{{\sum }_{j}{P}_{{jjk}}{W}_{j}^{(r)}}.\end{eqnarray} \tag{ 4 }$$

The only difference from the RL_sv algorithm, Equation (2), is the G(N_k )/A_k term. N (counts) is the map of detector counts. A (photons cm⁻² s⁻¹) is the Hadamard product of effective area and exposure time. G(N_k ) is a random number generator following a Poisson distribution with a count in the kth pixel of N_k ; i.e., G(N_k )/A_k is a flux in units of photons cm⁻² s⁻¹. The reason for normalizing G(N) by dividing it with A is to account for the slight variations in effective area and exposure time among pixels.

The performance of the RL_sv algorithm using Equation (4) is compared to that using Equation (2). The convergence is evaluated using the mean squared error (MSE) of the two images: one step before and after iterations. Figure 4 shows the history of the MSE during the iteration. The curve obtained by the RL_sv algorithm, Equation (4), saturates at a certain level, while the other continues to decrease. The saturation level is caused by the injection of Poisson fluctuation at each step, which is considered an indicator of stopping. In Figure 4, the iteration number of ∼30 seems appropriate.

Zoom In Zoom Out Reset image size

We then designed a simplified method for evaluating a certain amount of confidence. The RL_sv-deconvolved image should have a similar amount of fluctuation accompanying the observation image. The principle of the method is to propagate errors of the converged RL_sv image into the next step. Our choice of using the last step for the error propagation is just to simplify the task. Here, each error in the observed image is considered statistically independent. Assuming only uncertainties on H_k , using the law of error propagation, the image uncertainty can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}{\sigma }_{{W}_{i}^{{\prime} }} & = & \sqrt{\displaystyle \sum _{k}{\left[\displaystyle \frac{\partial }{\partial {H}_{k}}\left({W}_{i}\displaystyle \sum _{k}\displaystyle \frac{{P}_{{iik}}{H}_{k}}{\displaystyle \sum _{j}{P}_{{jjk}}{W}_{j}}\right){\sigma }_{{H}_{k}}\right]}^{2}}\\ & = & {W}_{i}\sqrt{\displaystyle \sum _{k}{\left(\displaystyle \frac{{P}_{{iik}}}{\displaystyle \sum _{j}{P}_{{jjk}}{W}_{j}}\displaystyle \frac{\sqrt{{N}_{k}}}{{A}_{k}}\right)}^{2}},\end{array}\end{eqnarray} \tag{ 5 }$$

where ${W}^{{\prime} }$ is the image of the next iteration number of any estimated true image of W, and $\sqrt{{N}_{k}}$ is the statistical error of N_k .

We compared the off-axis RL image with the error of Equation (5) to an on-axis observation in Figure 5. Figures 5(a) and (b) are the off-axis southeastern images from Figures 3(a) and (b), respectively. Figure 5(c) shows a southeastern on-axis image of Obs. ID = 5196. The exposure time for the on-axis observation is ∼50 ks, resulting in a larger statistical error compared to the off-axis observation of ∼400 ks. The fan-shaped regions in Figures 5(a)–(c), along the filament, are chosen to create the radial profiles, which is the one-dimensional profile of the photons in each region extending from the central compact object (CCO) of Cas A toward the outer regions. These radial profiles were created using dmextract in CIAO. In Figure 5(d), the framed regions from Figures 5(a)–(c) are color-coded as blue, orange, and black, respectively. The error bars of the radial profiles in Figure 5(d) correspond to statistical errors, represented by blue and black, and the result obtained by applying Equation (5) to the 199th iteration of the RL_sv image, is indicated by orange. From Figure 5(d), the profile of the off-axis RL_sv image agreed with that of the on-axis image within the statistical errors. This method gives a guideline for a certain level of confidence associated with the RL_sv method.

Zoom In Zoom Out Reset image size

In this section, further enhancements to the RL_sv method are discussed. The first is to reduce the loss of down-sampling PSFs. The positional dependence of the PSF does not contain high-frequency components, so the decimation of PSF sampling should work to some extent. For a small image such as a core plus jet structure in an active galactic nucleus, keeping a high sampling rate of the PSFs might be possible. However, for a largely extended source such as a supernova remnant or galaxy cluster, to minimize the sampling rate is critically important for practical use. The higher the decimation, the more emphasized the boundary of the segment. Taking Cas A as an example, the edges of specific segments clearly appear when the sampling interval is 35 × 35 pixels. To smooth out the edges, we propose that the PSFs' boundaries be randomly selected from nearby PSFs (see more details in the Appendix).

Second, this method can be developed by incorporating several regularization methods. We implemented an RL_sv method incorporating the TV regularization expressed in Equation (3). Finally, the RL_sv method is naturally applied to color images. By decomposing observed images into several colors (or energy bands) and generating PSFs for an appropriate energy in each band, an energy-dependent RL_sv method can be realized.

We compare the RL_sv method with and without the TV regularization. We use the same image as in Section 3.1. The PSF sampling is 35 × 35 pixels and the number of iterations is 200. The PSFs' boundaries are randomly selected following the Appendix. Figure 6 presents the enlarged eastern image after applying the RL_sv method to the entire region. Figures 6(a) and (b) show the results of the RL_sv method of Equation (2) and the regularization version of Equation (3), respectively. The TV regularization preserves the sharp structure and smoothes out statistical errors. In this way, regularization can be added to the RL_sv method.

Zoom In Zoom Out Reset image size

We implement the RL_sv method including these enhancements and adapt it to the Cas A observational data described in Section 3.1. Figure 2(a) is the observed image of Section 3.3 divided into three energies in RGB: 0.2–1.2 keV (red), 1.2–2.0 keV (green), and 2.0–7.0 keV (blue). Cas A is dominated by the thermal radiation in ≤4 keV and the nonthermal radiation in ≥4 keV. We applied the RL_sv method with the TV regularization Equation (3) to each energy image of Figure 2(a) using the appropriate energy of the PSF (red is 0.92 keV, green is 1.56 keV, and blue is 3.8 keV) based on the official CIAO page. ⁴ The sampling interval of PSFs is 35 × 35 pixels. The PSF is randomly selected at the sampling boundaries according to the Appendix. The number of iterations is 30, according to Section 4.2. The result of the RL_sv method is presented in Figure 2(b). The energy dependence in Figures 2(b-1) and (b-2) are clearly visible by this method.

We propose two complementary ways to obtain a guideline on the stop condition of the RL_sv method. One is to obtain a minimum of residuals by inserting statistical uncertainties into each update during the iteration process. This gives a rough estimate of the limit of the iterations. The other is to include the statistical uncertainties in the last step of the iteration. By combining the two methods, it is possible to derive uncertainties using the errors obtained by the latter method at the optimal iteration number estimated by the former. This is a quick and convenient way to derive the systematic uncertainties associated with the RL_sv method.

The systematic uncertainty in the former method is a way of defining the optimal number of iterations. One easy way is to use the same level of residuals as shown in Section 4.1. It is intrinsically difficult to distinguish the signal of the celestial objects from the statistical noise. This difficulty needs to be overcome by comparing the deconvolved images without errors and the error-estimated images around the optimal iteration number recommended by the former method. This method is based on a compromise between the computational cost and the simplicity of use while keeping a reasonable statistical error.