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);