Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14312)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14313)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-19 02:03:42 $
+ *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-19 02:33:33 $
  *
  *  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?
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -125,4 +125,25 @@
             } \
             break; \
+          } \
+          case PM_SUBTRACTION_KERNEL_GUNK: { \
+              if (i < (KERNELS)->inner) { \
+                  /* Using pre-calculated function */ \
+                  psKernel *preCalc = (KERNELS)->preCalc->data[i]; /* Precalculated values */ \
+                  for (int v = -size; v <= size; v++) { \
+                      for (int u = -size; u <= size; u++) { \
+                          (TARGET)->kernel[v][u] += value * FUNC(preCalc->kernel[v][u]); \
+                      } \
+                  } \
+              } else { \
+                  /* Using delta function */ \
+                  int u = (KERNELS)->u->data.S32[i]; /* Offset in x */ \
+                  int v = (KERNELS)->v->data.S32[i]; /* Offset in y */ \
+                  (TARGET)->kernel[v][u] += value; \
+              } \
+              /* The (0,0) kernel is subtracted from other kernels to preserve photometric scaling */ \
+              if ((KERNELS)->spatialOrder > 0 && i != (KERNELS)->subIndex) { \
+                  (TARGET)->kernel[0][0] += subValue; \
+              } \
+              break; \
           } \
           case PM_SUBTRACTION_KERNEL_ISIS: { \
@@ -309,4 +330,29 @@
           break; \
       } \
+      case PM_SUBTRACTION_KERNEL_GUNK: { \
+          double value;                 /* The value to return */ \
+          if ((INDEX) < (KERNELS)->inner) { \
+              /* Using pre-calculated function */ \
+              psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
+              int size = (KERNELS)->size;     /* Kernel half-size */ \
+              double sum = 0.0;             /* Accumulated sum from convolution */ \
+              for (int v = -size; v <= size; v++) { \
+                  for (int u = -size; u <= size; u++) { \
+                      sum += FUNC(kernel->kernel[v][u]) * (IMAGE)->data.F32[(Y) + v][(X) + u]; \
+                  } \
+              } \
+              value = FUNC(polyValue) * sum; \
+          } else { \
+              /* Using delta function */ \
+              int u = (KERNELS)->u->data.S32[(INDEX)]; /* Offset in x */ \
+              int v = (KERNELS)->v->data.S32[(INDEX)]; /* Offset in y */ \
+              value = FUNC(polyValue) * (IMAGE)->data.F32[(Y) + v][(X) + u]; /* Value of convolution */ \
+          } \
+          /* The (0,0) delta function is subtracted from most kernels to preserve photometric scaling */ \
+          if ((KERNELS)->spatialOrder > 0 && (INDEX) != (KERNELS)->subIndex) { \
+              value += FUNC(-1.0) * (IMAGE)->data.F32[(Y)][(X)]; \
+          } \
+          (TARGET) = value; \
+      } \
       case PM_SUBTRACTION_KERNEL_ISIS: { \
           psKernel *kernel = (KERNELS)->preCalc->data[(INDEX)]; /* The convolution kernel */ \
@@ -383,5 +429,4 @@
               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++) {
@@ -389,5 +434,5 @@
                   }
               }
-              value = weightFunc(polyValue) * sum + weightFunc(-1.0) * sub;
+              value = weightFunc(polyValue) * sum;
           } else {
               // Using delta function