Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14716)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14738)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.54 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-30 23:09:52 $
+ *  @version $Revision: 1.55 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-04 23:56:07 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -548,4 +548,32 @@
 }
 
+bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint)
+{
+    PS_ASSERT_PTR_NON_NULL(stamp, false);
+    PS_ASSERT_KERNEL_NON_NULL(stamp->reference, false);
+    PS_ASSERT_PTR_NON_NULL(kernels, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(kernels->u, kernels->v, false);
+    PS_ASSERT_INT_NONNEGATIVE(footprint, false);
+
+    if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Stamp not marked for calculation.");
+        return false;
+    }
+
+    if (stamp->convolutions) {
+        // Already done
+        return true;
+    }
+
+    stamp->convolutions = psArrayAlloc(kernels->num);
+    psKernel *reference = stamp->reference; // Reference postage stamp
+
+    for (int i = 0; i < kernels->num; i++) {
+        stamp->convolutions->data[i] = convolveRef(kernels, i, reference, footprint);
+    }
+
+    return true;
+}
+
 
 bool pmSubtractionCalculateEquation(psArray *stamps, const pmSubtractionKernels *kernels, int footprint)
@@ -591,31 +619,14 @@
         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
+        if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i);
+            return NULL;
+        }
+
         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
 
@@ -900,4 +911,7 @@
             stamp->y = 0;
             stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+            // Recalculate convolutions
+            psFree(stamp->convolutions);
+            stamp->convolutions = NULL;
         }
     }