Image Subtraction in Fourier Space

As the core engine of transient detection, image subtraction is often the most computation-intensive individual task in the data processing pipeline. The ongoing trend of very massive data flow from time-domain surveys is making real-time data reduction increasingly challenging. It motivated us to develop a new tool to relieve the computational bottleneck while reconciling the subtraction performance and the computing cost.

The image subtraction algorithms emanated from Alard & Lupton (1998) have been broadly applied in astronomical transient surveys. Our proposed method SFFT is a new member of this category. The main purpose of the approach aims to perform PSF-matching via image convolution with pixelized kernels formulated as linear combinations of a pre-defined set of basis functions. The fact that astronomical images from wide-field surveys generally possess non-constant PSF across the entire field of view leads to spatial variations of the convolution kernels. Moreover, the spatial varying photometric scale also needs to be taken into account in order to match the images from wide-field surveys.

There have been some attempts to accelerate the calculations involved in image subtraction, either aimed at faster kernel determination or speed-up of the subtraction afterward. For image subtraction, constructing the least-squares matrix to solve the convolution kernels is computationally expensive, especially for spatial variant kernels. Alard (2000) proposed to fit the space-varying kernel on a set of sub-areas of the field, with a simple assumption that the kernel spatial variation within each sub-area can be negligible. The useful strategy has been adopted in the software HOTPANTS, a widely-used implementation of Alard & Lupton (1998). Apart from the efficient simplification in computing, the flexibility of using sub-areas has some additional benefits in practice. One may pre-select an optimal set of sub-areas to exclude the sources in observational data that cannot be perfectly modeled by the image subtraction algorithm. Doing so can also avoid the fitting being strongly misled by specific regions of the field (e.g., the brightest and densest regions).

Very recently, Hitchcock et al. (2021) provided an innovative way (i.e., PyTorchDIA) to bypass the construction of the analytical least-squares matrix. PyTorchDIA finds the kernel solution by optimizing a loss function using the automatic differentiation in PyTorch, which brings a considerable gain in computing cost compared with Bramich (2008). However, PyTorchDIA is incompatible with kernel spatial variation and can only solve a constant kernel for image subtraction. In addition to the efforts on kernel determination, Hartung et al. (2012) successfully implemented fast image convolution with spatial-varying kernel on GPUs, which helps accelerate the subtraction after the convolution kernel has been computed.

The strategy of SFFT is different from the methods mentioned above: we present the least-squares question in Fourier space. By Parseval’s theorem, a least-squares minimization in real space is equivalent to a least-squares minimization in Fourier’s domain. Given that $\delta$-functions have higher flexibility and can maintain simple forms after Fourier transform, here we adopt the DBFs as the kernel basis functions following Bramich (2008) and Miller et al. (2008). Finally, SFFT allows for spatial variations across the field, whether from PSF, photometric scaling, or sky background. The following section will show how the calculations for image subtraction can be reduced to be FFTs and element-wise matrix operations.

Given a reference image $R$ and a science image $S$, PSF-homogenization is carried out by convolving $R$ or $S$ with a spatially-varying kernel $K_{x,y}$. Considering that the sky background between the two epochs may change, we introduce an additional term $B$ to model the sky background difference between the two images. The image subtraction problem can be written in real space as the minimization of the difference image $D$, defined as

where the input images $R$ and $S$ are image with dimensions $(N_{0},N_{1})$, and $x$ and $y$ are indices in the ranges of $[0,N_{0}-1]$ and $[0,N_{1}-1]$, respectively. Note that the kernel $K_{x,y}$ is spatially varying, so the integral in the equation above, strictly speaking, is not a convolution, such an integral is denoted by $R\circledast K$. The differential background map $B$ is modeled as a polynomial form following Alard & Lupton (1998), that is,

where the polynomial power indices $p$ and $q$ are in the range of $[0,D_{B}]$ and $[0,D_{B}-p]$, respectively. We follow Miller et al. (2008); Hartung et al. (2012) to decompose the kernel $K_{x,y}$ into the “shape-free” $\delta$-function basis $\mathcal{K}$. The kernel dimension is assumed to be $(L_{0},L_{1})$, with $L_{0}$ and $L_{1}$ being odd integers, such that $L_{0}=2w_{0}+1$ and $L_{1}=2w_{1}+1$. We assume

and,

where $\mathcal{D}_{\alpha\beta}(u,v)=\vec{\delta}(u-\alpha,v-\beta)$ is the standard Cartesian orthonormal basis, which consists of $\delta$ functions at each kernel pixel, $\alpha$ and $\beta$ are indices in the range of $[-w_{0},w_{0}+1)$ and $[-w_{1},w_{1}+1)$, respectively, u and v are kernel coordinate indices in the range of $[-w_{0},w_{0}+1)$ and $[-w_{1},w_{1}+1)$, respectively, and $\vec{\delta}$ is a binary function on integers: $\vec{\delta}(\rho,\epsilon)=1$ if $\rho=\epsilon=0$, otherwise $\vec{\delta}(\rho,\epsilon)=0$ with $\rho$ and $\epsilon$ being any integers. With such a construction, the sum of the convolution kernel is uniquely determined by the coefficient of the basis vector $\mathcal{K}_{00}$, which simplifies the way we control the photometric scaling through convolution.

The spatial variation of the kernel can be fitted by the sum of the kernel function given above, modulated by a two-dimensional polynomial function to account for their varying contributions across the image field

where the polynomial power indices $i$ and $j$ are in the range of $[0,D_{K}]$ and $[0,D_{K}-i]$, respectively. The spatial variation information of the kernel is encoded in the coefficients $\mathring{a}_{ij\alpha\beta}$, which will be eventually solved by the subsequent linear system during the minimization of the residuals. More specifically, once $\mathring{a}_{ij\alpha\beta}$ is known, we can calculate $\mathring{A}_{xy\alpha\beta}$ for any given pixel coordinate $(x,y)$ following Equation (6), and further derive the certain kernel associated with the pixel via the expansion on $\delta$-function basis by Equation (5).

Equations (2), (5) and (6) can be replaced by more complex functions. The code we have developed includes an option for using B-splines to model the PSF and background level variations. For the B-splines cases, the space-varying kernel $K_{x,y}$ and differential background B(x, y) are modeled as

and

where $B_{p;k,t}$ ($B_{p;k^{\prime},t^{\prime}}$) are B-spline basis functions of degree k ($k^{\prime}$) and knots t ($t^{\prime}$). For simplicity, in this paper we focus only on the performance with polynomial models.

The pixel-to-pixel flux variations of the two images are accounted for by the coefficients $\mathring{A}_{xy00}=\sum_{ij}\mathring{a}_{ij00}x^{i}y^{j}$. If we consider the flux level of the image-pair to be well-calibrated, the constant flux scaling between images requires a constant kernel integral, that is, $\mathring{A}_{xy00}=\mathring{a}_{0000}$. Note that a constant flux scaling was first presented in Alard & Lupton (1998). Having a constant photometric ratio across the entire field is optional in our program. Like Bramich et al. (2013), SFFT allows for space-varying flux scaling with a polynomial form to accommodate the effect of imperfect flat-field correction or cirrus cloud attenuation. Note that the current SFFT does not disentangle the polynomial degrees of spatial variations accounting for the convolution kernel and the photometric sensitivity.

With abbreviation $T^{\rho\epsilon}(x,y)=x^{\rho}y^{\epsilon}$ and $R^{\rho\epsilon}=T^{\rho\epsilon}R$ ,where $\rho$ and $\epsilon$ are any integers, Equation (1) is rewritten as,

The convolution kernel is typically very small in size, at such scale its spatial variation is expected to be negligible locally. An approximation (see Appendix A) can be made by moving the polynomial functions describing the spatial variation outside the integral in Equation (1). This leads to

and,

Note the notation $\circ$ in Equation $(\ref{eqn:sfft_eq12})$ indicates circular convolution and $\mathcal{K}_{\alpha\beta}$ can be estimated to generate the convoluted source image that best matches the PSF of the reference image across the entire field.

The $\delta$-function basis we have adopted for $\mathcal{K}$ as shown in equations (3) and (4) allows for a simple Fourier space representation of the image difference procedure. In Fourier space Equation (12) becomes

where the symbols with a hat denote the Fourier transform of the images, $a_{ij\alpha\beta}=N_{0}N_{1}\mathring{a}_{ij\alpha\beta}$ with $N_{0}$ and $N_{1}$ being the dimensional of the images in $x$ and $y$, respectively. The Discrete Fourier Transform (DFT) of the basis function has the simple form

here $W_{N_{0}}=e^{-i2\pi/N_{0}}$ and $W_{N_{1}}=e^{-i2\pi/N_{1}}$, with $i$ being the unitary imaginary number.

Now define $G=\widehat{D}\widehat{D}^{*}$, where ^∗ stands for complex conjugate, we find,

Our goal of optimal subtraction is to minimize the power of the difference image in Fourier space, as given in Equation (15). Let $F=\sum_{l,m}G(l,m)$, the minimization in Fourier space for image subtraction has a simple least-square solution given by $\nabla F=0$ where the gradients are calculated with respect to all the parameters for image matching,

Equations $(\ref{eqn:sfft_eq16})$ and $(\ref{eqn:sfft_eq17})$ form a linear system with the array elements shown in Equation $(\ref{eqn:sfft_eq18})$.

The elements of Equation (18) are given explicitly as the following,

These equations can be further simplified as shown in Appendix B. Finally, the difference image $D$ is calculated from Equation (13) with applying the linear system solution $\{...,a_{ijab},...,b_{pq},...\}$ of Equation (18). Given that the difference image will possess correlated noise on the background, we present a decorrelation process to whiten the background noise in Appendix C.

Applying SFFT is straightforward for signal-dominated crowded fields. Mostly, saturation is the predominant factor in all distraction-sources for crowded fields. The crowded-flavor SFFT will automatically mask the pixels affected by the saturated sources in the field and replace them with the neighboring background. To further reduce image matching errors, it is allowed to use a more elaborate customized mask instead, e.g., to include the pixels affected by known variables and transients in the field.

On the flip side, the sparse-flavor SFFT is relatively more sophisticated on image masking. We select a set of sources of astronomical significance in the field and mask all other irrelevant regions so that the selected objects can dominate the parameters-solving process. It is reminiscent of the similar strategy in HOTPANTS, which fits on the rectangle stamps (i.e., “sub-stamps” in HOTPANTS terminology) encompassing selected astronomical objects. The analogous intention of both methods is to restrict the calculation on a properly pre-selected set of sources to eliminate the effect of distraction-sources on the solution of convolution kernels. In this scheme, the regions masked by SFFT have also included abundant background pixels. This is a reasonable consideration in our framework. Pixel uncertainty is not considered as weights in SFFT subtraction. Accordingly, the overwhelming noise-dominated background pixels are more likely to degrade the construction of accurate PSF-matching kernels rather than to contribute (see the same consideration in Kochanski et al. (1996)).

In practice, the input image-pair of sparse-flavor SFFT is required to be sky subtracted. This requirement is to simplify the image-masking process so that all the pixels enclosed in masked regions can be replaced by a constant of zero. As a result, the differential background term in the SFFT algorithm becomes trivial as the background offset between the input image-pair has been minimized by sky subtraction. That is to say, we pass on the function of matching differential background, which was embedded in numerical calculations of SFFT, to a customized sky subtraction as a preliminary process. Fortunately, modeling sky background is usually feasible in sparse fields (e.g., using an interpolation-based method). Note this requirement is not a prerequisite for crowded-flavor SFFT. This is because properly modeling sky background can be tricky for crowded fields: the modeled sky is more susceptible to being biased by the signal harboring in the pixels misidentified as background (e.g., the outskirts of nearby galaxies, see Akhlaghi & Ichikawa (2015)).

We developed a morphological classifier to carry out the source selection in sparse-flavor SFFT. The classifier was initially proposed for PSFEx (Bertin, 2011), which enables PSFEx to select a subset of point sources for constructing the PSF model. Given a photometry catalog generated by SEXTRACTOR (Bertin & Arnouts, 1996), one can draw a figure of instrumental magnitudes (SEXTRACTOR catalog value MAG_AUTO) against flux radius (SEXTRACTOR catalog value FLUX_RADIUS) for all detected photometry objects. Generally, bright point sources tend to stay around a nearly vertical straight line on this figure. PSFEx leverages the statistical feature to select appropriate samples to establish the PSF model. Although the goal of source selection for sparse-flavor SFFT is not entirely aligned to that for PSFEx, it inspires us to make the selection based on SEXTRACTOR parameters and their statistical characters.

To demonstrate the source selection criterion in sparse-flavor SFFT, we show an example of an individual DECam (Honscheid & DePoy, 2008) image in Figure 1. By contrast, the SEXTRACTOR photometry catalog of this image is cross-matched with the Tractor catalog from Legacy Survey (Blum et al., 2016), which offers their own fitted morphological types. As shown in Figure 1, the point sources and extended sources form two conspicuous branches, which intersect at the faint end but become well-separated towards the bright side. According to the Tractor types, the branch with a nearly vertical orientation primarily comprises point-like sources with Tractor type PSF, while the extended sources with Tractor types DEV, EXP, REX are likely to be found in the other branch that sprawls out horizontally. Moreover, the discrete gray dots represent the objects that cannot be found in the Tractor catalog, which are casual detections such as cosmic rays. SFFT uses Hough Transformation (Hough, 1959) to identify the straight line surrounded by the branch of point sources. The figure is first pixelized by counting the number of objects in each small grid. Hough Transformation is applied so that we can search for the strongest line feature with vertical orientation. The thin belt-like region around the detected straight line with a fixed width is referred to as point-source-belt. In SFFT, the parameter -BeltHW describes the belt half-width (default value is 0.2). With this terminology, it is easier to describe the specific selection criteria.

For the source selection, sparse-flavor SFFT will first generate a SEXTRACTOR photometry catalog for reference image and science image, respectively. A source is selected if it lies in point-source-belt or out of point-source-belt from the right side for both images. The extended sources in the field are not excluded, as SFFT can use both stars and galaxies to derive the solution for the matching kernel. However, typical cosmic rays and the faintest astronomical sources will be discarded. Given that each pixel of selected samples is equally weighted no matter bright or faint in SFFT, rejecting the noisiest subset should be appropriate in our framework. Furthermore, SEXTRACTOR has been configured to only output the objects with SEXTRACTOR catalog value FLAG being zero. Namely, this guarantees that the saturated sources and blending objects (relatively rare in sparse fields) are not in the selection.

Significant variable sources in the field should be further rejected. SFFT cross-matches the photometry catalogs of the input image-pair, then calculates the difference in instrumental magnitude for each matched object. The outliers of the distribution of magnitude difference indicating violent brightness change will be discarded from the initial set of selection. In SFFT, the outliers are identified by the threshold parameter -MAGD$\_$THRESH (default value is 0.12 mag). The survived selection set is applied for the image-masking in sparse-flavor SFFT. The outlier rejection is usually sufficient enough for a successful image subtraction, but one is also allowed to further refine the automatic selection by removing the known but not recognized variables and transients. Unlike HOTPANTS, SFFT does not enclose each source of the selection within a fixed-size rectangle stamp. Instead, we employ the detection mask of SEXTRACTOR (check-image SEGMENTATION generated by SEXTRACTOR) to define the pixel domain of each selected source. Finally, the pixels not in any pixel domain will be masked by a constant of zero.

The next step is to perform subtraction with the unmasked and the masked image-pairs, in the same way for the two flavors. We implemented the subtraction algorithm in CUDA following the mathematical reasoning in Section 2.2.

There are a few free parameters in this process. In Section 2.2, the spatial variation of the convolution kernel was modeled as a polynomial form of degree $D_{K}$. Similarly, the differential background is assumed to be a polynomial with degree $D_{B}$. In SFFT, the two free parameters correspond to -KerPolyOrder and -BGPolyOrder, respectively. The default values of the two parameters are 2. Another parameter that should be specified is the pixel size of the convolution kernels. The default way is to determine kernel size as seeing-related, i.e., the kernel half-width is calculated as $\textbf{\text{-KerHWRatio}}\times\text{FWHM}_{\text{L}}$, where the ratio -KerHWRatio is a free parameter of SFFT and $\text{FWHM}_{\text{L}}$ is the worse seeing of the input image-pair. The last free parameter -ConstPhotRatio is boolean, which controls whether SFFT subtraction is subject to a constant photometric ratio or not. The default value of this parameter is True that assumes the input image-pair has been well-calibrated. Setting the parameter to be False means the photometric ratio will also become a polynomial form of degree -KerPolyOrder.

By design, SFFT does not require isolated sources to solve the image matching. In crowded fields where few or no isolated stars are present, SFFT can perform exceptionally well. The ZTF test images are downloaded from ZTF Data Release 3 ²²2https://www.ztf.caltech.edu/page/dr3, covering the nearby galaxy M31. Figure 2 and Figure 3 show the subtraction performance of both the SFFT and ZOGY method, with a wider view and a close-up view on the galaxy core, respectively.

The difference image generated by SFFT is much cleaner than the ZOGY approach, as indicated by the lower panels of the two figures. One can also notice that the difference image by the ZOGY method appears to have a bias seemingly due to the photometric mismatch of image subtraction. This guess is confirmed in Figure 4 that analyses the flux residuals of five field stars near the galaxy core. ZOGY subtraction uniformly remains negative net flux for these samples, suggesting an improperly-estimated photometric ratio between reference and science image. By contrast, SFFT does not have the same problem. Unlike ZOGY, the photometric ratio (a constant in our default configuration, i.e., $\mathring{a}_{0000}$) in SFFT is straightforwardly solved from the linear system, rather than an extra value calculated by some other photometric process.

Figure 5 and Figure 6 show SFFT subtraction performance on LMC observations of the Antarctica telescope AST3-II and test data from space-based TESS 30min-cadence Full Frame Images (FFIs), respectively. A known variable star of RRAB type is easily identified at the image center for both cases. It illustrates how the flux variability can be effectively captured by SFFT subtraction from the densely packed fields.

It is also essential to investigate the performance of SFFT subtraction on the ordinary sparse sky fields that are more common in extragalactic transient surveys. For comparison, we additionally employed the widely-used image subtraction software HOTPANTS as a baseline method. As one of specific implementations of Alard & Lupton (1998), HOTPANTS adopts Gaussian basis functions to construct the convolution kernel. In real space, the use of DBFs has been shown to yield better performance than Gaussian basis functions in terms of photometric accuracy (Albrow et al., 2009). Here we present the comparisons to HOTPANTS to confirm that the SFFT, though working in Fourier space, will still inherit the advantage of DBFs over Gaussian basis functions. The parameters of HOTPANTS used for this section are listed in Table A2. Recall that sparse-flavor SFFT restricts the calculations of convolution kernels on a pre-selected set of sources, this source list used by SFFT will be shared by HOTPANTS for the sake of fairness. Besides, both methods always apply the consistent size of convolution kernels and use the same polynomial degree of spatial variation for image subtraction. We selected 25 DECam observations with the same pointing but diverse seeing conditions from our DECamERON survey (Mould et al., 2017) (see Table B). With the set of DECam images, we conduct two types of tests to check the subtraction cleanness and photometric accuracy, respectively.

For the subtraction cleanness test, we use the image DECam-SREF as the shared reference, having the best seeing among all of them. The remaining 23 images aliased by initial DECam-OBS are regarded as science images in this test. Image subtractions are performed for the sequence of science images with a shared reference using sparse-flavor SFFT and HOTPANTS, respectively. An example of the subtraction results is shown in Figure 7 from the observation DECam-OBS02f. For both methods, the calculations involved in parameter-solving have been restricted to the same pre-selected sources, and the minimization aims at optimal subtraction on these sources as clean as possible. It is useful to give a close-up view of subtraction performance specifically for the selection. For each image in Figure 7, we cut a small thumbnail image around each selected source and combine the cutouts into a grid to be a new synthetic image, as shown in Figure 8. Note that the source selection, in this case, is more than 1100, so we only exhibit a subset of size 64 from the bright end of the complete set for a clear display. One can notice that the most conspicuous subtraction-induced artifacts have a dipole-like pattern which is especially prominent for the bright sources. The pattern can be found in both difference images, but much less profound in SFFT subtractions than in HOTPANTS subtractions. It is an expected improvement as the SFFT kernel solution with DBFs allows more degrees of freedom than HOTPANTS using Gaussian basis functions.

To further verify such an improvement is real in a statistical sense, here we defined a naive quantifiable metric to describe the subtraction cleanness and show the performance comparisons with the test data. Note that the subtraction-induced artifacts on difference images are usually detectable by SEXTRACTOR like other real transient and variable objects. Generally, a huge number of artifacts in transient surveys will survive until an AI-based stamp classifier recognizes them. This fact makes it possible to use a simple dichotomous metric to describe the subtraction cleanness. A given source is either clean-subtracted or not according to the existence of SEXTRACTOR detection of its subtraction residuals. Now we can investigate how the probability of clean-subtraction for a given set of sources is related to the employed image subtraction method. For transient surveys, it is essential to probe the subtraction performance over the galaxies in the field, which are potential hosts of transient events. When a transient emerges at a position very close to its host galaxy core, the subtraction-induced artifacts from the galaxy can severely hinder the discovery and photometric measurement of the transient.

Our deep DECam observations in the test (3 square degrees, down to  23.5 mag in i-band) have covered a vast number of galaxies with various morphological types. It is interesting to make statistics for the probability of clean-subtraction for those sources identified as extended galaxies (i.e., Tractor morphological types DEV, EXP, and REX) in the Legacy Survey (Blum et al., 2016). The check radius from the galaxy core for searching subtraction-induced artifacts using SEXTRACTOR is specified as 11 pixels ($\sim$ 2 times median FWHM for DECam). Figure 9 shows the rate of clean-subtraction as a function of the magnitude of the examined galaxy. As expected, the brighter galaxies have a higher probability of leaving detectable subtraction-induced artifacts. However, the SFFT method is more likely to achieve clean-subtraction than HOTPANTS.

For the photometric accuracy test, we use the image DECam-PTAR as the PSF target, which has a median level seeing among all selected DECam observations. The remaining 23 images aliased by initial DECam-OBS are convolved to achieve homogeneous PSF using sparse-flavor SFFT and HOTPANTS, respectively. We derive the light curves of the field stars by conducting photometry on the sequence of convolved images with a fixed aperture. The photometric accuracy measured by the light curves can work as a metric to assess the PSF matching quality. We investigate the aperture photometry of the point sources classified as Tractor morphological type PSF. Figure 10 shows the RMS error of the derived light curves as a function of the star brightness. As the photometric uncertainty from SFFT is systematically lower than that from HOTPANTS especially at the bright end, the improvement in PSF matching is again confirmed.

Figure 11 shows an example with TMTS data. The subtraction artifacts generated by SFFT and HOTPANTS are not identical. HOTPANTS produces relatively circular and symmetric patterns, while SFFT shows more random residuals. This obviously originates from the different base kernel construction in SFFT and HOTPANTS. SFFT uses a $\delta$-basis kernel which is more flexible, while HOTPANTS utilizes Gaussian functions which possess a higher intrinsic symmetry about the kernel center.