Image Co-Addition with Temporally Varying Kernels

In astronomy today, the standard technique for acquiring deep astronomical images is through a series of short, often dithered, exposures. Applying this procedure enables the removal of cosmic rays, the filtering of bad pixels, and the masking of device defects. It also increases the dynamic range of the resulting image (as bright sources in single, long exposures saturate) and provides better control of the underlying point-spread function (PSF) as long as the PSF is fully sampled.

To achieve the depth of a long exposure from a series of exposures requires that we combine the individual images, accounting for the variation in seeing, sky transparency, and background variability. This process is typically referred to as "image co-addition" and is at the heart of many image processing systems. The criteria that we optimize in the image co-addition depend on the science in question at hand. For example, in faint galaxy surveys (e.g., the Hubble Ultra Deep Field; Beckwith 2006) photometric accuracy might be paramount, but for weak-lensing surveys (e.g., the Canada-France-Hawaii Telescope Legacy Survey [CFHTLS]; Parker 2007), the preservation of the underlying image resolution may be the primary concern. Given these differing objectives, two open questions remain: how do we optimize the co-addition of astronomical images, and is there is an optimal criterion that is relevant to all scientific goals?

With a new generation of deep imaging surveys coming online throughout this decade, including the Panoramic Survey Telescope & Rapid Response System (Pan-STARRS¹), the Dark Energy Survey (DES²), and the Large Synoptic Sky Survey (LSST³), each producing petabytes of data, image co-addition cannot just be assessed by the statistical analysis; it is also a question of the computational efficiency of the different approaches. In this article, we compare several different approaches for image co-addition, taking into account their statistical properties and their computational cost. We consider in detail the Fourier-domain methods proposed by N. Kaiser (2004, private communication) along with several straightforward methods, such as running means or medians. We compare these methods using two performance criteria: flux conservation (i.e., pixelwise photometric accuracy) and image resolution (using a metric we call image quality). Using these comparisons we can make statements about whether the additional complexity outweighs the extra costs they incur. Our goal is to provide a framework for assessing the magnitude of the practical advantages that more complex methods might provide.

In § 2, we set up the statistical model underlying our analysis and outline the properties of the various techniques. In § 3, we describe the simulations and define the performance criteria, and in § 4, we present the results of our simulations.

1.1. Mathematical Setup

Consider a single image of a particular patch of sky. The object of interest is the true scene, which we represent as a function g that takes a two-dimensional argument corresponding to position in the image and returns an intensity. For ground-based viewing, however, we observe only a blurred, discretized, and noisy version of g. We represent the blurring (caused by atmospheric seeing and instrument artifacts) by an operator K. In this case, K is associated with a PSF k such that the action of K on g is to integrate g against k.

Now, suppose that instead of getting one image we get a stack of L images, all of the same patch of sky with scene g. We assume for this discussion that the images have been registered.⁴ In those cases where differential chromatic refraction plays a role in the quality of the co-added images, such as data taken in bluer passbands and at high air mass, users may want to consider creating both high and low air-mass co-adds or adjusting the positions of stars based upon their known colors and the air mass of observation. However, we do not pursue these ideas in this article.

Then, due to atmospheric conditions changing over time, for each image i there is a possibly different operator K_i with associated PSF k_i. Then, if we define a regular grid of pixel centers x₁,x₂,..., the observation on the i th image at the j th pixel is⁵

where the _ij are noise random variables with suitable distributions. Note that equation (1) also defines the operator K_i as integrating g against the PSF k_i. The statistical problem is to estimate g given observation of Y = (Y₁,Y₂,...).

What makes the problem challenging is that the operators K_i are ill-conditioned, meaning that large changes in g can make relatively small changes in Y_ij, making recovery of g from Y_ij unstable. This is so because blurring dampens high-frequency structure, so the singular values of this integral operator decay quickly to zero. Thus, even small levels of noise can critically obscure structure in g.

To see this in a simple example, assume that the PSF does not vary in shape across the image, which makes k(x,y) = k(x - y) a convolution kernel. It follows that by using the Fourier transform , we can recover g exactly in the absence of noise by

because . While the inverse exists in the absence of noise, the decay of at high frequencies can produce catastrophically high variances in the presence of noise, which produces nontrivial and random contribution at high frequencies (so the numerator is nonzero, but the denominator goes to zero). Figure 1 demonstrates this effect in a one-dimensional analog of the problem, comparing noise-free reconstruction via equation (2) with reconstruction under noise with a 10e9:1 signal-to-noise ratio (S/N). The estimate in the latter case is wholly untenable.

Zoom In Zoom Out Reset image size

1.2. Statistical Model

We assume that each observation is drawn from a distribution with mean

Because the observation process in the telescope involves counting photons, the distribution of Y_ij is well modeled by a Poisson distribution (Hu 2007; Scully 1969). Thus, we assume that the Y_ij are independent Poisson 〈θ_ij〉 random variables. If the mean counts θ_ij are large enough, the Gaussian approximation to the Poisson distribution is accurate, and we can model the Y_ij as independent normal 〈θ_ij,θ_ij〉 random variables.

The statistical problem is to form an estimator of g given observations Y_ij. In the present context, is the template or co-add we are attempting to create. For each we need ways to measure its misfit relative to the true scene g. In this article we consider two. The first criterion corresponds to flux conservation, which we measure through mean integrated squared error (MISE), and the second criterion we refer to as image quality. We give a brief description of these here. See § 3 for a more thorough discussion.

Flux conservation gives an indication about how well a method maintains the spatial location of a given amount flux. MISE, which corresponds to the expected value of the integrated squared difference between any recovered image and the true image g, is the measure of flux conservation we use. Specifically, suppose we form a template . Then

See Silverman (1985) for a thorough overview of MISE. Note that in practice, we compute MISE on the pixelized g and hence use

where (x_i) are the locations of each pixel. The technical implications of this approximation for this article are not relevant. Without strong assumptions, we can only hope to recover the pixelized g. Therefore, we can treat it as the object we wish to recover and take the MISE to be defined by (4).

While MISE measures the degree to which maintains flux in the correct location, it does not distinguish between different scatterings of the remaining misspecified flux. We define a second comparison, which we call image quality, that is designed to measure how spread out the source is, as this is not captured by MISE. For image quality, we do a very rough PSF estimation on a source in the co-added image . The relative size of these PSFs give us an indication about how spread out an object gets by the co-addition. See Figure 2 for an example of procedures with combinations of good and bad MISE and image quality. Note that image quality introduced here for stars is perhaps different from an analogous metric for extended objects such as galaxies. Though interesting, we do not consider a separate criterion for this article.

Zoom In Zoom Out Reset image size

One significant complication in practice is that the PSFs are unknown and must be estimated as part of the procedure. While it is ideal to simultaneously estimate the PSF and g, it is common, and typically effective, to first estimate the PSF and then use the estimate as a known PSF in the co-addition procedure. Because our focus is on the relative performance of techniques for estimating g, we will follow standard practices and assume that the PSFs are known. This produces a comparison of the techniques at their best, but does not account for how well the procedures tolerate misestimation of the kernels. Note that the straightforward comparison methods do not require PSF estimates. So any method that does require knowledge of the PSF would need to be substantially better, given that it has access to much more information.

In this section, we describe a variety of methods that have been applied to the co-addition problem, with particular focus on the Fourier-domain approach we are analyzing. Note that we restrict our attention to noniterative methods.

Lucky Imaging.—In lucky imaging, a large number of images are observed and only the best, according to some criterion such as seeing, are retained. Fried (1978) and Tubbs (2003) describe such implementations in detail. An advantage of this approach is that the reconstruction of the true scene is based entirely on high-quality data. A disadvantage is that the method requires storing many images to determine which are best. Moreover, the images that are discarded can contain useful information about the scene that is, in effect, wasted.

Pixelwise Statistics.—If a stack of images is aligned, the values at a particular pixel in each image represent a random sample from a common distribution. The aligned pixel stacks can thus be used to estimate the parameters of these distributions and, in turn, the true scene. Basic choices would be a mean, median, or perhaps a trimmed mean to account for heavy-tailed noise, but more sophisticated estimators tailored to particular distributional assumptions can be constructed. The advantages of these approaches are computational and conceptual simplicity. Moreover, the mean can be computed sequentially with a fixed amount of storage, although this is not true of the median or trimmed mean, for which the entire image stack needs to be maintained. A big disadvantage is that pixelwise methods tend to create rough, discontinuous images, as no information sharing is permitted between pixels.

Fourier Deconvolution.—In the spatially constant seeing case, equation (1) has a convolutional structure k_i(x,y) ≡ k_i(x - y). The advantage of Fourier-based methods is that the operator K_i decomposes nicely in the Fourier domain. An inherent disadvantage to the Fourier approach is that great care must be used to avoid outcomes such as in Figure 1.

The approach outlined by N. Kaiser (2004, private communication), building on the previous work of Hook (1992) and Lucy (1992), is to Fourier transform each image and estimate the u th Fourier coefficient of g by a weighted average of the u th Fourier coefficient of each image. The weighting is accomplished so that the images with better seeing are weighted more heavily in the average. This method has some associated optimality properties. However, in practice, these properties turn out to not be useful.

Note that, as presented, this method makes three nontrivial assumptions. First, it assumes that each k_i(x,y) ≡ k_i(x - y), that is, spatially constant seeing. As we describe in § 4, while this assumption could have significant consequences when seeing varies spatially, in practice, it does not appear to cause much problem. The second assumption is that the θ_ij are large enough that the Gaussian approximation to the Poisson is accurate across the image. The third assumption is that the variance of the Gaussian is a known value that depends only on the image i. These last two assumptions make the analysis easier, and while they are often reasonable, they need not hold in practice.

To define the Fourier deconvolution estimator, first expand g into the Fourier basis, which we write as (ϕ_u),

In general, we use the notation for the Fourier transform of the function f.

Using the preceding assumptions and the resulting form of the θ_ij, N. Kaiser (2004, private communication) proposes two estimators of . See Table 1 for these estimators and some of their cumulants. The first estimator in the table, , has the smallest variance of all unbiased estimators of . Also, the resulting estimator of g, defined to be

is the best linear unbiased estimator of g.

However, as alluded to previously, these optimality properties are not useful. To see this, observe that the variance of is the sum of the variances of . Also, since the K 's are smoothing operators, is small for large |u|. In fact, the worse the seeing, the smaller becomes, and it can be smaller than machine error in many cases. Hence, as depends on the reciprocal of , the variance of can be arbitrarily large. Note that we have seen a case of this already in Figure 1, where this large variance property causes unusable results.

Generally speaking, demanding unbiased estimators is less effective than choosing an estimator with some bias but smaller variance. In this case, the distinction is crucial because the unbiased estimator is unusable. For the interested reader, Wasserman (2006) provides a good introduction to this surprisingly broad and deep issue.⁶ See Figure 3 for an example of equation (6) in a simple simulated situation. The ability of to reconstruct even a simple g decays rapidly, as the full width at half-maximum⁷ (FWHM) of the PSF increases. In particular, even the small amount of seeing (FWHM of 1.1 pixels) in the rightmost column makes the estimator unusable.⁸

Zoom In Zoom Out Reset image size

N. Kaiser (2004, private communication) provides a correction that prevents the variance from exploding for large |u| by multiplying by a term that is proportional to the reciprocal of its standard deviation. See the second row of Table 1. This multiplication biases the estimator but bounds the variance. We define

and we spend the balance of the article investigating its properties. We refer to as the Fourier deconvolution (FD) method. As an aside, since the correction corresponds to multiplication in Fourier space, we see that the bias of is given by (k^∗ ∗ g) - g, where the Fourier transform of k^∗ is defined in the caption of Table 1 and (a ∗ b) is the convolution of a and b.

2.1. Computational and Space Complexities

For our evaluation metrics, we need to know the true sky before any distortions. Hence, we simulate an idealized view of the sky (above the atmosphere) using a catalog of point and extended sources. Extended sources are represented by Sérsic profiles, and the densities of point and extended sources are designed to match observations to a depth of r ∼ 28. Figure 4 shows a representation of one such image. From these true scenes, we apply a blurring operator K and noise with variance σ² to represent the effect of the atmosphere, telescope, and instrument.

Zoom In Zoom Out Reset image size

Under the assumptions of the FD method, we assume that the PSFs are convolutional (by that we mean spatially constant seeing), but this leaves open the specific functional form for the kernel k. We choose to set k to be the two-dimensional Gaussian probability density function. Our results are, however, robust against more complicated parametrization of k (e.g., a mixture of Gaussians).

For evaluating the effectiveness of the FD method, we compare it with the pixelwise mean and pixelwise median. Our goal is to assess the value added by the FD method and evaluate whether the improved performance overcomes the added cost and complexity of the method. The mean and median procedures are intended as extreme comparisons rather than serious proposals and we are not endorsing their use in this article. However, if the value added of more sophisticated methods relative to even these simplistic methods is not great, then it suggests that simple methods may suffice in practice. Note that the computational complexity of the methods is not a purely academic concern. In large imaging surveys such as the LSST, a near-constant data stream requires efficient, near-real-time processing, and the cost of generating template images is an important factor.

For this comparison, we use the two general criteria introduced in § 1.2 corresponding to MISE and image quality. In principle, at least two different errors can be made when comparing an estimated template with the true scene g. These errors can most easily be thought of in the context of a point source such as a star. Our estimated image can remove some of the flux from the point source and/or it can smooth the image (degrading the photometric accuracy and the resolution of the resulting template). Figure 2 shows examples of a pixelated point source that has been reconstructed by four made-up methods.

In the MISE comparison we compute the expected value of the integrated squared difference between the true image and our estimate, assuming zero background. Under the non-Gaussian distributions, a larger variance results in more total flux in the image. This, in turn, results in a large MISE. To compensate for this increase, and to make all the cases more directly comparable, we rescale the MISE by the total signal of each image to compensate for this additional flux.

For the evaluation of image quality we fit a two-dimensional, spherical¹¹ Gaussian density to a source using least squares. The covariance matrix of this fitted Gaussian is then used as a metric for the width of the point source.

In this section, we present the results of our simulations. The section is divided into two parts, corresponding to the two evaluation criteria outlined previously.

4.1. MISE

When looking at equation (1) we see two major influences on the quality of the observation. One is the type and severity of the seeing, determined by the PSF k. The other is the distribution of the noise . In this section, we look at several different scenarios for comparing the mean or median with the FD method by changing these two factors.

4.1.1. Increasing FWHM Comparison

Overall, we wish to understand the impact of seeing and noise on the performance of the methods. To do this, we simulate a stack of images with each image having a FWHM drawn from a distribution with mean, μ_FWHM. We do this procedure for μ_FWHM ranging from 2 pixels to 7 pixels (0.4'' to 1.4''). This interval is chosen so that it contains the expected seeing width for the LSST. This procedure creates a sequence of stacks of images.

We apply this increasing FWHM to six different noise parameterizations. For both high variance (5:1 S/N) and low variance (20:1 S/N) we consider noise models that are Gaussian, heavy-tailed shot, and inhomogeneous Poisson. In all six combinations, a lower MISE indicates that the method conserved the flux better, on average. See Figure 5 for the results.

Zoom In Zoom Out Reset image size

In all cases, the mean severity of seeing, μ_FWHM, does not impact the ranking of the methods. In the Gaussian case, the methods are ranked, from best to worst, as FD, mean, and then median. In both the low- and high-noise cases, absolute difference between the mean and the median is constant for all levels of μ_FWHM, due to the well-known efficiency gains of the mean over the median in the Gaussian case. However, in the high-noise case, the difference in MISE between the FD and mean methods increases 11% as μ_FWHM ranges from 2 to 7 pixels. As the FD method uses the additional information of a known seeing kernel, this nonconstant difference is expected.

In the inhomogeneous Poisson or heavy-tailed shot-noise cases, the median outperforms both other methods. This is due to the particular and well-known feature of the median to be more robust against heavier-tailed distributions. In the shot-noise case, the median ignores the random noise spikes with overwhelming probability and we get a noise-free version of the true scene g with the median amount of seeing. The very small advantage of the FD over the mean method is due to the FD's slightly narrower impulse response function, which we discuss in § 4.2. In the inhomogeneous Poisson noise case, the observed image is created by simulating independently from a Poisson〈θ_ij + σ²〉 from equation (3). This creates a very noisy image, especially near sources, and hence the median's ability to ignore large transient spikes results in the large MISE advantage.

Under all three noise distributions, the FD method displays a small but increasing benefit in flux conservation over the mean as μ_FWHM goes from 2 to 7. It should be kept in mind that the FD method uses the known kernel assumption. Hence, as the seeing becomes more severe, this informational advantage should become more pronounced. In reality, the kernel would need to be estimated and this advantage would most likely decay.

4.2. Image Quality

To do this comparison we use an image where the function g is a δ-function (i.e., a point source). We use the expected distribution of seeing at the LSST stack to draw 10 independent and identically distributed Gaussian PSFs and apply them to the point source. This mimics the forward process of an image getting blurred by the atmosphere and corresponds to operating on g with K_i for i = 1,2,...,10. For this stack of images we apply the three methods described earlier, FD, mean, and median, to derive a template image. We then fit a spherical Gaussian kernel to the template via least squares and use the diagonal element of the covariance matrix as a measure of the width of the method's PSF.

We make 1000 draws from the seeing distribution and compute this statistic for each method and each draw. The results are summarized in Figure 6 with box plots of the FWHM of the fitted Gaussian kernel. box plots are a graphical tool for quickly conveying the distribution of observed data and are composed of a "box" and "whiskers." The box is the main rectangle, which has horizontal lines, increasing with FWHM, at the 25th, 50th,¹² and 75th percentiles of the data. The whiskers correspond to the dotted vertical lines with horizontal lines at the end and the plus signs. These horizontal lines are at the first and 99th percentiles. The pluses are at extreme data points.

Zoom In Zoom Out Reset image size

The box plots are very similar and show little difference in the FWHM of the fitted kernels. The FD method has longer whiskers than the other methods, indicating a wider range of fitted FWHM that can be expected. It also has a very slightly lower box than the other two methods, indicating that the majority of draws result in a slightly better fitted FWHM. For instance, the 50th percentile of the FD method's fitted FWHM is 1.2% lower than the mean method's fitted FWHM.

Interestingly, one of the draws resulted in six very good seeing images and four poor images. This resulted in the median method having the smallest fitted FWHM recorded for any method on any draw. This is part of a more general property of the median to behave noncontinuously with the composition of the stack. For instance, if there were instead four very good seeing images and six poor seeing images, the result would be a large fitted FWHM.

4.2.1. Differences in Expectation

The mean and FD methods are very similar in that they are both based around summing the images and hence are both linear filters on the observed images. Let us denote the estimated co-adds of the two methods as and , respectively. We can understand the filtering that we are applying to the true image g by considering and . Note that these are

where the Fourier transform of the function k^∗ is defined in the footnote of Table 1.

By examining the impulse response functions (IRFs) of these filters we can find the average width of the PSFs that will result from these methods. Note that in this case this corresponds to the bias of each method. See Figure 7 for a one dimensional representation of the IRFs of the expected filters. The mean IRF and the Fourier deconvolution IRF are virtually identical and would create no visible difference in the image quality.

Zoom In Zoom Out Reset image size

A variety of template generation/image co-addition techniques have been applied in current imaging surveys, including the CFHTLS (Gwyn 2008; mean technique), Pan-STARRS (Price 2007; mean technique), Model Parameter Estimation Experiment (MOPEX) for Spitzer images (Makovoz 2005; mean and median techniques), and the VIRMOS survey (Radovich 2004; median technique). In this article, we address the broad question of whether sophisticated methods of template generation and image co-addition are worth the added cost and complexity, in the context of modern image surveys with large and nearly constant data streams. To this end, we compare Fourier-domain reconstruction techniques, which have a variety of nice theoretical properties to straightforward pixelwise statistics. We consider two cost metrics to evaluate the image reconstructions: image quality by way of the resulting Gaussian FWHM of an image and the conservation of flux, measured by MISE.

Under the assumptions that motivate the FD method, namely, Gaussian noise and spatially constant seeing, the FD and mean methods both have lower MISE than the median. However, in the heavy-tailed shot and inhomogeneous Poisson noise cases, the median outperforms the FD method.

In situations when there is a small amount of seeing in an image stack, corresponding to μ_FWHM ≈ 2 pixels, we see the mean and FD methods have nearly the same MISE properties (less than 0.7% difference across all considered noise distributions). However, as the severity of the seeing increases, the disparity between the methods increases as well. As the FD method incorporates knowledge of the seeing kernel, this property is to be expected. However, the MISE for the high-noise Gaussian, shot, and inhomogeneous Poisson cases are only 5.4%, 1.2%, and 0.13% lower, respectively, at the most severe levels of seeing.

These benefits are rather small, particularly in the more realistic heavy-tailed and Poisson noise cases. This is especially so considering the extra computational complexities involved in the FD method over the mean and median methods.

For comparing image quality of the templates, we find that the resolution of the resulting images is very similar, as the 50th percentile of the fitted FWHM of the FD method is less than 2% smaller than the 50th percentile of the fitted FWHM of the mean method.

We find that the value added by the FD method over the mean or median procedures does not overcome its added cost and complexity. This leads us to the conclusion that there is room for improvement over all three methods. We are currently exploring various possible improvements and theoretical justifications for some of the approaches outlined in this article.