Changeset 15756

Timestamp:

Dec 6, 2007, 3:57:15 PM (18 years ago)

Author:

Paul Price

Message:

Implementing dual-convolution. This required reworking those
functions involved with calculating and solving the equations; moved
these into a separate file and made various other cleanups (e.g.,
assertions for pmSubtractionKernels). Removed the normalisation
(central pixel) term from all kernels, because this shouldn't be in
the second kernel (there's only one normalisation term between the two
kernels); the normalisation is treated explicitly in the equations,
along with the background (still only a constant background is
supported for now, but there's a lot of support for a polynomial
background for when I get around to putting it in the equations).
Tested the single convolution, and it's working; same results as
before, apparently. Haven't tested dual convolution yet.

Location:

trunk/psModules/src

Files:

• camera/pmReadoutFake.c (modified) (4 diffs)
• camera/pmReadoutFake.h (modified) (1 diff)
• imcombine/Makefile.am (modified) (2 diffs)
• imcombine/pmStackReject.c (modified) (2 diffs)
• imcombine/pmSubtraction.c (modified) (35 diffs)
• imcombine/pmSubtraction.h (modified) (6 diffs)
• imcombine/pmSubtractionEquation.c (added)
• imcombine/pmSubtractionEquation.h (added)
• imcombine/pmSubtractionKernels.c (modified) (31 diffs)
• imcombine/pmSubtractionKernels.h (modified) (10 diffs)
• imcombine/pmSubtractionMask.c (added)
• imcombine/pmSubtractionMask.h (added)
• imcombine/pmSubtractionMatch.c (modified) (19 diffs)
• imcombine/pmSubtractionMatch.h (modified) (1 diff)
• imcombine/pmSubtractionParams.c (modified) (3 diffs)
• imcombine/pmSubtractionStamps.c (modified) (20 diffs)
• imcombine/pmSubtractionStamps.h (modified) (9 diffs)
• psmodules.h (modified) (1 diff)

trunk/psModules/src/camera/pmReadoutFake.c

@@ -20 +20 @@
 #include "pmSource.h"
 #include "pmSourceUtils.h"
+#include "pmModelUtils.h"
 
 #define MODEL_TYPE "PS_MODEL_RGAUSS"     // Type of model to use
 
 
-pmReadout *pmReadoutFakeFromSources(int numCols, int numRows, const psArray *sources,
-                                    float fwhm, float minFlux)
+bool pmReadoutFakeFromSources(pmReadout *readout, int numCols, int numRows, const psArray *sources,
+                              const psVector *xOffset, const psVector *yOffset, pmPSF *psf, float minFlux,
+                              int radius)
 {
-    PS_ASSERT_INT_LARGER_THAN(numCols, 0, NULL);
-    PS_ASSERT_INT_LARGER_THAN(numRows, 0, NULL);
-    PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, NULL);
-    PS_ASSERT_FLOAT_LARGER_THAN(minFlux, 0.0, NULL);
+    PS_ASSERT_PTR_NON_NULL(readout, false);
+    PS_ASSERT_INT_LARGER_THAN(numCols, 0, false);
+    PS_ASSERT_INT_LARGER_THAN(numRows, 0, false);
+    PS_ASSERT_ARRAY_NON_NULL(sources, false);
+    if (xOffset || yOffset) {
+        PS_ASSERT_VECTOR_NON_NULL(xOffset, false);
+        PS_ASSERT_VECTOR_NON_NULL(yOffset, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(xOffset, yOffset, false);
+        PS_ASSERT_VECTOR_TYPE(xOffset, PS_TYPE_S32, false);
+        PS_ASSERT_VECTOR_TYPE_EQUAL(xOffset, yOffset, false);
+        if (xOffset->n != sources->n) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                    "Number of offset vectors (%ld) and sources (%ld) doesn't match",
+                    xOffset->n, sources->n);
+            return false;
+        }
+    }
+    PS_ASSERT_PTR_NON_NULL(psf, false);
+    if (radius > 0 && isfinite(minFlux) && minFlux > 0.0) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot define both minimum flux and fixed radius.");
+        return false;
+    }
 
-    pmReadout *readout = pmReadoutAlloc(NULL); // Output readout
-    readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+
+    readout->image = psImageRecycle(readout->image, numCols, numRows, PS_TYPE_F32);
     psImageInit(readout->image, 0);
 
     int numSources = sources->n;          // Number of stars
 
+#if 0
     pmModelType modelType = pmModelClassGetType(MODEL_TYPE); // Type of PSF model
     assert(modelType >= 0);
@@ -65 +85 @@
         psAbort("Unsupported model type: %d", modelType);
     }
+#endif
 
+    pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0,
+                                             1.0); // Fake model, with central intensity of 1.0
     float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0
-
 
     for (int i = 0; i < numSources; i++) {
@@ -93 +115 @@
         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
         fakeSource->peak = pmPeakAlloc(x, y, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
-        float radius = fakeModel->modelRadius(fakeModel->params, minFlux); // Radius of interest for source
+        float fakeRadius = radius > 0 ? radius : fakeModel->modelRadius(fakeModel->params, minFlux); // Radius
 
-        if (!pmSourceDefinePixels(fakeSource, readout, x, y, radius)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
-            psFree(readout);
-            psFree(fakeModel);
-            return NULL;
-        }
-
-        if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to add model of source to image.");
-            psFree(readout);
-            psFree(fakeModel);
-            return NULL;
+        if (xOffset) {
+            if (!pmSourceDefinePixels(fakeSource, readout, x + xOffset->data.S32[i],
+                                      y + yOffset->data.S32[i], fakeRadius)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
+                psFree(readout);
+                psFree(fakeModel);
+                return false;
+            }
+            if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0,
+                                      - xOffset->data.S32[i], - yOffset->data.S32[i])) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add model of source to image.");
+                psFree(readout);
+                psFree(fakeModel);
+                return false;
+            }
+        } else {
+            if (!pmSourceDefinePixels(fakeSource, readout, x, y, fakeRadius)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source.");
+                psFree(readout);
+                psFree(fakeModel);
+                return false;
+            }
+            if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add model of source to image.");
+                psFree(readout);
+                psFree(fakeModel);
+                return false;
+            }
         }
         psFree(fakeSource);
@@ -113 +151 @@
     psFree(fakeModel);
 
-    return readout;
+    return true;
 }

trunk/psModules/src/camera/pmReadoutFake.h

@@ -6 +6 @@
 #include <pmFPA.h>
 
-// Generate a fake readout from an array of sources
-pmReadout *pmReadoutFakeFromSources(int numCols, int numRows, ///< Dimension of image
-                                    const psArray *sources, ///< Array of pmSource
-                                    float fwhm, ///< FWHM for sources
-                                    float minFlux ///< Minimum flux to bother about; for setting object size
+#include <pmMoments.h>
+#include <pmResiduals.h>
+#include <pmGrowthCurve.h>
+#include <pmTrend2D.h>
+#include <pmPSF.h>
+
+
+/// Generate a fake readout from an array of sources
+bool pmReadoutFakeFromSources(pmReadout *readout, ///< Output readout, or NULL
+                              int numCols, int numRows, ///< Dimension of image
+                              const psArray *sources, ///< Array of pmSource
+                              const psVector *xOffset, ///< x offsets for sources (source -> img), or NULL
+                              const psVector *yOffset, ///< y offsets for sources (source -> img), or NULL
+                              pmPSF *psf, ///< PSF for sources
+                              float minFlux, ///< Minimum flux to bother about; for setting source radius
+                              int radius ///< Fixed radius for sources
     );
 
-
-
-
 #endif

trunk/psModules/src/imcombine/Makefile.am

@@ -9 +9 @@
         pmStackReject.c         \
         pmSubtraction.c         \
+        pmSubtractionEquation.c \
         pmSubtractionKernels.c  \
+        pmSubtractionMask.c     \
+        pmSubtractionMatch.c    \
+        pmSubtractionParams.c   \
         pmSubtractionStamps.c   \
-        pmSubtractionMatch.c    \
-        pmSubtractionParams.c
+        pmPSFEnvelope.c
 
 pkginclude_HEADERS = \
@@ -20 +23 @@
         pmStackReject.h         \
         pmSubtraction.h         \
+        pmSubtractionEquation.h \
         pmSubtractionKernels.h  \
+        pmSubtractionMask.h     \
+        pmSubtractionMatch.h    \
+        pmSubtractionParams.h   \
         pmSubtractionStamps.h   \
-        pmSubtractionMatch.h    \
-        pmSubtractionParams.h
+        pmPSFEnvelope.h
 
 CLEANFILES = *~

trunk/psModules/src/imcombine/pmStackReject.c

@@ -57 +57 @@
     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(convRO, inRO, NULL, NULL, 0, region, solution, kernels,
-                                   PM_SUBTRACTION_MODE_1, true)) {
+        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -72 +70 @@
         float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols;
         float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows;
-        psImage *kernel = pmSubtractionKernelImage(solution, kernels, xNorm, yNorm);
+        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
         if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.73 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-11-05 23:43:38 $
+ *  @version $Revision: 1.74 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-07 01:57:15 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -23 +23 @@
 #include "pmFPA.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 
 #include "pmSubtraction.h"
@@ -34 +35 @@
 // Private (file-static) functions
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-// Return the number of polynomial terms there are for a given order
-static inline long polyTerms(int order  // Order of polynomial
-                             )
-{
-    return (order + 1) * (order + 2) / 2;
-}
-
-// Given a stamp coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
-static psImage *spatialPolyValues(int spatialOrder, // Maximum spatial polynomial order
-                                  float x, float y // Normalised position of interest, [-1,1]
-    )
-{
-    assert(spatialOrder >= 0);
-    assert(x >= -1 && x <= 1);
-    assert(y >= -1 && y <= 1);
-
-    psImage *polyValues = psImageAlloc(spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); // Matrix to return
-    polyValues->data.F64[0][0] = 1.0;
-
-    double value = 1.0;
-    for (int i = 1; i <= spatialOrder; i++) {
-        value *= x;
-        polyValues->data.F64[0][i] = value;
-    }
-
-    value = 1.0;
-    for (int j = 1; j <= spatialOrder; j++) {
-        value *= y;
-        polyValues->data.F64[j][0] = value;
-    }
-
-    for (int j = 1; j <= spatialOrder; j++) {
-        for (int i = 1; i <= spatialOrder - j; i++) {
-            polyValues->data.F64[j][i] = polyValues->data.F64[j][0] * polyValues->data.F64[0][i];
-        }
-    }
-
-    return polyValues;
-}
 
 // Generate the kernel to apply to the variance from the normal kernel
@@ -108 +69 @@
 
 // Generate an image of the solved kernel
-static inline psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
-                                     const psVector *solution, // Solution to the least-squares problem
-                                     const pmSubtractionKernels *kernels, // Kernel basis functions
-                                     const psImage *polyValues // Spatial polynomial values
-                                     )
-{
+static psKernel *solvedKernel(psKernel *kernel, // Kernel, to return
+                              const pmSubtractionKernels *kernels, // Kernel basis functions
+                              const psImage *polyValues, // Spatial polynomial values
+                              bool wantDual // Want the dual (second) kernel?
+                              )
+{
+    assert(kernels);
+    assert(polyValues);
+
     int numKernels = kernels->num;      // Number of kernel basis functions
-    assert(kernels->u->n == numKernels);
-    assert(kernels->v->n == numKernels);
-    int spatialOrder = kernels->spatialOrder; // Maximum spatial polynomial order
-    assert(solution->n == numKernels * polyTerms(spatialOrder) + polyTerms(kernels->bgOrder));
-
-    // Ensure the subIndex for POIS kernels is what is expected
-    assert((kernels->type != PM_SUBTRACTION_KERNEL_POIS &&
-            kernels->type != PM_SUBTRACTION_KERNEL_SPAM &&
-            kernels->type != PM_SUBTRACTION_KERNEL_FRIES &&
-            kernels->type != PM_SUBTRACTION_KERNEL_GUNK &&
-            kernels->type != PM_SUBTRACTION_KERNEL_RINGS) ||
-           (kernels->u->data.S32[kernels->subIndex] == 0 && kernels->v->data.S32[kernels->subIndex] == 0));
-
     int size = kernels->size;           // Kernel half-size
     if (!kernel) {
@@ -135 +86 @@
 
     for (int i = 0; i < numKernels; i++) {
-        double value = 0.0;             // Value of this component, from the sum of multiple polynomial terms
-        for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) {
-            for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
-                value += polyValues->data.F64[yOrder][xOrder] * solution->data.F64[index];
-            }
-        }
+        double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value
 
         switch (kernels->type) {
@@ -147 +93 @@
               int v = kernels->v->data.S32[i]; // Offset in y
               kernel->kernel[v][u] += value;
-              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
-                  // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                  kernel->kernel[0][0] -= value;
-              }
+              kernel->kernel[0][0] -= value;
               break;
           }
@@ -167 +110 @@
                   for (int u = uStart; u <= uStop; u++) {
                       kernel->kernel[v][u] += norm * value;
-                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
-                          // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                          kernel->kernel[0][0] -= value;
-                      }
+                      kernel->kernel[0][0] -= value;
                   }
               }
@@ -188 +128 @@
               } else {
                   // Using delta function
-                  bool subtract = (kernels->spatialOrder > 0 && i != kernels->subIndex); // Subtract (0,0)?
                   int u = kernels->u->data.S32[i]; // Offset in x
                   int v = kernels->v->data.S32[i]; // Offset in y
                   kernel->kernel[v][u] += value;
-                  if (subtract) {
-                      // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                      kernel->kernel[0][0] -= value;
-                  }
+                  kernel->kernel[0][0] -= value;
               }
               break;
@@ -211 +147 @@
           }
           case PM_SUBTRACTION_KERNEL_RINGS: {
-              if (i == kernels->subIndex) {
-                  kernel->kernel[0][0] += value;
-                  break;
-              }
               psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data
               psVector *uCoords = preCalc->data[0]; // u coordinates
@@ -232 +164 @@
         }
     }
+
+    if (!wantDual) {
+        // Put in the normalisation component
+        kernel->kernel[0][0] += p_pmSubtractionSolutionNorm(kernels);
+    }
+
 
     return kernel;
@@ -328 +266 @@
 }
 
-// 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 a region of an image
+static inline void convolveRegion(psImage *convImage, // Convolved image (output)
+                                  psImage *convWeight, // Convolved weight map (output), or NULL
+                                  psKernel **kernelImage, // Convolution kernel for the image
+                                  psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL
+                                  const psImage *image, // Image to convolve
+                                  const psImage *weight, // Weight map to convolve, or NULL
+                                  const psImage *subMask, // Subtraction mask
+                                  psMaskType maskVal, // Mask value to apply in convolution
+                                  const pmSubtractionKernels *kernels, // Kernels
+                                  psImage *polyValues, // Polynomial values
+                                  float background, // Background value to apply
+                                  psRegion region, // Region to convolve
+                                  bool useFFT  // Use FFT to convolve?
     )
 {
-    image->data.F32[y][x] = NAN;
-    if (mask) {
-        mask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
-    }
+    *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, false);
     if (weight) {
-        weight->data.F32[y][x] = NAN;
-    }
+        *kernelWeight = varianceKernel(*kernelWeight, *kernelImage);
+    }
+
+    if (useFFT) {
+        // Use Fast Fourier Transform to do the convolution
+        // This provides a big speed-up for large kernels
+        convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size);
+        if (weight) {
+            convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size);
+        }
+    } else {
+        convolveDirect(convImage, image, *kernelImage, region, background, kernels->size);
+        if (weight) {
+            convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size);
+        }
+    }
+
     return;
 }
@@ -350 +308 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-// Generate the convolution given a precalculated kernel
 psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel)
 {
@@ -363 +320 @@
 }
 
+psImage *p_pmSubtractionPolynomial(psImage *output, int spatialOrder, float x, float y)
+{
+    assert(spatialOrder >= 0);
+    assert(x >= -1 && x <= 1);
+    assert(y >= -1 && y <= 1);
+
+    output = psImageRecycle(output, spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64);
+    output->data.F64[0][0] = 1.0;
+
+    double value = 1.0;
+    for (int i = 1; i <= spatialOrder; i++) {
+        value *= x;
+        output->data.F64[0][i] = value;
+    }
+
+    value = 1.0;
+    for (int j = 1; j <= spatialOrder; j++) {
+        value *= y;
+        output->data.F64[j][0] = value;
+    }
+
+    for (int j = 1; j <= spatialOrder; j++) {
+        for (int i = 1; i <= spatialOrder - j; i++) {
+            output->data.F64[j][i] = output->data.F64[j][0] * output->data.F64[0][i];
+        }
+    }
+
+    return output;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -368 +354 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-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
@@ -538 +362 @@
     )
 {
+#if 0
+    if (index == kernels->num) {
+        // Normalisation component
+        return convolveOffset(image, 0, 0, footprint);
+    }
+#endif
+
     switch (kernels->type) {
       case PM_SUBTRACTION_KERNEL_POIS: {
@@ -543 +374 @@
           int v = kernels->v->data.S32[index]; // Offset in y
           psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
-          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              convolveSub(convolved, image, footprint);
-          }
+          convolveSub(convolved, image, footprint);
           return convolved;
       }
@@ -569 +398 @@
               }
           }
-          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              convolveSub(convolved, image, footprint);
-          }
+          convolveSub(convolved, image, footprint);
           return convolved;
       }
@@ -583 +410 @@
           int v = kernels->v->data.S32[index]; // Offset in y
           psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
-          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              convolveSub(convolved, image, footprint);
-          }
+          convolveSub(convolved, image, footprint);
           return convolved;
       }
@@ -593 +418 @@
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
-          if (index == kernels->subIndex) {
-              // Simply copying over the image
-              return convolveOffset(image, 0, 0, footprint);
-          }
           psKernel *convolved = psKernelAlloc(-footprint, footprint,
                                               -footprint, footprint); // Convolved image
@@ -650 +471 @@
 
 
-bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint,
-                                pmSubtractionMode mode)
+bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
 {
     PS_ASSERT_PTR_NON_NULL(stamp, false);
-    PS_ASSERT_PTR_NON_NULL(kernels, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PS_ASSERT_INT_NONNEGATIVE(footprint, false);
 
@@ -663 +482 @@
     }
 
-    switch (mode) {
+    switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_TARGET:
       case PM_SUBTRACTION_MODE_1:
@@ -677 +496 @@
         break;
       default:
-        psAbort("Unsupported subtraction mode: %x", mode);
+        psAbort("Unsupported subtraction mode: %x", kernels->mode);
     }
 
@@ -684 +503 @@
 
 
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels,
-                                    int footprint, pmSubtractionMode mode)
-{
-    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
-    PS_ASSERT_PTR_NON_NULL(kernels, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
-    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
-
-    int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation
-    int numKernels = kernels->num;      // Number of kernel basis functions
-    int numSpatial = polyTerms(spatialOrder); // Number of spatial variations
-    int numBackground = polyTerms(kernels->bgOrder); // Number of background terms
-
-    // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the
-    // number of coefficients for the spatial polynomial, and a constant background offset.
-    int numParams = numKernels * numSpatial + numBackground;
-    int bgIndex = numParams - numBackground; // Index in matrix for the background
-
-    // We iterate over each stamp, allocate the matrix and vectors if
-    // necessary, and then calculate those matrix/vectors.
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
-            continue;
-        }
-
-        psImage *stampMatrix = stamp->matrix; // Least-squares matrix for this stamp
-        psVector *stampVector = stamp->vector; // Least-squares vector for this stamp
-
-        if (!stampMatrix) {
-            stamp->matrix = stampMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64);
-        }
-        assert(stampMatrix->type.type == PS_TYPE_F64);
-        assert(stampMatrix->numCols == numParams && stampMatrix->numRows == numParams);
-        psImageInit(stampMatrix, 0.0);
-
-        if (!stampVector) {
-            stamp->vector = stampVector = psVectorAlloc(numParams, PS_TYPE_F64);
-        }
-        assert(stampVector->type.type == PS_TYPE_F64);
-        assert(stampVector->n == numParams);
-        psVectorInit(stampVector, 0.0);
-
-        // Generate convolutions
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
-            return NULL;
-        }
-
-        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
-
-#ifdef TESTING
-        for (int j = 0; j < numKernels; j++) {
-            psKernel *conv = convolutions->data[j]; // Convolution of interest
-            psString filename = NULL; // Filename for output
-            psStringAppend(&filename, "conv_%03d_%03d.fits", i, j);
-            psFits *fits = psFitsOpen(filename, "w"); // FITS file pointer
-            psFree(filename);
-            psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
-            psFitsClose(fits);
-        }
-#endif
-
-        // Generate least-squares vector and matrix
-        bool bad = false;               // Are there bad values?
-        psString badString = NULL;      // Details on bad values
-        for (int j = 0; j < numKernels && !bad; j++) {
-            psKernel *jConv = convolutions->data[j]; // Convolution for j-th element
-
-            // Generate matrix
-            for (int k = j; k < numKernels && !bad; k++) {
-                psKernel *kConv = convolutions->data[k]; // Convolution for k-th element
-                double sumCC = 0.0; // Sum of the convolution products
-                for (int y = - footprint; y <= footprint; y++) {
-                    for (int x = - footprint; x <= footprint; x++) {
-                        sumCC += jConv->kernel[y][x] * kConv->kernel[y][x] / weight->kernel[y][x];
-                    }
-                }
-                if (!isfinite(sumCC)) {
-                    psStringAppend(&badString, "\nBad sumCC at %d, %d", j, k);
-                    bad = true;
-                    continue;
-                }
-                // Generate the pseudo-convolutions from the spatial polynomial terms
-                for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
-                    for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
-                         jxOrder++, jIndex += numKernels) {
-                        for (int kyOrder = 0, kIndex = k; kyOrder <= spatialOrder; kyOrder++) {
-                            for (int kxOrder = 0; kxOrder <= spatialOrder - kyOrder;
-                                 kxOrder++, kIndex += numKernels) {
-                                double convPoly = sumCC * polyValues->data.F64[jyOrder][jxOrder] *
-                                    polyValues->data.F64[kyOrder][kxOrder]; // Temporary value
-                                stampMatrix->data.F64[jIndex][kIndex] = convPoly;
-                                stampMatrix->data.F64[kIndex][jIndex] = convPoly;
-                            }
-                        }
-                    }
-                }
-            }
-
-            // Vector and background term for matrix
-            double sumC = 0.0;      // Sum of the convolution
-            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;
-                    sumTC += target->kernel[y][x] * convProduct;
-                }
-            }
-            if (!isfinite(sumC)) {
-                psStringAppend(&badString, "\nBad sumC at %d", j);
-                bad = true;
-                continue;
-            }
-            if (!isfinite(sumTC)) {
-                psStringAppend(&badString, "\nBad sumIC detected at %d", j);
-                bad = true;
-                continue;
-            }
-
-            for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
-                for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
-                     jxOrder++, jIndex += numKernels) {
-                    stampVector->data.F64[jIndex] = sumTC * polyValues->data.F64[jyOrder][jxOrder];
-
-                    double convPoly = sumC * polyValues->data.F64[jyOrder][jxOrder]; // Value
-                    stampMatrix->data.F64[jIndex][bgIndex] = convPoly;
-                    stampMatrix->data.F64[bgIndex][jIndex] = convPoly;
-                }
-            }
-
-        }
-        psFree(polyValues);
-
-        // Background only terms
-        double sum1 = 0.0;              // Sum of the weighting
-        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;
-                sumT += target->kernel[y][x] * invNoise2;
-            }
-        }
-        if (!isfinite(sum1)) {
-            psStringAppend(&badString, "\nBad sum1 detected");
-            bad = true;
-        } else if (!isfinite(sumT)) {
-            psStringAppend(&badString, "\nBad sumI detected");
-            bad = true;
-        }
-
-        stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
-        stampVector->data.F64[bgIndex] = sumT;
-
-        if (bad) {
-            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-            psWarning("Rejecting stamp %d (%d,%d) because of bad equation:%s",
-                      i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), badString);
-            psFree(badString);
-        } else {
-            stamp->status = PM_SUBTRACTION_STAMP_USED;
-        }
-
-        if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
-            psString matrixName = NULL;
-            psStringAppend(&matrixName, "matrix%d.fits", i);
-            psFits *matrixFile = psFitsOpen(matrixName, "w");
-            psFree(matrixName);
-            psFitsWriteImage(matrixFile, NULL, stampMatrix, 0, NULL);
-            psFitsClose(matrixFile);
-        }
-
-    }
-
-    return true;
-}
-
-
-psVector *pmSubtractionSolveEquation(psVector *solution, const pmSubtractionStampList *stamps)
-{
-    PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL);
-
-    // Check inputs, while summing the stamp matrices and vectors
-    long numParams = -1;                // Number of parameters
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        PS_ASSERT_PTR_NON_NULL(stamp, NULL);
-        if (stamp->status != PM_SUBTRACTION_STAMP_USED) {
-            continue;
-        }
-        if (numParams == -1) {
-            numParams = stamp->vector->n;
-        } else {
-            PS_ASSERT_VECTOR_SIZE(stamp->vector, numParams, NULL);
-        }
-        PS_ASSERT_VECTOR_TYPE(stamp->vector, PS_TYPE_F64, NULL);
-        PS_ASSERT_IMAGE_TYPE(stamp->matrix, PS_TYPE_F64, NULL);
-        PS_ASSERT_IMAGE_SIZE(stamp->matrix, numParams, numParams, NULL);
-    }
-    if (numParams == -1) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found.");
-        return NULL;
-    }
-
-    if (solution) {
-        PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, NULL);
-        PS_ASSERT_VECTOR_SIZE(solution, numParams, NULL);
-    }
-
-    // Accumulate the least-squares matricies and vectors
-    psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix for all stamps
-    psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector for all stamps
-    psVectorInit(sumVector, 0.0);
-    psImageInit(sumMatrix, 0.0);
-    for (int i = 0; i < stamps->num; i++) {
-        pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_USED) {
-            (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix);
-            (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector);
-        }
-    }
-
-    psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
-    psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
-    psFree(sumMatrix);
-    if (!luMatrix) {
-        psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n");
-        psFree(sumVector);
-        psFree(luMatrix);
-        psFree(permutation);
-        return NULL;
-    }
-
-    solution = psMatrixLUSolve(solution, luMatrix, sumVector, permutation);
-    psFree(sumVector);
-    psFree(luMatrix);
-    psFree(permutation);
-    if (!solution) {
-        psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n");
-        return NULL;
-    }
-
-    if (psTraceGetLevel("psModules.imcombine") >= 7) {
-        for (int i = 0; i < solution->n; i++) {
-            psTrace("psModules.imcombine", 7, "Solution %d: %f\n", i, solution->data.F64[i]);
-        }
-    }
-
-    return solution;
-}
-
-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;
-}
 
 
@@ -1148 +589 @@
 }
 
-psImage *pmSubtractionKernelImage(const psVector *solution, const pmSubtractionKernels *kernels,
-                                  float x, float y
-                                 )
-{
-    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
-    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
-    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
-                          polyTerms(kernels->bgOrder), NULL);
+psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
 
     // Precalulate polynomial values
-    psImage *polyValues = spatialPolyValues(kernels->spatialOrder, x, y);
+    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y);
 
     // The appropriate kernel
-    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues);
+    psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
 
     psFree(polyValues);
@@ -1175 +610 @@
 }
 
-psArray *pmSubtractionKernelSolutions(const psVector *solution, const pmSubtractionKernels *kernels,
-                                      float x, float y
-                                      )
-{
-    PS_ASSERT_VECTOR_NON_NULL(solution, NULL);
-    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
-    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
-                          polyTerms(kernels->bgOrder), NULL);
+#if 0
+psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual)
+{
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL);
     PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL);
     PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, NULL);
 
     psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return
@@ -1194 +624 @@
     for (int i = 0; i < solution->n - 1; i++) {
         fakeSolution->data.F64[i] = solution->data.F64[i];
-        images->data[i] = pmSubtractionKernelImage(fakeSolution, kernels, x, y);
+        images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual);
         fakeSolution->data.F64[i] = 0.0;
     }
@@ -1202 +632 @@
     return images;
 }
-
-bool pmSubtractionConvolve(pmReadout *out, const pmReadout *ro1, const pmReadout *ro2,
+#endif
+
+
+
+bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, 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);
+                           const pmSubtractionKernels *kernels, bool useFFT)
+{
+    PS_ASSERT_PTR_NON_NULL(out1, false);
     PS_ASSERT_PTR_NON_NULL(ro1, false);
     PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
@@ -1217 +649 @@
         PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false);
     }
-    if (mode != PM_SUBTRACTION_MODE_1) {
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false);
+    if (kernels->mode != PM_SUBTRACTION_MODE_1) {
+        PS_ASSERT_PTR_NON_NULL(out2, false);
         PS_ASSERT_PTR_NON_NULL(ro2, false);
         PS_ASSERT_IMAGE_NON_NULL(ro2->image, false);
@@ -1232 +667 @@
         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);
+    if (out1->image) {
+        PS_ASSERT_IMAGE_NON_NULL(out1->image, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out1->image, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out1->image, PS_TYPE_F32, false);
+    }
+    if (out1->mask) {
+        PS_ASSERT_IMAGE_NON_NULL(out1->mask, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out1->mask, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out1->mask, PS_TYPE_MASK, false);
         PS_ASSERT_IMAGE_NON_NULL(subMask, 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);
-    PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, false);
-    PS_ASSERT_PTR_NON_NULL(kernels, false);
-    PS_ASSERT_VECTOR_SIZE(solution, kernels->num * polyTerms(kernels->spatialOrder) +
-                          polyTerms(kernels->bgOrder), -1);
-    PS_ASSERT_VECTOR_SIZE(kernels->u, kernels->num, false);
-    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
+    if (out1->weight) {
+        PS_ASSERT_IMAGE_NON_NULL(out1->weight, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(out1->weight, ro1->image, false);
+        PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false);
+    }
     if (region) {
         if (psRegionIsNaN(*region)) {
@@ -1273 +701 @@
 
     psImage *image, *weight;            // Image and weight map to convolve
-    switch (mode) {
+    switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_TARGET:
       case PM_SUBTRACTION_MODE_1:
+      case PM_SUBTRACTION_MODE_DUAL:
         image = ro1->image;
         weight = ro1->weight;
@@ -1281 +710 @@
       case PM_SUBTRACTION_MODE_2:
         image = ro2->image;
-        weight = ro2->image;
+        weight = ro2->weight;
         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
+        psAbort("Unsupported subtraction mode: %x", kernels->mode);
+    }
+
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
 
-    psImage *convImage = out->image;   // Convolved image
-    if (!convImage) {
-        out->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        convImage = out->image;
+    psImage *convImage1 = out1->image;   // Convolved image
+    if (!convImage1) {
+        convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
     }
     psImage *convMask = NULL;           // Convolved mask image
     if (subMask) {
-        if (!out->mask) {
-            out->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
-        }
-        convMask = out->mask;
+        if (!out1->mask) {
+            out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
+        }
+        convMask = out1->mask;
         psImageInit(convMask, 0);
     }
-    psImage *convWeight = NULL;         // Convolved weight (variance) image
+    psImage *convWeight1 = NULL;         // Convolved weight (variance) image
     if (weight) {
-        if (!out->weight) {
-            out->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        convWeight = out->weight;
-        psImageInit(convWeight, 0.0);
+        if (!out1->weight) {
+            out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        convWeight1 = out1->weight;
+        psImageInit(convWeight1, 0.0);
+    }
+
+    psImage *convImage2, *convWeight2;   // Convolved products for dual mode
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        convImage2 = out2->image;
+        if (!convImage2) {
+            convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        convWeight2 = NULL;         // Convolved weight (variance) image
+        if (ro2->weight) {
+            if (!out2->weight) {
+                out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            convWeight2 = out2->weight;
+            psImageInit(convWeight2, 0.0);
+        }
+        if (subMask) {
+            // Copying mask --- they're the same!
+            psFree(out2->mask);
+            out2->mask = psMemIncrRefCounter(convMask);
+        }
     }
 
     int size = kernels->size;           // Half-size of kernel
     int fullSize = 2 * size + 1;        // Full size of kernel
-    psImage *polyValues = NULL;         // Pre-calculated polynomial values
-
-    psKernel *kernelImage = NULL;       // Kernel for the image
-    psKernel *kernelWeight = NULL;      // Kernel for the weight map
 
     // Get region for convolution: [xMin:xMax,yMin:yMax]
@@ -1330 +773 @@
     }
 
-    psMaskType maskSource, maskTarget;  // Mask values for source and target
-    switch (mode) {
+    psMaskType maskSource;              // Mask these pixels when convolving
+    psMaskType maskTarget;              // Mark these pixels as bad when propagating the subtractionMask
+    switch (kernels->mode) {
       case PM_SUBTRACTION_MODE_TARGET:
       case PM_SUBTRACTION_MODE_1:
@@ -1341 +785 @@
         maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
         break;
+      case PM_SUBTRACTION_MODE_DUAL:
+        maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2;
+        maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2;
+        break;
       default:
-        psAbort("Unsupported subtraction mode: %x", mode);
-    }
+        psAbort("Unsupported subtraction mode: %x", kernels->mode);
+    }
+
+    psImage *polyValues = NULL;         // Pre-calculated polynomial values
+    psKernel *kernelImage = NULL;       // Kernel for the images
+    psKernel *kernelWeight = NULL;      // Kernel for the weight maps
 
     for (int j = yMin; j < yMax; j += fullSize) {
         int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
+        float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate
         for (int i = xMin; i < xMax; i += fullSize) {
             int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
+            float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
             // (with the stamps) that it does not vary rapidly on this scale.
-            psFree(polyValues);
-            polyValues = spatialPolyValues(kernels->spatialOrder,
-                                           2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols,
-                                           2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows);
-
-            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues);
-            if (weight) {
-                kernelWeight = varianceKernel(kernelWeight, kernelImage);
-            }
+            polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm);
+            float background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background term
 
             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, image, subMask, maskSource, kernelImage, subRegion,
-                            background, size);
-                if (weight) {
-                    convolveFFT(convWeight, weight, subMask, maskSource, kernelWeight, subRegion,
-                                0.0, size);
-                }
-            } else {
-                convolveDirect(convImage, image, kernelImage, subRegion, background, size);
-                if (weight) {
-                    convolveDirect(convWeight, weight, kernelWeight, subRegion, 0.0, size);
-                }
+
+            convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight,
+                           subMask, maskSource, kernels, polyValues, background, subRegion, useFFT);
+
+            if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+                convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight,
+                               subMask, maskSource, kernels, polyValues, background, subRegion, useFFT);
             }
 
@@ -1385 +825 @@
                         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;
+                            convImage1->data.F32[y][x] = NAN;
                             if (weight) {
-                                convWeight->data.F32[y][x] = NAN;
+                                convWeight1->data.F32[y][x] = NAN;
                             }
                         }
@@ -1402 +842 @@
     return true;
 }
-
-bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask,
-                         int size, psMaskType blank)
-{
-    PS_ASSERT_IMAGE_NON_NULL(image, false);
-    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
-    if (mask) {
-        PS_ASSERT_IMAGE_NON_NULL(mask, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false);
-        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false);
-    }
-    if (weight) {
-        PS_ASSERT_IMAGE_NON_NULL(weight, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false);
-        PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false);
-    }
-
-    int numCols = image->numCols, numRows = image->numRows; // Image dimensions
-
-    for (int y = size; y < numRows - size; y++) {
-        for (int x = 0; x < size; x++) {
-            markBlank(image, mask, weight, x, y, blank);
-        }
-        for (int x = numCols - size; x < numCols; x++) {
-            markBlank(image, mask, weight, x, y, blank);
-        }
-    }
-    for (int y = 0; y < size; y++) {
-        for (int x = 0; x < numCols; x++) {
-            markBlank(image, mask, weight, x, y, blank);
-        }
-    }
-    for (int y = numRows - size; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
-            markBlank(image, mask, weight, x, y, blank);
-        }
-    }
-
-    return true;
-}
-
-

trunk/psModules/src/imcombine/pmSubtraction.h

@@ -6 +6 @@
  * @author GLG, MHPCC
  *
- * @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- * @date $Date: 2007-11-10 01:09:20 $
+ * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ * @date $Date: 2007-12-07 01:57:15 $
  * Copyright 2004-207 Institute for Astronomy, University of Hawaii
  */
@@ -37 +37 @@
 } pmSubtractionMasks;
 
-/// Generate a mask for use in the subtraction process
-psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved)
-                           const psImage *inMask, ///< Mask for the input image, or NULL
-                           psMaskType maskVal, ///< Value to mask out
-                           int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
-                           int footprint, ///< Half-size of the kernel footprint
-                           float badFrac, ///< Maximum fraction of bad input pixels to accept
-                           bool useFFT  ///< Use FFT to do convolution?
-    );
+
+/// Number of terms in a polynomial
+#define PM_SUBTRACTION_POLYTERMS(ORDER) (((ORDER) + 1) * ((ORDER) + 2) / 2)
+
+/// Set the indices for the normalisation and background terms
+#define PM_SUBTRACTION_INDICES(NORM,BG,KERNELS) { \
+    int numSpatial = PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder); /* Number of spatial terms */ \
+    NORM = (KERNELS)->num * numSpatial; \
+    BG = NORM + 1; \
+}
+
+/// Return the index for the start of the normalisation terms
+#define PM_SUBTRACTION_INDEX_NORM(KERNELS) \
+    ((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder))
+
+/// Return the index for the start of the background terms
+#define PM_SUBTRACTION_INDEX_BG(KERNELS) \
+    (((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder)) + 1)
+
 
 /// Convolve the reference stamp with the kernel components
 bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve
                                 const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                int footprint, ///< Half-size of region over which to calculate equation
-                                pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
-    );
-
-/// Calculate the least-squares equation to match the image quality
-bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
-                                    const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                    int footprint, ///< Half-size of region over which to calculate equation
-                                    pmSubtractionMode mode ///< Mode for subtraction (which to convolve)
-                                    );
-
-/// Solve the least-squares equation to match the image quality
-psVector *pmSubtractionSolveEquation(psVector *solution, ///< Solution vector, or NULL
-                                     const pmSubtractionStampList *stamps ///< Stamps
-                                     );
-
-/// 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
+                                int footprint ///< Half-size of region over which to calculate equation
     );
 
@@ -83 +72 @@
 
 /// Generate an image of the convolution kernel
-psImage *pmSubtractionKernelImage(const psVector *solution, ///< Solution vector
-                                  const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                  float x, float y ///< Normalised position [-1,1] for which to generate image
+psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, ///< Kernel parameters
+                                  float x, float y,///< Normalised position [-1,1] for which to generate image
+                                  bool wantDual ///< Calculate for the dual kernel?
                                   );
 
@@ -95 +84 @@
 
 /// Convolve image in preparation for subtraction
-bool pmSubtractionConvolve(pmReadout *out, ///< Output image
+bool pmSubtractionConvolve(pmReadout *out1, ///< Output image 1
+                           pmReadout *out2, ///< Output image 2 (DUAL mode only)
                            const pmReadout *ro1, // Input image 1
                            const pmReadout *ro2, // Input image 2
@@ -101 +91 @@
                            psMaskType blank, ///< Mask value for blank regions
                            const psRegion *region, ///< Region to convolve (or NULL)
-                           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?
-    );
-
-/// Mark the non-convolved part of the image as blank
-bool pmSubtractionBorder(psImage *image,///< Image
-                         psImage *weight, ///< Weight map (or NULL)
-                         psImage *mask, ///< Mask (or NULL)
-                         int size,      ///< Kernel half-size
-                         psMaskType blank ///< Mask value for blank regions
     );
 
@@ -122 +102 @@
     );
 
+/// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
+psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL
+                                   int spatialOrder, ///< Maximum spatial polynomial order
+                                   float x, float y ///< Normalised position of interest, [-1,1]
+    );
+
 /// @}
 #endif

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <string.h>
 #include <strings.h>
 #include <pslib.h>
 
+#include "pmSubtraction.h"
 #include "pmSubtractionKernels.h"
 
@@ -22 +24 @@
     psFree(kernels->vStop);
     psFree(kernels->preCalc);
-}
-
-// Raise a number to an integer power
+    psFree(kernels->solution1);
+    psFree(kernels->solution2);
+}
+
+// Raise an integer to an integer power
 static inline long power(int value,     // Value
                          int exp        // Exponent
@@ -45 +49 @@
 bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size)
 {
-    PS_ASSERT_PTR_NON_NULL(kernels, false);
-    PS_ASSERT_VECTOR_NON_NULL(kernels->u, false);
-    PS_ASSERT_VECTOR_NON_NULL(kernels->v, false);
-    PS_ASSERT_VECTOR_NON_NULL(kernels->widths, false);
-    PS_ASSERT_VECTOR_TYPE(kernels->u, PS_TYPE_S32, false);
-    PS_ASSERT_VECTOR_TYPE(kernels->v, PS_TYPE_S32, false);
-    PS_ASSERT_VECTOR_TYPE(kernels->widths, PS_TYPE_F32, false);
-    PS_ASSERT_ARRAY_NON_NULL(kernels->preCalc, false);
+    PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false);
     PS_ASSERT_INT_NONNEGATIVE(start, false);
     PS_ASSERT_INT_NONNEGATIVE(size, false);
 
-    int numNew = PS_SQR(2 * size + 1);  // Number of new kernel parameters to add
+    int numNew = PS_SQR(2 * size + 1) - 1;  // Number of new kernel parameters to add
 
     // Ensure the sizes match
@@ -68 +65 @@
     for (int v = - size, index = start; v <= size; v++) {
         for (int u = - size; u <= size; u++, index++) {
+            if (v == 0 && u == 0) {
+                // Skip normalisation component: added explicitly
+                index--;
+                continue;
+            }
             kernels->widths->data.F32[index] = NAN;
             kernels->u->data.S32[index] = u;
@@ -81 +83 @@
 
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
-                                                  const psVector *fwhms, const psVector *orders)
+                                                    const psVector *fwhms, const psVector *orders,
+                                                    pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL);
@@ -97 +100 @@
     for (int i = 0; i < numGaussians; i++) {
         int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian
-        psStringAppend(&params, "(%.2f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
+        psStringAppend(&params, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]);
         num += (gaussOrder + 1) * (gaussOrder + 2) / 2;
     }
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS,
-                                                              size, spatialOrder); // The kernels
+                                                              size, spatialOrder, mode); // The kernels
     psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder);
 
@@ -149 +152 @@
     }
 
-    kernels->subIndex = -1;
-
     return kernels;
 }
@@ -159 +160 @@
 
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
-                                                int size, int spatialOrder)
+                                                int size, int spatialOrder, pmSubtractionMode mode)
 {
     pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return
@@ -173 +174 @@
     kernels->uStop = NULL;
     kernels->vStop = NULL;
-    kernels->subIndex = 0;
     kernels->size = size;
     kernels->inner = 0;
     kernels->spatialOrder = spatialOrder;
     kernels->bgOrder = 0;
-
-    return kernels;
-}
-
-pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder)
+    kernels->mode = mode;
+    kernels->solution1 = NULL;
+    kernels->solution2 = NULL;
+
+    return kernels;
+}
+
+pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL);
 
-    int num = PS_SQR(2 * size + 1); // Number of basis functions
+    int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS,
-                                                              size, spatialOrder); // The kernels
+                                                              size, spatialOrder, mode); // The kernels
     psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
@@ -199 +202 @@
     }
 
-    kernels->subIndex = num / 2;
-    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
-           kernels->v->data.S32[kernels->subIndex] == 0);
-
     return kernels;
 }
@@ -208 +207 @@
 
 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
-                                               const psVector *fwhms, const psVector *orders)
+                                               const psVector *fwhms, const psVector *orders,
+                                               pmSubtractionMode mode)
 {
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders); // Kernels
+                                                                  fwhms, orders, mode); // Kernels
     if (!kernels) {
         return NULL;
     }
 
-    // Subtract a reference kernel from all the others to maintain flux scaling
+    // Subtract normalisation from all the others to maintain flux scaling
     if (spatialOrder > 0) {
-        int subIndex = 0;                   // Index of kernel to subtract from others
-        psKernel *subtract = kernels->preCalc->data[subIndex]; // Kernel to subtract from others
         int numGaussians = fwhms->n;       // Number of Gaussians
         for (int i = 0, index = 0; i < numGaussians; i++) {
             for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
                 for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                    if (uOrder % 2 == 0 && vOrder % 2 == 0 && index != subIndex) {
+                    if (uOrder % 2 == 0 && vOrder % 2 == 0) {
                         psKernel *kernel = kernels->preCalc->data[index]; // Kernel to subtract
-                        psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+                        kernel->kernel[0][0] -= 1.0;
                     }
                 }
@@ -242 +240 @@
             psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
             psFitsClose(kernelFile);
-        }
-    }
-
-    return kernels;
-}
-
-/// Generate SPAM kernels
-pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel
-                                               int spatialOrder, ///< Order of spatial variations
-                                               int inner, ///< Inner radius to preserve unbinned
-                                               int binning ///< Kernel binning factor
-    )
+            double sum = 0.0;
+            for (int y = 0; y < kernel->image->numRows; y++) {
+                for (int x = 0; x < kernel->image->numCols; x++) {
+                    sum += kernel->image->data.F32[y][x];
+                }
+            }
+            psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %lf\n", i, sum);
+        }
+    }
+
+    return kernels;
+}
+
+pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
+                                               pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -269 +270 @@
     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
 
-    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
+    int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions
 
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM,
-                                                              size, spatialOrder); // The kernels
+                                                              size, spatialOrder, mode); // The kernels
     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder);
 
@@ -280 +281 @@
              size, inner, binning, spatialOrder, num);
 
-    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
-    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
+    kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
+    kernels->vStop = psVectorAlloc(num, PS_TYPE_S32);
 
     psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element
@@ -314 +315 @@
 
         for (int j = - numTotal; j <= numTotal; j++, index++) {
+            if (i == 0 && j == 0) {
+                // Skip normalisation component: added explicitly
+                index--;
+                continue;
+            }
             int v = locations->data.S32[numTotal + j]; // Location of pixel
             int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel
@@ -327 +333 @@
     }
 
-    kernels->subIndex = num / 2;
-    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
-           kernels->v->data.S32[kernels->subIndex] == 0 &&
-           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->vStop->data.S32[kernels->subIndex] == 0);
-
     psFree(locations);
     psFree(widths);
@@ -340 +340 @@
 
 
-/// Generate FRIES kernels
-pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
-                                                int spatialOrder, ///< Order of spatial variations
-                                                int inner ///< Inner radius to preserve unbinned
-    )
+pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -366 +362 @@
     psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter);
 
-    int num = PS_SQR(2 * numTotal + 1); // Number of basis functions
+    int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions
 
     psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num);
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES,
-                                                              size, spatialOrder); // The kernels
+                                                              size, spatialOrder, mode); // The kernels
     psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder);
 
@@ -377 +373 @@
              size, inner, spatialOrder, num);
 
-    kernels->uStop = psVectorAlloc(num, PS_TYPE_F32);
-    kernels->vStop = psVectorAlloc(num, PS_TYPE_F32);
+    kernels->uStop = psVectorAlloc(num, PS_TYPE_S32);
+    kernels->vStop = psVectorAlloc(num, PS_TYPE_S32);
 
     psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32);
@@ -409 +405 @@
         int uStop = stop->data.S32[numTotal + i]; // Width of pixel
         for (int j = - numTotal; j <= numTotal; j++, index++) {
+            if (i == 0 && j == 0) {
+                // Skip normalisation component: added explicitly
+                index--;
+                continue;
+            }
             int v = start->data.S32[numTotal + j]; // Location of pixel
             int vStop = stop->data.S32[numTotal + j]; // Width of pixel
@@ -422 +423 @@
     }
 
-    kernels->subIndex = num / 2;
-    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
-           kernels->v->data.S32[kernels->subIndex] == 0 &&
-           kernels->uStop->data.S32[kernels->subIndex] == 0 &&
-           kernels->vStop->data.S32[kernels->subIndex] == 0);
-
     psFree(start);
     psFree(stop);
@@ -436 +431 @@
 // Grid United with Normal Kernel
 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
-                                               const psVector *orders, int inner)
+                                               const psVector *orders, int inner, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -449 +444 @@
 
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders); // Kernels
+                                                                  fwhms, orders, mode); // Kernels
     psStringPrepend(&kernels->description, "GUNK=");
     psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder);
@@ -469 +464 @@
 
     int numISIS = kernels->num;         // Number of ISIS kernels
-    int numInner = PS_SQR(2 * inner + 1); // Number of inner kernel elements
 
     if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) {
@@ -475 +469 @@
     }
 
-    kernels->subIndex = numInner/2 + numISIS;
-    assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
-           kernels->v->data.S32[kernels->subIndex] == 0);
-
     return kernels;
 }
 
 // RINGS --- just what it says
-pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder)
+pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
+                                                pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -505 +496 @@
     }
 
-    int numInner = inner;               // Number of pixels in the inner region
+    int numInner = inner - 1;           // Number of pixels in the inner region
     int numOuter = fibNum;              // Number of summed pixels in the outer region
 
     int numRings = numOuter + numInner; // Number of rings (not including the central pixel)
-    int numPoly = (ringsOrder + 1) * (ringsOrder + 2) / 2; // Number of polynomial variants of each ring
-
-    int num = numRings * numPoly + 1; // Total number of basis functions
+    int numPoly = PM_SUBTRACTION_POLYTERMS(ringsOrder); // Number of polynomial variants of each ring
+
+    int num = numRings * numPoly; // Total number of basis functions
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS,
-                                                              size, spatialOrder); // The kernels
+                                                              size, spatialOrder, mode); // The kernels
     psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder);
 
@@ -522 +513 @@
     // Set the Gaussian kernel parameters
     int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters
-    int radiusLast = 0;                 // Last radius
-    for (int i = 0, index = 0; i < numRings + 1; i++) {
-
+    int radiusLast = 1;                 // Last radius
+    for (int i = 1, index = 0; i < numRings + 1; i++) {
         float lower2;                   // Lower limit of radius^2
         float upper2;                   // Upper limit of radius^2
-        if (i == 0) {
-            lower2 = 0;
-            upper2 = PS_SQR(0.5);
-        } else if (i <= numInner) {
+        if (i <= inner) {
             // A ring every pixel width
             float radius = i;
@@ -572 +559 @@
                     for (int v = -size; v <= size; v++) {
                         int v2 = PS_SQR(v);   // Square of v
-//                        float vPoly = power(v, vOrder); // Value of v^vOrder
                         float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder
 
@@ -579 +565 @@
                             int distance2 = u2 + v2; // Distance from the centre
                             if (distance2 > lower2 && distance2 < upper2) {
-//                                float uPoly = power(u, uOrder); // Value of u^uOrder
                                 float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder
 
@@ -627 +612 @@
     }
 
-    kernels->subIndex = 0;
-
     return kernels;
 }
@@ -634 +617 @@
 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
                                                    const psVector *fwhms, const psVector *orders, int inner,
-                                                   int binning, int ringsOrder)
+                                                   int binning, int ringsOrder, pmSubtractionMode mode)
 {
     switch (type) {
       case PM_SUBTRACTION_KERNEL_POIS:
-        return pmSubtractionKernelsPOIS(size, spatialOrder);
+        return pmSubtractionKernelsPOIS(size, spatialOrder, mode);
       case PM_SUBTRACTION_KERNEL_ISIS:
-        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders);
+        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode);
       case PM_SUBTRACTION_KERNEL_SPAM:
-        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning);
+        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode);
       case PM_SUBTRACTION_KERNEL_FRIES:
-        return pmSubtractionKernelsFRIES(size, spatialOrder, inner);
+        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode);
       case PM_SUBTRACTION_KERNEL_GUNK:
-        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner);
+        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode);
       case PM_SUBTRACTION_KERNEL_RINGS:
-        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder);
+        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);

trunk/psModules/src/imcombine/pmSubtractionKernels.h

@@ -31 +31 @@
     long num;                           ///< Number of kernel components (not including the spatial ones)
     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS)
-    psVector *widths;                   ///< Gaussian FWHMs (ISIS) or ring radius (RINGS)
+    psVector *widths;                   ///< Gaussian FWHMs (ISIS)
     psVector *uStop, *vStop;            ///< Width of kernel element (SPAM,FRIES only)
-    int subIndex;                       ///< Index of kernel to be subtracted (to maintain flux conservation)
     psArray *preCalc;                   ///< Array of images containing pre-calculated kernel (for ISIS)
     int size;                           ///< The half-size of the kernel
@@ -39 +38 @@
     int spatialOrder;                   ///< The spatial order of the kernels
     int bgOrder;                        ///< The order for the background fitting
+    pmSubtractionMode mode;             ///< Mode for subtraction
+    psVector *solution1, *solution2;    ///< Solution for the PSF matching
 } pmSubtractionKernels;
+
+// Assertion to check pmSubtractionKernels
+#define PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(KERNELS, RETURNVALUE) { \
+    PS_ASSERT_PTR_NON_NULL(KERNELS, RETURNVALUE); \
+    PS_ASSERT_STRING_NON_EMPTY((KERNELS)->description, RETURNVALUE); \
+    PS_ASSERT_INT_POSITIVE((KERNELS)->num, RETURNVALUE); \
+    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->u, RETURNVALUE); \
+    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->v, RETURNVALUE); \
+    PS_ASSERT_VECTOR_TYPE((KERNELS)->u, PS_TYPE_S32, RETURNVALUE); \
+    PS_ASSERT_VECTOR_TYPE((KERNELS)->v, PS_TYPE_S32, RETURNVALUE); \
+    PS_ASSERT_VECTOR_SIZE((KERNELS)->u, (KERNELS)->num, RETURNVALUE); \
+    PS_ASSERT_VECTOR_SIZE((KERNELS)->v, (KERNELS)->num, RETURNVALUE); \
+    if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \
+    } \
+    if ((KERNELS)->uStop || (KERNELS)->vStop) { \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \
+        PS_ASSERT_VECTOR_NON_NULL((KERNELS)->vStop, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->uStop, PS_TYPE_S32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->vStop, PS_TYPE_S32, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->uStop, (KERNELS)->num, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->vStop, (KERNELS)->num, RETURNVALUE); \
+    } \
+    if ((KERNELS)->preCalc) { \
+        PS_ASSERT_ARRAY_NON_NULL((KERNELS)->preCalc, RETURNVALUE); \
+        PS_ASSERT_ARRAY_SIZE((KERNELS)->preCalc, (KERNELS)->num, RETURNVALUE); \
+    } \
+    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->size, RETURNVALUE); \
+    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->inner, RETURNVALUE); \
+    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->spatialOrder, RETURNVALUE); \
+    PS_ASSERT_INT_NONNEGATIVE((KERNELS)->bgOrder, RETURNVALUE); \
+}
+
+// Assertion to check that the solution is attached
+#define PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(KERNELS, RETURNVALUE) { \
+    PS_ASSERT_VECTOR_NON_NULL((KERNELS)->solution1, RETURNVALUE); \
+    PS_ASSERT_VECTOR_TYPE((KERNELS)->solution1, PS_TYPE_F64, RETURNVALUE); \
+    PS_ASSERT_VECTOR_SIZE((KERNELS)->solution1, \
+                          (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder) + 1 + \
+                              PM_SUBTRACTION_POLYTERMS((KERNELS)->bgOrder), \
+                          RETURNVALUE); \
+    if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { \
+        PS_ASSERT_VECTOR_NON_NULL(kernels->solution2, RETURNVALUE); \
+        PS_ASSERT_VECTOR_TYPE((KERNELS)->solution2, PS_TYPE_F64, RETURNVALUE); \
+        PS_ASSERT_VECTOR_SIZE((KERNELS)->solution2, \
+                              (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder), \
+                               RETURNVALUE); \
+    } \
+}
 
 /// Generate a delta-function grid for subtraction kernels (like the POIS kernel)
@@ -54 +106 @@
                                                 pmSubtractionKernelsType type, ///< Kernel type
                                                 int size, ///< Half-size of kernel
-                                                int spatialOrder ///< Order of spatial variations
+                                                int spatialOrder, ///< Order of spatial variations
+                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 
 /// Generate POIS kernels
 pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims)
-                                               int spatialOrder ///< Order of spatial variations
+                                               int spatialOrder, ///< Order of spatial variations
+                                               pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -66 +120 @@
                                                     int spatialOrder, ///< Order of spatial variations
                                                     const psVector *fwhms, ///< Gaussian FWHMs
-                                                    const psVector *orders ///< Polynomial order of gaussians
+                                                    const psVector *orders, ///< Polynomial order of gaussians
+                                                    pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -73 +128 @@
                                                int spatialOrder, ///< Order of spatial variations
                                                const psVector *fwhms, ///< Gaussian FWHMs
-                                               const psVector *orders ///< Polynomial order of gaussians
+                                               const psVector *orders, ///< Polynomial order of gaussians
+                                               pmSubtractionMode mode ///< Mode for subtraction
                                                );
 
@@ -80 +136 @@
                                                int spatialOrder, ///< Order of spatial variations
                                                int inner, ///< Inner radius to preserve unbinned
-                                               int binning ///< Kernel binning factor
+                                               int binning, ///< Kernel binning factor
+                                               pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -86 +143 @@
 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel
                                                 int spatialOrder, ///< Order of spatial variations
-                                                int inner ///< Inner radius to preserve unbinned
+                                                int inner, ///< Inner radius to preserve unbinned
+                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -94 +152 @@
                                                const psVector *fwhms, ///< Gaussian FWHMs
                                                const psVector *orders, ///< Polynomial order of gaussians
-                                               int inner ///< Inner radius containing grid of delta functions
+                                               int inner, ///< Inner radius containing grid of delta functions
+                                               pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -101 +160 @@
                                                 int spatialOrder, ///< Order of spatial variations
                                                 int inner, ///< Inner radius to preserve unbinned
-                                                int ringsOrder ///< Polynomial order
+                                                int ringsOrder, ///< Polynomial order
+                                                pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -113 +173 @@
                                                    int inner, ///< Inner radius to preserve unbinned
                                                    int binning, ///< Kernel binning factor
-                                                   int ringsOrder ///< Polynomial order for RINGS
+                                                   int ringsOrder, ///< Polynomial order for RINGS
+                                                   pmSubtractionMode mode ///< Mode for subtraction
     );
 

trunk/psModules/src/imcombine/pmSubtractionMatch.c

@@ -13 +13 @@
 #include "pmSubtractionKernels.h"
 #include "pmSubtractionStamps.h"
+#include "pmSubtractionEquation.h"
 #include "pmSubtraction.h"
+#include "pmSubtractionMask.h"
 #include "pmSubtractionMatch.h"
 
+#define TOL 1.0e-3                      // Tolerance for subtraction solution
+#define KERNEL_MOSAIC 2                 // Half-number of kernel instances in the mosaic image
 
 static bool useFFT = true;              // Do convolutions using FFT
@@ -63 +67 @@
 {
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
-    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, stampSpacing, mode);
+    *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, footprint,
+                                      stampSpacing, mode);
     if (!*stamps) {
         psError(PS_ERR_UNKNOWN, false, "Unable to find stamps.");
@@ -72 +77 @@
 
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, footprint, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -86 +91 @@
 
 
-bool pmSubtractionMatch(pmReadout *convolved, const pmReadout *ro1, const pmReadout *ro2,
+bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2,
                         int footprint, float regionSize, float stampSpacing, float threshold,
                         const psArray *sources, const char *stampsName,
@@ -95 +100 @@
                         psMaskType maskBlank, float badFrac, pmSubtractionMode mode)
 {
-    PS_ASSERT_PTR_NON_NULL(convolved, false);
+    PS_ASSERT_PTR_NON_NULL(conv1, false);
+    if (mode == PM_SUBTRACTION_MODE_DUAL) {
+        PS_ASSERT_PTR_NON_NULL(conv2, false);
+    }
     PS_ASSERT_PTR_NON_NULL(ro1, false);
     PS_ASSERT_IMAGE_NON_NULL(ro1->image, false);
@@ -174 +182 @@
 
     // Reset the output readout, just in case
-    if (convolved->image) {
-        psFree(convolved->image);
-        convolved->image = NULL;
-    }
-    if (convolved->mask) {
-        psFree(convolved->mask);
-        convolved->mask = NULL;
-    }
-    if (convolved->weight) {
-        psFree(convolved->weight);
-        convolved->weight = NULL;
+    if (conv1->image) {
+        psFree(conv1->image);
+        conv1->image = NULL;
+    }
+    if (conv1->mask) {
+        psFree(conv1->mask);
+        conv1->mask = NULL;
+    }
+    if (conv1->weight) {
+        psFree(conv1->weight);
+        conv1->weight = NULL;
     }
 
@@ -250 +258 @@
 
             if (sources) {
-                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, stampSpacing, mode);
+                stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint,
+                                                           stampSpacing, mode);
             } else if (stampsName && strlen(stampsName) > 0) {
-                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, stampSpacing, mode);
+                stamps = pmSubtractionStampsSetFromFile(stampsName, subMask, region, footprint,
+                                                        stampSpacing, mode);
             }
 
@@ -295 +305 @@
                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder);
-            }
-            psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
+                                                       inner, binning, ringsOrder, mode);
+            }
+            psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL",
                              PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels);
 
@@ -312 +322 @@
 
                 psTrace("psModules.imcombine", 3, "Calculating equation...\n");
-                if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) {
+                if (!pmSubtractionCalculateEquation(stamps, kernels)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -320 +330 @@
 
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                solution = pmSubtractionSolveEquation(solution, stamps);
-                if (!solution) {
+
+                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
@@ -328 +338 @@
                 memCheck("  solve equation");
 
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
-                                                                        mode); // Deviations for each stamp
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
@@ -351 +360 @@
             if (numRejected > 0) {
                 psTrace("psModules.imcombine", 3, "Solving equation...\n");
-                solution = pmSubtractionSolveEquation(solution, stamps);
-                if (!solution) {
+                if (!pmSubtractionSolveEquation(kernels, stamps, TOL)) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation.");
                     goto MATCH_ERROR;
                 }
-                psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels,
-                                                                        mode); // Deviations for each stamp
+                psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations
                 if (!deviations) {
                     psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations.");
@@ -376 +383 @@
                 psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32);
                 psImageInit(convKernels, NAN);
-                for (int j = -2; j <= 2; j++) {
-                    for (int i = -2; i <= 2; i++) {
-                        psImage *kernel = pmSubtractionKernelImage(solution, kernels, (float)i / 2.0,
-                                                                   (float)j / 2.0); // Image of the kernel
+                for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) {
+                    for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) {
+                        psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC,
+                                                                   (float)j / (float)KERNEL_MOSAIC,
+                                                                   false); // Image of the kernel
                         if (!kernel) {
                             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
@@ -386 +394 @@
                         }
 
-                        if (psImageOverlaySection(convKernels, kernel, (i + 2) * fullSize,
-                                                  (j + 2) * fullSize, "=") == 0) {
+                        if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize,
+                                                  (j + KERNEL_MOSAIC) * fullSize, "=") == 0) {
                             psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image.");
                             psFree(kernel);
@@ -399 +407 @@
                 psString comment = NULL; // Comment for metadata
                 psStringAppend(&comment, "Subtraction kernel for region %s", regionString);
-                psMetadataAddImage(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
+                psMetadataAddImage(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE",
                                    PS_META_DUPLICATE_OK, comment, convKernels);
                 psFree(comment);
@@ -427 +435 @@
 
             psTrace("psModules.imcombine", 2, "Convolving...\n");
-            if (!pmSubtractionConvolve(convolved, ro1, ro2, subMask, maskBlank, region,
-                                       solution, kernels, mode, useFFT)) {
+            if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to convolve image.");
                 goto MATCH_ERROR;
@@ -439 +446 @@
                 psString comment = NULL; // Comment for metadata
                 psStringAppend(&comment, "Subtraction solution for region %s", regionString);
-                psMetadataAddVector(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
+                psMetadataAddVector(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.SOLUTION",
                                     PS_META_DUPLICATE_OK, comment, solution);
                 psFree(comment);
                 if (region) {
-                    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
+                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
                 } else {
                     region = psRegionAlloc(0, numCols, 0, numRows);
-                    psMetadataAddPtr(convolved->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
+                    psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.REGION",
                                      PS_META_DUPLICATE_OK | PS_DATA_REGION, comment, region);
                     psFree(region);
@@ -458 +465 @@
 
             // There is data in the readout now
-            convolved->data_exists = true;
-            if (convolved->parent) {
-                convolved->parent->data_exists = true;
-                convolved->parent->parent->data_exists = true;
+            conv1->data_exists = true;
+            if (conv1->parent) {
+                conv1->parent->data_exists = true;
+                conv1->parent->parent->data_exists = true;
+            }
+            if (mode == PM_SUBTRACTION_MODE_DUAL) {
+                conv2->data_exists = true;
+                if (conv2->parent) {
+                    conv2->parent->data_exists = true;
+                    conv2->parent->parent->data_exists = true;
+                }
             }
         }
@@ -474 +488 @@
     weight = NULL;
 
-    if (!pmSubtractionBorder(convolved->image, convolved->weight, convolved->mask, size, maskBlank)) {
+    if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image.");
         goto MATCH_ERROR;

trunk/psModules/src/imcombine/pmSubtractionMatch.h

@@ -10 +10 @@
 
 /// Match two images
-bool pmSubtractionMatch(pmReadout *convolved, ///< Output convolved data
+bool pmSubtractionMatch(pmReadout *conv1, ///< Output convolved data for image 1
+                        pmReadout *conv2, ///< Output convolved data for image 2
                         const pmReadout *ro1, ///< Image 1
                         const pmReadout *ro2, ///< Image 2

trunk/psModules/src/imcombine/pmSubtractionParams.c

@@ -232 +232 @@
     psVectorInit(orders, maxOrder);
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder,
-                                                                  fwhms, orders); // Kernels
+                                                                  fwhms, orders, mode); // Kernels
     psFree(orders);
     psFree(kernels->description);
@@ -279 +279 @@
             continue;
         }
-        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint, mode)) {
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
             psFree(targets);
@@ -505 +505 @@
 
         kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
-
-        kernels->subIndex = PS_SQR(2 * inner + 1) / 2 + numGaussians;
-        assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
-               kernels->v->data.S32[kernels->subIndex] == 0);
     }
 

trunk/psModules/src/imcombine/pmSubtractionStamps.c

@@ -44 +44 @@
                                  )
 {
-    psFree(stamp->matrix);
-    psFree(stamp->vector);
     psFree(stamp->image1);
     psFree(stamp->image2);
@@ -51 +49 @@
     psFree(stamp->convolutions1);
     psFree(stamp->convolutions2);
+
+    psFree(stamp->matrix1);
+    psFree(stamp->matrix2);
+    psFree(stamp->matrixX);
+    psFree(stamp->vector);
+
 }
 
@@ -103 +107 @@
 
 pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region,
-                                                    float spacing)
+                                                    int footprint, float spacing)
 {
     pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return
@@ -124 +128 @@
 
     for (int y = 0, index = 0; y < yStamps; y++) {
-        int yStart = yMin + y * ySize / (yStamps + 1); // Subregion starts here
-        int yStop = yMin + (y + 1) * ySize / (yStamps + 1) - 1; // Subregion stops here
+        int yStart = yMin + y * ((float)ySize / (float)(yStamps + 1)); // Subregion starts here
+        int yStop = yMin + (y + 1) * ((float)ySize / (float)(yStamps + 1)) - 1; // Subregion stops here
         assert(yStart >= yMin && yStop < yMax);
 
         for (int x = 0; x < xStamps; x++, index++) {
-            int xStart = xMin + x * xSize / (xStamps + 1); // Subregion starts here
-            int xStop = xMin + (x + 1) * xSize / (xStamps + 1) - 1; // Subregion stops here
+            int xStart = xMin + x * ((float)xSize / (float)(xStamps + 1)); // Subregion starts here
+            int xStop = xMin + (x + 1) * ((float)xSize / (float)(xStamps + 1)) - 1; // Subregion stops here
             assert(xStart >= xMin && xStop < xMax);
 
             list->stamps->data[index] = pmSubtractionStampAlloc();
+            psTrace("psModules.imcombine", 6, "Stamp region %d: [%d:%d,%d:%d]\n",
+                    index, xStart, xStop, yStart, yStop);
             list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop);
         }
@@ -141 +147 @@
     list->y = NULL;
     list->flux = NULL;
+    list->footprint = footprint;
 
     return list;
@@ -155 +162 @@
     stamp->xNorm = NAN;
     stamp->yNorm = NAN;
-    stamp->matrix = NULL;
-    stamp->vector = NULL;
     stamp->status = PM_SUBTRACTION_STAMP_INIT;
 
@@ -165 +170 @@
     stamp->convolutions2 = NULL;
 
+    stamp->matrix1 = NULL;
+    stamp->matrix2 = NULL;
+    stamp->matrixX = NULL;
+    stamp->vector = NULL;
+
     return stamp;
 }
@@ -171 +181 @@
 pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image,
                                                 const psImage *subMask, const psRegion *region,
-                                                float threshold, float spacing, pmSubtractionMode mode)
+                                                float threshold, int footprint, float spacing,
+                                                pmSubtractionMode mode)
 {
     PS_ASSERT_IMAGE_NON_NULL(image, NULL);
@@ -180 +191 @@
         PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL);
     }
+    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
     PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL);
     if (region) {
@@ -201 +213 @@
 
     if (!stamps) {
-        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, spacing);
+        stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing);
     }
 
@@ -244 +256 @@
             fluxStamp = threshold;
             psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest
-            for (int y = subRegion->y0; y <= subRegion->y1 ; y++) {
-                for (int x = subRegion->x0; x <= subRegion->x1 ; x++) {
+            for (int y = subRegion->y0; y <= subRegion->y1; y++) {
+                for (int x = subRegion->x0; x <= subRegion->x1; x++) {
                     if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) {
                         fluxStamp = image->data.F32[y][x];
@@ -289 +301 @@
 pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux,
                                                const psImage *image, const psImage *subMask,
-                                               const psRegion *region, float spacing, pmSubtractionMode mode)
+                                               const psRegion *region, int footprint, float spacing,
+                                               pmSubtractionMode mode)
 
 {
@@ -316 +329 @@
     int numStars = x->n;                // Number of stars
     pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows,
-                                                                 region, spacing); // Stamp list
+                                                                 region, footprint, spacing); // Stamp list
     int numStamps = stamps->num;        // Number of stamps
 
@@ -401 +414 @@
 
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
-                                psImage *weight, int footprint, int kernelSize)
+                                psImage *weight, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -414 +427 @@
     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 = image1->numCols, numRows = image1->numRows; // Size of images
-    int size = kernelSize + footprint; // Size of postage stamps
+    int size = kernelSize + stamps->footprint; // Size of postage stamps
 
     for (int i = 0; i < stamps->num; i++) {
@@ -467 +479 @@
 
 #if 0
-bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int footprint, int kernelSize)
+bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
     PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false);
-    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
     PS_ASSERT_INT_NONNEGATIVE(kernelSize, false);
 
-    int size = kernelSize + footprint; // Size of postage stamps
+    int size = kernelSize + stamps->footprint; // Size of postage stamps
     int num = stamps->num;              // Number of stamps
     float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma
@@ -514 +525 @@
 
 pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask,
-                                                          const psRegion *region, float spacing,
-                                                          pmSubtractionMode mode)
+                                                          const psRegion *region, int footprint,
+                                                          float spacing, pmSubtractionMode mode)
 {
     PS_ASSERT_ARRAY_NON_NULL(sources, NULL);
@@ -539 +550 @@
 
     pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region,
-                                                            spacing, mode); // Stamps to return
+                                                            footprint, spacing, mode); // Stamps to return
     psFree(x);
     psFree(y);
@@ -553 +564 @@
 
 pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *subMask,
-                                                       const psRegion *region, float spacing,
+                                                       const psRegion *region, int footprint, float spacing,
                                                        pmSubtractionMode mode)
 {
@@ -572 +583 @@
     psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32));
 
-    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, spacing,
-                                                            mode);
+    pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, footprint,
+                                                            spacing, mode);
     psFree(data);
 

trunk/psModules/src/imcombine/pmSubtractionStamps.h

@@ -23 +23 @@
     psArray *x, *y;                     ///< Coordinates for possible stamps (or NULL)
     psArray *flux;                      ///< Fluxes for possible stamps (or NULL)
+    int footprint;                      ///< Half-size of stamps
 } pmSubtractionStampList;
 
@@ -29 +30 @@
                                                     int numRows, // Number of rows in image
                                                     const psRegion *region, // Region for stamps, or NULL
+                                                    int footprint, // Half-size of stamps
                                                     float spacing // Rough average spacing between stamps
     );
@@ -40 +42 @@
     PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \
     PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \
+    PS_ASSERT_INT_NONNEGATIVE((LIST)->footprint, RETURNVALUE); \
     if ((LIST)->x || (LIST)->y || (LIST)->flux) { \
         PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \
@@ -57 +60 @@
     psKernel *image1;                   ///< Reference image postage stamp
     psKernel *image2;                   ///< Input image postage stamp
-    psKernel *weight;                   ///< Weight image postage stamp
-    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
+    psKernel *weight;                   ///< Weight image postage stamp, or NULL
+    psArray *convolutions1;             ///< Convolutions of image 1 for each kernel component, or NULL
+    psArray *convolutions2;             ///< Convolutions of image 2 for each kernel component, or NULL
+    psImage *matrix1, *matrix2;         ///< Matrices for each image, or NULL
+    psImage *matrixX;                   ///< Cross-matrix (for mode DUAL), or NULL
+    psVector *vector;                   ///< Associated vector (when mode not DUAL), or NULL
     pmSubtractionStampStatus status;    ///< Status of stamp
 } pmSubtractionStamp;
@@ -74 +78 @@
                                                 const psRegion *region, ///< Region to search, or NULL
                                                 float threshold, ///< Threshold for stamps in the image
+                                                int footprint, ///< Half-size for stamps
                                                 float spacing, ///< Rough spacing for stamps
                                                 pmSubtractionMode mode ///< Mode for subtraction
@@ -85 +90 @@
                                                const psImage *mask, ///< Mask, or NULL
                                                const psRegion *region, ///< Region to search, or NULL
+                                               int footprint, ///< Half-size for stamps
                                                float spacing, ///< Rough spacing for stamps
                                                pmSubtractionMode mode ///< Mode for subtraction
@@ -94 +100 @@
     const psImage *subMask,             ///< Mask, or NULL
     const psRegion *region,             ///< Region to search, or NULL
+    int footprint,                      ///< Half-size for stamps
     float spacing,                      ///< Rough spacing for stamps
     pmSubtractionMode mode              ///< Mode for subtraction
@@ -103 +110 @@
     const psImage *subMask,             ///< Mask, or NULL
     const psRegion *region,             ///< Region to search, or NULL
+    int footprint,                      ///< Half-size for stamps
     float spacing,                      ///< Rough spacing for stamps
     pmSubtractionMode mode              ///< Mode for subtraction
@@ -109 +117 @@
 /// Extract stamps from the images
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps
-                                psImage *reference, ///< Reference image
-                                psImage *input, ///< Input image (or NULL)
+                                psImage *image1, ///< Reference image
+                                psImage *image2, ///< Input image (or NULL)
                                 psImage *weight, ///< Weight (variance) map
-                                int footprint, ///< Stamp footprint size
                                 int kernelSize ///< Kernel half-size
     );

trunk/psModules/src/psmodules.h

@@ -82 +82 @@
 #include <pmSubtractionMatch.h>
 #include <pmSubtractionParams.h>
+#include <pmSubtractionMask.h>
+#include <pmSubtractionEquation.h>
 #include <pmReadoutCombine.h>