Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 13389)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 13390)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-05-16 00:01:03 $
+ *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-05-16 01:04:43 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -58,5 +58,5 @@
                               const pmSubtractionKernels *kernels, // Kernel basis functions
                               const psImage *polyValues, // Spatial polynomial values
-                              float (*weightFunc)(float value) // Function for weighting
+                              double (*weightFunc)(double value) // Function for weighting
     )
 {
@@ -88,7 +88,7 @@
               int v = kernels->v->data.S32[i]; // Offset in y
               kernel->kernel[v][u] += weightFunc(polyValue * solution->data.F64[i]);
-              if (i != kernels->subIndex) {
+              if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
                   // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-                  kernel->kernel[0][0] -= weightFunc(solution->data.F64[i]);
+                  kernel->kernel[0][0] += weightFunc(-solution->data.F64[i]);
               }
               break;
@@ -96,7 +96,14 @@
           case PM_SUBTRACTION_KERNEL_ISIS: {
               psKernel *preCalc = kernels->preCalc->data[i];// Precalculated values
+              psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
               for (int v = -size; v <= size; v++) {
                   for (int u = -size; u <= size; u++) {
-                      kernel->kernel[v][u] += weightFunc(polyValue * preCalc->kernel[v][u]);
+                      kernel->kernel[v][u] += weightFunc(solution->data.F64[i] * polyValue *
+                                                         preCalc->kernel[v][u]);
+                      // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
+                      if (kernels->spatialOrder > 0 && i != kernels->subIndex) {
+                          kernel->kernel[v][u] += weightFunc(-solution->data.F64[i] *
+                                                             subKernel->kernel[v][u]);
+                      }
                   }
               }
@@ -117,5 +124,5 @@
                                    const psImage *image, // Image to convolve
                                    const psImage *polyValues, // Spatial polynomial values
-                                   float (*weightFunc)(float value) // Function for weighting
+                                   double (*weightFunc)(double value) // Function for weighting
     )
 {
@@ -130,7 +137,7 @@
           int v = kernels->v->data.S32[index]; // Offset in y
           double value = weightFunc(polyValue) * image->data.F32[y + v][x + u]; // Value of convolution
-          if (index != kernels->subIndex) {
+          if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
               // The (0,0) element is subtracted from most kernels to preserve photometric scaling
-              value -= image->data.F32[y][x];
+              value += weightFunc(-1.0) * image->data.F32[y][x];
           }
           return value;
@@ -138,12 +145,18 @@
       case PM_SUBTRACTION_KERNEL_ISIS: {
           psKernel *kernel = kernels->preCalc->data[index]; // The convolution kernel
+          psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
           int size = kernels->size;     // Kernel half-size
           double sum = 0.0;             // Accumulated sum from convolution
+          double sub = 0.0;             // Accumulated sum to subtract
           for (int v = -size; v <= size; v++) {
               for (int u = -size; u <= size; u++) {
                   sum += weightFunc(kernel->kernel[v][u]) * image->data.F32[y + v][x + u];
+                  // The (0,0) kernel is subtracted from other kernels to preserve photometric scaling
+                  if (kernels->spatialOrder > 0 && index != kernels->subIndex) {
+                      sub += weightFunc(subKernel->kernel[v][u]) * image->data.F32[y + v][x + u];
+                  }
               }
           }
-          return weightFunc(polyValue) * sum;
+          return weightFunc(polyValue) * sum + weightFunc(-1.0) * sub;
       }
       default:
@@ -154,5 +167,5 @@
 
 // Weighting function for use with convolvePixel: no weighting applied, suitable for combining image pixels
-static inline float imageWeighting(float value)
+static inline double imageWeighting(double value)
 {
     return value;
@@ -160,5 +173,5 @@
 
 // Weighting function for use with convolvePixel: weighting suitable for combining variances
-static inline float varianceWeighting(float value)
+static inline double varianceWeighting(double value)
 {
     return PS_SQR(value);
Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 13389)
+++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 13390)
@@ -137,4 +137,12 @@
                 }
 
+#if 0
+                // Subtract a particular kernel in order to preserve photometric calibration across image
+                if (spatialOrder > 0 && index != kernels->subIndex) {
+                    psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
+                    (void)psBinaryOp(preCalc->image, preCalc->image, "-", subKernel->image);
+                }
+#endif
+
                 // Iterate over spatial order.  This loop creates the terms for
                 // x^xOrder * y^yOrder  such that (xOrder+yOrder) <= spatialOrder.
@@ -153,14 +161,4 @@
                 }
 
-
-                if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
-                    psString kernelName = NULL;
-                    psStringAppend(&kernelName, "kernel%d.fits", index);
-                    psFits *kernelFile = psFitsOpen(kernelName, "w");
-                    psFree(kernelName);
-                    psFitsWriteImage(kernelFile, NULL, preCalc->image, 0, NULL);
-                    psFitsClose(kernelFile);
-                }
-
                 psFree(preCalc);        // Drop reference
             }
@@ -168,5 +166,4 @@
     }
 
-    // Subtract a particular kernel in order to preserve photometric calibration across image
     kernels->subIndex = 0;
     assert(kernels->u->data.S32[kernels->subIndex] == 0 &&
@@ -175,14 +172,15 @@
            kernels->yOrder->data.S32[kernels->subIndex] == 0);
 
-    psKernel *subKernel = kernels->preCalc->data[kernels->subIndex]; // Kernel to subtract
-    for (int i = 0; i < kernels->subIndex; i++) {
-        psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-        (void)psBinaryOp(kernel->image, kernel->image, "-", subKernel->image);
-    }
-    for (int i = kernels->subIndex + 1; i < num; i++) {
-        psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
-        (void)psBinaryOp(kernel->image, kernel->image, "-", subKernel->image);
-    }
-
+    if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) {
+        for (int i = 0; i < num; i++) {
+            psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest
+            psString kernelName = NULL;
+            psStringAppend(&kernelName, "kernel%03d.fits", i);
+            psFits *kernelFile = psFitsOpen(kernelName, "w");
+            psFree(kernelName);
+            psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL);
+            psFitsClose(kernelFile);
+        }
+    }
 
     return kernels;