## 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.