Index: trunk/psModules/src/imcombine/pmSubtractionParams.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionParams.c	(revision 14671)
+++ trunk/psModules/src/imcombine/pmSubtractionParams.c	(revision 14709)
@@ -9,7 +9,10 @@
 
 #include "pmSubtractionStamps.h"
+#include "pmSubtraction.h"
+#include "pmSubtractionParams.h"
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+#if 0
 // Convolve the reference stamp by the kernel
 static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved
@@ -45,5 +48,5 @@
     return convolution;
 }
-
+#endif
 
 // Accumulate cross-term sums for a stamp
@@ -53,20 +56,16 @@
                             const pmSubtractionStamp *stamp, // Stamp with weight
                             const psKernel *input, // Input image, I(x)
-                            const psKernel *convolution // Convolution, conv(x)
-    )
-{
-    // Range of convolution
-    int xMin = convolution->xMin;
-    int xMax = convolution->xMax;
-    int yMin = convolution->yMin;
-    int yMax = convolution->yMax;
-
+                            int kernelIndex, // Index for kernel component
+                            int footprint // Size of region of interest
+    )
+{
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-
-    for (int y = yMin; y <= yMax; y++) {
-        psF32 *in = &input->kernel[y][xMin]; // Dereference input
-        psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
-        psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
-        for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) {
+    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+
+    for (int y = -footprint; y <= footprint; y++) {
+        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
+        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
+        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
             double temp = *in / *wt; // Temporary product
             *sumI += temp;
@@ -82,19 +81,15 @@
                                    double *sumCC, // Sum of conv(x)^2/sigma(x)^2
                                    const pmSubtractionStamp *stamp, // Stamp with input and weight
-                                   const psKernel *convolution // Convolution, conv(x)
-    )
-{
-    // Range of convolution
-    int xMin = convolution->xMin;
-    int xMax = convolution->xMax;
-    int yMin = convolution->yMin;
-    int yMax = convolution->yMax;
-
+                                   int kernelIndex, // Index for kernel component
+                                   int footprint // Size of region of interest
+    )
+{
     psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
-
-    for (int y = yMin; y <= yMax; y++) {
-        psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
-        psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
-        for (int x = xMin; x <= xMax; x++, wt++, conv++) {
+    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+
+    for (int y = -footprint; y <= footprint; y++) {
+        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
+        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
+        for (int x = -footprint; x <= footprint; x++, wt++, conv++) {
             double convNoise = *conv / *wt; // Temporary product
             *sumC += convNoise;
@@ -106,24 +101,20 @@
 
 static double accumulateChi2(psKernel *input, // Input stamp
-                             psKernel *convolution, // Convolution, R(x)*k(u)
                              pmSubtractionStamp *stamp, // Stamp with weight
+                             int kernelIndex, // Index for kernel component
                              double coeff, // Coefficient of convolution
-                             double bg  // Background term
-    )
-{
-    // Range of convolution
-    int xMin = convolution->xMin;
-    int xMax = convolution->xMax;
-    int yMin = convolution->yMin;
-    int yMax = convolution->yMax;
-
+                             double bg,  // Background term
+                             int footprint // Size of region of interest
+    )
+{
     double chi2 = 0.0;
-    psKernel *weight = stamp->weight;
-
-    for (int y = yMin; y <= yMax; y++) {
-        psF32 *in = &input->kernel[y][xMin]; // Dereference input
-        psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
-        psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
-        for (int x = xMin; x <= xMax; x++, in++, wt++, conv++) {
+    psKernel *weight = stamp->weight;   // Weight, sigma(x)^2
+    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+
+    for (int y = -footprint; y <= footprint; y++) {
+        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight
+        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
+        for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) {
             chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt;
         }
@@ -159,22 +150,16 @@
 static void subtractConvolution(psKernel *input, // Input stamp
                                 const pmSubtractionStamp *stamp, // Stamp with weight
-                                const psKernel *convolution, // Convolution, R(x)*k(u)
+                                int kernelIndex, // Index for kernel component
                                 float coeff, // Coefficient of subtraction
-                                float bg // Background term
-    )
-{
-    // Range of convolution
-    int xMin = convolution->xMin;
-    int xMax = convolution->xMax;
-    int yMin = convolution->yMin;
-    int yMax = convolution->yMax;
-
-    psKernel *weight = stamp->weight;   // Weight map
-
-    for (int y = yMin; y <= yMax; y++) {
-        psF32 *in = &input->kernel[y][xMin]; // Dereference input
-        psF32 *conv = &convolution->kernel[y][xMin]; // Dereference convolution
-        psF32 *wt = &weight->kernel[y][xMin]; // Dereference weight
-        for (int x = xMin; x <= xMax; x++, in++, conv++, wt++) {
+                                float bg, // Background term
+                                int footprint // Size of region of interest
+    )
+{
+    psKernel *convolution = stamp->convolutions->data[kernelIndex]; // Convolution of interest
+
+    for (int y = -footprint; y <= footprint; y++) {
+        psF32 *in = &input->kernel[y][-footprint]; // Dereference input
+        psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution
+        for (int x = -footprint; x <= footprint; x++, in++, conv++) {
             *in -= *conv * coeff + bg;
         }
@@ -229,18 +214,18 @@
         psKernel *input = stamp->input; // Input image of interest
         psImage *copy = psImageCopy(NULL, input->image, PS_TYPE_F32); // Copy of the image
-        inputs->data[i] = psKernelAllocFromImage(copy, size + footprint + 1, size + footprint + 1);
+        inputs->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint);
         psFree(copy);                   // Drop reference
     }
 
     // Generate the convolutions
-    psArray *convolutions = psArrayAlloc(numKernels); // Convolutions per kernel component and stamp
-    for (int i = 0; i < numKernels; i++) {
-        psKernel *kernel = kernels->preCalc->data[i]; // Precalculated kernel
-
-        psArray *stampConv = psArrayAlloc(numStamps); // Convolutions for each stamp
-        convolutions->data[i] = stampConv;
-
-        for (int j = 0; j < numStamps; j++) {
-            stampConv->data[j] = convolveStamp(stamps->data[j], kernel, footprint);
+    for (int i = 0; i < numStamps; i++) {
+        pmSubtractionStamp *stamp = stamps->data[i]; // Stamp of interest
+        if (!stamp->convolutions) {
+            stamp->convolutions = psArrayAlloc(numKernels);
+        }
+
+        for (int j = 0; j < numKernels; j++) {
+            psKernel *kernel = kernels->preCalc->data[j]; // Precalculated kernel
+            stamp->convolutions->data[j] = p_pmSubtractionConvolveStampPrecalc(stamp->reference, kernel);
         }
     }
@@ -264,8 +249,6 @@
     psVectorInit(sumCC, 0.0);
     for (int i = 0; i < numKernels; i++) {
-        psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp
         for (int j = 0; j < numStamps; j++) {
-            accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j],
-                                   stampsConv->data[j]);
+            accumulateConvolutions(&sumC->data.F64[i], &sumCC->data.F64[i], stamps->data[j], i, footprint);
         }
     }
@@ -298,8 +281,6 @@
             double sumIC = 0.0;         // sum of I(x)C(x)/sigma(x)^2
 
-            psArray *stampsConv = convolutions->data[i]; // Convolutions per stamp
-
             for (int j = 0; j < numStamps; j++) {
-                accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], stampsConv->data[j]);
+                accumulateCross(&sumI, &sumII, &sumIC, stamps->data[j], inputs->data[j], i, footprint);
             }
 
@@ -310,5 +291,5 @@
             double chi2 = 0.0;          // Chi^2
             for (int j = 0; j < numStamps; j++) {
-                chi2 += accumulateChi2(inputs->data[j], stampsConv->data[j], stamps->data[j], coeff, bg);
+                chi2 += accumulateChi2(inputs->data[j], stamps->data[j], i, coeff, bg, footprint);
             }
 
@@ -329,8 +310,6 @@
         ranking->data.S32[bestIndex] = iter;
         // Remove its contribution, and don't include it in the future.
-        psArray *stampsConv = convolutions->data[bestIndex]; // Convolutions per stamp
         for (int j = 0; j < numStamps; j++) {
-            subtractConvolution(inputs->data[j], stamps->data[j], stampsConv->data[j],
-                                           bestCoeff, bestBG);
+            subtractConvolution(inputs->data[j], stamps->data[j], bestIndex, bestCoeff, bestBG, footprint);
         }
 
@@ -349,5 +328,4 @@
     }
     psFree(inputs);
-    psFree(convolutions);
     psFree(sumC);
     psFree(sumCC);
@@ -392,6 +370,8 @@
         psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
         for (int i = 1; i < newSize; i++) {
-            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-            psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
+                psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
+                psBinaryOp(kernel->image, kernel->image, "-", subtract->image);
+            }
         }
     } else if (type == PM_SUBTRACTION_KERNEL_GUNK) {
@@ -400,6 +380,8 @@
 
         for (int i = 0; i < newSize; i++) {
-            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-            kernel->kernel[0][0] -= 1.0;
+            if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
+                psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
+                kernel->kernel[0][0] -= 1.0;
+            }
         }