High Order Transform of PSF ANd Template Subtraction
This package is designed to photometrically align one input image with another, after they have been astrometrically aligned. This is an implementation of the Alard 1999 algorithm for image subtraction. This is very similar to the ISIS package, but with a variety of changes and (dare I say it again) improvements. The software has been applied in a variety of real-time data pipelines with great success, and much work has gone into making this a robust package.
Description
The crux of the difference imaging problem is to find a convolution kernel K that matches the PSFs of two astronomical images, I (referred to as the image) and T (referred to as the template). These images in general are taken under different conditions, including atmospheric transparency, atmospheric seeing, or exposure times. One may even difference data taken from different sites and equipment entirely. However, this technique does not work well for differencing data taken in different filters - to the extent that instrumental response as a function of wavelength is different on different equipment, differencing data taken from different sites does require some care.
Mathematically, the equation you want to solve is the minimization of the function
These defaults are not optimal in the global sense that you really
want to feed in the relative sizes of your PSFs so that you can
estimate the scales at which they are expected to be different, and
thus the size of the gaussians that make up your convolution kernel.
HOTPANTS allows you to enter on the command
line the sizes of the kernel gaussians and their spatial orders, as
well as the size of the kernel. In practice, HOTPANTS is called by Perl or
Python scripts that scale these as functions of the
relative PSF sizes. However, the way that HOTPANTS is implemented, these cannot be fit
for.
Method
HOTPANTS divides the image up into multiple regions. You fit for one convolution kernel within each region. An image might have multiple regions if e.g. the aspect ratio is large (4 x 1) and you can't expect a single function to model its spatial variation. Each region is divided up into stamps, generally 1 per every 100 pixels (a 2k x 4k image has 20 x 40 stamps). Within each stamp, you choose multiple sub-stamps that are centered on individual astronomical objects. The idea is that you want constraints on the kernel in each stamp, and you have multiple options within each stamp in case particular objects are sigma-clipped from the process. Since you are fitting for a function that is spatially varying, you want to ensure constraints across the entire image (especially in the corners).
Input masks may be fed to the code, which will allow it to avoid bad pixels while looking for sub-stamps. This allows it to distinguish between completely bad pixels which should be ignored (e.g. saturated pixels), pixels which you might trust but don't want to use while fitting for the kernel (e.g. interpolated pixels), and pristine untouched pixels. These masks are convolved and may be output from the software to be used in subsequent stages of the processing (e.g. diffim photometry) that are designed to take advantage of this information. In addition, noise maps may be fed to, and delivered from, the code. Because of the convolution, some of the noise will get lost in off-diagonal covariance terms. This is unavoidable, the best you can do is track the image noise so that you know the significance of detections in your difference image.
In practice, the program determines (or is fed a list of) objects to fill the sub-stamps. These objects are extracted from the template and image, and fed to the the matrix equations
For each object, this kernel is determined by solving for a above. Because some objects may be varying in brightness or position, the kernel sum is used as an initial way to sigma clip outliers from this distribution. Another figure of merit used to reject bad stamps is the average pixel residual (after convolution and template subtraction) divided by the estimated variance in that pixel (the true noise is not known until the final convolution kernel is determined). Another figure of merit that may be used is the width of the histogram of residuals in this difference stamp. HOTPANTS will choose which image to convolve based upon these figures of merit (deconvolution tends to leave larger residuals), unless the user demands otherwise.
Since the PSF varies spatially in all astronomical images, we must model K as a spatially varying function also. In practice, we use
Since we have decomposed the spatial variation of the kernel using basis functions, we are able to solve for the spatial coefficients of the kernel with another set of linear equations. We have constraints on the kernel at the location x, y of each sub-stamp. We solve for the spatially varying coefficients using the set of all sub-stamp solutions. Outliers from the fit are sigma-clipped out, and neighboring sub-stamps from the stamp are brought into the fit in their place. This process repeats until no sub-stamps are clipped out.
Metrics For Success
The photometric scaling of the two images comes out in the convolution
kernel sum. Ideally, this should relate to the
relative zero-points in the image through -2.5 *
\log({\rm ksum}) = {\rm Zpt}_I - {\rm Zpt}_T. This plot is
shown here for an ensemble of SuperMACHO data. The \chi^2 / dof for this fit is very nearly 1.
The distribution of residuals ({\rm diff} /
\sigma)^2 in each sub-stamp also provides an estimate of how
good the overall fit is. A complication is that, because of
convolution, some of the noise is lost in covariance, and thus the
noise per pixel in the final image is an underestimate of the true
noise. This biases this metric to values larger than 1.0. To
demonstrate this, two sets of figures are shown below. For each, the
plot on the left shows a histogram of FSIG, the average value
across sub-stamps of the average sigma per pixel; on the
right is how many sigma this value is from the expectation value (1.0)
using the RMS of the distribution across sub-stamps (FSCAT).
5-sigma outliers are clipped. At the top is the plot for 100k
SuperMACHO difference images; at the bottom is the same for 14k
ESSENCE difference images.
