Changeset 15443

Timestamp:

Nov 2, 2007, 4:28:24 PM (19 years ago)

Author:

Paul Price

Message:

Adding function pmSubtractionOrder (choice of name is probably not the best) to determine which of the two images under consideration should be convolved to match the other. Originally was doing this by solving for a RING kernel of width 1, going both ways, and comparing the deviations for each. However, when doing stacks this didn't work (convolving the wider fake Gaussian image gave smaller deviations than convolving the narrower input image), so have settled on measuring the second moments for each stamp, and using the ratio of moments between the two images to determine which is wider; this seems to work for both subtracting and stacking images. In the process, plugged a few memory leaks, and added code to support convolving either way (or both ways; this will be useful when it comes time to code the dual convolution algorithm).

Location:

trunk/psModules/src/imcombine

Files:

• pmStackReject.c (modified) (3 diffs)
• pmSubtraction.c (modified) (25 diffs)
• pmSubtraction.h (modified) (7 diffs)
• pmSubtractionKernels.h (modified) (1 diff)
• pmSubtractionMatch.c (modified) (17 diffs)
• pmSubtractionMatch.h (modified) (3 diffs)
• pmSubtractionParams.c (modified) (19 diffs)
• pmSubtractionParams.h (modified) (1 diff)
• pmSubtractionStamps.c (modified) (21 diffs)
• pmSubtractionStamps.h (modified) (5 diffs)

trunk/psModules/src/imcombine/pmStackReject.c

@@ -52 +52 @@
 
     // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding
-    psImage *convolved = NULL;          // Convolved image
+    pmReadout *convRO = pmReadoutAlloc(NULL); // Readout with convolved image
+    pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
+    inRO->image = image;
     for (int i = 0; i < numRegions; i++) {
         psRegion *region = regions->data[i]; // Region of interest
         psVector *solution = solutions->data[i]; // Solution of interest
-        if (!pmSubtractionConvolve(&convolved, NULL, NULL, image, NULL, NULL, 0,
-                                   region, solution, kernels, true)) {
+        if (!pmSubtractionConvolve(convRO, inRO, NULL, NULL, 0, region, solution, kernels, true,
+                                   PM_SUBTRACTION_MODE_1)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
-            psFree(convolved);
-            psFree(image);
+            psFree(convRO);
+            psFree(inRO);
             return NULL;
         }
@@ -73 +75 @@
         if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
-            psFree(convolved);
-            psFree(image);
+            psFree(convRO);
+            psFree(inRO);
             return NULL;
         }
@@ -85 +87 @@
         psFree(kernel);
 
-        psImage *subConv = psImageSubset(convolved, *region); // Sub-image of convolved image
+        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
         psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
-    psFree(image);
+    psFree(inRO);
+    psImage *convolved = psMemIncrRefCounter(convRO->image);
+    psFree(convRO);
 
     // Threshold the convolved image

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.71 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-11-01 01:11:03 $
+ *  @version $Revision: 1.72 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-11-03 02:28:24 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -76 +76 @@
 
 // Generate the kernel to apply to the variance from the normal kernel
-static psKernel *varianceKernel(psKernel *normalKernel // Normal kernel
+static psKernel *varianceKernel(psKernel *out, // Output kernel
+                                psKernel *normalKernel // Normal kernel
                                 )
 {
@@ -83 +84 @@
     int yMin = normalKernel->yMin, yMax = normalKernel->yMax;
 
-    psKernel *kernel = psKernelAlloc(xMin, xMax, yMin, yMax); // Kernel to return
+    if (!out) {
+        out = psKernelAlloc(xMin, xMax, yMin, yMax);
+    }
 
     // Take the square of the normal kernel
@@ -93 +96 @@
             sumNormal += value;
             sumVariance += value2;
-            kernel->kernel[v][u] = value2;
+            out->kernel[v][u] = value2;
         }
     }
@@ -99 +102 @@
     // Normalise so that the sum of the variance kernel is the square of the sum of the normal kernel
     // This is required to keep the relative scaling between the image and the weight map
-    psBinaryOp(kernel->image, kernel->image, "*",
-               psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
-
-    return kernel;
+    psBinaryOp(out->image, out->image, "*", psScalarAlloc(PS_SQR(sumNormal) / sumVariance, PS_TYPE_F32));
+
+    return out;
 }
 
@@ -264 +266 @@
 }
 
-// Generate the convolved reference image
-static psKernel *convolveRef(const pmSubtractionKernels *kernels, // Kernel basis functions
-                             int index, // Kernel basis function index
-                             const psKernel *image, // Image to convolve (a kernel for convenience)
-                             int footprint // Size of region of interest
-                             )
+
+// Convolve an image using FFT
+static void convolveFFT(psImage *target,// Place the result in here
+                        const psImage *image, // Image to convolve
+                        const psImage *mask, // Mask image
+                        psMaskType maskVal, // Value to mask
+                        const psKernel *kernel, // Kernel by which to convolve
+                        psRegion region,// Region of interest
+                        float background, // Background to add
+                        int size        // Size of (square) kernel
+                        )
+{
+    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
+                                  region.y0 - size, region.y1 + size); // Add a border
+
+    // Casting away const so psImageSubset can add the child
+    psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
+    psImage *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask
+    psImage *convolved = psImageConvolveFFT(subImage, subMask, maskVal, kernel, 0.0); // Convolution
+    psFree(subImage);
+    psFree(subMask);
+
+    // Now, we have to chop off the borders, and stick it in where it belongs
+    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
+    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
+    if (background != 0.0) {
+        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
+    } else {
+        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+        for (int y = 0; y < subTarget->numRows; y++) {
+            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
+        }
+    }
+    psFree(subConv);
+    psFree(convolved);
+    psFree(subTarget);
+    return;
+}
+
+// Convolve an image directly
+static void convolveDirect(psImage *target, // Put the result here
+                           const psImage *image, // Image to convolve
+                           const psKernel *kernel, // Kernel by which to convolve
+                           psRegion region,// Region of interest
+                           float background, // Background to add
+                           int size        // Size of (square) kernel
+                           )
+{
+    for (int y = region.y0; y < region.y1; y++) {
+        for (int x = region.x0; x < region.x1; x++) {
+            target->data.F32[y][x] = background;
+            for (int v = -size; v <= size; v++) {
+                for (int u = -size; u <= size; u++) {
+                    target->data.F32[y][x] += kernel->kernel[v][u] *
+                        image->data.F32[y - v][x - u];
+                }
+            }
+        }
+    }
+    return;
+}
+
+// Mark a pixel as blank in the image, mask and weight
+static inline void markBlank(psImage *image, // Image to mark as blank
+                             psImage *mask, // Mask to mark as blank (or NULL)
+                             psImage *weight, // Weight map to mark as blank (or NULL)
+                             int x, int y, // Coordinates to mark blank
+                             psMaskType blank // Blank mask value
+    )
+{
+    image->data.F32[y][x] = NAN;
+    if (mask) {
+        mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+    }
+    if (weight) {
+        weight->data.F32[y][x] = NAN;
+    }
+    return;
+}
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Semi-public functions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+// Generate the convolution given a precalculated kernel
+psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel)
+{
+    PS_ASSERT_KERNEL_NON_NULL(image, NULL);
+    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
+
+    psImage *conv = psImageConvolveFFT(image->image, NULL, 0, kernel, 0.0); // Convolved image
+    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
+    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
+    psFree(conv);
+    return convolved;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+// Public functions
+//////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psMaskType maskVal,
+                           int size, int footprint, float badFrac, bool useFFT)
+{
+    PS_ASSERT_IMAGE_NON_NULL(mask1, NULL);
+    PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_MASK, NULL);
+    if (mask2) {
+        PS_ASSERT_IMAGE_NON_NULL(mask2, NULL);
+        PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_MASK, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL);
+    }
+    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
+    if (isfinite(badFrac)) {
+        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL);
+        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL);
+    }
+
+    int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images
+
+    // Dereference inputs for convenience
+    psMaskType **data1 = mask1->data.PS_TYPE_MASK_DATA;
+    psMaskType **data2 = NULL;
+    if (mask2) {
+        data2 = mask2->data.PS_TYPE_MASK_DATA;
+    }
+
+    // First, a pass through to determine the fraction of bad pixels
+    if (isfinite(badFrac) && badFrac != 1.0) {
+        int numBad = 0;                 // Number of bad pixels
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if (data1[y][x] & maskVal) {
+                    numBad++;
+                    continue;
+                }
+                if (data2 && data2[y][x] & maskVal) {
+                    numBad++;
+                }
+            }
+        }
+        if (numBad > badFrac * numCols * numRows) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
+                    numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac);
+            return NULL;
+        }
+    }
+
+    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
+    psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask
+    psImageInit(mask, 0);
+    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience
+
+    // Block out a border around the edge of the image
+
+    // Bottom stripe
+    for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) {
+        for (int x = 0; x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+    // Either side
+    for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) {
+        for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+        for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+    // Top stripe
+    for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+
+    // XXX Could do something smarter here --- we will get images that are predominantly masked (where the
+    // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad
+    // pixel is wasting time.  As a first cut, I've put in a check on the fraction of bad pixels, but we could
+    // imagine looking for the edge of big regions and convolving just at the edge.  As a second cut, allow
+    // use of FFT convolution.
+
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            if (data1[y][x] & maskVal) {
+                maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1;
+            }
+            if (data2 && data2[y][x] & maskVal) {
+                maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2;
+            }
+        }
+    }
+
+    // Block out the entire stamp footprint around bad input pixels.
+
+    // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad
+    // reference pixel (within 'size').  Then we want to block out with the FOOTPRINT mask everything within a
+    // footprint's distance of those (within 'footprint').
+
+    if (useFFT) {
+        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2,
+                                    PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1,
+                                    -size, size, -size, size, 0.5)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_2,
+                                    PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2,
+                                    -size, size, -size, size, 0.5)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1,
+                                    PM_SUBTRACTION_MASK_FOOTPRINT_1,
+                                    -footprint, footprint, -footprint, footprint, 0.5)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2,
+                                    PM_SUBTRACTION_MASK_FOOTPRINT_2,
+                                    -footprint, footprint, -footprint, footprint, 0.5)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2.");
+            psFree(mask);
+            return NULL;
+        }
+    } else {
+        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_1,
+                                       PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1,
+                                       -size, size, -size, size)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_2,
+                                       PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2,
+                                       -size, size, -size, size)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1,
+                                       PM_SUBTRACTION_MASK_FOOTPRINT_1,
+                                       -footprint, footprint, -footprint, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1.");
+            psFree(mask);
+            return NULL;
+        }
+        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2,
+                                       PM_SUBTRACTION_MASK_FOOTPRINT_2,
+                                       -footprint, footprint, -footprint, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2.");
+            psFree(mask);
+            return NULL;
+        }
+    }
+
+    return mask;
+}
+
+// Convolve a stamp by a single kernel basis function
+static psKernel *convolveStampSingle(const pmSubtractionKernels *kernels, // Kernel basis functions
+                                     int index, // Kernel basis function index
+                                     const psKernel *image, // Image to convolve (a kernel for convenience)
+                                     int footprint // Size of region of interest
+    )
 {
     switch (kernels->type) {
@@ -326 +593 @@
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
-          psKernel *convolved = psKernelAlloc(-footprint, footprint,
-                                              -footprint, footprint); // Convolved image
           if (index == kernels->subIndex) {
               // Simply copying over the image
               return convolveOffset(image, 0, 0, footprint);
           }
-
+          psKernel *convolved = psKernelAlloc(-footprint, footprint,
+                                              -footprint, footprint); // Convolved image
           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
           psVector *uCoords = preCalc->data[0]; // u coordinates
@@ -357 +623 @@
 }
 
-// Convolve an image using FFT
-static void convolveFFT(psImage *target,// Place the result in here
-                        const psImage *image, // Image to convolve
-                        const psImage *mask, // Mask image
-                        const psKernel *kernel, // Kernel by which to convolve
-                        psRegion region,// Region of interest
-                        float background, // Background to add
-                        int size        // Size of (square) kernel
-                        )
-{
-    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
-                                  region.y0 - size, region.y1 + size); // Add a border
-
-    // Casting away const so psImageSubset can add the child
-    psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
-    psImage *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask
-    psImage *convolved = psImageConvolveFFT(subImage, subMask, PM_SUBTRACTION_MASK_REF,
-                                            kernel, 0.0); // Convolution
-    psFree(subImage);
-    psFree(subMask);
-
-    // Now, we have to chop off the borders, and stick it in where it belongs
-    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges
-    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
-    if (background != 0.0) {
-        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
-    } else {
-        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
-        for (int y = 0; y < subTarget->numRows; y++) {
-            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
-        }
-    }
-    psFree(subConv);
-    psFree(convolved);
-    psFree(subTarget);
-    return;
-}
-
-// Convolve an image directly
-static void convolveDirect(psImage *target, // Put the result here
-                           const psImage *image, // Image to convolve
-                           const psKernel *kernel, // Kernel by which to convolve
-                           psRegion region,// Region of interest
-                           float background, // Background to add
-                           int size        // Size of (square) kernel
-                           )
-{
-    for (int y = region.y0; y < region.y1; y++) {
-        for (int x = region.x0; x < region.x1; x++) {
-            target->data.F32[y][x] = background;
-            for (int v = -size; v <= size; v++) {
-                for (int u = -size; u <= size; u++) {
-                    target->data.F32[y][x] += kernel->kernel[v][u] *
-                        image->data.F32[y - v][x - u];
-                }
-            }
-        }
-    }
-    return;
-}
-
-// Mark a pixel as blank in the image, mask and weight
-static inline void markBlank(psImage *image, // Image to mark as blank
-                             psImage *mask, // Mask to mark as blank (or NULL)
-                             psImage *weight, // Weight map to mark as blank (or NULL)
-                             int x, int y, // Coordinates to mark blank
-                             psMaskType blank // Blank mask value
+// Convolve the stamp by each of the kernel basis functions
+static psArray *convolveStamp(psArray *convolutions, // The convolutions
+                              const psKernel *image, // Image to convolve
+                              const pmSubtractionKernels *kernels, // Kernel basis functions
+                              int footprint // Stamp half-size
     )
 {
-    image->data.F32[y][x] = NAN;
-    if (mask) {
-        mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-    }
-    if (weight) {
-        weight->data.F32[y][x] = NAN;
-    }
-    return;
-}
-
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Semi-public functions
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-// Generate the convolution given a precalculated kernel
-psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel)
-{
-    PS_ASSERT_KERNEL_NON_NULL(image, NULL);
-    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
-
-    psImage *conv = psImageConvolveFFT(image->image, NULL, 0, kernel, 0.0); // Convolved image
-    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
-    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
-    psFree(conv);
-    return convolved;
-}
-
-
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-// Public functions
-//////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-psImage *pmSubtractionMask(const psImage *refMask, const psImage *inMask, psMaskType maskVal,
-                           int size, int footprint, float badFrac, bool useFFT)
-{
-    PS_ASSERT_IMAGE_NON_NULL(refMask, NULL);
-    PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL);
-    if (inMask) {
-        PS_ASSERT_IMAGE_NON_NULL(inMask, NULL);
-        PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL);
-    }
-    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
-    if (isfinite(badFrac)) {
-        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL);
-        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL);
-    }
-
-    // Size of the images
-    int numCols = refMask->numCols;
-    int numRows = refMask->numRows;
-
-    // Dereference inputs for convenience
-    psMaskType **inData = NULL;
-    if (inMask) {
-        inData = inMask->data.PS_TYPE_MASK_DATA;
-    }
-    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;
-
-    // First, a pass through to determine the fraction of bad pixels
-    if (isfinite(badFrac) && badFrac != 1.0) {
-        int numBad = 0;                 // Number of bad pixels
-        for (int y = 0; y < numRows; y++) {
-            for (int x = 0; x < numCols; x++) {
-                if (inData && inData[y][x] & maskVal) {
-                    numBad++;
-                    continue;
-                }
-                if (refData[y][x] & maskVal) {
-                    numBad++;
-                }
-            }
-        }
-        if (numBad > badFrac * numCols * numRows) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                    "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
-                    numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac);
-            return NULL;
-        }
-    }
-
-    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
-    psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask
-    psImageInit(mask, 0);
-    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience
-
-    // Block out a border around the edge of the image
-
-    // Bottom stripe
-    for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) {
-        for (int x = 0; x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-    // Either side
-    for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) {
-        for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-        for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-    // Top stripe
-    for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
-            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
-        }
-    }
-
-    // XXX Could do something smarter here --- we will get images that are predominantly masked (where the
-    // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad
-    // pixel is wasting time.  As a first cut, I've put in a check on the fraction of bad pixels, but we could
-    // imagine looking for the edge of big regions and convolving just at the edge.  As a second cut, allow
-    // use of FFT convolution.
-
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
-            if (inData && inData[y][x] & maskVal) {
-                maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT;
-            }
-            if (refData[y][x] & maskVal) {
-                maskData[y][x] |= PM_SUBTRACTION_MASK_REF;
-            }
-        }
-    }
-
-    // Block out the entire stamp footprint around bad input pixels.
-
-    // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad
-    // reference pixel (within 'size').  Then we want to block out with the FOOTPRINT mask everything within a
-    // footprint's distance of those (within 'footprint').
-
-    if (useFFT) {
-        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_INPUT, PM_SUBTRACTION_MASK_FOOTPRINT,
-                                    -footprint, footprint, -footprint, footprint, 0.5)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad input pixels.");
-            psFree(mask);
-            return NULL;
-        }
-        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_REF, PM_SUBTRACTION_MASK_CONVOLVE,
-                                    -size, size, -size, size, 0.5)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad reference pixels.");
-            psFree(mask);
-            return NULL;
-        }
-        if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE, PM_SUBTRACTION_MASK_FOOTPRINT,
-                                    -footprint, footprint, -footprint, footprint, 0.5)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad reference pixels.");
-            psFree(mask);
-            return NULL;
-        }
-    } else {
-        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_INPUT, PM_SUBTRACTION_MASK_FOOTPRINT,
-                                       -footprint, footprint, -footprint, footprint)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad input pixels.");
-            psFree(mask);
-            return NULL;
-        }
-        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_REF, PM_SUBTRACTION_MASK_CONVOLVE,
-                                   -size, size, -size, size)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad reference pixels.");
-            psFree(mask);
-            return NULL;
-        }
-        if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE,
-                                       PM_SUBTRACTION_MASK_FOOTPRINT,
-                                       -footprint, footprint, -footprint, footprint)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad reference pixels.");
-            psFree(mask);
-            return NULL;
-        }
-    }
-
-    return mask;
-}
-
-bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
+    assert(image);
+    assert(kernels);
+    assert(footprint >= 0);
+
+    if (convolutions) {
+        // Already done
+        return convolutions;
+    }
+
+    int numKernels = kernels->num;      // Number of kernels
+    convolutions = psArrayAlloc(numKernels);
+
+    for (int i = 0; i < numKernels; i++) {
+        convolutions->data[i] = convolveStampSingle(kernels, i, image, footprint);
+    }
+
+    return convolutions;
+}
+
+
+bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint,
+                                pmSubtractionMode mode)
 {
     PS_ASSERT_PTR_NON_NULL(stamp, false);
-    PS_ASSERT_KERNEL_NON_NULL(stamp->reference, false);
     PS_ASSERT_PTR_NON_NULL(kernels, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
@@ -617 +663 @@
     }
 
-    if (stamp->convolutions) {
-        // Already done
-        return true;
-    }
-
-    stamp->convolutions = psArrayAlloc(kernels->num);
-    psKernel *reference = stamp->reference; // Reference postage stamp
-
-    for (int i = 0; i < kernels->num; i++) {
-        stamp->convolutions->data[i] = convolveRef(kernels, i, reference, footprint);
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_TARGET:
+      case PM_SUBTRACTION_MODE_1:
+        stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint);
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint);
+        break;
+      case PM_SUBTRACTION_MODE_UNSURE:
+      case PM_SUBTRACTION_MODE_DUAL:
+        stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint);
+        stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint);
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
     }
 
@@ -634 +685 @@
 
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
-                                    int footprint)
+                                    int footprint, pmSubtractionMode mode)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -676 +727 @@
         psVectorInit(stampVector, 0.0);
 
-        psKernel *input = stamp->input; // Input postage stamp
-        psKernel *weight = stamp->weight; // Weight map postage stamp
-
-        // Generate convolutions of the reference
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+        // Generate convolutions
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
             return NULL;
         }
 
-        psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel
+        psKernel *weight = stamp->weight; // Weight map postage stamp
+        psKernel *target;               // Target postage stamp
+        psArray *convolutions;          // Convolutions for each kernel
+
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_TARGET:
+          case PM_SUBTRACTION_MODE_1:
+            target = stamp->image2;
+            convolutions = stamp->convolutions1;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            target = stamp->image1;
+            convolutions = stamp->convolutions2;
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
+
         psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
 
@@ -739 +804 @@
             // Vector and background term for matrix
             double sumC = 0.0;      // Sum of the convolution
-            double sumIC = 0.0;     // Sum of the convolution/input product
+            double sumTC = 0.0;     // Sum of the convolution/target product
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
                     double convProduct = jConv->kernel[y][x] / weight->kernel[y][x]; // Convolution / noise^2
                     sumC += convProduct;
-                    sumIC += input->kernel[y][x] * convProduct;
+                    sumTC += target->kernel[y][x] * convProduct;
                 }
             }
@@ -752 +817 @@
                 continue;
             }
-            if (!isfinite(sumIC)) {
+            if (!isfinite(sumTC)) {
                 psStringAppend(&badString, "\nBad sumIC detected at %d", j);
                 bad = true;
@@ -761 +826 @@
                 for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
                      jxOrder++, jIndex += numKernels) {
-                    stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder];
+                    stampVector->data.F64[jIndex] = sumTC * polyValues->data.F64[jyOrder][jxOrder];
 
                     double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value
@@ -774 +839 @@
         // Background only terms
         double sum1 = 0.0;              // Sum of the weighting
-        double sumI = 0.0;              // Sum of the input
+        double sumT = 0.0;              // Sum of the target
         for (int y = - footprint; y <= footprint; y++) {
             for (int x = - footprint; x <= footprint; x++) {
                 double invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse noise, squared
                 sum1 += invNoise2;
-                sumI += input->kernel[y][x] * invNoise2;
+                sumT += target->kernel[y][x] * invNoise2;
             }
         }
@@ -785 +850 @@
             psStringAppend(&badString, "\nBad sum1 detected");
             bad = true;
-        } else if (!isfinite(sumI)) {
+        } else if (!isfinite(sumT)) {
             psStringAppend(&badString, "\nBad sumI detected");
             bad = true;
@@ -791 +856 @@
 
         stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
-        stampVector->data.F64[bgIndex] = sumI;
+        stampVector->data.F64[bgIndex] = sumT;
 
         if (bad) {
@@ -890 +955 @@
 }
 
-
-int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, psImage *subMask, const psVector *solution,
-                              int footprint, float sigmaRej, const pmSubtractionKernels *kernels)
+psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, const psVector *solution,
+                                           int footprint, const pmSubtractionKernels *kernels,
+                                           pmSubtractionMode mode)
+{
+    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
+    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
+    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
+                          polyTerms(kernels->bgOrder), NULL);
+    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, NULL);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
+
+    psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
+    double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations
+    int numKernels = kernels->num;      // Number of kernels
+    int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations
+    float background = solution->data.F64[solution->n-1]; // The difference in background
+
+    psVector *coeff = psVectorAlloc(numKernels, PS_TYPE_F64); // Coefficients for the kernel basis functions
+
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            deviations->data.F32[i] = NAN;
+            continue;
+        }
+
+        // Calculate coefficients of the kernel basis functions
+        psImage *polyValues = spatialPolyValues(kernels->spatialOrder,
+                                                stamp->xNorm, stamp->yNorm); // Polynomial terms
+        for (int j = 0; j < numKernels; j++) {
+            double polynomial = 0.0; // Value of the polynomial
+            for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
+                for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
+                    polynomial += solution->data.F64[index] * polyValues->data.F64[yOrder][xOrder];
+                }
+            }
+            coeff->data.F64[j] = polynomial;
+        }
+        psFree(polyValues);
+
+        // Calculate residuals
+        psKernel *weight = stamp->weight; // Weight postage stamp
+        psKernel *target;               // Target postage stamp
+        psArray *convolutions;          // Convolution postage stamps for each kernel basis function
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_TARGET:
+          case PM_SUBTRACTION_MODE_1:
+            target = stamp->image2;
+            convolutions = stamp->convolutions1;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            target = stamp->image1;
+            convolutions = stamp->convolutions2;
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
+
+#ifdef TESTING
+        psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint);
+#endif
+        float deviation = 0.0;      // Deviation for this stamp
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                float conv = background; // The value of the convolution
+                for (int j = 0; j < numKernels; j++) {
+                    psKernel *convolution = convolutions->data[j]; // Convolution
+                    conv += convolution->kernel[y][x] * coeff->data.F64[j];
+                }
+                float diff = target->kernel[y][x] - conv;
+                deviation += PS_SQR(diff) / weight->kernel[y][x];
+#ifdef TESTING
+                residual->kernel[y][x] = diff;
+#endif
+            }
+        }
+
+        deviations->data.F32[i] = sqrtf(devNorm * deviation);
+        psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
+                i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
+        if (!isfinite(deviations->data.F32[i])) {
+            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            psTrace("psModules.imcombine", 5,
+                    "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
+                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+            continue;
+        }
+
+#ifdef TESTING
+        {
+            psString filename = NULL;
+            psStringAppend(&filename, "resid_%03d.fits", i);
+            psFits *fits = psFitsOpen(filename, "w");
+            psFree(filename);
+            psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
+            psFitsClose(fits);
+        }
+        psFree(residual);
+#endif
+
+    }
+    psFree(coeff);
+
+    return deviations;
+}
+
+
+int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, const psVector *deviations,
+                              psImage *subMask, float sigmaRej, int footprint)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1);
+    PS_ASSERT_VECTOR_NON_NULL(deviations, -1);
+    PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1);
     PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
     PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
-    PS_ASSERT_VECTOR_NON_NULL(solution, -1);
-    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1);
-    PS_ASSERT_PTR_NON_NULL(kernels, -1);
-    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
-                          polyTerms(kernels->bgOrder), -1);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, -1);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, -1);
-
-    // Measure deviations
-    psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps
+
     double totalSquareDev = 0.0;        // Total square deviation from zero
     int numStamps = 0;                  // Number of used stamps
-    int numKernels = kernels->num;      // Number of kernels
-    int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations
-    double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations
-    {
-        float background = solution->data.F64[solution->n-1]; // The difference in background
-
-        for (int i = 0; i < stamps->num; i++) {
-            pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest
-            if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-                continue;
-            }
-
-            psImage *polyValues = spatialPolyValues(kernels->spatialOrder, stamp->xNorm,
-                                                    stamp->yNorm); // Polynomial terms
-
-            psKernel *input = stamp->input; // Input postage stamp
-            psKernel *weight = stamp->weight; // Weight postage stamp
-            psArray *convolutions = stamp->convolutions; // Convolutions of reference image for each kernel
-#ifdef TESTING
-            psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint);
-#endif
-            float deviation = 0.0;      // Deviation for this stamp
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    float conv = background; // The value of the convolution
-                    for (int j = 0; j < numKernels; j++) {
-                        psKernel *convolution = convolutions->data[j]; // Convolution of reference
-
-                        // XXX Precalculate these values, so that we're not calculating the same thing for
-                        // every pixel.
-                        double polynomial = 0.0; // Value of the polynomial
-                        for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
-                            for (int xOrder = 0; xOrder <= spatialOrder - yOrder;
-                                 xOrder++, index += numKernels) {
-                                polynomial += solution->data.F64[index] *
-                                    polyValues->data.F64[yOrder][xOrder];
-                            }
-                        }
-                        conv += convolution->kernel[y][x] * polynomial;
-                    }
-                    float diff = input->kernel[y][x] - conv;
-                    deviation += PS_SQR(diff) / weight->kernel[y][x];
-#ifdef TESTING
-                    residual->kernel[y][x] = diff;
-#endif
-                }
-            }
-            psFree(polyValues);
-
-            deviations->data.F32[i] = devNorm * deviation;
-            if (!isfinite(deviations->data.F32[i])) {
-                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                psTrace("psModules.imcombine", 5,
-                        "Rejecting stamp %d (%d,%d) because of non-finite deviation\n",
-                        i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
-                continue;
-            }
-            psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
-                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
-            totalSquareDev += PS_SQR(deviations->data.F32[i]);
-            numStamps++;
-
-#ifdef TESTING
-            {
-                psString filename = NULL;
-                psStringAppend(&filename, "resid_%03d.fits", i);
-                psFits *fits = psFitsOpen(filename, "w");
-                psFree(filename);
-                psFitsWriteImage(fits, NULL, residual->image, 0, NULL);
-                psFitsClose(fits);
-            }
-            psFree(residual);
-#endif
-
-        }
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            continue;
+        }
+        totalSquareDev += PS_SQR(deviations->data.F32[i]);
+        numStamps++;
     }
 
     if (numStamps == 0) {
         psError(PS_ERR_UNKNOWN, true, "No good stamps found.");
-        psFree(deviations);
         return -1;
     }
 
     if (!isfinite(sigmaRej) || sigmaRej <= 0.0) {
+        // User just wanted to calculate and record the RMS for posterity
         psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation: %f", sqrt(totalSquareDev / numStamps));
-        psFree(deviations);
         return 0;
     }
+
+    float limit = sigmaRej * sqrt(totalSquareDev / (double)numStamps); // Limit on maximum deviation
+    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
 
     int numRejected = 0;                // Number of stamps rejected
     int numGood = 0;                    // Number of good stamps
     double newSquareDev = 0.0;          // New square deviation
-
-    float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation
-    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
-
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
@@ -1024 +1120 @@
                 stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
                 // Recalculate convolutions
-                psFree(stamp->convolutions);
-                stamp->convolutions = NULL;
+                psFree(stamp->convolutions1);
+                psFree(stamp->convolutions2);
+                stamp->convolutions1 = stamp->convolutions2 = NULL;
+                psFree(stamp->image1);
+                psFree(stamp->image2);
+                psFree(stamp->weight);
+                stamp->image1 = stamp->image2 = stamp->weight = NULL;
             } else {
                 numGood++;
@@ -1033 +1134 @@
     }
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
-             numGood, numRejected, sqrt(totalSquareDev / numStamps), sqrt(newSquareDev / numGood));
-
-    psFree(deviations);
+    if (numRejected > 0) {
+        psLogMsg("psModules.imcombine", PS_LOG_INFO,
+                 "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n",
+                 numGood, numRejected, sqrt(totalSquareDev / numStamps),
+                 sqrt(newSquareDev / (double)numGood));
+    } else {
+        psLogMsg("psModules.imcombine", PS_LOG_INFO,
+                 "%d good stamps; 0 rejected.\nRMS deviation: %f\n",
+                 numGood, sqrt(totalSquareDev / numStamps));
+    }
 
     return numRejected;
@@ -1096 +1203 @@
 }
 
-bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
-                           const psImage *inImage, const psImage *inWeight, const psImage *subMask,
-                           psMaskType blank, const psRegion *region, const psVector *solution,
-                           const pmSubtractionKernels *kernels, bool useFFT)
-{
-    PS_ASSERT_IMAGE_NON_NULL(inImage, false);
-    PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false);
+bool pmSubtractionConvolve(pmReadout *out, const pmReadout *ro1, const pmReadout *ro2,
+                           const psImage *subMask, psMaskType blank, const psRegion *region,
+                           const psVector *solution, const pmSubtractionKernels *kernels,
+                           pmSubtractionMode mode, bool useFFT)
+{
+    PS_ASSERT_PTR_NON_NULL(out, false);
+    PS_ASSERT_PTR_NON_NULL(ro1, false);
+    PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
+    PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false);
+    if (ro1->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
+    }
+    PS_ASSERT_PTR_NON_NULL(ro2, false);
+    PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
+    PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false);
+    if (ro2->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false);
+        PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false);
+    }
     if (subMask) {
         PS_ASSERT_IMAGE_NON_NULL(subMask, false);
         PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false);
-    }
-    if (inWeight) {
-        PS_ASSERT_IMAGE_NON_NULL(inWeight, false);
-        PS_ASSERT_IMAGE_TYPE(inWeight, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(inWeight, inImage, false);
-    }
-    PS_ASSERT_PTR_NON_NULL(outImage, false);
-    if (*outImage) {
-        PS_ASSERT_IMAGE_NON_NULL(*outImage, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(*outImage, inImage, false);
-        PS_ASSERT_IMAGE_TYPE(*outImage, PS_TYPE_F32, false);
-    }
-    if (outMask && *outMask) {
-        PS_ASSERT_IMAGE_NON_NULL(*outMask, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false);
-        PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, ro1->image, false);
+    }
+    if (out->image) {
+        PS_ASSERT_IMAGE_NON_NULL(out->image, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out->image, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out->image, PS_TYPE_F32, false);
+    }
+    if (out->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(out->mask, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out->mask, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out->mask, PS_TYPE_MASK, false);
         PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     }
-    if (outWeight && *outWeight) {
-        PS_ASSERT_IMAGE_NON_NULL(*outWeight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(*outWeight, inImage, false);
-        PS_ASSERT_IMAGE_TYPE(*outWeight, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGE_NON_NULL(inWeight, false);
+    if (out->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(out->weight, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out->weight, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out->weight, PS_TYPE_F32, false);
     }
     PS_ASSERT_VECTOR_NON_NULL(solution, false);
@@ -1145 +1260 @@
             return false;
         }
-        if (region->x0 < 0 || region->x1 > inImage->numCols ||
-            region->y0 < 0 || region->y1 > inImage->numRows) {
+        if (region->x0 < 0 || region->x1 > ro1->image->numCols ||
+            region->y0 < 0 || region->y1 > ro1->image->numRows) {
             psString string = psRegionToString(*region);
             psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)",
-                    string, inImage->numCols, inImage->numRows);
+                    string, ro1->image->numCols, ro1->image->numRows);
             psFree(string);
             return false;
@@ -1155 +1270 @@
     }
 
+    psImage *image, *weight;            // Image and weight map to convolve
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_TARGET:
+      case PM_SUBTRACTION_MODE_1:
+        image = ro1->image;
+        weight = ro1->weight;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        image = ro2->image;
+        weight = ro2->image;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+
     int numBackground = polyTerms(kernels->bgOrder); // Number of background terms
     float background = solution->data.F64[solution->n - numBackground]; // The difference in background
-    int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions
-
-    psImage *convImage = *outImage;    // Convolved image
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
+
+    psImage *convImage = out->image;   // Convolved image
     if (!convImage) {
-        *outImage = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        convImage = *outImage;
+        out->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        convImage = out->image;
     }
     psImage *convMask = NULL;           // Convolved mask image
-    if (outMask && subMask) {
-        if (!*outMask) {
-            *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
-        }
-        convMask = *outMask;
+    if (subMask) {
+        if (!out->mask) {
+            out->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
+        }
+        convMask = out->mask;
         psImageInit(convMask, 0);
     }
     psImage *convWeight = NULL;         // Convolved weight (variance) image
-    if (outWeight && inWeight) {
-        if (!*outWeight) {
-            *outWeight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        convWeight = *outWeight;
+    if (weight) {
+        if (!out->weight) {
+            out->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        convWeight = out->weight;
         psImageInit(convWeight, 0.0);
     }
@@ -1198 +1328 @@
     }
 
+    psMaskType maskSource, maskTarget;  // Mask values for source and target
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_TARGET:
+      case PM_SUBTRACTION_MODE_1:
+        maskSource = PM_SUBTRACTION_MASK_BAD_1;
+        maskTarget = PM_SUBTRACTION_MASK_BAD_2 | PM_SUBTRACTION_MASK_CONVOLVE_1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        maskSource = PM_SUBTRACTION_MASK_BAD_2;
+        maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+
     for (int j = yMin; j < yMax; j += fullSize) {
         int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
@@ -1211 +1356 @@
 
             kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues);
-            if (inWeight) {
-                kernelWeight = varianceKernel(kernelImage);
-            }
-
-
+            if (weight) {
+                kernelWeight = varianceKernel(kernelWeight, kernelImage);
+            }
+
+            psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve
             if (useFFT) {
                 // Use Fast Fourier Transform to do the convolution
                 // This provides a big speed-up for large kernels
-                convolveFFT(convImage, inImage, subMask, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
+                convolveFFT(convImage, image, subMask, maskSource, kernelImage, subRegion,
                             background, size);
-                if (inWeight) {
-                    convolveFFT(convWeight, inWeight, subMask, kernelWeight,
-                                psRegionSet(i, xSubMax, j, ySubMax), 0.0, size);
+                if (weight) {
+                    convolveFFT(convWeight, weight, subMask, maskSource, kernelWeight, subRegion,
+                                0.0, size);
                 }
             } else {
-                convolveDirect(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
-                               background, size);
-                if (inWeight) {
-                    convolveDirect(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax),
-                                   0.0, size);
+                convolveDirect(convImage, image, kernelImage, subRegion, background, size);
+                if (weight) {
+                    convolveDirect(convWeight, weight, kernelWeight, subRegion, 0.0, size);
                 }
             }
 
             // Propagate the mask
-            for (int y = j; y < ySubMax; y++) {
-                for (int x = i; x < xSubMax; x++) {
-                    if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
-                                    (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {
-                        convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-                        convImage->data.F32[y][x] = NAN;
-                        if (inWeight) {
-                            convWeight->data.F32[y][x] = NAN;
+            if (subMask) {
+                for (int y = j; y < ySubMax; y++) {
+                    for (int x = i; x < xSubMax; x++) {
+                        if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) {
+                            convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                            convImage->data.F32[y][x] = NAN;
+                            if (weight) {
+                                convWeight->data.F32[y][x] = NAN;
+                            }
                         }
                     }

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -6 +6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-10-17 02:45:40 $
+ * @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-11-03 02:28:24 $
  *
  * Copyright 2004-207 Institute for Astronomy, University of Hawaii
@@ -16 +16 @@
 
 #include <pslib.h>
-#include "pmSubtractionKernels.h"
-#include "pmSubtractionStamps.h"
+
+#include <pmHDU.h>
+#include <pmFPA.h>
+#include <pmSubtractionKernels.h>
+#include <pmSubtractionStamps.h>
 
 /// @addtogroup imcombine Image Combinations
 /// @{
 
+/// Mask values for the subtraction mask
 typedef enum {
-    PM_SUBTRACTION_MASK_CLEAR     = 0x00, // No masking
-    PM_SUBTRACTION_MASK_REF       = 0x01, // Reference image is bad
-    PM_SUBTRACTION_MASK_INPUT     = 0x02, // Input image is bad
-    PM_SUBTRACTION_MASK_CONVOLVE  = 0x04, // If convolved, would be bad
-    PM_SUBTRACTION_MASK_FOOTPRINT = 0x08, // Bad pixel within the stamp footprint
-    PM_SUBTRACTION_MASK_BORDER    = 0x10, // Image border
-    PM_SUBTRACTION_MASK_REJ       = 0x20, // Previously tried as a stamp, and rejected
+    PM_SUBTRACTION_MASK_CLEAR       = 0x00, // No masking
+    PM_SUBTRACTION_MASK_BAD_1       = 0x01, // Image 1 is bad
+    PM_SUBTRACTION_MASK_BAD_2       = 0x02, // Image 2 is bad
+    PM_SUBTRACTION_MASK_CONVOLVE_1  = 0x04, // If image 1 is convolved, would be bad
+    PM_SUBTRACTION_MASK_CONVOLVE_2  = 0x08, // If image 2 is convolved, would be bad
+    PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1
+    PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2
+    PM_SUBTRACTION_MASK_BORDER      = 0x40, // Image border
+    PM_SUBTRACTION_MASK_REJ         = 0x80, // Previously tried as a stamp, and rejected
 } pmSubtractionMasks;
 
@@ -45 +51 @@
 bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve
                                 const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                int footprint ///< Half-size of region over which to calculate equation
+                                int footprint, ///< Half-size of region over which to calculate equation
+                                pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
     );
 
@@ -51 +58 @@
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
                                     const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    int footprint ///< Half-size of region over which to calculate equation
+                                    int footprint, ///< Half-size of region over which to calculate equation
+                                    pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
                                     );
 
@@ -59 +67 @@
                                      );
 
+/// Calculate deviations
+psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps
+                                           const psVector *solution, ///< Solution vector
+                                           int footprint, ///< Half-size of stamp
+                                           const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                           pmSubtractionMode mode ///< Mode for subtraction
+    );
+
 /// Reject stamps
 int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps
+                              const psVector *deviations, ///< Deviations for each stamp
                               psImage *subMask, ///< Subtraction mask
-                              const psVector *solution, ///< Solution vector
-                              int footprint, ///< Region to mask if stamp is bad
                               float sigmaRej, ///< Number of RMS deviations above zero at which to reject
-                              const pmSubtractionKernels *kernels ///< Kernel parameters
+                              int footprint ///< Half-size of stamp
     );
 
@@ -81 +96 @@
 
 /// Convolve image in preparation for subtraction
-bool pmSubtractionConvolve(psImage **outImage, ///< Output image
-                           psImage **outWeight, ///< Output weight map (or NULL)
-                           psImage **outMask, ///< Output mask (or NULL)
-                           const psImage *inImage, ///< Input image
-                           const psImage *inWeight, ///< Input weight map (or NULL)
+bool pmSubtractionConvolve(pmReadout *out, ///< Output image
+                           const pmReadout *ro1, // Input image 1
+                           const pmReadout *ro2, // Input image 2
                            const psImage *subMask, ///< Subtraction mask (or NULL)
                            psMaskType blank, ///< Mask value for blank regions
@@ -91 +104 @@
                            const psVector *solution, ///< The solution vector
                            const pmSubtractionKernels *kernels, ///< Kernel parameters
+                           pmSubtractionMode mode, ///< Mode for subtraction
                            bool useFFT  ///< Use Fast Fourier Transform for the convolution?
     );

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -14 +14 @@
     PM_SUBTRACTION_KERNEL_RINGS,        ///< Rings Instead of the Normal Gaussian Subtraction
 } pmSubtractionKernelsType;
+
+/// Modes --- specifies which image to convolve
+typedef enum {
+    PM_SUBTRACTION_MODE_ERR,            // Error in the mode
+    PM_SUBTRACTION_MODE_TARGET,         // Convolve image 1 to match target PSF
+    PM_SUBTRACTION_MODE_1,              // Convolve image 1
+    PM_SUBTRACTION_MODE_2,              // Convolve image 2
+    PM_SUBTRACTION_MODE_UNSURE,         // Not sure yet which image to convolve so try to satisfy both
+    PM_SUBTRACTION_MODE_DUAL,           // Dual convolution
+} pmSubtractionMode;
 
 /// Kernels specification

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -50 +50 @@
 
 static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read
-                      const pmReadout *reference, // Reference readout
-                      const pmReadout *input, // Input readout, or NULL to generate fake stamps
+                      const pmReadout *ro1, // Readout 1
+                      const pmReadout *ro2, // Readout 2
                       const psImage *subMask, // Mask for subtraction, or NULL
                       psImage *weight,  // Weight map
@@ -57 +57 @@
                       float threshold,  // Threshold for stamp finding
                       float stampSpacing, // Spacing between stamps
-                      float targetWidth,// Target width for fake stamps
                       int size,         // Kernel half-size
-                      int footprint     // Convolution footprint for stamps
+                      int footprint,     // Convolution footprint for stamps
+                      pmSubtractionMode mode // Mode for subtraction
     )
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, reference->image, subMask, region, threshold, stampSpacing);
+    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
     if (!*stamps) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
@@ -71 +71 @@
     memCheck("  find stamps");
 
-    if (!input && !pmSubtractionStampsGenerate(*stamps, targetWidth, footprint, size)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to generate target stamps.");
-        return false;
-    }
-
-    memCheck("   generate stamps");
-
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, reference->image, input ? input->image : NULL,
-                                    weight, footprint, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -94 +86 @@
 
 
-bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *reference, const pmReadout *input,
+bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
-                        const psArray *sources, const char *stampsName, float targetWidth,
+                        const psArray *sources, const char *stampsName,
                         pmSubtractionKernelsType type, int size, int spatialOrder,
                         const psVector *isisWidths, const psVector *isisOrders,
                         int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs,
                         int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad,
-                        psMaskType maskBlank, float badFrac)
+                        psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
 {
     PS_ASSERT_PTR_NON_NULL(convolved, false);
-    PS_ASSERT_PTR_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference->image, false);
-    PS_ASSERT_IMAGE_TYPE(reference->image, PS_TYPE_F32, false);
-    if (reference->mask) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->mask, false);
-        PS_ASSERT_IMAGE_TYPE(reference->mask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->mask, reference->image, false);
-    }
-    if (reference->weight) {
-        PS_ASSERT_IMAGE_NON_NULL(reference->weight, false);
-        PS_ASSERT_IMAGE_TYPE(reference->weight, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(reference->weight, reference->image, false);
-    }
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input->image, false);
-        PS_ASSERT_IMAGE_TYPE(input->image, PS_TYPE_F32, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input->image, reference->image, false);
-        if (input->mask) {
-            PS_ASSERT_IMAGE_NON_NULL(input->mask, false);
-            PS_ASSERT_IMAGE_TYPE(input->mask, PS_TYPE_MASK, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->mask, reference->image, false);
-        }
-        if (input->weight) {
-            PS_ASSERT_IMAGE_NON_NULL(input->weight, false);
-            PS_ASSERT_IMAGE_TYPE(input->weight, PS_TYPE_F32, false);
-            PS_ASSERT_IMAGES_SIZE_EQUAL(input->weight, reference->image, false);
+    PS_ASSERT_PTR_NON_NULL(ro1, false);
+    PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
+    PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false);
+    if (ro1->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false);
+    }
+    if (ro1->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false);
+        PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
+    }
+    if (ro2) {
+        PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
+        PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false);
+        if (ro2->mask) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false);
+        }
+        if (ro2->weight) {
+            PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false);
+            PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false);
+            PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false);
         }
     } else if (!stampsName && !sources) {
@@ -144 +136 @@
     }
     // stampsName may be anything
-    // targetWidth can be just about anything (except maybe negative, but it can be NAN)
     // We'll check kernel type when we allocate the kernels
     PS_ASSERT_INT_POSITIVE(size, false);
@@ -196 +187 @@
     }
 
-    psImage *inImage = NULL, *inMask = NULL; // Input image, mask, weight
-    if (input) {
-        inImage = input->image;
-        inMask = input->mask;
-    }
-
     // Where does our weight map come from?
     psImage *weight = NULL;             // Weight image to use
-    if (input && input->weight) {
-        weight = input->weight;
-    } else if (reference->weight) {
-        weight = reference->weight;
-    } else if (input) {
-        weight = input->image;
+    if (ro1->weight && ro2 && ro2->weight) {
+        weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight);
+    } else if (ro1->weight) {
+        weight = psMemIncrRefCounter(ro1->weight);
+    } else if (ro2) {
+        if (ro2->weight) {
+            weight = psMemIncrRefCounter(ro2->weight);
+        } else {
+            weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image);
+        }
     } else {
-        weight = reference->image;
+        weight = psMemIncrRefCounter(ro1->image);
     }
 
@@ -221 +210 @@
     psVector *solution = NULL;          // Solution to match PSF
     pmSubtractionKernels *kernels = NULL; // Kernel basis functions
-
-    int numCols = reference->image->numCols, numRows = reference->image->numRows; // Image dimensions
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
 
     memCheck("start");
 
-    subMask = pmSubtractionMask(reference->mask, inMask, maskBad, size, footprint, badFrac, useFFT);
+    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint,
+                                badFrac, useFFT);
     if (!subMask) {
         psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask.");
+        psFree(weight);
         return false;
     }
@@ -260 +250 @@
 
             if (sources) {
-                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing,
-                                                           input ? 0 : 2 * footprint);
+                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                // Read stamps from file
-                psTrace("psModules.imcombine", 3, "Reading stamps from %s...\n", stampsName);
-                const char *stampFormat = input ? "%f %f" : "%f %f %f"; // Format for reading stamp file
-                psArray *stampsData = stampsData = psVectorsReadFromFile(stampsName, stampFormat);
-                psVector *xStamp = NULL, *yStamp = NULL, *fluxStamp = NULL; // Stamp positions and fluxes
-                if (!stampsData) {
-                    psError(PS_ERR_IO, false, "Unable to read stamps file %s", stampsName);
-                    goto ERROR;
-                }
-                xStamp = stampsData->data[0];
-                yStamp = stampsData->data[1];
-                if (!input) {
-                    fluxStamp = stampsData->data[2];
-                }
-
-                // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
-                psBinaryOp(xStamp, xStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-                psBinaryOp(yStamp, yStamp, "-", psScalarAlloc(1.0, PS_TYPE_F32));
-
-                stamps = pmSubtractionStampsSet(xStamp, yStamp, fluxStamp, reference->image, subMask,
-                                                region, stampSpacing, input ? 0 : footprint + size);
-                psFree(stampsData);
+                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
+            }
+
+            // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
+            // doesn't matter.
+            if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing,
+                           size, footprint, mode)) {
+                goto MATCH_ERROR;
+            }
+
+            if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) {
+                pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode
+                switch (newMode) {
+                  case PM_SUBTRACTION_MODE_1:
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2.");
+                    break;
+                  case PM_SUBTRACTION_MODE_2:
+                    if (mode == PM_SUBTRACTION_MODE_TARGET) {
+                        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                                "Input PSF is larger than target PSF --- can't match image.");
+                        goto MATCH_ERROR;
+                    }
+                    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1.");
+                    break;
+                  default:
+                    psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order.");
+                    goto MATCH_ERROR;
+                }
+                mode = newMode;
             }
 
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                if (!getStamps(&stamps, reference, input, subMask, weight, NULL,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
-                }
                 kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold);
+                                                          stamps, footprint, optThreshold, mode);
                 if (!kernels) {
                     psErrorClear();
@@ -314 +306 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-                if (!getStamps(&stamps, reference, input, subMask, weight, region,
-                               threshold, stampSpacing, targetWidth, size, footprint)) {
-                    goto ERROR;
+                if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing,
+                               size, footprint, mode)) {
+                    goto MATCH_ERROR;
                 }
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
@@ -331 +323 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
+                    goto MATCH_ERROR;
                 }
 
                 memCheck("  solve equation");
 
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+
+                memCheck("   calculate deviations");
+
                 psTrace("psModules.imcombine", 3, "Rejecting stamps...\n");
-                numRejected = pmSubtractionRejectStamps(stamps, subMask, solution, footprint, rej, kernels);
+                numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint);
                 if (numRejected < 0) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps.");
-                    goto ERROR;
-                }
+                    psFree(deviations);
+                    goto MATCH_ERROR;
+                }
+                psFree(deviations);
+
                 memCheck("  reject stamps");
             }
@@ -350 +354 @@
                 if (!solution) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
-                    goto ERROR;
-                }
-                (void)pmSubtractionRejectStamps(stamps, subMask, solution, footprint, NAN, kernels);
+                    goto MATCH_ERROR;
+                }
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
+                                                                        mode); // Deviations for each stamp
+                if (!deviations) {
+                    psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+                    goto MATCH_ERROR;
+                }
+                (void)pmSubtractionRejectStamps(stamps, deviations, subMask, footprint, NAN);
+                psFree(deviations);
             }
             psFree(stamps);
@@ -372 +383 @@
                             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
 
@@ -380 +391 @@
                             psFree(kernel);
                             psFree(convKernels);
-                            goto ERROR;
+                            goto MATCH_ERROR;
                         }
                         psFree(kernel);
@@ -416 +427 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(&convolved->image, &convolved->weight, &convolved->mask,
-                                       reference->image, reference->weight, subMask, maskBlank, region,
-                                       solution, kernels, useFFT)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to convolve reference image.");
-                goto ERROR;
+            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
+                                       solution, kernels, mode, useFFT)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
+                goto MATCH_ERROR;
             }
             psFree(kernels);
@@ -461 +471 @@
     psFree(subMask);
     subMask = NULL;
+    psFree(weight);
+    weight = NULL;
 
     if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
-        goto ERROR;
+        goto MATCH_ERROR;
     }
 
@@ -472 +484 @@
     return true;
 
-ERROR:
+MATCH_ERROR:
     psFree(region);
     psFree(regionString);
@@ -479 +491 @@
     psFree(stamps);
     psFree(solution);
+    psFree(weight);
     return false;
 }
+
+// Calculate the second order moments for an image
+static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments
+                                    int radius   // Maximum radius
+                                    )
+{
+    assert(kernel && kernel->kernel);
+
+    int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x
+    int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y
+
+    float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment)
+    float sum = 0.0;       // Sum (zero-th moment)
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            xCentroid += kernel->kernel[y][x] * x;
+            yCentroid += kernel->kernel[y][x] * y;
+            sum += kernel->kernel[y][x];
+        }
+    }
+    xCentroid /= sum;
+    yCentroid /= sum;
+
+    float eta20 = 0.0, eta02 = 0.0;     // Second moments
+    for (int y = yMin; y <= yMax; y++) {
+        float yDiff = y - yCentroid;
+        for (int x = xMin; x <= xMax; x++) {
+            float xDiff = x - xCentroid;
+            eta20 += PS_SQR(xDiff) * kernel->kernel[y][x];
+            eta02 += PS_SQR(yDiff) * kernel->kernel[y][x];
+        }
+    }
+
+    // Normalise to calculate the scale-invariant
+    float sum2 = PS_SQR(sum);
+    eta20 /= sum2;
+    eta02 /= sum2;
+    // eta11 /= sum2;
+
+    return eta20 + eta02;
+}
+
+#if 0
+// Calculate the deviations for a particular subtraction order
+static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel
+                                           pmSubtractionStampList *stamps, // Stamps to convolve
+                                           const pmSubtractionKernels *kernels, // Kernel basis functions
+                                           int footprint, // Stamp footprint
+                                           pmSubtractionMode mode // Mode of subtraction
+                                           )
+{
+    assert(stamps);
+    assert(footprint >= 0);
+    assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2);
+
+    if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    psVector *solution = pmSubtractionSolveEquation(NULL, stamps);
+    if (!solution) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
+        return NULL;
+    }
+
+    if (sumKernel) {
+        float sum = 0.0;                // Sum of the kernel
+        psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel
+        for (int y = 0; y < image->numRows; y++) {
+            for (int x = 0; x < image->numCols; x++) {
+                sum += image->data.F32[y][x];
+            }
+        }
+        psFree(image);
+        *sumKernel = sum;
+    }
+
+    psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode);
+    psFree(solution);
+    if (!deviations) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
+        return NULL;
+    }
+
+    return deviations;
+}
+#endif
+
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius)
+{
+    PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR);
+
+    psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps
+    psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments
+    for (int i = 0; i < stamps->num; i++) {
+        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
+        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) {
+            mask->data.PS_TYPE_MASK_DATA[i] = 0xff;
+            continue;
+        }
+        mask->data.PS_TYPE_MASK_DATA[i] = 0;
+        moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) /
+            subtractionOrderMoment(stamp->image2, radius);
+    }
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
+    if (!psVectorStats(stats, moments, NULL, mask, 0xff)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio.");
+        psFree(mask);
+        psFree(moments);
+        psFree(stats);
+        return PM_SUBTRACTION_MODE_ERR;
+    }
+    psFree(moments);
+    psFree(mask);
+
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian);
+    pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2);
+    psFree(stats);
+
+    return mode;
+}
+
+
+

trunk/psModules/src/imcombine/pmSubtractionMatch.h

@@ -4 +4 @@
 #include <pslib.h>
 
-#include "pmHDU.h"
-#include "pmFPA.h"
-#include "pmSubtractionKernels.h"
+#include <pmHDU.h>
+#include <pmFPA.h>
+#include <pmSubtractionKernels.h>
+#include <pmSubtractionStamps.h>
 
-
+/// Match two images
 bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data
-                        const pmReadout *reference, ///< Reference data
-                        const pmReadout *input, ///< Input data
+                        const pmReadout *ro1, ///< Image 1
+                        const pmReadout *ro2, ///< Image 2
                         // Stamp parameters
                         int footprint,  ///< Stamp half-size
@@ -19 +20 @@
                         const psArray *sources, ///< Sources for stamps
                         const char *stampsName, ///< Filename for stamps
-                        float targetWidth, ///< Width of PSF for simulated target
                         // Kernel parameters
                         pmSubtractionKernelsType type, ///< Kernel type
@@ -38 +38 @@
                         psMaskType maskBad, ///< Value to mask
                         psMaskType maskBlank, ///< Mask for blank region
-                        float badFrac   ///< Maximum fraction of bad input pixels to accept
+                        float badFrac,   ///< Maximum fraction of bad input pixels to accept
+                        pmSubtractionMode mode ///< Mode of subtraction
     );
 
+/// Determine which image to convolve
+pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted
+                                     int footprint ///< Stamp half-size
+    );
+
+
 #endif

trunk/psModules/src/imcombine/pmSubtractionParams.c

@@ -50 +50 @@
 #endif
 
+/// Select the appropriate convolution, given the kernel basis function and subtraction mode
+static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp
+                                          int kernelIndex, // Index for kernel component
+                                          pmSubtractionMode mode // Mode of subtraction
+    )
+{
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+        return stamp->convolutions1->data[kernelIndex];
+      case PM_SUBTRACTION_MODE_2:
+        return stamp->convolutions2->data[kernelIndex];
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+    return NULL;                        // Unreached
+}
+
 // Accumulate cross-term sums for a stamp
 static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2
@@ -55 +72 @@
                             double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2
                             const pmSubtractionStamp *stamp, // Stamp with weight
-                            const psKernel *input, // Input image, I(x)
+                            const psKernel *target, // Target stamp
                             int kernelIndex, // Index for kernel component
-                            int footprint // Size of region of interest
+                            int footprint, // Size of region of interest
+                            pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
@@ -82 +100 @@
                                    const pmSubtractionStamp *stamp, // Stamp with input and weight
                                    int kernelIndex, // Index for kernel component
-                                   int footprint // Size of region of interest
+                                   int footprint, // Size of region of interest
+                                   pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
@@ -100 +119 @@
 }
 
-static double accumulateChi2(psKernel *input, // Input stamp
+static double accumulateChi2(const psKernel *target, // Target stamp
                              pmSubtractionStamp *stamp, // Stamp with weight
                              int kernelIndex, // Index for kernel component
                              double coeff, // Coefficient of convolution
                              double bg,  // Background term
-                             int footprint // Size of region of interest
+                             int footprint, // Size of region of interest
+                             pmSubtractionMode mode // Mode of subtraction
     )
 {
     double chi2 = 0.0;
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
 
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
@@ -125 +145 @@
 
 // Return the initial value of chi^2
-static double initialChi2(psKernel *input, // Input stamp
+static double initialChi2(const psKernel *target, // Target stamp
                           const pmSubtractionStamp *stamp, // Stamp with weight
-                          int footprint // Size of convolution
+                          int footprint, // Size of convolution
+                          pmSubtractionMode mode // Mode of subtraction
     )
 {
     psKernel *weight = stamp->weight;   // Weight map
-    psKernel *reference = stamp->reference; // Reference stamp
+    psKernel *source;                   // Source stamp
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+        source = stamp->image1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        source = stamp->image2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
 
     double chi2 = 0.0;                  // Chi^2
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
-        psF32 *ref = &reference->kernel[y][-footprint]; // Derference reference
+        psF32 *ref = &source->kernel[y][-footprint]; // Derference reference
         for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) {
             float diff = *in - *ref;    // Temporary value
@@ -148 +179 @@
 
 // Subtract a convolution from the input
-static void subtractConvolution(psKernel *input, // Input stamp
+static void subtractConvolution(psKernel *target, // Target stamp
                                 const pmSubtractionStamp *stamp, // Stamp with weight
                                 int kernelIndex, // Index for kernel component
                                 float coeff, // Coefficient of subtraction
                                 float bg, // Background term
-                                int footprint // Size of region of interest
-    )
-{
-    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
-
+                                int footprint, // Size of region of interest
+                                pmSubtractionMode mode // Mode of subtraction
+    )
+{
+    psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest
     for (int y = -footprint; y <= footprint; y++) {
-        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *in = &target->kernel[y][-footprint]; // Dereference input
         psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
         for (int x = -footprint; x <= footprint; x++, in++, conv++) {
@@ -173 +204 @@
                                                       int spatialOrder, const psVector *fwhms, int maxOrder,
                                                       const pmSubtractionStampList *stamps, int footprint,
-                                                      float tolerance)
+                                                      float tolerance, pmSubtractionMode mode)
 {
     if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
@@ -210 +241 @@
     // Need to save the stamp inputs --- we're changing the values!
     int numStamps = stamps->num;        // Number of stamps
-    psArray *inputs = psArrayAlloc(numStamps); // Deep copies of the inputs
+    psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets
     psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps
     psVectorInit(badStamps, 0);
@@ -219 +250 @@
             continue;
         }
-        psKernel *input = stamp->input; // Input image of interest
-        psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image
-        inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
+        psKernel *target;               // Target image of interest
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_1:
+            target = stamp->image2;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            target = stamp->image1;
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
+        psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image
+        targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
         psFree(copy);                   // Drop reference
     }
@@ -238 +279 @@
             continue;
         }
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
-            psFree(inputs);
+            psFree(targets);
             psFree(kernels);
             psFree(badStamps);
@@ -258 +299 @@
                     "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n",
                     i, (int)stamp->x, (int)stamp->y);
-            psFree(inputs);
+            psFree(targets);
             psFree(kernels);
             psFree(badStamps);
@@ -265 +306 @@
 
         for (int j = 0; j < numKernels; j++) {
-            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint);
-        }
-
-        lastChi2 += initialChi2(inputs->data[i], stamp, footprint);
+            accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode);
+        }
+
+        lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode);
         numPixels += PS_SQR(2 * footprint + 1);
     }
@@ -297 +338 @@
                 }
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-                accumulateCross(&sumI, &sumII, &sumIC, stamp, inputs->data[j], i, footprint);
+                accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode);
             }
 
@@ -310 +351 @@
                 }
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-                chi2 += accumulateChi2(inputs->data[j], stamp, i, coeff, bg, footprint);
+                chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode);
             }
 
@@ -328 +369 @@
         if (bestIndex == -1) {
             psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter);
-            psFree(inputs);
+            psFree(targets);
             psFree(sumC);
             psFree(sumCC);
@@ -345 +386 @@
             }
             pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
-            subtractConvolution(inputs->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint);
+            subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode);
         }
 
@@ -361 +402 @@
         lastChi2 = bestChi2;
     }
-    psFree(inputs);
+    psFree(targets);
     psFree(sumC);
     psFree(sumCC);
@@ -401 +442 @@
                 pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest
                 psArray *convolutions = convNew->data[j]; // Convolutions for this stamp
-                convolutions->data[rank] = psMemIncrRefCounter(stamp->convolutions->data[i]);
+                convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode));
             }
         }
@@ -420 +461 @@
         }
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        psFree(stamp->convolutions);
-        stamp->convolutions = convNew->data[i];
+        psFree(stamp->convolutions1);
+        psFree(stamp->convolutions2);
+        switch (mode) {
+          case PM_SUBTRACTION_MODE_1:
+            stamp->convolutions1 = convNew->data[i];
+            stamp->convolutions2 = NULL;
+            break;
+          case PM_SUBTRACTION_MODE_2:
+            stamp->convolutions1 = NULL;
+            stamp->convolutions2 = convNew->data[i];
+            break;
+          default:
+            psAbort("Unsupported subtraction mode: %x", mode);
+        }
     }
 

trunk/psModules/src/imcombine/pmSubtractionParams.h

@@ -15 +15 @@
                                                       const pmSubtractionStampList *stamps, ///< Stamps
                                                       int footprint, ///< Convolution footprint for stamps
-                                                      float tolerance ///< Maximum difference in chi^2
+                                                      float tolerance, ///< Maximum difference in chi^2
+                                                      pmSubtractionMode mode // Mode for subtraction
     );
 

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -45 +46 @@
     psFree(stamp->matrix);
     psFree(stamp->vector);
-    psFree(stamp->reference);
-    psFree(stamp->input);
+    psFree(stamp->image1);
+    psFree(stamp->image2);
     psFree(stamp->weight);
-    psFree(stamp->convolutions);
+    psFree(stamp->convolutions1);
+    psFree(stamp->convolutions2);
 }
 
@@ -65 +67 @@
 // Is this position unmasked?
 static bool checkStampMask(int x, int y, // Coordinates of stamp
-                           const psImage *mask
+                           const psImage *mask, // Mask
+                           pmSubtractionMode mode // Mode for subtraction
                            )
 {
@@ -74 +77 @@
         return false;
     }
-    return (mask->data.PS_TYPE_MASK_DATA[y][x] & (PM_SUBTRACTION_MASK_BORDER |
-                                                  PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_REJ)) ?
-        false : true;
+
+    psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value
+    switch (mode) {
+      case PM_SUBTRACTION_MODE_1:
+      case PM_SUBTRACTION_MODE_TARGET:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1;
+        break;
+      case PM_SUBTRACTION_MODE_2:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      case PM_SUBTRACTION_MODE_UNSURE:
+      case PM_SUBTRACTION_MODE_DUAL:
+        maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2;
+        break;
+      default:
+        psAbort("Unsupported subtraction mode: %x", mode);
+    }
+
+    return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true;
 }
 
@@ -140 +159 @@
     stamp->status = PM_SUBTRACTION_STAMP_INIT;
 
-    stamp->reference = NULL;
-    stamp->input = NULL;
+    stamp->image1 = NULL;
+    stamp->image2 = NULL;
     stamp->weight = NULL;
-    stamp->convolutions = NULL;
+    stamp->convolutions1 = NULL;
+    stamp->convolutions2 = NULL;
 
     return stamp;
@@ -151 +171 @@
 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
                                                 const psImage *subMask, const psRegion *region,
-                                                float threshold, float spacing)
+                                                float threshold, float spacing, pmSubtractionMode mode)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -226 +246 @@
             for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
                 for (int x = subRegion->x0; x <= subRegion->y1 ; x++) {
-                    if (checkStampMask(x, y, subMask) && image->data.F32[y][x] > fluxStamp) {
+                    if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
                         fluxStamp = image->data.F32[y][x];
                         xStamp = x;
@@ -242 +262 @@
 
             // Reset the postage stamps since we're making a new stamp
-            psFree(stamp->reference);
-            psFree(stamp->input);
+            psFree(stamp->image1);
+            psFree(stamp->image2);
             psFree(stamp->weight);
-            stamp->reference = stamp->input = stamp->weight = NULL;
-
-            stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
+            psFree(stamp->convolutions1);
+            psFree(stamp->convolutions2);
+            stamp->image1 = stamp->image2 = stamp->weight = NULL;
+            stamp->convolutions1 = stamp->convolutions2 = NULL;
+
+            stamp->status = PM_SUBTRACTION_STAMP_FOUND;
             numFound++;
             psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n",
@@ -266 +289 @@
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
                                                const psImage *image, const psImage *subMask,
-                                               const psRegion *region, float spacing, int exclusionZone)
+                                               const psRegion *region, float spacing, pmSubtractionMode mode)
+
 {
     PS_ASSERT_VECTOR_NON_NULL(x, NULL);
@@ -289 +313 @@
     }
     PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(exclusionZone, NULL);
 
     int numStars = x->n;                // Number of stars
-    psVector *exclude = psVectorAlloc(numStars, PS_TYPE_U8); // Exclude a star?
-    psVectorInit(exclude, 0);
-
-    // Apply exclusion zone around each star --- important when we're convolving to a specified PSF; less so
-    // when we're trying to get a good subtraction.
-    if (exclusionZone > 0) {
-        psTrace("psModules.imcombine", 2, "Applying exclusion zone of %d pixels for stamps\n", exclusionZone);
-        // We use something resembling a binary search --- coordinates are sorted in the x dimension, and then
-        // to exclude stars within a nominated distance, we need only examine (i.e., calculate the
-        // 2-dimensional distance to) other stars in the list that are within that distance of the x
-        // coordinate of the star of interest.  By marking both stars when we find a violation of the
-        // exclusion zone, we need only search upwards from the x coordinate of the star of interest.
-
-        int minDistance2 = PS_SQR(exclusionZone); // Minimum square distance for other stars
-        psVector *indexes = psVectorSortIndex(NULL, x); // Indices for sorting in x
-        for (int i = 0; i < numStars - 1; i++) {
-            int iIndex = indexes->data.S32[i]; // Index for star i
-            if (exclude->data.U8[iIndex]) {
-                continue;
-            }
-            float ix = x->data.F32[iIndex], iy = y->data.F32[iIndex]; // Coordinates for star i
-            float jx, jy;                   // Coordinates for star j
-            for (int j = i + 1, jIndex = indexes->data.S32[j];
-                 j < numStars && (jx = x->data.F32[jIndex]) < ix + exclusionZone;
-                 j++, jIndex = indexes->data.S32[j]) {
-                jy = y->data.F32[jIndex];
-                if (PS_SQR(ix - jx) + PS_SQR(iy - jy) < minDistance2) {
-                    exclude->data.U8[iIndex] = 0xff;
-                    exclude->data.U8[jIndex] = 0xff;
-                    psTrace("psModules.imcombine", 5, "Excluding stamps %d,%d and %d,%d\n",
-                            (int)x->data.F32[iIndex], (int)y->data.F32[iIndex],
-                            (int)x->data.F32[jIndex], (int)y->data.F32[jIndex]);
-                    // Would 'break' here, but there might be more stamps within the exclusion zone.
-                }
-            }
-        }
-        psFree(indexes);
-    }
-
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
                                                                  region, spacing); // Stamp list
@@ -347 +331 @@
     // Put the stars into their appropriate subregions
     for (int i = 0; i < numStars; i++) {
-        if (exclude->data.U8[i]) {
-            // Star has been excluded
-            continue;
-        }
         float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp
         int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp
@@ -357 +337 @@
             continue;
         }
-        if (!checkStampMask(xPix, yPix, subMask)) {
+        if (!checkStampMask(xPix, yPix, subMask, mode)) {
             // Not a good stamp
             continue;
@@ -386 +366 @@
         }
     }
-    psFree(exclude);
 
     // Sort the list by flux, with the brightest last
@@ -421 +400 @@
 
 
-bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *reference, psImage *input,
+bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
                                 psImage *weight, int footprint, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_IMAGE_NON_NULL(reference, false);
-    PS_ASSERT_IMAGE_TYPE(reference, PS_TYPE_F32, false);
-    if (input) {
-        PS_ASSERT_IMAGE_NON_NULL(input, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(input, reference, false);
-        PS_ASSERT_IMAGE_TYPE(input, PS_TYPE_F32, false);
-    }
-    if (weight) {
-        PS_ASSERT_IMAGE_NON_NULL(weight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, reference, false);
-        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
-    }
+    PS_ASSERT_IMAGE_NON_NULL(image1, false);
+    PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false);
+    if (image2) {
+        PS_ASSERT_IMAGE_NON_NULL(image2, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false);
+        PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false);
+    }
+    PS_ASSERT_IMAGE_NON_NULL(weight, false);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false);
+    PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
-    int numCols = reference->numCols, numRows = reference->numRows; // Size of images
+    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + footprint; // Size of postage stamps
-
-    if (!weight) {
-        // Use the input (or if I must, the reference) as a rough approximation to the variance map, and HOPE
-        // that it's positive.
-        weight = input ? input : reference;
-    }
 
     for (int i = 0; i < stamps->num; i++) {
         pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+        if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) {
             continue;
         }
@@ -468 +439 @@
         }
 
+        // Catch memory leaks --- these should have been freed and NULLed before
+        assert(stamp->image1 == NULL);
+        assert(stamp->image2 == NULL);
+        assert(stamp->weight == NULL);
+
         psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest
 
-        psImage *refSub = psImageSubset(reference, region); // Subimage with stamp
-        stamp->reference = psKernelAllocFromImage(refSub, size, size);
-        psFree(refSub);                 // Drop reference
-
-        if (input) {
-            psImage *inSub = psImageSubset(input, region); // Subimage with stamp
-            stamp->input = psKernelAllocFromImage(inSub, size, size);
-            psFree(inSub);              // Drop reference
+        psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp
+        stamp->image1 = psKernelAllocFromImage(sub1, size, size);
+        psFree(sub1);                   // Drop reference
+
+        if (image2) {
+            psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp
+            stamp->image2 = psKernelAllocFromImage(sub2, size, size);
+            psFree(sub2);               // Drop reference
         }
 
         psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp
         stamp->weight = psKernelAllocFromImage(wtSub, size, size);
-        psFree(wtSub);                 // Drop reference
+        psFree(wtSub);                  // Drop reference
+
+        stamp->status = PM_SUBTRACTION_STAMP_CALCULATE;
     }
 
@@ -488 +466 @@
 }
 
-
+#if 0
 bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
 {
@@ -517 +495 @@
         float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame
 
-        psFree(stamp->input);
-        stamp->input = psKernelAlloc(-size, size, -size, size);
-        psKernel *input = stamp->input; // Target stamp
+        psFree(stamp->image2);
+        stamp->image2 = psKernelAlloc(-size, size, -size, size);
+        psKernel *target = stamp->image2; // Target stamp
 
         // Put in a Waussian, just for fun!
@@ -525 +503 @@
             for (int u = -size; u <= size; u++) {
                 float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma));
-                input->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
+                target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z));
             }
         }
@@ -533 +511 @@
     return true;
 }
-
+#endif
 
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
                                                           const psRegion *region, float spacing,
-                                                          int exclusionZone)
+                                                          pmSubtractionMode mode)
 {
     PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
@@ -561 +539 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
-                                                            spacing, exclusionZone); // Stamps to return
+                                                            spacing, mode); // Stamps to return
     psFree(x);
     psFree(y);
@@ -572 +550 @@
     return stamps;
 }
+
+
+pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask,
+                                                       const psRegion *region, float spacing,
+                                                       pmSubtractionMode mode)
+{
+    PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
+    // Let pmSubtractionStampsSet take care of the rest of the assertions
+
+    const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file
+    psArray *data = psVectorsReadFromFile(filename, format);
+    if (!data) {
+        psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename);
+        return NULL;
+    }
+    psVector *x = data->data[0], *y = data->data[1]; // Stamp positions
+    psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes
+
+    // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions
+    psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+    psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
+
+    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing,
+                                                            mode);
+    psFree(data);
+
+    return stamps;
+
+}

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -9 +9 @@
 typedef enum {
     PM_SUBTRACTION_STAMP_INIT,          ///< Initial state
+    PM_SUBTRACTION_STAMP_FOUND,         ///< Found a suitable source for this stamp
+    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
     PM_SUBTRACTION_STAMP_USED,          ///< Use this stamp
     PM_SUBTRACTION_STAMP_REJECTED,      ///< This stamp has been rejected
-    PM_SUBTRACTION_STAMP_CALCULATE,     ///< Calculate matrix and vector values for this stamp
     PM_SUBTRACTION_STAMP_NONE           ///< No stamp in this region
 } pmSubtractionStampStatus;
@@ -54 +55 @@
     float flux;                         ///< Flux
     float xNorm, yNorm;                 ///< Normalised position
-    psKernel *reference;                ///< Reference image postage stamp
-    psKernel *input;                    ///< Input image postage stamp
+    psKernel *image1;                   ///< Reference image postage stamp
+    psKernel *image2;                   ///< Input image postage stamp
     psKernel *weight;                   ///< Weight image postage stamp
-    psArray *convolutions;              ///< Convolutions of the reference for each kernel component
+    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component
+    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component
     psImage *matrix;                    ///< Associated matrix
     psVector *vector;                   ///< Assoicated vector
@@ -72 +74 @@
                                                 const psRegion *region, ///< Region to search, or NULL
                                                 float threshold, ///< Threshold for stamps in the image
-                                                float spacing ///< Rough spacing for stamps
+                                                float spacing, ///< Rough spacing for stamps
+                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 
 /// Set stamps based on a list of x,y
-///
-/// We may optionally apply an exclusion zone around each star --- this is important when we're convolving to
-/// a specified PSF; less so when we're only trying to get a good subtraction.
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp
                                                const psVector *y, ///< y coordinates for each stamp
@@ -86 +86 @@
                                                const psRegion *region, ///< Region to search, or NULL
                                                float spacing, ///< Rough spacing for stamps
-                                               int exclusionZone ///< Exclusion zone around each star
+                                               pmSubtractionMode mode ///< Mode for subtraction
     );
 
-///
+/// Set stamps based on a list of sources
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(
-    const psArray *sources, ///< Sources for each stamp
-    const psImage *subMask, ///< Mask, or NULL
-    const psRegion *region, ///< Region to search, or NULL
-    float spacing, ///< Rough spacing for stamps
-    int exclusionZone ///< Exclusion zone around each star
+    const psArray *sources,             ///< Sources for each stamp
+    const psImage *subMask,             ///< Mask, or NULL
+    const psRegion *region,             ///< Region to search, or NULL
+    float spacing,                      ///< Rough spacing for stamps
+    pmSubtractionMode mode              ///< Mode for subtraction
+    );
+
+/// Set stamps based on values in a file
+pmSubtractionStampList *pmSubtractionStampsSetFromFile(
+    const char *filename,               ///< Filename of file containing x,y (or x,y,flux) on each line
+    const psImage *subMask,             ///< Mask, or NULL
+    const psRegion *region,             ///< Region to search, or NULL
+    float spacing,                      ///< Rough spacing for stamps
+    pmSubtractionMode mode              ///< Mode for subtraction
     );
 
@@ -107 +116 @@
     );
 
-/// Generate target for stamps based on list
-bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, ///< Stamps
-                                 float fwhm, ///< FWHM for each stamp
-                                 int footprint, ///< Stamp footprint size
-                                 int size ///< Kernel half-size
-    );
-
 #endif