Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14317)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14318)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-19 03:36:06 $
+ *  @version $Revision: 1.28 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-19 21:26:12 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -25,5 +25,5 @@
 #include "pmSubtraction.h"
 
-//#define USE_FUNCTIONS_INSTEAD_OF_MACROS          // Can I pass the address of a static inline function?
+#define USE_FUNCTIONS_INSTEAD_OF_MACROS          // Can I pass the address of a static inline function?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -851,7 +851,5 @@
     }
 
-    int numKernels = solution->n - 1;   // Number of kernel basis functions
-
-    // Measure
+    // Measure deviations
     psVector *deviations = psVectorAlloc(stamps->n, PS_TYPE_F32); // Mean deviation for stamps
     double totalSquareDev = 0.0;        // Total square deviation from zero
@@ -859,9 +857,10 @@
     {
         int numCols = refImage->numCols, numRows = refImage->numRows; // Image dimensions
-
         psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // Statistics
-
-        // Convolved image of stamp
-        psImage *convolvedStamp = psImageAlloc(2 * footprint + 1, 2 * footprint + 1, PS_TYPE_F32);
+        psImage *convolvedStamp = psImageAlloc(2 * footprint + 1, 2 * footprint + 1,
+                                               PS_TYPE_F32); // Convolved image of stamp
+        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++) {
@@ -871,5 +870,4 @@
             }
             int xStamp = stamp->x, yStamp = stamp->y; // Coordinates of stamp
-            psImageInit(convolvedStamp, 0.0);
 
             // Precalulate polynomial values
@@ -878,18 +876,21 @@
                                                     2.0 * (float)(stamp->y - numRows/2.0) / (float)numRows);
 
-            for (int y = - footprint; y <= footprint; y++) {
-                for (int x = footprint; x <= footprint; x++) {
-                    for (int k = 0; k < numKernels; k++) {
 #ifdef USE_FUNCTIONS_INSTEAD_OF_MACROS
-                        convolvedStamp->data.F32[y + footprint][x + footprint] +=
-                            convolvePixel(kernels, k, xStamp + x, yStamp + y, refImage, polyValues,
-                                          imageWeighting);
+            kernelImage = solvedKernel(kernelImage, solution, kernels, polyValues, imageWeighting);
 #else
-                        CONVOLVE_PIXEL(convolvedStamp->data.F32[y + footprint][x + footprint], kernels,
-                                       k, xStamp + x, yStamp + y, refImage, polyValues, imageWeighting);
+            SOLVED_KERNEL(kernelImage, solution, kernels, polyValues, imageWeighting);
 #endif
+            for (int y = yStamp - footprint, j = 0; y <= yStamp + footprint; y++, j++) {
+                for (int x = xStamp - footprint, i = 0; x <= xStamp + footprint; x++, i++) {
+                    convolvedStamp->data.F32[j][i] = background;
+                    for (int v = -size; v <= size; v++) {
+                        for (int u = -size; u <= size; u++) {
+                            convolvedStamp->data.F32[j][i] += kernelImage->kernel[v][u] *
+                                refImage->data.F32[y + v][x + u];
+                        }
                     }
                 }
             }
+
             psFree(polyValues);
 
@@ -900,7 +901,7 @@
                    inStamp->numRows == convolvedStamp->numRows);
             (void)psBinaryOp(convolvedStamp, inStamp, "-", convolvedStamp);
+            (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp);
             (void)psBinaryOp(convolvedStamp, convolvedStamp, "/", inStamp);
             psFree(inStamp);
-            (void)psBinaryOp(convolvedStamp, convolvedStamp, "*", convolvedStamp);
             if (!psImageStats(stats, convolvedStamp, NULL, 0)) {
                 psError(PS_ERR_UNKNOWN, false,
@@ -911,14 +912,17 @@
                 return -1;
             }
-            deviations->data.F32[i] = stats->sampleMean;
+            deviations->data.F32[i] = sqrt(stats->sampleMean / 2.0);
+            psTrace("psModules.imcombine", 1, "Deviation for stamp %d: %f\n", i, deviations->data.F32[i]);
             totalSquareDev += PS_SQR(deviations->data.F32[i]);
             numStamps++;
         }
 
+        psFree(kernelImage);
         psFree(convolvedStamp);
         psFree(stats);
     }
 
-    float limit = sigmaRej * sqrt(totalSquareDev / (numStamps - 1)); // Limit on maximum deviation
+    float limit = sigmaRej * sqrt(totalSquareDev / numStamps); // Limit on maximum deviation
+    psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit);
 
     int numRejected = 0;