Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 19764)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 19765)
@@ -30,4 +30,5 @@
 #define MIN_SAMPLE_STATS    7           // Minimum number to use sample statistics; otherwise use quartiles
 
+#define SYS_ERROR 0.05                  // Relative error in derived convolution kernel pixels
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -135,11 +136,16 @@
           }
           case PM_SUBTRACTION_KERNEL_ISIS: {
-              psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values
+              psArray *preCalc = kernels->preCalc->data[i]; // Precalculated values
+              psVector *xKernel = preCalc->data[0]; // Kernel in x
+              psVector *yKernel = preCalc->data[1]; // Kernel in y
               // Iterating over the kernel
-              for (int v = -size; v <= size; v++) {
-                  for (int u = -size; u <= size; u++) {
-                      kernel->kernel[v][u] += preCalc->kernel[v][u] * value;
-                      // Photometric scaling is built into the preCalc kernel --- no subtraction!
+              for (int y = 0, v = -size; v <= size; y++, v++) {
+                  for (int x = 0, u = -size; u <= size; x++, u++) {
+                      kernel->kernel[v][u] +=  value * xKernel->data.F32[x] * yKernel->data.F32[y];
                   }
+              }
+              // Photometric scaling for even kernels only
+              if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) {
+                  kernel->kernel[0][0] -= value;
               }
               break;
@@ -200,4 +206,41 @@
 }
 
+// Take a subset of an image
+static inline psImage *convolveSubsetAlloc(psImage *image, // Image to be convolved
+                                           psRegion region, // Region of interest (with border)
+                                           bool threaded // Are we running threaded?
+                                           )
+{
+    if (!image) {
+        return NULL;
+    }
+    if (threaded) {
+        psMutexLock(image);
+    }
+    psImage *subset = psImageSubset(image, region); // Subset image, to return
+    if (threaded) {
+        psMutexUnlock(image);
+    }
+    return subset;
+}
+
+// Free the subset of an image
+static inline void convolveSubsetFree(psImage *parent, // Parent image
+                                      psImage *child, // Child (subset) image
+                                      bool threaded // Are we running threaded?
+                                      )
+{
+    if (!child || !parent) {
+        return;
+    }
+    if (threaded) {
+        psMutexLock(parent);
+    }
+    psFree(child);
+    if (threaded) {
+        psMutexUnlock(parent);
+    }
+    return;
+}
 
 // Convolve an image using FFT
@@ -217,44 +260,11 @@
     bool threaded = pmSubtractionThreaded(); // Are we running threaded?
 
-    if (threaded) {
-        psMutexLock(image);
-    }
-    psImage *subImage = psImageSubset(image, border); // Subimage to convolve
-    if (threaded) {
-        psMutexUnlock(image);
-    }
-
-    psImage *subMask;                   // Subimage mask
-    if (mask) {
-        if (threaded) {
-            psMutexLock(mask);
-        }
-        subMask = psImageSubset(mask, border);
-        if (threaded) {
-            psMutexUnlock(mask);
-        }
-    } else {
-        subMask = NULL;
-    }
+    psImage *subImage = convolveSubsetAlloc(image, border, threaded); // Subimage to convolve
+    psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Subimage mask
 
     psImage *convolved = psImageConvolveFFT(NULL, subImage, subMask, maskVal, kernel); // Convolution
 
-    if (threaded) {
-        psMutexLock(image);
-    }
-    psFree(subImage);
-    if (threaded) {
-        psMutexUnlock(image);
-    }
-
-    if (mask) {
-        if (threaded) {
-            psMutexLock(mask);
-        }
-        psFree(subMask);
-        if (threaded) {
-            psMutexUnlock(mask);
-        }
-    }
+    convolveSubsetFree(image, subImage, threaded);
+    convolveSubsetFree(mask, subMask, threaded);
 
     // Now, we have to stick it in where it belongs
@@ -276,4 +286,56 @@
     return;
 }
+
+
+// Convolve an image using FFT
+static void convolveWeightFFT(psImage *target,// Place the result in here
+                              psImage *weight, // Weight map to convolve
+                              psImage *sys, // Systematic error image
+                              psImage *mask, // Mask image
+                              psMaskType maskVal, // Value to mask
+                              const psKernel *kernel, // Kernel by which to convolve
+                              psRegion region,// Region of interest
+                              int size        // Size of (square) kernel
+                              )
+{
+    psRegion border = psRegionSet(region.x0 - size, region.x1 + size,
+                                  region.y0 - size, region.y1 + size); // Add a border
+
+    bool threaded = pmSubtractionThreaded(); // Are we running threaded?
+
+    psImage *subWeight = convolveSubsetAlloc(weight, border, threaded); // Weight map
+    psImage *subSys = convolveSubsetAlloc(sys, border, threaded); // Systematic error image
+    psImage *subMask = convolveSubsetAlloc(mask, border, threaded); // Mask
+
+    // XXX Can trim this a little by combining the convolution: only have to take the FFT of the kernel once
+    psImage *convWeight = psImageConvolveFFT(NULL, subWeight, subMask, maskVal, kernel); // Convolved weight
+    psImage *convSys = subSys ? psImageConvolveFFT(NULL, subSys, subMask, maskVal, kernel) : NULL; // Conv sys
+
+    convolveSubsetFree(weight, subWeight, threaded);
+    convolveSubsetFree(sys, subSys, threaded);
+    convolveSubsetFree(mask, subMask, threaded);
+
+    // Now, we have to stick it in where it belongs
+    int xMin = region.x0, xMax = region.x1, yMin = region.y0, yMax = region.y1; // Bounds of region
+    if (convSys) {
+        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
+            for (int xTarget = xMin, xSource = size; xTarget < xMax; xTarget++, xSource++) {
+                target->data.F32[yTarget][xTarget] = convWeight->data.F32[ySource][xSource] +
+                    convSys->data.F32[ySource][xSource];
+            }
+        }
+    } else {
+        int numBytes = (xMax - xMin) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy
+        for (int yTarget = yMin, ySource = size; yTarget < yMax; yTarget++, ySource++) {
+            memcpy(&target->data.F32[yTarget][xMin], &convWeight->data.F32[ySource][size], numBytes);
+        }
+    }
+
+    psFree(convWeight);
+    psFree(convSys);
+
+    return;
+}
+
 
 // Convolve an image directly
@@ -307,4 +369,5 @@
                                   psImage *image, // Image to convolve
                                   psImage *weight, // Weight map to convolve, or NULL
+                                  psImage *sys, // Systematic error image, or NULL
                                   psImage *subMask, // Subtraction mask
                                   const pmSubtractionKernels *kernels, // Kernels
@@ -343,5 +406,5 @@
         convolveFFT(convImage, image, subMask, subBad, *kernelImage, region, background, kernels->size);
         if (weight) {
-            convolveFFT(convWeight, weight, subMask, subBad, *kernelWeight, region, 0.0, kernels->size);
+            convolveWeightFFT(convWeight, weight, sys, subMask, subBad, *kernelWeight, region, kernels->size);
         }
     } else {
@@ -361,23 +424,13 @@
             bool threaded = pmSubtractionThreaded(); // Are we running threaded?
 
-            if (threaded) {
-                psMutexLock(subMask);
-            }
-            psImage *image = psImageSubset(subMask, psRegionSet(colMin - box, colMax + box,
-                                                                rowMin - box, rowMax + box)); // Mask
-            if (threaded) {
-                psMutexUnlock(subMask);
-            }
+            psRegion region = psRegionSet(colMin - box, colMax + box,
+                                          rowMin - box, rowMax + box); // Region to convolve
+
+            psImage *image = convolveSubsetAlloc(subMask, region, threaded); // Mask to convolve
 
             psImage *convolved = psImageConvolveMask(NULL, image, subBad, subConvBad,
                                                      -box, box, -box, box); // Convolved subtraction mask
 
-            if (threaded) {
-                psMutexLock(subMask);
-            }
-            psFree(image);
-            if (threaded) {
-                psMutexUnlock(subMask);
-            }
+            convolveSubsetFree(subMask, image, threaded);
 
             psAssert(convolved->numCols - 2 * box == colMax - colMin, "Bad number of columns");
@@ -405,5 +458,19 @@
 }
 
-
+// Generate an image that can be used to track systematic errors
+static psImage *subtractionSysErrImage(const psImage *image     // Image from which to make sys err image
+                                       )
+{
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+    psImage *sys = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Systematic error image
+
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            sys->data.F32[y][x] = PS_SQR(image->data.F32[y][x]) * PS_SQR(SYS_ERROR);
+        }
+    }
+
+    return sys;
+}
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -575,6 +642,52 @@
       }
       case PM_SUBTRACTION_KERNEL_ISIS: {
+          psArray *preCalc = kernels->preCalc->data[index]; // Precalculated values
+          psVector *xKernel = preCalc->data[0]; // Kernel in x
+          psVector *yKernel = preCalc->data[1]; // Kernel in y
+          int size = kernels->size;     // Size of kernel
+
+          // Convolve in x
+          // Need to convolve a bit more than the footprint, for the y convolution
+          int yMin = -size - footprint, yMax = size + footprint; // Range for y
+          psKernel *temp = psKernelAlloc(yMin, yMax,
+                                         -footprint, footprint); // Temporary convolution; NOTE: wrong way!
+          for (int y = yMin; y <= yMax; y++) {
+              for (int x = -footprint; x <= footprint; x++) {
+                  float value = 0.0;    // Value of convolved pixel
+                  int uMin = x - size, uMax = x + size; // Range for u
+                  psF32 *xKernelData = xKernel->data.F32; // Kernel values
+                  psF32 *imageData = &image->kernel[y][uMin]; // Image values
+                  for (int u = uMin; u <= uMax; u++, xKernelData++, imageData++) {
+                      value += *xKernelData * *imageData;
+                  }
+                  temp->kernel[x][y] = value; // NOTE: putting in wrong way!
+              }
+          }
+
+          // Convolve in y
+          psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint);// Convolved image
+          for (int x = -footprint; x <= footprint; x++) {
+              for (int y = -footprint; y <= footprint; y++) {
+                  float value = 0.0;    // Value of convolved pixel
+                  int vMin = y - size, vMax = y + size; // Range for v
+                  psF32 *yKernelData = yKernel->data.F32; // Kernel values
+                  psF32 *imageData = &temp->kernel[x][vMin]; // Image values; NOTE: wrong way!
+                  for (int v = vMin; v <= vMax; v++, yKernelData++, imageData++) {
+                      value += *yKernelData * *imageData;
+                  }
+                  convolved->kernel[y][x] = value;
+              }
+          }
+          psFree(temp);
+
+          // Photometric scaling for even kernels only
+          if (kernels->u->data.S32[index] % 2 == 0 && kernels->v->data.S32[index] % 2 == 0) {
+              convolveSub(convolved, image, footprint);
+          }
+          return convolved;
+#if 0
           // Photometric scaling is already built in to the precalculated kernel
           return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]);
+#endif
       }
       case PM_SUBTRACTION_KERNEL_RINGS: {
@@ -889,4 +1002,5 @@
                                      psImage *convMask, // Output convolved mask
                                      const pmReadout *ro1, const pmReadout *ro2, // Input readouts
+                                     psImage *sys1, psImage *sys2, // Systematic error images
                                      psImage *subMask, // Input subtraction mask
                                      psMaskType maskBad, // Mask value to give bad pixels
@@ -914,10 +1028,10 @@
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out1->image, out1->weight, convMask, &kernelImage, &kernelWeight,
-                       ro1->image, ro1->weight, subMask, kernels, polyValues, background, *region,
+                       ro1->image, ro1->weight, sys1, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, false);
     }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
         convolveRegion(out2->image, out2->weight, convMask, &kernelImage, &kernelWeight,
-                       ro2->image, ro2->weight, subMask, kernels, polyValues, background, *region,
+                       ro2->image, ro2->weight, sys2, subMask, kernels, polyValues, background, *region,
                        maskBad, maskPoor, poorFrac, useFFT, kernels->mode == PM_SUBTRACTION_MODE_DUAL);
     }
@@ -965,15 +1079,17 @@
     const pmReadout *ro1 = args->data[7]; // Input readout 1
     const pmReadout *ro2 = args->data[8]; // Input readout 2
-    psImage *subMask = args->data[9]; // Subtraction mask
-    psMaskType maskBad = PS_SCALAR_VALUE(args->data[10], U8); // Output mask value for bad pixels
-    psMaskType maskPoor = PS_SCALAR_VALUE(args->data[11], U8); // Output mask value for poor pixels
-    float poorFrac = PS_SCALAR_VALUE(args->data[12], F32); // Fraction for "poor"
-    const psRegion *region = args->data[13]; // Region to convolve
-    const pmSubtractionKernels *kernels = args->data[14]; // Kernels
-    bool doBG = PS_SCALAR_VALUE(args->data[15], U8); // Do background subtraction?
-    bool useFFT = PS_SCALAR_VALUE(args->data[16], U8); // Use FFT for convolution?
-
-    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
-                                    maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
+    psImage *sys1 = args->data[9]; // Systematic error image 1
+    psImage *sys2 = args->data[10]; // Systematic error image 2
+    psImage *subMask = args->data[11]; // Subtraction mask
+    psMaskType maskBad = PS_SCALAR_VALUE(args->data[12], U8); // Output mask value for bad pixels
+    psMaskType maskPoor = PS_SCALAR_VALUE(args->data[13], U8); // Output mask value for poor pixels
+    float poorFrac = PS_SCALAR_VALUE(args->data[14], F32); // Fraction for "poor"
+    const psRegion *region = args->data[15]; // Region to convolve
+    const pmSubtractionKernels *kernels = args->data[16]; // Kernels
+    bool doBG = PS_SCALAR_VALUE(args->data[17], U8); // Do background subtraction?
+    bool useFFT = PS_SCALAR_VALUE(args->data[18], U8); // Use FFT for convolution?
+
+    return subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, sys1, sys2,
+                                    subMask, maskBad, maskPoor, poorFrac, region, kernels, doBG, useFFT);
 }
 
@@ -1022,29 +1138,15 @@
     }
 
-    // Inputs
-    psImage *image1 = NULL, *weight1 = NULL; // Image and weight map from input 1
-    if (ro1) {
-        image1 = ro1->image;
-        weight1 = ro1->weight;
-    }
-    psImage *image2 = NULL, *weight2 = NULL; // Image and weight map from input 2
-    if (ro2) {
-        image2 = ro2->image;
-        weight2 = ro2->weight;
-    }
-
     bool threaded = pmSubtractionThreaded(); // Running threaded?
 
     // Outputs
-    psImage *convImage1 = NULL, *convWeight1 = NULL; // Convolved image and weight from input 1
     if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convImage1 = out1->image;
-        if (!convImage1) {
-            convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        if (!out1->image) {
+            out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
             if (threaded) {
-                psMutexInit(convImage1);
+                psMutexInit(out1->image);
             }
         }
-        if (weight1) {
+        if (ro1->weight) {
             if (!out1->weight) {
                 out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -1053,18 +1155,15 @@
                 }
             }
-            convWeight1 = out1->weight;
-            psImageInit(convWeight1, 0.0);
-        }
-    }
-    psImage *convImage2 = NULL, *convWeight2 = NULL; // Convolved image and weight from input 2
+            psImageInit(out1->weight, 0.0);
+        }
+    }
     if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        convImage2 = out2->image;
-        if (!convImage2) {
-            convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        if (!out2->image) {
+            out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
             if (threaded) {
-                psMutexInit(convImage2);
+                psMutexInit(out2->image);
             }
         }
-        if (weight2) {
+        if (ro2->weight) {
             if (!out2->weight) {
                 out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32);
@@ -1073,6 +1172,5 @@
                 }
             }
-            convWeight2 = out2->weight;
-            psImageInit(convWeight2, 0.0);
+            psImageInit(out2->weight, 0.0);
         }
     }
@@ -1104,4 +1202,17 @@
     }
 
+    psImage *sys1 = NULL, *sys2 = NULL; // Systematic error images
+    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        sys1 = subtractionSysErrImage(ro1->image);
+        if (threaded) {
+            psMutexInit(sys1);
+        }
+    }
+    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
+        sys2 = subtractionSysErrImage(ro2->image);
+        if (threaded) {
+            psMutexInit(sys2);
+        }
+    }
 
     int size = kernels->size;           // Half-size of kernel
@@ -1144,4 +1255,18 @@
                 psArrayAdd(args, 1, (pmReadout*)ro2); // Casting away const
                 // Since adding to the array can impact the reference count, we need to lock
+                if (sys1) {
+                    psMutexLock(sys1);
+                }
+                psArrayAdd(args, 1, sys1);
+                if (sys1) {
+                    psMutexUnlock(sys1);
+                }
+                if (sys2) {
+                    psMutexLock(sys2);
+                }
+                psArrayAdd(args, 1, sys2);
+                if (sys2) {
+                    psMutexUnlock(sys2);
+                }
                 if (subMask) {
                     psMutexLock(subMask);
@@ -1165,6 +1290,7 @@
                 psFree(job);
             } else {
-                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2, subMask,
-                                         maskBad, maskPoor, poorFrac, subRegion, kernels, doBG, useFFT);
+                subtractionConvolvePatch(numCols, numRows, x0, y0, out1, out2, convMask, ro1, ro2,
+                                         sys1, sys2, subMask, maskBad, maskPoor, poorFrac, subRegion,
+                                         kernels, doBG, useFFT);
             }
             psFree(subRegion);
@@ -1190,4 +1316,10 @@
             psMutexDestroy(subMask);
         }
+        if (sys1) {
+            psMutexDestroy(sys1);
+        }
+        if (sys2) {
+            psMutexDestroy(sys2);
+        }
     }
 
Index: /trunk/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 19764)
+++ /trunk/psModules/src/imcombine/pmSubtractionEquation.c	(revision 19765)
@@ -758,4 +758,25 @@
         calculatePenalty(sumVector, kernels);
 
+#ifdef TESTING
+        {
+            psImage *inverse = psMatrixInvert(NULL, sumMatrix, NULL);
+            psFits *fits = psFitsOpen("matrixInv.fits", "w");
+            psFitsWriteImage(fits, NULL, inverse, 0, NULL);
+            psFitsClose(fits);
+            psFree(inverse);
+        }
+        {
+            psImage *X = psMatrixInvert(NULL, sumMatrix, NULL);
+            psImage *Xt = psMatrixTranspose(NULL, X);
+            psImage *XtX = psMatrixMultiply(NULL, Xt, X);
+            psFits *fits = psFitsOpen("matrixErr.fits", "w");
+            psFitsWriteImage(fits, NULL, XtX, 0, NULL);
+            psFitsClose(fits);
+            psFree(X);
+            psFree(Xt);
+            psFree(XtX);
+        }
+#endif
+
         psVector *permutation = NULL;       // Permutation vector, required for LU decomposition
         psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix);
@@ -769,4 +790,5 @@
         }
         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
+
         psFree(sumVector);
         psFree(luMatrix);
@@ -806,13 +828,4 @@
         calculatePenalty(sumVector2, kernels);
 
-#if 0
-        // Apply weighting to maximise the difference between solution coefficients for the same functions
-        double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor
-        for (int i = 0; i < kernels->num; i++) {
-            sumMatrix1->data.F64[i][i] -= fudge;
-            sumMatrix2->data.F64[i][i] += fudge;
-        }
-#endif
-
         // Pure matrix operations
 
@@ -941,4 +954,5 @@
 #endif
 
+#ifdef TESTING
         {
             psFits *fits = psFitsOpen("sumMatrix1.fits", "w");
@@ -961,5 +975,5 @@
             psFitsClose(fits);
         }
-
+#endif
 
         kernels->solution1 = a;
Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 19764)
+++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 19765)
@@ -42,4 +42,22 @@
     }
     return result;
+}
+
+// Generate 1D convolution kernel for ISIS
+static psVector *subtractionKernelISIS(float sigma, // Gaussian width
+                                       int order, // Polynomial order
+                                       int size // Kernel half-size
+    )
+{
+    int fullSize = 2 * size + 1;        // Full size of kernel
+    psVector *kernel = psVectorAlloc(fullSize, PS_TYPE_F32); // Kernel to return
+
+    float expNorm = -0.5 / PS_SQR(sigma); // Normalisation for exponential
+    float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
+    for (int i = 0, x = -size; x <= size; i++, x++) {
+        kernel->data.F32[i] = norm * power(x, order) * expf(expNorm * PS_SQR(x));
+    }
+
+    return kernel;
 }
 
@@ -116,30 +134,38 @@
 
     // Set the kernel parameters
+    int fullSize = 2 * size + 1;        // Full size of kernels
     for (int i = 0, index = 0; i < numGaussians; i++) {
         float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma
-        float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian
         // Iterate over (u,v) order
         for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) {
             for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) {
-                // Set the pre-calculated kernel
-                psKernel *preCalc = psKernelAlloc(-size, size, -size, size);
-                double sum = 0.0;       // Normalisation
+                psArray *preCalc = psArrayAlloc(2); // Array to hold precalculated values
+                psVector *xKernel = preCalc->data[0] = subtractionKernelISIS(sigma, uOrder, size); // x Kernel
+                psVector *yKernel = preCalc->data[1] = subtractionKernelISIS(sigma, vOrder, size); // y Kernel
+
+                // Calculate moments
+                double sum = 0.0;       // Sum of kernel component, for normalisation
                 double moment = 0.0;    // Moment, for penalty
-                for (int v = -size; v <= size; v++) {
-                    for (int u = -size; u <= size; u++) {
-                        sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) *
-                            expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma));
-                        moment += preCalc->kernel[v][u] * (PS_SQR(u) + PS_SQR(v));
+                for (int v = -size, y = 0; v <= size; v++, y++) {
+                    for (int u = -size, x = 0; u <= size; u++, x++) {
+                        double value = xKernel->data.F32[x] * yKernel->data.F32[y]; // Value of kernel
+                        sum += value;
+                        moment += value * (PS_SQR(u) + PS_SQR(v));
                     }
                 }
+                moment = 0.0;
 
                 // Normalise sum of kernel component to unity for even functions
                 if (uOrder % 2 == 0 && vOrder % 2 == 0) {
-                    for (int v = -size; v <= size; v++) {
-                        for (int u = -size; u <= size; u++) {
-                            preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum;
+                    double sum = 0.0;   // Sum of kernel component
+                    for (int v = 0; v < fullSize; v++) {
+                        for (int u = 0; u < fullSize; u++) {
+                            sum += xKernel->data.F32[u] * yKernel->data.F32[v];
                         }
                     }
-                    preCalc->kernel[0][0] -= 1.0;
+                    sum = 1.0 / sum;
+                    psBinaryOp(xKernel, xKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    psBinaryOp(yKernel, yKernel, "*", psScalarAlloc(sum, PS_TYPE_F32));
+                    moment *= sum;
                 }
 
@@ -702,5 +728,5 @@
             // Count the number of Gaussians
             int numGauss = 0;
-            for (char *string = ptr; string; string = strchr(string, '(')) {
+            for (char *string = ptr; string; string = strchr(string + 1, '(')) {
                 numGauss++;
             }
Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionThreads.c	(revision 19764)
+++ /trunk/psModules/src/imcombine/pmSubtractionThreads.c	(revision 19765)
@@ -76,5 +76,5 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE", 17);
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE", 19);
         task->function = &pmSubtractionConvolveThread;
         psThreadTaskAdd(task);