Changeset 15443 for trunk/psModules/src/imcombine/pmSubtraction.c

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

File:

• trunk/psModules/src/imcombine/pmSubtraction.c (modified) (25 diffs)

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;
+                            }
                         }
                     }