%%% $Id: imageSubtraction.tex,v 1.2 2005-03-04 02:48:52 price Exp $
\documentclass[panstarrs]{panstarrs}

% basic document variables
\title{Image Subtraction Algorithm}
\subtitle{Exploring an Alternative}
\shorttitle{Modules SDRS}
\author{Paul Price}
\audience{Pan-STARRS PMO}
\group{Pan-STARRS Algorithm Group}
\project{Pan-STARRS Image Processing Pipeline}
\organization{Institute for Astronomy}
\version{DR}
\docnumber{PSDC-430-???}

\setlength{\topsep}{-2pt}
  
\begin{document}
\maketitle
\sloppy

% -- Revision History --
% provide explicit values for the old versions
% use '\theversion' for the current version (set above)
% use \hline between each table row
\RevisionsStart
% version  Date            Description
DR & 2005 Feb 24 & Draft \\ \hline
\RevisionsEnd

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\DocumentsInternal
PSDC-430-011  &   Pan-STARRS PS-1 IPP System/Subsystem Design Description \\ \hline
\DocumentsExternal
Posix Standard & Open Group Based Specifications Issue 6, IEEE Std 1003.1, 2003 \\
\DocumentsEnd

\tableofcontents
\pagebreak 
\pagenumbering{arabic}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

The goals of the \PS{} system can be (overly?) simply expressed
as the following:
\begin{enumerate}
\item To make deep images of the static sky; and
\item To identify variable objects.
\end{enumerate}

We have explored the first of these goals in a companion document
on image combination.  Here, we explore the second goal.

By ``identification of variable objects'', we mean recognising a
source which has been detected in an image but was not detected at the
same position in a template image.  There are two methods which are
used to identify variable objects.

The first method simply generates a catalog of all sources on each
image and compares the two catalogs.  This is a fast and simple method
for identifying variable sources, but suffers from a strong drawback:
it does not account for differences between the point spread functions
(PSFs) of the two images.  For example, say that the template image is
taken in poor seeing, and the new image has been taken in excellent
seeing.  What appears to be a single star in the template image might
be resolved as a double star in the input image, leading to the false
identification of a variable object.  This problem has led to the
general rejection of this method for identifying variable sources in
major synoptic surveys, apart from under-sampled systems such as
RAPTOR and ROTSE.

The second method overcomes the differences in PSFs between the
images, but it is more computationally expensive.  It involves
matching the PSFs between the two images by solving for the
appropriate convolution kernel.  This method has been in use for quite
a while, but the most famous and generally used implementation is that
of Alard (2000), known as ISIS
(\href{http://www2.iap.fr/users/alard/package.html}{website}).

ISIS provides a fast method of solving for a convolution kernel that
matches the PSFs in each image.  The resultant subtracted image has
deviations close to that expected for a Poisson distribution.  ISIS
uses a linear combination of Gaussians modified by polynomials as a
basis set for the convolution kernel, avoiding the tortuously slow
approach of non-linear fitting using a more physically motivated
kernel.  One obvious question, then, is why is ISIS so successful?  In
some cases, especially around brighter stars, ISIS leaves systematic
residuals visible in the subtracted image.

In this document, we explore the possibility of using a more general
image subtraction algorithm, in the expectation that adding generality
would allow smaller systematic residuals in the subtracted images.
One important issue, which is how to robustly determine which
residuals in the subtracted image are real, and which are simply part
of the afore-mentioned systematic residuals, is not dealt with here,
but deferred in the expectation that a clean subtracted image (with a
minimum of systematic residuals) will simplify this matter.  We are,
for now, solely concerned with achieving the cleanest possible image
subtraction, believing that this is the best way to robustly identify
variable sources.

\section{Convolution Kernels}

PSF matching is achieved by solving for the convolution kernel:
\begin{equation}
I(x,y) = R(x,y) \otimes K(u,v)
\end{equation}
where $I(x,y)$ is the input image, $R(x,y)$ is the reference image,
$K(u,v)$ is the convolution kernel, and $\otimes$ denotes convolution.

The speed of ISIS comes from choosing to represent the convolution
kernel as a linear combination of basis functions:
\begin{equation}
K(u,v) = \sum a_i B_i(u,v)
\end{equation}

In this case, the problem of minimising $\chi^2$ reduces to solving a
matrix equation:
\begin{equation}
b_i = \sum_j M_{ij} a_j
\end{equation}
where
\begin{equation}
b_i = \sum_{x,y} I(x,y) [ R(x,y) \otimes B_i(u,v) ] / \sigma(x,y)^2
\end{equation}
and
\begin{equation}
M_{ij} = \sum_{x,y} \left[ R(x,y) \otimes B_i(u,v) \right] \  \left[ R(x,y) \otimes B_j(u,v) \right] / \sigma(x,y)^2
\end{equation}
Here the sum over $x,y$ indicates summing over select regions (known
as ``stamps'') that are known not to be variable (or, in practise,
assumed not to be variable, with iterative rejection of outliers).
$\sigma(x,y)$ is the error map (which can easily be constructed from
the input images by assuming that the noise is Poisson).

In this section, we consider two choices for the basis functions,
$B_i(u,v)$.

\subsection{ISIS}

ISIS chooses for its basis functions a set of Gaussians, each modified
by a polynomial:
\begin{equation}
B_{ijk}(u,v) = e^{(u^2 + v^2)/2\sigma_i^2} u^j v^k
\end{equation}

Alard \& Lupton (1998) claim that, ``In practice, it seems that three
Gaussian components with associated polynomial degrees in the range of
2 to 6 can give subtracted images with residuals comparable to
$2^{1/2} \times$ photon noise.''

It is not immediately clear why such a solution would work well ---
there is no physical motivation for this choice of basis functions.
Nevertheless, it does seem to work well, and is widely used.

\subsection{POIS}

Pricey's Optimal Image Subtraction (POIS) tries a different choice of
basis functions --- a set of delta functions:
\begin{equation}
B_{ij}(u,v) = \delta(u - i)\ \delta(v - j)
\end{equation}
In essence, this choice of basis functions amounts to measuring (a
discrete version of) the kernel directly, and is therefore
intrinsically more flexible than ISIS' choice of basis functions.
While this could conceivably yield better subtractions, it may also
introduce a couple of problems:
\begin{itemize}
\item Speed: the number parameters goes as the number of pixels in the
  kernel.  Especially when adding spatial dependence to the kernel,
  this can lead to a very large number of matrix elements.
\item Noise: neighbouring pixels in the kernel are not correlated,
  which raises the possibility of the noise being increased relative
  to a kernel which is smooth.
\end{itemize}

\section{Test cases}

We will place the two kernels in a head-to-head competition to produce
the smallest stellar residuals.  May the best kernel win.

In comparing the stellar residuals, there are two quantities of
interest:
\begin{enumerate}
\item The deviation, $\sum_{x,y} s(x,y)$; and
\item $\chi^2 = \sum_{x,y} [ s(x,y) / I(x,y) ]^2$.
\end{enumerate}
where $s(x,y) = I(x,y) - R(x,y) \otimes B(u,v)$ is the subtracted
image.  The first simply looks at the absolute deviation from what is
expected (zero), while the second is a measure of the deviation
relative to the noise.  We will consider these quantities, or rather,
the distribution of these quantities, for both kernel choices.
Obviously, lower values indicate a better subtraction.

\subsection{Test data}

Our test data consists of images taken with MegaCam of a photometric
standard field.  We will examine two chips --- one near the centre of
the mosaic, and another near the corner.  The former is expected to be
relatively free from spatial variations of the PSF, while the former
may suffer strongly from these effects.  It is worth noting, however,
that the pointing of the observations were not very different, and so
spatial variation may not be an issue.

\subsection{Results}

We are not using any spatial variation of the kernel.  The quantities
are being measured within a 5 pixel radius aperture, centered on stars
detected in the image.  We use clipping (2 iterations at 3.5$\sigma$)
to obtain statistics on the population unaffected by small issues
(such as very bright stars which haven't been properly masked).

\begin{table}[h]
\begin{tabular}{|l|r|r|c|r|r|c|} \hline
Dataset & \multicolumn{3}{c|}{Deviation} & \multicolumn{3}{c|}{$\chi^2$} \\
        & Mean & Std.\ Dev.\  & Clipped & Mean & Std.\ Dev.\ & Clipped \\ \hline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Centre, ISIS & -58 & 221 & 32/884 & 0.6 & 3.2 & 9/884 \\
Centre, POIS & -7 & 238 & 31/884 & 0.5 & 2.3 & 9/884 \\
Edge, ISIS & -224 & 269 & 102/2968 & 0.05 & 0.19 & 67/2968 \\
Edge, ISIS & -281 & 278 & 107/2968 & 0.04 & 0.14 & 67/2968 \\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\hline
\end{tabular}
\end{table}

These test cases indicate that the POIS kernel is able to achieve
similar deviations as ISIS, but it has a better distribution of
$\chi^2$.  This is clearly not due to clipping more values when
calculating the statistics, but indicates a real improvement in the
quality of the subtraction.

\begin{figure}[h]
\psfig{file=subtraction.ps,width=15cm,angle=0}
\caption{Input image (left), ISIS subtraction (centre) and POIS
subtraction (right).  The stretch in the subtractions are the same.}
\end{figure}

\subsection{Timing}

An optimised (\code{-O2}) version version of POIS performs the
subtraction on the MegaCam images ($2112\times 4644$ pixels) on
\code{alala} (dual Opteron 2.2~GHz; note that the code is not
multithreaded) in 15 sec (including I/O), using a $9\times 9$ kernel.

On \code{mithrandir} (dual Xeon 2.4~GHz), POIS runs in 14~sec,
compared to ISIS in 41~sec.  It is worth noting that the ISIS
implementation uses strictly integer arithmetic, which should be
faster than the floating point arithmetic of POIS.

\section{Conclusion}

Our head-to-head comparison indicates that POIS performs a subtraction
of comparable quality to ISIS, but is faster by a factor of a few.  We
are yet to compare the quality of subtractions where there is a low
density of stars with good signal-to-noise; in this regime, ISIS might
provide a better subtraction than POIS, simply because the POIS kernel
pixels are not correlated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliographystyle{plain} \bibliography{panstarrs}

\end{document}

