Index: trunk/doc/modules/ModulesSDRS.tex
===================================================================
--- trunk/doc/modules/ModulesSDRS.tex	(revision 3300)
+++ trunk/doc/modules/ModulesSDRS.tex	(revision 3425)
@@ -1,3 +1,3 @@
-%%% $Id: ModulesSDRS.tex,v 1.32 2005-02-22 19:30:45 eugene Exp $
+%%% $Id: ModulesSDRS.tex,v 1.33 2005-03-15 23:13:16 price Exp $
 \documentclass[panstarrs]{panstarrs}
 
@@ -1615,4 +1615,577 @@
 \end{verbatim}
 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Image Combination}
+
+The image combination for \PS{} will employ an iterative approach, in
+order to identify cosmic rays.  The first pass involves transforming
+and combining the input images, and noting pixels which are apparently
+deviant.  These pixels are examined in further detail, before a subset
+of them are declared to be bad, whereupon these pixels are
+re-transformed, and the images are combined properly.  Here we
+introduce two functions which will perform the combination and
+examination steps.  Prototype code exists for each of these functions.
+\tbd{For further details, see the document about image combination for
+\PS{}.}
+
+\subsection{Combining images}
+
+\begin{verbatim}
+psImage *pmCombineImages(psImage *combined, // Combined image
+			 psArray **questionablePixels, // Array of rejection masks
+			 const psArray *images, // Array of input images
+			 const psArray *errors, // Array of input error images
+			 const psArray *masks,// Array of input masks
+			 unsigned int maskVal, // Mask value
+			 const psArray *pixels,	// Pixels to combine
+			 int numIter,	// Number of rejection iterations
+			 float sigmaClip, // Number of standard deviations at which to reject
+			 const psStats *stats // Statistics to use in the combination
+			 );
+\end{verbatim}
+
+\code{pmCombineImages} shall combine the input \code{images},
+returning the \code{combined} image and a list of
+\code{questionablePixels} in each input image.  The array of error
+images, \code{errors}, shall be used to calculate the value in the
+combined image and the list of questionable pixels, if
+non-\code{NULL}.  Pixels whose corresponding value in the array of
+mask images, \code{masks}, matches \code{maskVal} shall be masked from
+the combination.  The \code{images}, \code{errors} and \code{masks}
+arrays, if non-\code{NULL}, shall all carry the same number of images;
+otherwise the function shall generate an error and return \code{NULL}.
+The sizes of all images in the \code{images}, \code{errors} and
+\code{masks} arrays shall be identical; otherwise the function shall
+generate an error and return \code{NULL}.
+
+If \code{pixels} is non-\code{NULL}, only those pixels specified shall
+be combined.  The combination consists of \code{numIter} iterations in
+which a stack of pixels is combined using the specified \code{stats}.
+In each iteration, questionable pixels are identified as lying more
+than \code{sigmaClip} standard deviations from the combined value;
+these pixels are excluded from the stack for the next iteration.  The
+value for the combined image is that produced by the \textit{first}
+iteration (i.e., with no pixels excluded except those which have their
+corresponding mask match the \code{maskVal}); this allows subsequent
+calls to the function to only act on a small fraction of the pixels,
+since questionable pixels identified in the first call of the function
+will be properly rejected at a later point (see the example, below).
+
+In the event that \code{images} or \code{stats} are \code{NULL}, the
+function shall generate an error and return \code{NULL}.
+
+\subsection{Rejecting pixels}
+
+\begin{verbatim}
+psArray *pmRejectPixels(const psArray *images, // Array of input images
+			const psArray *pixels, // These are the pixels which were rejected in the combination
+			const psArray *inToOut, // Transformations from input to output system
+			const psArray *outToIn, // Transformations from output to input system
+			float rejThreshold, // Rejection threshold
+			float gradLimit // Gradient limit
+			);
+\end{verbatim}
+
+\code{pmRejectPixels} inspects those questionable \code{pixels}
+identified by \code{pmCombineImages} to determine if they are truly
+discrepant.  This inspection is performed in the coordinate frame of
+the detector, where the pixels haven't been smeared by transformation.
+Two tests are applied to each of the \code{images}:
+\begin{enumerate}
+\item The list of questionable pixels for an image is converted to an
+  image which is transformed back to the coordinate frame of the
+  detector.  Those pixels in the detector frame which have a value
+  exceeding \code{rejThreshold} are suspected cosmic rays and
+  subjected to the next test.  Depending on the value of the
+  \code{rejThreshold}, this test basically amounts to demanding that
+  questionable pixels neighbor each other in the transformed image.
+\item The cores of point sources may mimic a cosmic ray, especially in
+  under-sampled images.  To minimise flagging stars as cosmic rays, we
+  determine the gradient around the pixel of interest; if the gradient
+  is large, then the pixel is likely the core of a point source.  In
+  order to reliably measure the gradient in the presence of a
+  suspected cosmic ray, we use the companion images --- the gradient
+  is the mean gradient at the corresponding position on the other
+  images.  In order to calculate the corresponding positions, the
+  \code{inToOut} and \code{outToIn} transformations are required.  If
+  the gradient is less than \code{gradLimit}, then the pixel is
+  identified as a cosmic ray.
+\end{enumerate}
+
+The function shall return an array of \code{psPixels}, one for each of
+the input \code{images}, containing pixels that have been identified
+as cosmic rays according to the above criteria.
+
+If any of the input pointers are \code{NULL}, then the function shall
+generate an error and return \code{NULL}.
+
+\subsection{Example}
+
+Here is an example of what the image combination routine looks like,
+demonstrating how the various pieces fit together.  The inputs are:
+\begin{itemize}
+\item \code{psArray *inputs}: Input detector images, each a
+  \code{psImage} of type \code{psF32}
+\item \code{psArray *inputMask}: Input mask images, each a
+  \code{psImage} of type \code{psU8}
+\item \code{psArray *inputsErr}: Input error images, each a
+  \code{psImage} of type \code{psF32}
+\item \code{psPlaneTransform *skyToDetector}: Maps from sky
+  coordinates to detector coordinates, each a \code{psPlaneTransform}
+\item \code{psRegion *combineRegion}: Sky coordinate pixels to combine
+\item \code{int numIter}: Number of iterations in combination
+\item \code{float rejThreshold}: Threshold for rejection
+\item \code{float gradLimit}: Limit for gradient
+\end{itemize}
+
+The output is the combined image.
+
+\begin{verbatim}
+    psArray *transformed = psArrayAlloc(nImages); // Array of transformed images
+    psArray *transformedErr = psArrayAlloc(nImages); // Array of transformed error images
+    psArray *transformedMask = psArrayAlloc(nImages); // Array of masks for transformed images
+
+    for (int i = 0; i < nImages; i++) {
+	psPixels *blanks = NULL;	// List of blank pixels
+	transformed->data[i] = psImageTransform(NULL, &blanks, inputs->data[i],
+						inputMask->data[i], inputMaskVal, NAN, skyToDetector,
+						combineRegion, NULL, PS_INTERPOLATE_BILINEAR);
+	transformedErr->data[i] = psImageTransform(NULL, NULL, inputsErr->data[i], inputMask->data[i],
+						   inputMaskVal, NAN, skyToDetector, combineRegion, NULL,
+						   PS_INTERPOLATE_BILINEAR_VARIANCE);
+	psImage *skyImage = transformed->data[i]; // Dereference the transformed image
+	psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
+											   // transformed
+											   // image
+	transformedMask->data[i] = psPixelsToMask(NULL, blanks, blankRegion, PS_MASK_BLANK);
+	psFree(blankRegion);
+	psFree(blanks);
+    }
+
+    psArray *rejected = NULL;		// Array of rejected pixel lists
+    psStats *combineStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistic to use in doing the combination
+    psImage *combined = pmCombineImages(NULL, &rejected, transformed, transformedErr, transformedMask, 0,
+					NULL, numIter, sigmaClip, combineStats); // Combined image
+    psArray *bad = pmRejectPixels(inputs, rejected, NULL, skyToDetector, rejThreshold, gradLimit); // Bad pix
+    psPixels *combinePixels = NULL;	// Pixels to combine
+    for (int i = 0; i < nImages; i++) {
+	psPixels *badSource = psPixelsTransform(NULL, bad->data[i], skyToDetector); // Bad pixels on the input
+	psImage *badMask = psPixelsToMask(NULL, badSource, PS_MASK_COSMICRAY); // Mask image for the input
+	(void)psBinaryOp(inputMask->data[i], inputMask->data[i], "|", badMask);	// Put CRs into original mask
+	psFree(badSource);
+	psFree(badMask);
+
+	combinePixels = psPixelsConcatenate(redo, bad->data[i]);
+
+	// Update transformed image
+	psPixels *blanks = NULL;	// List of blank pixels
+	transformed->data[i] = psImageTransform(transformed->data[i], &blanks, inputs->data[i],
+						inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY, NAN,
+						skyToDetector, combineRegion, bad->data[i],
+						PS_INTERPOLATE_BILINEAR);
+	transformedErr->data[i] = psImageTransform(transformedErr->data[i], NULL, inputsErr->data[i],
+						   inputMask->data[i], inputMaskVal | PS_MASK_COSMICRAY,
+						   NAN, skyToDetector, combineRegion, bad->data[i],
+						   PS_INTERPOLATE_BILINEAR_VARIANCE);
+	psImage *skyImage = transformed->data[i]; // Dereference the transformed image
+	psRegion *blankRegion = psRegionAlloc(0, 0, skyImage->numCols, skyImage->numRows); // Size of
+											   // transformed
+											   // image
+	transformedMask->data[i] = psPixelsToMask(transformedMask->data[i], blanks, blankRegion,
+						  PS_MASK_BLANK);
+	psFree(blankRegion);
+	psFree(blanks);
+    }
+    psFree(bad);
+
+    // Combine with no rejection
+    combined = pmCombineImages(combined, NULL, transformed, transformedErr, transformedMask,
+			       PS_MASK_BLANK, combinePixels, 0, 0.0, combineStats);
+    psFree(combineStats);
+    psFree(combinePixels);
+    psFree(transformed);
+    psFree(transformedErr);
+    psFree(transformedMask);
+\end{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\section{Image Subtraction}
+
+Image subtraction is arguably the best method of identifying faint
+variable sources in images with different point-spread functions.  It
+relies on fitting for a convolution kernel that minimises the
+residuals in subtracting small regions of the image.  The use of a
+convolution kernel consisting of a linear combination of basis
+functions allows the problem to be solved with only modest computing
+power.
+
+\subsection{The kernels}
+
+We will allow for the use of two convolution kernels.  The first is
+that employed by the popular image subtraction program,
+\href{http://www2.iap.fr/users/alard/package.html}{ISIS}, consisting
+of Gaussians modified by polynomials:
+\begin{equation}
+B_{ijk}(u,v) = e^{(u^2 + v^2)/2\sigma_i^2} u^j v^k
+\end{equation}
+The second simply consists of delta functions, which we refer to as
+POIS (Pan-STARRS Optimal Image Subtraction):
+\begin{equation}
+B_{ij}(u,v) = \delta(u - i)\ \delta(v - j)
+\end{equation}
+\tbd{For further details, see the document about image subtraction for
+\PS{}.}  The former is widely used, while the second appears to be
+equally useful and faster, though not as tried and proven.
+
+\begin{verbatim}
+typedef enum {
+    PM_SUBTRACTION_KERNEL_POIS,		// POIS kernel --- delta functions
+    PM_SUBTRACTION_KERNEL_ISIS		// ISIS kernel --- gaussians modified by polynomials
+} pmSubtractionKernelsType;
+\end{verbatim}
+
+In order to simplify the book-keeping for the kernels, we will define
+a \code{pmSubtractionKernels}, which keeps track of the details of the
+each of the kernel basis functions:
+
+\begin{verbatim}
+typedef struct {
+    pmSubtractionKernelType type;	// Type of kernels --- allowing the use of multiple kernels
+    psVector *u, *v;			// Offset (for POIS) or polynomial order (for ISIS)
+    psVector *sigma;			// Width of Gaussian (for ISIS)
+    psVector *xOrder, *yOrder;		// Spatial polynomial order (for all)
+    int subIndex;                       // Index of kernel to be subtracted to maintain flux conservation
+    psArray *preCalc;			// Array of images containing pre-calculated kernel (to
+                                        // accelerate ISIS; don't use for POIS)
+} pmSubtractionKernels;
+\end{verbatim}
+
+This structure caters for both choices of kernel type.  For a POIS
+kernel, the \code{u} and \code{v} vectors shall be set to the
+coordinates for the delta functions for the corresponding kernel.  For
+an ISIS kernel, the \code{sigma} vector shall be set to the Gaussian
+widths and the \code{u} and \code{v} vectors shall be set to the
+orders of the modifying polynomials for the corresponding kernel.  For
+both choices of kernel, the \code{xOrder} and \code{yOrder} vectors
+specify the order of the spatial variation.
+
+In order to maintain flux conservation when the kernel is spatially
+variable, we need to treat one kernel in the set differently.  The
+convolutions for this kernel, identified by the \code{subIndex}, are
+calculated in the usual way, while all others have the \code{subIndex}
+kernel subtracted from them.  For details, see the
+\href{http://www.edpsciences.org/journal/index.cfm?v_url=aas/full/2000/11/ds8706/ds8706.html}{paper
+by Alard (2000, A\&AS, 144, 363)}.
+
+Since the ISIS kernels are continuous functions, it is worth
+pre-calculating them instead of calculating them each time they are
+required.  The \code{preCalc} array, consisting of \code{psImage}s is
+provided for this purpose.
+
+The \code{pmSubtractionKernels} are generated by the following functions:
+
+\begin{verbatim}
+pmSubtractionKernels *pmSubtractionKernelsAllocPOIS(int size, int spatialOrder);
+pmSubtractionKernels *pmSubtractionKernelsAllocISIS(const psVector *sigmas, const psVector *orders,
+                                                    int size, int spatialOrder);
+\end{verbatim}
+
+\code{pmSubtractionKernelsAllocPOIS} shall generate the
+\code{pmSubtractionKernels} suitable for the POIS kernel basis set.
+This involves setting the \code{u}, \code{v}, \code{xOrder} and
+\code{yOrder} to the appropriate values.  \code{size} is the half-size
+of the kernel, and \code{spatialOrder} is the maximum spatial order
+(the spatial variation is $x^i y^j$ with $i+j <$ \code{spatialOrder}).
+The \code{subIndex} is set to the kernel which has \code{u = 0},
+\code{v = 0}, \code{xOrder = 0} and \code{yOrder = 0}.  There should
+be \code{(2 * size + 1) * (2 * size + 1) * (spatialOrder + 1) *
+(spatialOrder + 2) / 2} kernels.
+
+\code{pmSubtractionKernelsAllocISIS} shall generate the
+\code{pmSubtractionKernels} suitable for the ISIS kernel basis set.
+This involves setting the \code{sigma}, \code{u}, \code{v},
+\code{xOrder} and \code{yOrder} to the appropriate values, as well as
+generating the \code{preCalc} images.  Note that the \code{sigma}
+vector contained within the \code{pmSubtractionKernels} is not the
+same as the input \code{sigmas} vector, but contains repeated entries.
+\code{size} is the half-size of the kernel, which specifies the size
+of the \code{preCalc} images.  The \code{spatialOrder} is the maximum
+spatial order (the spatial variation is $x^i y^j$ with $i+j <$
+\code{spatialOrder}).  The \code{subIndex} is set to the kernel which
+has \code{u = 0}, \code{v = 0}, \code{xOrder = 0} and \code{yOrder =
+0}, for the first of the gaussian widths in the \code{sigmas} vector.
+
+\subsection{Stamps}
+
+Sub-regions on an image which are used to derive the best-fit
+convolution kernel are referred to as ``stamps''.
+
+\begin{verbatim}
+typedef struct {
+    int x, y;				// Position
+    psImage *matrix;			// Associated matrix
+    psVector *vector;			// Associated vector
+    psStampStatus status;		// Status of stamp
+} pmStamp;
+\end{verbatim}
+
+A stamp is the region around a central pixel, \code{x,y}.  The
+\code{matrix} and \code{vector} are generated in the process of
+solving for the best-fit convolution kernel; each of these will likely
+be of type \code{psF64} in order to maintain the best possible
+precision (we will be summing squares).  In order to allow us to throw
+out stamps without having to laboriously recompute the total
+least-squares matrix and vector, we use a separate matrix and vector
+for each stamp.
+
+To allow iteration on the choice of stamps, a stamp contains a
+\code{status}, an enumerated type:
+
+\begin{verbatim}
+typedef enum {
+    PM_STAMP_USED,			// Use this stamp
+    PM_STAMP_REJECTED,			// This stamp has been rejected
+    PM_STAMP_RECALC,			// Having been reset, this stamp needs to be recalculated
+    PM_STAMP_NONE			// No stamp in this region
+} pmStampStatus;
+\end{verbatim}
+
+\begin{verbatim}
+psArary *pmSubtractionFindStamps(psArray *stamps, // Output stamps, or NULL
+				 const psImage *image, // Image for which to find stamps
+				 const psImage *mask, // Mask
+				 unsigned int maskVal, // Value for mask
+				 float threshold, // Threshold for stamps in the image
+				 int xNum, int yNum, // Number of stamps in x and y
+				 int border // Border around image to ignore (should be size of kernel)
+				 );
+\end{verbatim}
+
+\code{pmSubtractionFindStamps} returns an array of stamps on the
+\code{image} suitable for use in calculating the best-fit convolution
+kernel.  Except for a \code{border} all the way around, the
+\code{image} is broken into \code{xNum} $\times$ \code{yNum}
+rectangles; there will be a stamp within each rectangle.  If
+\code{stamps} is non-\code{NULL}, then the function shall only attempt
+to identify a new stamp in a particular rectangle if the corresponding
+stamp \code{status} is \code{PM_STAMP_REJECTED}.
+
+A stamp shall be recognised as the pixel with the greatest value that
+does not have the corresponding pixel in the \code{mask} matching
+\code{maskVal}.  If the value of the this pixel does not exceed
+\code{threshold}, then the stamp \code{status} shall be marked as
+\code{PM_STAMP_NONE}, which means that the stamp will be ignored in
+future iterations.  If a legitimate stamp is found within the region,
+then its status shall be changed to \code{PM_STAMP_RECALC}.
+
+
+\subsection{Solving for the kernel}
+
+Calculating the best-fit convolution kernel requires solving a matrix
+equation, the elements of which are obtained by applying the kernel
+basis functions to the stamps.  The final matrix and vector are the
+sum of the matrices and vectors obtained for each of the individual
+stamps.
+
+\begin{verbatim}
+bool pmSubtractionCalculateEquation(psArray *stamps, // The stamps for which to calculate the equation
+				    const psImage *reference, // Reference image
+				    const psImage *input, // Input image
+				    const psSubtractionKernels *kernels, // The kernel basis functions
+				    int footprint // Half-size of region over which to calculate equation
+				    );
+\end{verbatim}
+
+\code{pmSubtractionCalculateEquation} shall calculate the
+\code{matrix} and \code{vector} for each of the \code{stamps} which
+have \code{status} set to \code{PM_STAMP_RECALC}.  The calculation is
+made over a region with a half size of \code{footprint} on the
+\code{reference} and \code{input} images, using each of the
+\code{kernels}.  In the event that any of the input pointers are
+\code{NULL}, the function shall generate an error and return
+\code{false}; otherwise, the function shall return \code{true}.
+
+The vector is:
+\begin{equation}
+v_i = \sum_{x,y} I(x,y) [ R(x,y) \otimes B_i(u,v) ] / \sigma(x,y)^2
+\end{equation}
+and the matrix is:
+\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}
+where $I(x,y)$ is the input image, $R(x,y)$ is the reference image,
+$B_i(u,v)$ is the $i$-th kernel basis function, $\otimes$ denotes
+convolution, $\sigma(x,y) = R(x,y)^{1/2}$ is an estimate of the error,
+and the sum over $x,y$ indicates summing over the stamp regions.
+
+In addition to the each of the \code{kernels}, an additional parameter
+for which we must solve is the difference in the background level
+between the \code{reference} and \code{input} images.  The appropriate
+term shall be added to the \code{matrix} and \code{vector}.
+
+In order to maintain flux conservation when the kernel is spatially
+variable, for each of the kernel basis functions apart from the first,
+the kernel actually employed shall be the first kernel function
+subtracted from the original kernel function.
+
+Having calculated the matrix equation for a stamp, its \code{status}
+is set to \code{PM_STAMP_USED}.
+
+Since this step is one of the major rate-limiting factors in image
+subtraction, care should be taken with optimisation.
+
+\begin{verbatim}
+psVector *pmSubtractionSolveEquation(psVector *solution,	// Solution vector, or NULL
+				     const psArray *stamps // Array of stamps
+				     );
+\end{verbatim}
+
+\code{pmSubtractionSolveEquation} shall solve the matrix equation
+provided by each of the \code{stamps}, returning the \code{solution}
+vector.  This involves summing the \code{matrix} and \code{vector} of
+each of the stamps which have \code{status} set to
+\code{PM_STAMP_USED}, and multiplying the inverse of the matrix by the
+\code{vector}.  If the \code{solution} is \code{NULL}, then the
+function shall allocate and return a new vector; otherwise, the
+\code{solution} vector shall be modified in-place.  If \code{stamps}
+is \code{NULL}, then the function shall generate an error and return
+\code{NULL}.  The type of the \code{solution} vector should be
+\code{psF64}, since the matrix equation involves summing squares.
+
+
+\subsection{Rejection of stamps}
+
+\begin{verbatim}
+bool pmSubtractionRejectStamps(psArray *stamps, // Array of stamps to check for rejection
+			       psImage *mask, // Mask image
+			       unsigned int badStampMaskVal, // Value to use in mask for bad stamp
+			       int footprint, // Region to mask if stamp is bad
+			       float sigmaRej, // Number of RMS deviations above zero at which to reject
+			       const psImage *refImage, // Reference image
+			       const psImage *inImage, // Input image
+			       const psVector *solution, // Solution vector
+			       const pmSubtractionKernels *kernels // Array of kernel parameters
+			       );
+\end{verbatim}
+
+\code{pmSubtractionRejectStamps} shall apply the \code{solution} to
+the \code{stamps}, rejecting stamps for which the mean square
+residuals exceed \code{sigmaRej} RMS deviations from zero.
+\code{stamps} which are rejected have their \code{status} set to
+\code{PM_STAMP_REJECTED}, and have pixels within \code{footprint} of
+the corresponding position in the \code{mask} set to
+\code{badStampMaskVal} so they will not be used again.
+
+The deviations are calculated through extracting the stamps from the
+\code{refImage} and \code{inImage}, convolving the reference stamp by
+the best-fit kernel (derived from the \code{solutions} vector and the
+\code{kernels}), subtracting and then dividing by the stamp from the
+input image, and then squaring to obtain the mean square residual.
+
+\subsection{Visualisation of kernel}
+
+Having solved for the best-fit kernel, it is often useful to visualise
+it.
+
+\begin{verbatim}
+psImage *pmSubtractionKernelImage(psImage *out, const psVector *solution,
+                                  const pmSubtractionKernels *kernels, float x, float y);
+\end{verbatim}
+
+\code{pmSubtractionKernelImage} shall create an image of the kernel
+from the \code{solution} vector and the \code{kernels}.  The relative
+position (between -1 and +1) on the image at which to evaluate the
+kernel (important if the kernel is spatially variable) is specified by
+\code{x} and \code{y}.  If \code{out} is \code{NULL}, then the
+function shall allocate a new image of sufficient size, and return the
+result; otherwise, \code{out} shall be modified in-place.
+
+
+\subsection{Example}
+
+Here is an example of what the image subtraction routine looks like,
+demonstrating how the various pieces fit together.  The inputs are:
+\begin{itemize}
+\item \code{psImage *reference}: Reference image
+\item \code{psImage *refMask}: Mask for reference image
+\item \code{psImage *input}: Input image
+\item \code{psImage *inMask}: Mask for input image
+\item \code{unsigned int maskVal}: Value to be masked
+\item \code{pmSubtractionKernelType kernelType}: Type of kernel to use
+\item \code{int kernelHalfSize}: Half the kernel size (full size is \code{2*kernelHalfSize + 1})
+\item \code{psVector *sigmas}: Widths for the ISIS Gaussians
+\item \code{psVector *polyOrders}: Polynomial orders for ISIS Gaussians
+\item \code{int spatialOrder}: Maximum spatial order for spatially variable kernel
+\item \code{float stampThreshold}: Threshold for finding stamps
+\item \code{int nStampsX, nStampsY}: Number of stamps in x and y
+\item \code{int stampSize}: Half size of stamp footprint
+\item \code{int numIter}: Number of iterations on the stamps
+\item \code{float sigmaRej}: Rejection threshold for stamps
+\end{itemize}
+
+The output is the subtracted image and the corresponding mask.
+
+\begin{verbatim}
+    // Mask around bad pixels in the reference image.  There are two cases to worry about:
+    // 1. Bad pixels within the kernel, which will affect the subtracted image
+    // 2. Bad pixels within the stamp, which affects the calculation of the kernel
+    psImage *subMask = psImageGrowMask(NULL, refMask, maskVal, kernelHalfSize, PS_MASK_NEAR_BAD);
+    (void)psImageGrowMask(subMask, refMask, maskVal, stampSize, PS_MASK_BAD_STAMP);
+    // Add in the mask for the input image.  Don't need to grow this, since it isn't convolved.
+    (void)psBinaryOp(subMask, subMask, "|", inMask);
+
+    // Generate kernel basis functions
+    psArray *kernels = NULL;		// Array of kernel basis functions
+    switch (kernelType) {
+      case PM_SUBTRACTION_KERNEL_POIS:
+	// Create the kernel basis functions
+	kernels = pmSubtractionKernelsGeneratePOIS(kernelHalfSize, spatialOrder);
+	break;
+      case PM_SUBTRACTION_KERNEL_ISIS:
+	kernels = pmSubtractionKernelsGenerateISIS(sigmas, polyOrders, kernelHalfSize, spatialOrder);
+	break;
+      default:
+	barf();
+    }
+
+    psArray *stamps = NULL;		// Array of stamps
+    psVector *kernelCoeffs = NULL;	// Coefficients for the kernels
+    bool rejected = true;		// Did we reject a stamp in the last iteration?
+
+    // Iterate for a solution
+    for (int iter = 0; iter < numIter && rejected; iter++) {
+
+	// Find stamps
+	stamps = pmSubtractionFindStamps(stamps, reference, subMask, maskVal | PS_MASK_BAD_STAMP,
+					 stampThreshold, nStampsX, nStampsY, stampSize, kernelHalfSize);
+
+	// Generate and solve matrix equations
+	(void)pmSubtractionCalculateEquation(stamps, reference, input, kernels, stampSize);
+	kernelCoeffs = pmSubtractionSolveEquation(kernelCoeffs, stamps);
+
+	// Reject bad stamps
+	rejected = pmSubtractionRejectStamps(stamps, subMask, PS_MASK_BAD_STAMP, stampSize, sigmaRej,
+					     reference, input, kernelCoeffs, kernels);
+    }
+
+    // Convolve the reference image
+    psImage *referenceConvolved = pmSubtractionConvolveImage(NULL, reference, subMask, kernelCoeffs, kernels);
+    // Subtract
+    psImage *subtracted = (psImage*)psBinaryOp(NULL, input, "-", referenceConvolved);
+
+    // What does the kernel look like?
+    psImage *kernelImage = pmSubtractionKernelImage(NULL, kernelCoeffs, kernels, 0.0, 0.0);
+    // Check/save kernel image, print statistics....
+
+    psFree(referenceConvolved);
+    psFree(stamps);
+    psFree(kernels);
+    psFree(kernelCoeffs);
+\ebd{verbatim}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 \section{Revision Change Log}
 \input{ChangeLogSDRS.tex}