Index: trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 19209)
+++ 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++;
             }