Changeset 19765

Timestamp:

Sep 26, 2008, 4:03:24 PM (18 years ago)

Author:

Paul Price

Message:

Two major changes. 1. Optimising ISIS kernels: instead of doing an FFT convolution, they're separable. 2. Trying to make the convolved weight map account for systematic errors in the convolution kernel by assuming a relative error in each kernel pixel (currently hardwired with a #define, but we can turn it into a recipe value later, or get an estimate from the data) and propagating that into the variance map. The effect is to increase the variance near the bright stars, which is hopefully sufficient to decrease the number of false sources detected in bad subtractions.

Location:

trunk/psModules/src/imcombine

Files:

• pmSubtraction.c (modified) (20 diffs)
• pmSubtractionEquation.c (modified) (5 diffs)
• pmSubtractionKernels.c (modified) (3 diffs)
• pmSubtractionThreads.c (modified) (1 diff)

trunk/psModules/src/imcombine/pmSubtraction.c

@@ -30 +30 @@
 #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 +136 @@
           }
           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 +206 @@
 }
 
+// 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 +260 @@
     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 +286 @@
     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 +369 @@
                                   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 +406 @@
         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 +424 @@
             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 +458 @@
 }
 
-
+// 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 +642 @@
       }
       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 +1002 @@
                                      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 +1028 @@
     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 +1079 @@
     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 +1138 @@
     }
 
-    // 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 +1155 @@
                 }
             }
-            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 +1172 @@
                 }
             }
-            convWeight2 = out2->weight;
-            psImageInit(convWeight2, 0.0);
+            psImageInit(out2->weight, 0.0);
         }
     }
@@ -1104 +1202 @@
     }
 
+    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 +1255 @@
                 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 +1290 @@
                 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 +1316 @@
             psMutexDestroy(subMask);
         }
+        if (sys1) {
+            psMutexDestroy(sys1);
+        }
+        if (sys2) {
+            psMutexDestroy(sys2);
+        }
     }
 

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -758 +758 @@
         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 +790 @@
         }
         kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation);
+
         psFree(sumVector);
         psFree(luMatrix);
@@ -806 +828 @@
         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 +954 @@
 #endif
 
+#ifdef TESTING
         {
             psFits *fits = psFitsOpen("sumMatrix1.fits", "w");
@@ -961 +975 @@
             psFitsClose(fits);
         }
-
+#endif
 
         kernels->solution1 = a;

trunk/psModules/src/imcombine/pmSubtractionKernels.c

@@ -42 +42 @@
     }
     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 +134 @@
 
     // 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 +728 @@
             // 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++;
             }

trunk/psModules/src/imcombine/pmSubtractionThreads.c

@@ -76 +76 @@
 
     {
-        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE", 17);
+        psThreadTask *task = psThreadTaskAlloc("PSMODULES_SUBTRACTION_CONVOLVE", 19);
         task->function = &pmSubtractionConvolveThread;
         psThreadTaskAdd(task);