Index: /trunk/doc/misc/imageSubtraction.tex =================================================================== --- /trunk/doc/misc/imageSubtraction.tex (revision 3368) +++ /trunk/doc/misc/imageSubtraction.tex (revision 3368) @@ -0,0 +1,268 @@ +%%% $Id: imageSubtraction.tex,v 1.1 2005-03-03 00:38:39 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 the density of +high signal-to-noise stars is low; 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} +