Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 19271)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 19282)
@@ -300,5 +300,5 @@
                                   psImage *subMask, // Subtraction mask
                                   const pmSubtractionKernels *kernels, // Kernels
-                                  psImage *polyValues, // Polynomial values
+                                  const psImage *polyValues, // Polynomial values
                                   float background, // Background value to apply
                                   psRegion region, // Region to convolve
@@ -340,32 +340,6 @@
     // Convolve the mask for bad pixels
     if (subMask && convMask) {
-        psKernel *kernel = *kernelWeight; // Kernel of interest
-        int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
-
-        // Determine the thresholds
-        double sumKernel2 = 0.0;        // Sum of the kernel-squared
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                sumKernel2 += kernel->kernel[y][x];
-            }
-        }
-        float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad
-
-        // Get bounds of threshold region
-        // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold
-        int box = kernels->size;                    // Size of box with bad pixels
-        for (double sumBox = sumKernel2; box > 0; box--) {
-            for (int x = -box; x <= box; x++) {
-                sumBox -= kernel->kernel[-box][x] + kernel->kernel[box][x];
-            }
-            for (int y = -box + 1; y <= box - 1; y++) {
-                // Note: not doing corners
-                sumBox -= kernel->kernel[y][-box] + kernel->kernel[y][box];
-            }
-            if (sumBox < threshold) {
-                break;
-            }
-        }
-
+        int box = p_pmSubtractionBadRadius(*kernelImage, kernels, polyValues,
+                                           wantDual, poorFrac); // Size of bad box
         if (box > 0) {
             int colMin = region.x0, colMax = region.x1, rowMin = region.y0, rowMax = region.y1; // Bounds
@@ -415,4 +389,6 @@
     return;
 }
+
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -430,4 +406,65 @@
     psFree(conv);
     return convolved;
+}
+
+int p_pmSubtractionBadRadius(psKernel *preKernel, const pmSubtractionKernels *kernels,
+                             const psImage *polyValues, bool wantDual, float poorFrac)
+{
+    psKernel *kernel;                   // Kernel to use
+    if (!preKernel) {
+        kernel = solvedKernel(NULL, kernels, polyValues, wantDual);
+    } else {
+        kernel = psMemIncrRefCounter(preKernel);
+    }
+    PS_ASSERT_IMAGE_NON_NULL(polyValues, -1);
+
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+
+    // Determine the threshold between bad and poor
+    double sumKernel2 = 0.0;            // Sum of the kernel-squared
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sumKernel2 += PS_SQR(kernel->kernel[y][x]);
+        }
+    }
+    float threshold = sumKernel2 * poorFrac; // Threshold between poor and bad
+
+    // Get bounds of threshold region
+    // Start with the entire kernel, and keep reducing the size of the box until it drops below threshold
+    int box = kernels->size;                    // Size of box with bad pixels
+    for (double sumBox = sumKernel2; box > 0; box--) {
+        for (int x = -box; x <= box; x++) {
+            sumBox -= PS_SQR(kernel->kernel[-box][x]) + PS_SQR(kernel->kernel[box][x]);
+        }
+        for (int y = -box + 1; y <= box - 1; y++) {
+            // Note: not doing corners
+            sumBox -= PS_SQR(kernel->kernel[y][-box]) + PS_SQR(kernel->kernel[y][box]);
+        }
+        if (sumBox < threshold) {
+            break;
+        }
+    }
+
+    psFree(kernel);
+
+    return box;
+}
+
+psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
+                                             int numCols, int numRows, int x, int y)
+{
+    assert(kernels);
+    assert(numCols > 0 && numRows > 0);
+
+    // Size to use when calculating normalised coordinates (different from actual size when convolving
+    // subimage)
+    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
+    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
+
+    // Normalised coordinates
+    float yNorm = 2.0 * (float)(y - yNormSize/2.0) / (float)yNormSize;
+    float xNorm = 2.0 * (float)(x - xNormSize/2.0) / (float)xNormSize;
+
+    return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
 }
 
@@ -847,13 +884,4 @@
     int xMin = region->x0, xMax = region->x1, yMin = region->y0, yMax = region->y1; // Bounds of patch
 
-    // Size to use when calculating normalised coordinates (different from actual size when convolving
-    // subimage)
-    int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols);
-    int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows);
-
-    // Normalised coordinates
-    float yNorm = 2.0 * (float)(yMin + y0 + size + 1 - yNormSize/2.0) / (float)yNormSize; // Normalised coord
-    float xNorm = 2.0 * (float)(xMin + x0 + size + 1 - xNormSize/2.0) / (float)xNormSize; // Normalised coord
-
     psKernel *kernelImage = NULL;       // Kernel for the images
     psKernel *kernelWeight = NULL;      // Kernel for the weight maps
@@ -861,6 +889,7 @@
     // Only generate polynomial values every kernel footprint, since we have already assumed
     // (with the stamps) that it does not vary rapidly on this scale.
-    psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder,
-                                                    xNorm, yNorm); // Pre-calculated polynomial values
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
+                                                              xMin + x0 + size + 1,
+                                                              yMin + y0 + size + 1);
     float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term