Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 13379)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 13380)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-15 01:23:19 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-05-15 02:55:05 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -51,4 +51,52 @@
 
     return polyValues;
+}
+
+// Generate an image of the solved kernel
+static 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
+                              float (*weightFunc)(float value) // Function for weighting
+    )
+{
+    int numKernels = solution->n - 1;   // Number of kernel basis functions
+    assert(kernels->u->n == numKernels);
+    assert(kernels->v->n == numKernels);
+    assert(kernels->xOrder->n == numKernels);
+    assert(kernels->yOrder->n == numKernels);
+
+    int size = kernels->size;           // Kernel half-size
+    if (!kernel) {
+        kernel = psKernelAlloc(-size, size, -size, size);
+    }
+    psImageInit(kernel->image, 0.0);
+
+    for (int i = 0; i < numKernels; i++) {
+        int xOrder = kernels->xOrder->data.S32[i]; // Polynomial order in x
+        int yOrder = kernels->yOrder->data.S32[i]; // Polynomial order in y
+        double polyValue = polyValues->data.F64[yOrder][xOrder]; // Value of spatial polynomial
+        switch (kernels->type) {
+          case PM_SUBTRACTION_KERNEL_POIS: {
+              int u = kernels->u->data.S32[i]; // Offset in x
+              int v = kernels->v->data.S32[i]; // Offset in y
+              kernel->kernel[v][u] += weightFunc(polyValue * solution->data.F64[i]);
+              break;
+          }
+          case PM_SUBTRACTION_KERNEL_ISIS: {
+              psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values
+              for (int v = -size; v <= size; v++) {
+                  for (int u = -size; u <= size; u++) {
+                      kernel->kernel[v][u] += weightFunc(polyValue * preCalc->kernel[v][u]);
+                  }
+              }
+              break;
+          }
+          default:
+            psAbort("Should never get here.");
+        }
+    }
+
+    return kernel;
 }
 
@@ -463,31 +511,14 @@
     }
 
-    int numKernels = kernels->u->n;     // Number of kernel basis functions
-    int size = kernels->size;           // Size of kernel
-
-    // We could simply sum the kernel basis functions, but convolving a delta function gives a more 'real'
-    // picture of what's actually happening in the convolution, and the speed difference can't be too large,
-    // since the image is small.
-
-    psImage *image = psImageAlloc(2 * size + 1, 2 * size + 1, PS_TYPE_F32); // Output image
-    psImageInit(image, 0.0);
-    psImage *reference = psImageAlloc(4 * size + 1, 4 * size + 1, PS_TYPE_F32); // Reference image
-    psImageInit(reference, 0.0);
-    reference->data.F32[2 * size][2 * size] = 1.0;
-
     // Precalulate polynomial values
     psImage *polyValues = spatialPolyValues(kernels->spatialOrder, x, y);
 
-    for (int j = -size; j <= size; j++) {
-        for (int i = -size; i <= size; i++) {
-            for (int k = 0; k < numKernels; k++) {
-                image->data.F32[j + size][i + size] += solution->data.F64[k] *
-                    convolvePixel(kernels, k, 2 * size + i, 2 * size + j,
-                                  reference, polyValues, imageWeighting);
-            }
-        }
-    }
+    // The appropriate kernel
+    psKernel *kernel = solvedKernel(NULL, solution, kernels, polyValues, imageWeighting);
 
     psFree(polyValues);
+
+    psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel
+    psFree(kernel);
 
     return image;
@@ -541,5 +572,4 @@
     }
 
-    int numKernels = kernels->u->n;     // Number of kernel basis functions
     int size = kernels->size;           // Half-size of kernel
     float background = solution->data.F64[solution->n-1]; // The difference in background
@@ -570,4 +600,7 @@
     psImage *polyValues = NULL;         // Pre-calculated polynomial values
     int fullSize = 2 * size + 1;        // Full size of kernel
+
+    psKernel *kernelImage = NULL;       // Kernel for the image
+    psKernel *kernelWeight = NULL;      // Kernel for the weight map
 
     for (int y = size; y < inImage->numRows - size; y++) {
@@ -582,4 +615,9 @@
                                                2.0 * (float)(x - numCols/2.0) / (float)numCols,
                                                2.0 * (float)(y - numRows/2.0) / (float)numRows);
+                kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
+                if (inWeight) {
+                    kernelWeight = solvedKernel(kernelWeight, solution, kernels, polyValues,
+                                                varianceWeighting);
+                }
             }
 
@@ -602,18 +640,23 @@
 
             // Convolve the image
-            for (int i = 0; i < numKernels; i++) {
-                convImage->data.F32[y][x] += solution->data.F64[i] *
-                    convolvePixel(kernels, i, x, y, inImage, polyValues, imageWeighting);
+            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 i = 0; i < numKernels; i++) {
-                    convWeight->data.F32[y][x] += PS_SQR(solution->data.F64[i]) *
-                        convolvePixel(kernels, i, x, y, inWeight, polyValues, varianceWeighting);
-                }
-            }
-        }
-    }
+                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];
+                    }
+                }
+            }
+        }
+    }
+    psFree(kernelImage);
+    psFree(kernelWeight);
     psFree(polyValues);