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

Timestamp:

Aug 29, 2007, 5:50:28 PM (19 years ago)

Author:

Paul Price

Message:

Extensive changes to image subtraction codes. Discovered that I
wasn't (again) implementing the algorithm correctly --- I was
accumulating

sum_x,y C_i(x,y) sum_x,y C_j(x,y)

instead of

sum_x,y C_i(x,y) C_j(x,y)

i.e., I was doing the sums separately. I'm amazed that it was
actually getting fairly decent results before, but now it's doing it
properly. As part of this fix, added in, for each stamp, images
corresponding to the convolution of the reference with each of the
kernel components. This allows fast rejection (instead of generating
the kernel and convolving, can simply sum the individual convolutions
with the correct scaling).
Added use of FFT to do convolutions. This involved changing around
the convolvePixel function (now convolveRef, since it generates the
entire image for the convolution of the reference stamp with a kernel
component, rather than a single pixel as previously). Tested with
POIS and ISIS kernels, but no others yet. FFT is nice and fast for
the large (size = 10) kernels I've been using, so I also moved it in
to the convolveRef function for appropriate kernel types.
As part of the FFT move, discovered that my convolutions were around
the wrong way --- supposed to be convolved = image[y][x] * conv[y - v][x - u]
instead of conv[y + v][x + u]. However, this problem didn't manifest
itself previously because everything was consistently around the wrong
way. The FFT does it the Correct Way, so it revealed the problem. This
is now fixed throughout the image subtraction code.

File:

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

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -4 +4 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.48 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-23 23:43:12 $
+ *  @version $Revision: 1.49 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-30 03:50:28 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -17 +17 @@
 #include <stdio.h>
 #include <math.h>
+#include <string.h>
 #include <pslib.h>
 
@@ -24 +25 @@
 
 #include "pmSubtraction.h"
+
+//#define TESTING
 
 #define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
@@ -216 +219 @@
 }
 
-// Generate the convolved pixel value
-static inline double convolvePixel(const pmSubtractionKernels *kernels, // Kernel basis functions
-                                   int index, // Kernel basis function index
-                                   int x, int y, // Pixel around which to convolve
-                                   const psKernel *image // Image to convolve (a kernel for convenience)
-                                   )
+// Subtract the (0,0) element to preserve photometric scaling
+static void convolveSub(psKernel *convolved, // Convolved image
+                        const psKernel *image, // Image being convolved
+                        int footprint  // Size of region of interest
+                        )
+{
+    // Can't use psBinaryOp because the images are of different size
+    for (int y = -footprint; y <= footprint; y++) {
+        for (int x = -footprint; x <= footprint; x++) {
+            convolved->kernel[y][x] -= image->kernel[y][x];
+        }
+    }
+    return;
+}
+
+// Generate the convolution given some offset
+static psKernel *convolveOffset(const psKernel *image, // Image to convolve (a kernel for convenience)
+                                int u, int v, // Offset to apply
+                                int footprint // Size of region of interest
+                                )
+{
+    psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image
+    int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+    for (int y = -footprint; y <= footprint; y++) {
+        // Convolution with a delta function is just the value specified by the offset
+        memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes);
+    }
+    return convolved;
+}
+
+// Generate the convolution given a precalculated kernel
+static psKernel *convolvePrecalc(const psKernel *image, // Image to convolve (a kernel for convenience)
+                                 const psKernel *kernel, // Kernel by which to convolve
+                                 int footprint // Size of region of interest
+                                 )
+{
+    psImage *conv = psImageConvolveFFT(image->image, kernel, 0.0); // Convolved image
+    int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image
+    psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version
+    psFree(conv);
+    return convolved;
+}
+
+// Generate the convolved reference image
+static psKernel *convolveRef(const pmSubtractionKernels *kernels, // Kernel basis functions
+                             int index, // Kernel basis function index
+                             const psKernel *image, // Image to convolve (a kernel for convenience)
+                             int footprint // Size of region of interest
+                             )
 {
     switch (kernels->type) {
       case PM_SUBTRACTION_KERNEL_POIS: {
-          // Convolution with a delta function is just the value specified by the offset
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          float value = image->kernel[y + v][x + u]; // Value of convolution
+          psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-              value -= image->kernel[y][x];
+              convolveSub(convolved, image, footprint);
           }
-          return value;
+          return convolved;
       }
         // Method for SPAM and FRIES is the same
       case PM_SUBTRACTION_KERNEL_SPAM:
       case PM_SUBTRACTION_KERNEL_FRIES: {
+          psKernel *convolved = psKernelAlloc(-footprint, footprint,
+                                              -footprint, footprint); // Convolved image
           int uStart = kernels->u->data.S32[index];
           int uStop = kernels->uStop->data.S32[index];
           int vStart = kernels->v->data.S32[index];
           int vStop = kernels->vStop->data.S32[index];
-          double sum = 0.0;
-          for (int v = vStart; v <= vStop; v++) {
-              for (int u = uStart; u <= uStop; u++) {
-                  sum += image->kernel[y + v][x + u];
+          float norm = 1.0 / (uStop - uStart + 1) * (vStop - vStart + 1); // Normalisation
+          for (int y = -footprint; y <= footprint; y++) {
+              for (int x = -footprint; x <= footprint; x++) {
+                  double sum = 0.0;
+                  for (int v = vStart; v <= vStop; v++) {
+                      for (int u = uStart; u <= uStop; u++) {
+                          sum += image->kernel[y - v][x - u];
+                      }
+                  }
+                  convolved->kernel[y][x] = norm * sum;
               }
           }
-          sum /= (uStop - uStart + 1) * (vStop - vStart + 1); // Normalising sum of kernel component to unity
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-              sum -= image->kernel[y][x];
+              convolveSub(convolved, image, footprint);
           }
-          return sum;
+          return convolved;
       }
       case PM_SUBTRACTION_KERNEL_GUNK: {
           if (index < kernels->inner) {
-              // Using pre-calculated function
-              psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel
-              int size = kernels->size;     // Kernel half-size
-              double sum = 0.0;             // Accumulated sum from convolution
-              for (int v = -size; v <= size; v++) {
-                  for (int u = -size; u <= size; u++) {
-                      sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
-                  }
-              }
-              return sum;
+              // Photometric scaling is already built in to the precalculated kernel
+              return convolvePrecalc(image, kernels->preCalc->data[index], footprint);
           }
           // Using delta function
           int u = kernels->u->data.S32[index]; // Offset in x
           int v = kernels->v->data.S32[index]; // Offset in y
-          float value = image->kernel[y + v][x + u]; // Value of convolution
-          // The (0,0) delta function is subtracted from most kernels to preserve photometric scaling
+          psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image
           if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
-              value -= image->kernel[y][x];
+              convolveSub(convolved, image, footprint);
           }
-          return value;
+          return convolved;
       }
       case PM_SUBTRACTION_KERNEL_ISIS: {
-          psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel
-          int size = kernels->size;     // Kernel half-size
-          double sum = 0.0;             // Accumulated sum from convolution
-          for (int v = -size; v <= size; v++) {
-              for (int u = -size; u <= size; u++) {
-                  sum += kernel->kernel[v][u] * image->kernel[y + v][x + u];
-                  // Photometric scaling is already built in to the precalculated kernel
-              }
-          }
-          return sum;
+          // Photometric scaling is already built in to the precalculated kernel
+          return convolvePrecalc(image, kernels->preCalc->data[index], footprint);
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
+          psKernel *convolved = psKernelAlloc(-footprint, footprint,
+                                              -footprint, footprint); // Convolved image
           if (index == kernels->subIndex) {
-              return image->kernel[y][x];
+              // Simply copying over the image
+              return convolveOffset(image, 0, 0, footprint);
           }
+
           psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data
           psVector *uCoords = preCalc->data[0]; // u coordinates
@@ -299 +336 @@
           psVector *poly = preCalc->data[2]; // Polynomial values
           int num = uCoords->n;         // Number of pixels
-          double sum = 0.0;             // Accumulated sum from convolution
-          for (int j = 0; j < num; j++) {
-              int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
-              sum += image->kernel[y + v][x + u] * poly->data.F32[j];
+          for (int y = -footprint; y <= footprint; y++) {
+              for (int x = -footprint; x <= footprint; x++) {
+                  double sum = 0.0;             // Accumulated sum from convolution
+                  for (int j = 0; j < num; j++) {
+                      int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates
+                      sum += image->kernel[y - v][x - u] * poly->data.F32[j];
+                  }
+                  convolved->kernel[y][x] = sum;
+                  // Photometric scaling is built into the kernel --- no subtraction!
+              }
           }
-          // Photometric scaling is built into the kernel --- no subtraction!
-          return sum;
+          return convolved;
       }
       default:
         psAbort("Should never get here.");
     }
-    return NAN;
+    return NULL;
 }
 
@@ -460 +502 @@
     int numParams = numKernels * numSpatial + numBackground;
     int bgIndex = numParams - numBackground; // Index in matrix for the background
-    psVector *convolutions = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); // Convolutions
 
     // We iterate over each stamp, allocate the matrix and vectors if
@@ -466 +507 @@
     for (int i = 0; i < stamps->n; i++) {
         pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
-        if (stamp->status == PM_SUBTRACTION_STAMP_CALCULATE) {
-            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);
-
-            // Spatial polynomial terms
-            psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm);
-
-            psKernel *reference = stamp->reference; // Reference postage stamp
-            psKernel *input = stamp->input; // Input postage stamp
-            psKernel *weight = stamp->weight; // Weight map postage stamp
-
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = - footprint; x <= footprint; x++) {
-                    float invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse square noise
-
-                    // Generate the convolutions
-                    for (int j = 0; j < numKernels; j++) {
-                        double value = convolvePixel(kernels, j, x, y, reference); // Value from convolution
-                        // Generate the pseudo-convolutions from the spatial polynomial terms
-                        for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
-                            for (int xOrder = 0; xOrder <= spatialOrder - yOrder;
-                                 xOrder++, index += numKernels) {
-                                convolutions->data.F64[index] = value * polyValues->data.F64[yOrder][xOrder];
+        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);
+
+        psKernel *reference = stamp->reference; // Reference postage stamp
+        psKernel *input = stamp->input; // Input postage stamp
+        psKernel *weight = stamp->weight; // Weight map postage stamp
+
+        // Generate convolutions of the reference
+        psArray *convolutions = stamp->convolutions; // Convolutions of the reference for each kernel
+        if (!convolutions) {
+            stamp->convolutions = convolutions = psArrayAlloc(numKernels);
+        }
+        for (int j = 0; j < numKernels; j++) {
+            if (convolutions->data[j]) {
+                psFree(convolutions->data[j]);
+            }
+            convolutions->data[j] = convolveRef(kernels, j, reference, footprint);
+#ifdef TESTING
+            {
+                psKernel *conv = convolutions->data[j]; // Convolution of interest
+                psString filename = NULL;
+                psStringAppend(&filename, "conv_%03d_%03d.fits", i, j);
+                psFits *fits = psFitsOpen(filename, "w");
+                psFree(filename);
+                psFitsWriteImage(fits, NULL, conv->image, 0, NULL);
+                psFitsClose(fits);
+            }
+#endif
+       }
+
+        psImage *polyValues = spatialPolyValues(spatialOrder, stamp->xNorm, stamp->yNorm); // Polynomial terms
+
+        // Generate least-squares vector and matrix
+        for (int j = 0; j < numKernels; j++) {
+            psKernel *jConv = convolutions->data[j]; // Convolution for j-th element
+
+            // Generate upper diagonals of matrix
+            for (int k = 0; k < numKernels; 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];
+                    }
+                }
+                // 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) {
+                                stampMatrix->data.F64[jIndex][kIndex] = sumCC *
+                                    polyValues->data.F64[jyOrder][jxOrder] *
+                                    polyValues->data.F64[kyOrder][kxOrder];
                             }
                         }
                     }
-
-                    // Generate the least-squares matrix and vector
-                    // Upper diagonal only
-                    for (int i = 0; i < bgIndex; i++) {
-                        for (int j = i; j < bgIndex; j++) {
-                            stampMatrix->data.F64[i][j] += convolutions->data.F64[i] *
-                                convolutions->data.F64[j] * invNoise2;
+                }
+            }
+
+            // Vector and background term for matrix
+            double sumC = 0.0;      // Sum of the convolution
+            double sumIC = 0.0;     // Sum of the convolution/input product
+            for (int y = - footprint; y <= footprint; y++) {
+                for (int x = - footprint; x <= footprint; x++) {
+                    double convProduct = jConv->kernel[y][x] / weight->kernel[y][x]; // Convolution / noise^2
+                    sumC += convProduct;
+                    sumIC += input->kernel[y][x] * convProduct;
+                }
+            }
+            for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) {
+                for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder;
+                     jxOrder++, jIndex += numKernels) {
+                    stampVector->data.F64[jIndex] = sumIC * polyValues->data.F64[jyOrder][jxOrder];
+                    stampMatrix->data.F64[jIndex][bgIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder];
+                    stampMatrix->data.F64[bgIndex][jIndex] = sumC * polyValues->data.F64[jyOrder][jxOrder];
+                }
+            }
+
+        }
+        psFree(polyValues);
+
+        // Background only terms
+        double sum1 = 0.0;              // Sum of the weighting
+        double sumI = 0.0;              // Sum of the input
+        for (int y = - footprint; y <= footprint; y++) {
+            for (int x = - footprint; x <= footprint; x++) {
+                double invNoise2 = 1.0 / weight->kernel[y][x]; // Inverse noise, squared
+                sum1 += invNoise2;
+                sumI += input->kernel[y][x] * invNoise2;
+            }
+        }
+        stampMatrix->data.F64[bgIndex][bgIndex] = sum1;
+        stampVector->data.F64[bgIndex] = sumI;
+
+        // Fill in lower diagonals of symmetric matrix, while checking for bad values
+        // Note, there are two symmetries going on here --- the matrix is C_i C_j P_i P_j
+        bool bad = false;           // Are there bad values?
+#if 0
+        for (int j = 0; j < numKernels; j++) {
+            for (int jSpatial = 0, jIndex = j; jSpatial < numSpatial; jSpatial++, jIndex += numKernels) {
+                for (int k = 0; k < j; k++) {
+                    for (int kSpatial = 0, kIndex = k; kSpatial < numSpatial;
+                         kSpatial++, kIndex += numKernels) {
+                        double value = stampMatrix->data.F64[kIndex][jIndex]; // Value of matrix
+                        stampMatrix->data.F64[jIndex][kIndex] = value;
+                        if (!isfinite(value)) {
+                            bad = true;
                         }
-                        stampVector->data.F64[i] += input->kernel[y][x] * convolutions->data.F64[i] *
-                            invNoise2;
-
-                        // Background term
-                        stampMatrix->data.F64[i][bgIndex] += convolutions->data.F64[i] * invNoise2;
                     }
-                    // Background only terms
-                    stampMatrix->data.F64[bgIndex][bgIndex] += invNoise2;
-                    stampVector->data.F64[bgIndex] += input->kernel[y][x] * invNoise2;
-                }
-            }
-            psFree(polyValues);
-
-            // Fill in lower diagonal of symmetric matrix, while checking for bad values
-            bool bad = false;           // Are there bad values?
-            for (int i = 0; i < bgIndex; i++) {
-                for (int j = 0; j < i; j++) {
-                    stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i];
-                    if (!isfinite(stampMatrix->data.F64[j][i])) {
-                        bad = true;
-                    }
-                }
-                stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex];
-                if (!isfinite(stampMatrix->data.F64[i][bgIndex]) ||
-                    !isfinite(stampMatrix->data.F64[i][i]) ||
-                    !isfinite(stampVector->data.F64[i])) {
+                }
+
+                stampMatrix->data.F64[bgIndex][jIndex] = stampMatrix->data.F64[jIndex][bgIndex];
+                if (!isfinite(stampMatrix->data.F64[jIndex][bgIndex]) ||
+                    !isfinite(stampMatrix->data.F64[jIndex][jIndex]) ||
+                    !isfinite(stampVector->data.F64[jIndex])) {
                     bad = true;
                 }
             }
-            if (!isfinite(stampVector->data.F64[bgIndex])) {
-                bad = true;
-            }
-
-            if (bad) {
-                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
-                psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
-                        i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
-            } 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);
-            }
-
-        }
-    }
-
-    psFree(convolutions);
+        }
+#endif
+        if (!isfinite(stampVector->data.F64[bgIndex])) {
+            bad = true;
+        }
+
+        if (bad) {
+            stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d) because of bad equation\n",
+                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5));
+        } 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;
@@ -665 +765 @@
     double totalSquareDev = 0.0;        // Total square deviation from zero
     int numStamps = 0;                  // Number of used stamps
+    int numKernels = kernels->num;      // Number of kernels
+    int spatialOrder = kernels->spatialOrder; // Order of kernel spatial variations
+    double devNorm = 1.0 / PS_SQR(2 * footprint + 1); // Normalisation for deviations
     {
-        psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
-        psKernel *convolution = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolution
-        psKernel *kernelImage = NULL;       // The kernel, with which to convolve the stamps
         float background = solution->data.F64[solution->n-1]; // The difference in background
-        int size = kernels->size;       // Half-size of the kernel
 
         for (int i = 0; i < stamps->n; i++) {
@@ -681 +780 @@
                                                     stamp->yNorm); // Polynomial terms
 
-            psKernel *reference = stamp->reference; // Reference postage stamp
-
-            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, false);
+            psKernel *input = stamp->input; // Input postage stamp
+            psKernel *weight = stamp->weight; // Weight postage stamp
+            psArray *convolutions = stamp->convolutions; // Convolutions of reference image for each kernel
+#ifdef TESTING
+            psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint);
+#endif
+            float deviation = 0.0;      // Deviation for this stamp
             for (int y = - footprint; y <= footprint; y++) {
                 for (int x = - footprint; x <= footprint; x++) {
-                    convolution->kernel[y][x] = background;
-                    for (int v = -size; v <= size; v++) {
-                        for (int u = -size; u <= size; u++) {
-                            convolution->kernel[y][x] += kernelImage->kernel[v][u] *
-                                reference->kernel[y + v][x + u];
+                    float conv = background; // The value of the convolution
+                    for (int j = 0; j < numKernels; j++) {
+                        psKernel *convolution = convolutions->data[j]; // Convolution of reference
+                        double polynomial = 0.0; // Value of the polynomial
+                        for (int yOrder = 0, index = j; yOrder <= spatialOrder; yOrder++) {
+                            for (int xOrder = 0; xOrder <= spatialOrder - yOrder;
+                                 xOrder++, index += numKernels) {
+                                polynomial += solution->data.F64[index] *
+                                    polyValues->data.F64[yOrder][xOrder];
+                            }
                         }
+                        conv += convolution->kernel[y][x] * polynomial;
                     }
+                    float diff = input->kernel[y][x] - conv;
+                    deviation += PS_SQR(diff) / weight->kernel[y][x];
+#ifdef TESTING
+                    residual->kernel[y][x] = diff;
+#endif
                 }
             }
             psFree(polyValues);
 
-            psImage *convolvedStamp = convolution->image; // Image of the convolution
-            psImage *input = stamp->input->image; // Input image postage stamp
-            psImage *weight = stamp->weight->image; // Weight image postage stamp
-
-            // Region of interest
-            psRegion inRegion = psRegionSet(input->col0 + size, input->col0 + size + 2 * footprint + 1,
-                                          input->row0 + size, input->row0 + size + 2 * footprint + 1);
-            psRegion wtRegion = (input == weight ? inRegion :
-                                 psRegionSet(weight->col0 + size, weight->col0 + size + 2 * footprint + 1,
-                                             weight->row0 + size, weight->row0 + size + 2 * footprint + 1));
-
-            psImage *inStamp = psImageSubset(stamp->input->image, inRegion); // Image of stamp
-            psImage *wtStamp = psImageSubset(stamp->weight->image, wtRegion); // Image of stamp
-            assert(convolvedStamp->numCols == inStamp->numCols &&
-                   convolvedStamp->numRows == inStamp->numRows);
-            assert(convolvedStamp->numCols == wtStamp->numCols &&
-                   convolvedStamp->numRows == wtStamp->numRows);
-            (void)psBinaryOp(convolvedStamp, inStamp, "-", convolvedStamp);
-            (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp);
-            (void)psBinaryOp(convolvedStamp, convolvedStamp, "/", wtStamp);
-            psFree(inStamp);
-            psFree(wtStamp);
-            if (!psImageStats(stats, convolvedStamp, NULL, 0)) {
-                psError(PS_ERR_UNKNOWN, false,
-                        "Unable to calculate statistics on normalised residual image of stamp.");
-                psFree(deviations);
-                psFree(convolution);
-                psFree(stats);
-                return -1;
-            }
-            deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0);
-            psTrace("psModules.imcombine", 1, "Deviation for stamp %d (%d,%d): %f\n",
-                    i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]);
+            deviations->data.F32[i] = devNorm * deviation;
+            psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n",
+                    i, (int)stamp->x, (int)stamp->y, deviations->data.F32[i]);
             totalSquareDev += PS_SQR(deviations->data.F32[i]);
             numStamps++;
-        }
-
-        psFree(kernelImage);
-        psFree(convolution);
-        psFree(stats);
+
+#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
+
+        }
     }
 
@@ -823 +916 @@
 }
 
+// Convolve an image using FFT
+static void convolveFFT(psImage *target,// Place the result in here
+                        const psImage *image, // Image to convolve
+                        const psKernel *kernel, // Kernel by which to convolve
+                        psRegion region,// Region of interest
+                        float background, // Background to add
+                        int size        // Size of (square) kernel
+                        )
+{
+    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
+                                  region.y0 - size, region.y1 + size); // Add a border
+
+    // Casting away const so psImageSubset can add the child
+    psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve
+    psImage *convolved = psImageConvolveFFT(subImage, kernel, 0.0); // Convolution
+    psFree(subImage);
+
+    // Now, we have to chop off the borders, and stick it in where it belongs
+    psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size,
+                                                            size, -size)); // Cut off the edges
+    psImage *subTarget = psImageSubset(target, region); // Target for this subregion
+    if (background != 0.0) {
+        psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32));
+    } else {
+        int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+        for (int y = 0; y < subTarget->numRows; y++) {
+            memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes);
+        }
+    }
+    psFree(subConv);
+    psFree(convolved);
+    psFree(subTarget);
+    return;
+}
+
+// Convolve an image directly
+static void convolveDirect(psImage *target, // Put the result here
+                           const psImage *image, // Image to convolve
+                           const psKernel *kernel, // Kernel by which to convolve
+                           psRegion region,// Region of interest
+                           float background, // Background to add
+                           int size        // Size of (square) kernel
+                           )
+{
+    for (int y = region.y0; y < region.y1; y++) {
+        for (int x = region.x0; x < region.x1; x++) {
+            target->data.F32[y][x] = background;
+            for (int v = -size; v <= size; v++) {
+                for (int u = -size; u <= size; u++) {
+                    target->data.F32[y][x] += kernel->kernel[v][u] *
+                        image->data.F32[y - v][x - u];
+                }
+            }
+        }
+    }
+    return;
+}
 
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
                            const psImage *inImage, const psImage *inWeight, const psImage *subMask,
                            psMaskType blank, const psRegion *region, const psVector *solution,
-                           const pmSubtractionKernels *kernels)
+                           const pmSubtractionKernels *kernels, bool useFFT)
 {
     PS_ASSERT_IMAGE_NON_NULL(inImage, false);
@@ -927 +1077 @@
 
     for (int j = yMin; j < yMax; j += fullSize) {
+        int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
         for (int i = xMin; i < xMax; i += fullSize) {
+            int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
@@ -935 +1087 @@
                                            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, false);
             if (inWeight) {
@@ -940 +1093 @@
             }
 
-            for (int y = j; y < PS_MIN(j + fullSize, yMax); y++) {
-                for (int x = i; x < PS_MIN(i + fullSize, xMax); x++) {
-                    // Check and propagate the kernel footprint, if required
+
+            if (useFFT) {
+                // Use Fast Fourier Transform to do the convolution
+                // This provides a big speed-up for large kernels
+                convolveFFT(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
+                            background, size);
+                if (inWeight) {
+                    convolveFFT(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax),
+                                0.0, size);
+                }
+            } else {
+                convolveDirect(convImage, inImage, kernelImage, psRegionSet(i, xSubMax, j, ySubMax),
+                               background, size);
+                if (inWeight) {
+                    convolveDirect(convWeight, inWeight, kernelWeight, psRegionSet(i, xSubMax, j, ySubMax),
+                                   0.0, size);
+                }
+            }
+
+            // Propagate the mask
+            for (int y = j; y < ySubMax; y++) {
+                for (int x = i; x < xSubMax; x++) {
                     if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
                                     (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {
@@ -950 +1122 @@
                             convWeight->data.F32[y][x] = NAN;
                         }
-                    } else {
-                        // Convolve the image
-                        convImage->data.F32[y][x] = background;
-                        for (int v = -size; v <= size; v++) {
-                            for (int u = -size; u <= size; u++) {
-                                convImage->data.F32[y][x] += kernelImage->kernel[v][u] *
-                                    inImage->data.F32[y + v][x + u];
-                            }
-                        }
-
-                        // Convolve the weight (variance) map
-                        if (inWeight) {
-                            for (int v = -size; v <= size; v++) {
-                                for (int u = -size; u <= size; u++) {
-                                    convWeight->data.F32[y][x] += kernelWeight->kernel[v][u] *
-                                        inWeight->data.F32[y + v][x + u];
-                                }
-                            }
-                        }
                     }
                 }
             }
+
         }
     }