Index: /branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c
===================================================================
--- /branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c	(revision 26738)
+++ /branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c	(revision 26739)
@@ -27,4 +27,5 @@
     psImageMaskType maskIgnore,         // Ignore pixels with this mask
     psImageMaskType maskThresh,         // Give pixels this mask if below threshold
+    psRegion *region,                   // Region of interest
     const char *description             // Description of image
     )
@@ -37,9 +38,11 @@
     psImage *image = ro->image;         // Image
     psImage *mask = ro->mask;           // Mask
-    int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
+
+    psImage *subImage = psImageSubset(image, *region); // Image with region of interest
+    psImage *subMask = psImageSubset(mask, *region);  // Maks with region of interest
 
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);                               // Random number generator
-    if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) {
+    if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold.");
         psFree(rng);
@@ -49,10 +52,13 @@
     psFree(rng);
 
+    psFree(subImage);
+    psFree(subMask);
+
     float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip
     psFree(stats);
     psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description);
 
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
+    for (int y = region->y0; y < region->y1; y++) {
+        for (int x = region->x0; x < region->x1; x++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) {
                 continue;
@@ -94,8 +100,4 @@
         return false;
     }
-    if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
-        return false;
-    }
 
     pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image
@@ -104,8 +106,21 @@
         return false;
     }
-    if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
-        return false;
+
+    psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD,
+                                                          "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
+    psMetadataItem *regItem;        // Item with region
+    while ((regItem = psMetadataGetAndIncrement(regIter))) {
+        psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type");
+        psRegion *region = regItem->data.V; // Region of interest
+        if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+            return false;
+        }
+        if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+            return false;
+        }
     }
+    psFree(regIter);
 
     psFree(view);
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmStackReject.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmStackReject.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmStackReject.c	(revision 26739)
@@ -35,6 +35,6 @@
 {
     int size = kernels->size;           // Half-size of convolution kernel
-    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                              xMin + size + 1, yMin + size + 1); // Polynomial
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + size + 1,
+                                                              yMin + size + 1); // Polynomial
     int box = p_pmSubtractionBadRadius(NULL, kernels, polyValues, false, poorFrac); // Radius of bad box
     psTrace("psModules.imcombine", 10, "Growing by %d", box);
@@ -165,8 +165,8 @@
 
         // Image of the kernel at the centre of the region
-        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernels->numCols/2.0) /
-            (float)kernels->numCols;
-        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernels->numRows/2.0) /
-            (float)kernels->numRows;
+        float xNorm, yNorm;             // Normalised coordinates
+        p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, 0.5 * (region->x1 - region->x0),
+                                            0.5 * (region->y1 - region->y0),
+                                            kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
         psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
         if (!kernel) {
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c	(revision 26739)
@@ -424,4 +424,6 @@
                 } else if (*source & subConvPoor) {
                     *target |= maskPoor;
+                } else {
+                    *target &= ~maskBad & ~maskPoor;
                 }
             }
@@ -574,19 +576,21 @@
 }
 
+void p_pmSubtractionPolynomialNormCoords(float *xOut, float *yOut, float xIn, float yIn,
+                                         int xMin, int xMax, int yMin, int yMax)
+{
+    float xNormSize = xMax - xMin, yNormSize = yMax - yMin; // Size to use for normalisation
+    *xOut = 2.0 * (float)(xIn - xMin - xNormSize/2.0) / xNormSize;
+    *yOut = 2.0 * (float)(yIn - yMin - yNormSize/2.0) / yNormSize;
+    return;
+}
+
 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, const pmSubtractionKernels *kernels,
-                                             int numCols, int numRows, int x, int y)
+                                             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;
-
+
+    float xNorm, yNorm;                 // Normalised coordinates
+    p_pmSubtractionPolynomialNormCoords(&xNorm, &yNorm, x, y,
+                                        kernels->xMin, kernels->xMax, kernels->yMin, kernels->yMax);
     return p_pmSubtractionPolynomial(output, kernels->spatialOrder, xNorm, yNorm);
 }
@@ -1072,7 +1076,6 @@
     // 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_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                              xMin + x0 + size + 1,
-                                                              yMin + y0 + size + 1);
+    psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, xMin + x0 + size + 1,
+                                                              yMin + y0 + size + 1);        // Polynomial
     float background = doBG ? p_pmSubtractionSolutionBackground(kernels, polyValues) : 0.0; // Background term
 
@@ -1200,41 +1203,10 @@
     bool threaded = pmSubtractionThreaded(); // Running threaded?
 
-    // Outputs
-    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out1->image) {
-            out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        if (ro1->variance) {
-            if (!out1->variance) {
-                out1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            }
-            psImageInit(out1->variance, 0.0);
-        }
-    }
-    if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-        if (!out2->image) {
-            out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-        }
-        if (ro2->variance) {
-            if (!out2->variance) {
-                out2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
-            }
-            psImageInit(out2->variance, 0.0);
-        }
-    }
     psImage *convMask = NULL;           // Convolved mask image (common to inputs 1 and 2)
     if (subMask) {
         if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out1->mask) {
-                out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            }
-            psImageInit(out1->mask, 0);
             convMask = out1->mask;
         }
         if (kernels->mode == PM_SUBTRACTION_MODE_2 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-            if (!out2->mask) {
-                out2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
-            }
-            psImageInit(out2->mask, 0);
             if (!convMask) {
                 convMask = out2->mask;
@@ -1256,6 +1228,6 @@
 
     // Get region for convolution: [xMin:xMax,yMin:yMax]
-    int xMin = size, xMax = numCols - size;
-    int yMin = size, yMax = numRows - size;
+    int xMin = kernels->xMin + size, xMax = kernels->xMax - size;
+    int yMin = kernels->yMin + size, yMax = kernels->yMax - size;
     if (region) {
         xMin = PS_MAX(region->x0, xMin);
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.h	(revision 26739)
@@ -127,4 +127,11 @@
     );
 
+/// Return normalised coordinates
+void p_pmSubtractionPolynomialNormCoords(
+    float *xOut, float *yOut,           ///< Normalised coordinates, returned
+    float xIn, float yIn,               ///< Input coordinates
+    int xMin, int xMax, int yMin, int yMax ///< Bounds of validity
+    );
+
 /// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j
 psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL
@@ -138,5 +145,4 @@
 psImage *p_pmSubtractionPolynomialFromCoords(psImage *output, ///< Output matrix, or NULL
                                              const pmSubtractionKernels *kernels, ///< Kernel parameters
-                                             int numCols, int numRows, ///< Size of image of interest
                                              int x, int y ///< Position of interest
     );
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionAnalysis.c	(revision 26739)
@@ -118,5 +118,5 @@
 
     // sample difference images
-    { 
+    {
         psMetadataAddArray(analysis, PS_LIST_TAIL, "SUBTRACTION.SAMPLE.STAMP.SET", PS_META_DUPLICATE_OK, "Sample Difference Stamps", kernels->sampleStamps);
     }
@@ -273,6 +273,5 @@
     // Difference in background
     {
-        psImage *polyValues = p_pmSubtractionPolynomialFromCoords(NULL, kernels, numCols, numRows,
-                                                                  numCols / 2.0, numRows / 2.0); // Polynomial
+        psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, 0.0, 0.0); // Polynomial
         float bg = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background difference
 
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c	(revision 26739)
@@ -213,5 +213,4 @@
 
     *norm = normI2 / normI1;
-    fprintf(stderr, "Sums: %f %f %f, normWindow: %d\n", normI1, normI2, *norm, normWindow);
 
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
@@ -522,5 +521,4 @@
 
     *norm = normI2 / normI1;
-    fprintf(stderr, "Sums: %f %f %f - I1I1: %f\n", normI1, normI2, *norm, sumI1I1);
 
     if (mode & PM_SUBTRACTION_EQUATION_NORM) {
@@ -590,16 +588,6 @@
                 // Contribution to chi^2: a_i^2 P_i
                 psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
-                fprintf (stderr, "penalty main: %d %e vs %e x %f : %f\n",
-                         index,
-                         matrix->data.F64[index][index],
-                         norm, penalties->data.F32[i],
-                         (norm * penalties->data.F32[i]) / matrix->data.F64[index][index] );
                 matrix->data.F64[index][index] += norm * penalties->data.F32[i];
                 if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
-                    fprintf (stderr, "penalty dual: %d %e vs %e x %f : %f\n",
-                             index + numParams + 2,
-                             matrix->data.F64[index + numParams + 2][index + numParams + 2],
-                             norm, penalties->data.F32[i],
-                             (norm * penalties->data.F32[i]) / matrix->data.F64[index + numParams + 2][index + numParams + 2] );
                     matrix->data.F64[index + numParams + 2][index + numParams + 2] += norm * penalties->data.F32[i];
                     // matrix[i][i] is ~ (k_i * I_1)(k_i * I_1)
@@ -997,4 +985,5 @@
             // Solve normalisation and background separately
             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
 
             psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for norm
@@ -1014,34 +1003,28 @@
 
             fprintf(stderr, "Norm: %lf\n", normValue);
-            // fprintf(stderr, "BG: %lf\n", bgValue);
 
             // Solve kernel components
             for (int i = 0; i < numSolution1; i++) {
-                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
+                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
 
                 sumMatrix->data.F64[i][normIndex] = 0.0;
                 sumMatrix->data.F64[normIndex][i] = 0.0;
-
-                // sumMatrix->data.F64[i][bgIndex] = 0.0;
-                // sumMatrix->data.F64[bgIndex][i] = 0.0;
-            }
+            }
+            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
+            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
 
             sumMatrix->data.F64[normIndex][normIndex] = 1.0;
-            // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
-            // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
-
             sumVector->data.F64[normIndex] = 0.0;
-            // sumVector->data.F64[bgIndex] = 0.0;
 
             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
 
             solution->data.F64[normIndex] = normValue;
-            // solution->data.F64[bgIndex] = bgValue;
         }
 # endif
 
         if (!kernels->solution1) {
-            kernels->solution1 = psVectorAlloc (sumVector->n, PS_TYPE_F64);
-            psVectorInit (kernels->solution1, 0.0);
+            kernels->solution1 = psVectorAlloc(sumVector->n, PS_TYPE_F64);
+            psVectorInit(kernels->solution1, 0.0);
         }
 
@@ -1126,4 +1109,5 @@
             // Solve normalisation and background separately
             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
+            int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background
 
 #if 0
@@ -1159,38 +1143,31 @@
 
             double normValue = stats->robustMedian;
-            // double bgValue = 0.0;
 
             psFree(stats);
 
             fprintf(stderr, "Norm: %lf\n", normValue);
-            // fprintf(stderr, "BG: %lf\n", bgValue);
 
             // Solve kernel components
             for (int i = 0; i < numSolution2; i++) {
-                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i]; // + bgValue * sumMatrix->data.F64[bgIndex][i];
-                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1]; // + bgValue * sumMatrix->data.F64[bgIndex][i + numSolution1];
+                sumVector->data.F64[i] -= normValue * sumMatrix->data.F64[normIndex][i];
+                sumVector->data.F64[i + numSolution1] -= normValue * sumMatrix->data.F64[normIndex][i + numSolution1];
 
                 sumMatrix->data.F64[i][normIndex] = 0.0;
                 sumMatrix->data.F64[normIndex][i] = 0.0;
 
-                // sumMatrix->data.F64[i][bgIndex] = 0.0;
-                // sumMatrix->data.F64[bgIndex][i] = 0.0;
                 sumMatrix->data.F64[i + numSolution1][normIndex] = 0.0;
                 sumMatrix->data.F64[normIndex][i + numSolution1] = 0.0;
-                // sumMatrix->data.F64[i + numSolution1][bgIndex] = 0.0;
-                // sumMatrix->data.F64[bgIndex][i + numSolution1] = 0.0;
-            }
+            }
+            sumVector->data.F64[bgIndex] -= normValue * sumMatrix->data.F64[normIndex][bgIndex];
+            sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
+            sumMatrix->data.F64[normIndex][bgIndex] = 0.0;
 
             sumMatrix->data.F64[normIndex][normIndex] = 1.0;
-            // sumMatrix->data.F64[bgIndex][bgIndex] = 1.0;
-            // sumMatrix->data.F64[normIndex][bgIndex] = sumMatrix->data.F64[bgIndex][normIndex] = 0.0;
 
             sumVector->data.F64[normIndex] = 0.0;
-            // sumVector->data.F64[bgIndex] = 0.0;
 
             solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, NAN);
 
             solution->data.F64[normIndex] = normValue;
-            // solution->data.F64[bgIndex] = bgValue;
         }
 #endif
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.h	(revision 26739)
@@ -10,5 +10,5 @@
     PM_SUBTRACTION_EQUATION_BG      = 0x02,
     PM_SUBTRACTION_EQUATION_KERNELS = 0x04,
-    PM_SUBTRACTION_EQUATION_ALL     = 0x05, // value should be NORM | BG | KERNELS
+    PM_SUBTRACTION_EQUATION_ALL     = 0x07, // value should be NORM | BG | KERNELS
 } pmSubtractionEquationCalculationMode;
 
@@ -21,5 +21,5 @@
                                          const pmSubtractionKernels *kernels, ///< Kernel parameters
                                          int index, ///< Index of stamp
-					 const pmSubtractionEquationCalculationMode mode
+                                         const pmSubtractionEquationCalculationMode mode
     );
 
@@ -27,5 +27,5 @@
 bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps
                                     const pmSubtractionKernels *kernels, ///< Kernel parameters
-				    const pmSubtractionEquationCalculationMode mode
+                                    const pmSubtractionEquationCalculationMode mode
     );
 
@@ -33,5 +33,5 @@
 bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters
                                 const pmSubtractionStampList *stamps, ///< Stamps
-				const pmSubtractionEquationCalculationMode mode
+                                const pmSubtractionEquationCalculationMode mode
     );
 
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionIO.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionIO.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionIO.c	(revision 26739)
@@ -80,10 +80,5 @@
         while ((item = psMetadataGetAndIncrement(iter))) {
             assert(item->type == PS_DATA_UNKNOWN);
-            // Set the normalisation dimensions, since these will be otherwise unavailable when reading the
-            // images by scans.
             pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction
-            kernel->numCols = ro->image->numCols;
-            kernel->numRows = ro->image->numRows;
-
             kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel);
         }
@@ -126,6 +121,4 @@
                          kernel->bgOrder);
         psMetadataAddS32(row, PS_LIST_TAIL, NAME_MODE,  0, "Matching mode (enum)", kernel->mode);
-        psMetadataAddS32(row, PS_LIST_TAIL, NAME_COLS,  0, "Number of columns", kernel->numCols);
-        psMetadataAddS32(row, PS_LIST_TAIL, NAME_ROWS,  0, "Number of rows", kernel->numRows);
         if (kernel->mode == PM_SUBTRACTION_MODE_1 || kernel->mode == PM_SUBTRACTION_MODE_2) {
             psMetadataAddVector(row, PS_LIST_TAIL, NAME_SOL1, 0, "Solution vector 1", kernel->solution1);
@@ -318,6 +311,4 @@
 
         TABLE_LOOKUP(int, S32, bg,      NAME_BG);
-        TABLE_LOOKUP(int, S32, numCols, NAME_COLS);
-        TABLE_LOOKUP(int, S32, numRows, NAME_ROWS);
 
         TABLE_LOOKUP(float, F32, mean,      NAME_MEAN);
@@ -325,7 +316,5 @@
         TABLE_LOOKUP(int,   S32, numStamps, NAME_NUMSTAMPS);
 
-        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, mode);
-        kernels->numCols = numCols;
-        kernels->numRows = numRows;
+        pmSubtractionKernels *kernels = pmSubtractionKernelsFromDescription(description, bg, *region, mode);
         kernels->mean = mean;
         kernels->rms = rms;
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.c	(revision 26739)
@@ -197,5 +197,5 @@
     }
 
-#if 1
+#if 0
     fprintf(stderr, "%d raw: %lf, null: %f, min: %lf, max: %lf, moment: %lf\n", index, sum, preCalc->kernel->kernel[0][0], min, max, penalty);
 #endif
@@ -245,5 +245,5 @@
     }
 
-#if 1
+#if 0
     {
         double sum = 0.0;   // Sum of kernel component
@@ -276,5 +276,5 @@
 pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder,
                                                     const psVector *fwhmsIN, const psVector *ordersIN,
-                                                    float penalty, pmSubtractionMode mode)
+                                                    float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -305,5 +305,6 @@
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size, spatialOrder, penalty, mode); // The kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "ISIS(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 
@@ -331,5 +332,5 @@
 pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, int spatialOrder,
                                                       const psVector *fwhmsIN, const psVector *ordersIN,
-                                                      float penalty, pmSubtractionMode mode)
+                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -361,5 +362,6 @@
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size, spatialOrder, penalty, mode); // The kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS_RADIAL, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "ISIS_RADIAL(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 
@@ -388,5 +390,5 @@
 pmSubtractionKernels *pmSubtractionKernelsHERM(int size, int spatialOrder,
                                                const psVector *fwhmsIN, const psVector *ordersIN,
-                                               float penalty, pmSubtractionMode mode)
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -417,8 +419,10 @@
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size, spatialOrder, penalty, mode); // The kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_HERM, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 
-    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements", params, spatialOrder, num);
+    psLogMsg("psModules.imcombine", PS_LOG_INFO, "HERM kernel: %s,%d --> %d elements",
+             params, spatialOrder, num);
     psFree(params);
 
@@ -439,6 +443,6 @@
 
 pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, int spatialOrder,
-                                                     const psVector *fwhmsIN, const psVector *ordersIN,
-                                                     float penalty, pmSubtractionMode mode)
+                                                      const psVector *fwhmsIN, const psVector *ordersIN,
+                                                      float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_VECTOR_NON_NULL(fwhmsIN, NULL);
@@ -469,5 +473,6 @@
     }
 
-    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size, spatialOrder, penalty, mode); // The kernels
+    pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_DECONV_HERM, size,
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "DECONV_HERM(%d,%s,%d,%.2e)", size, params, spatialOrder, penalty);
 
@@ -544,5 +549,5 @@
 
 pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type,
-                                                int size, int spatialOrder, float penalty,
+                                                int size, int spatialOrder, float penalty, psRegion bounds,
                                                 pmSubtractionMode mode)
 {
@@ -558,4 +563,8 @@
     kernels->uStop = NULL;
     kernels->vStop = NULL;
+    kernels->xMin = bounds.x0;
+    kernels->xMax = bounds.x1;
+    kernels->yMin = bounds.y0;
+    kernels->yMax = bounds.y1;
     kernels->preCalc = psArrayAlloc(numBasisFunctions);
     kernels->penalty = penalty;
@@ -566,6 +575,4 @@
     kernels->bgOrder = 0;
     kernels->mode = mode;
-    kernels->numCols = 0;
-    kernels->numRows = 0;
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
@@ -640,5 +647,5 @@
 }
 
-pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty,
+pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, float penalty, psRegion bounds,
                                                pmSubtractionMode mode)
 {
@@ -649,5 +656,5 @@
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(0, PM_SUBTRACTION_KERNEL_POIS, size,
-                                                              spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     psStringAppend(&kernels->description, "POIS(%d,%d,%.2e)", size, spatialOrder, penalty);
     psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements",
@@ -664,8 +671,8 @@
 pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder,
                                                const psVector *fwhms, const psVector *orders,
-                                               float penalty, pmSubtractionMode mode)
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
-                                                                  penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     if (!kernels) {
         return NULL;
@@ -676,5 +683,5 @@
 
 pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning,
-                                               float penalty, pmSubtractionMode mode)
+                                               float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -697,5 +704,5 @@
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, size,
-                                                              spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d,%.2e)", size, inner, binning, spatialOrder,
@@ -768,5 +775,5 @@
 
 pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, float penalty,
-                                                pmSubtractionMode mode)
+                                                psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -795,5 +802,5 @@
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, size,
-                                                              spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "FRIES(%d,%d,%d,%.2e)", size, inner, spatialOrder, penalty);
@@ -864,5 +871,5 @@
 pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms,
                                                const psVector *orders, int inner, float penalty,
-                                               pmSubtractionMode mode)
+                                               psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -877,5 +884,5 @@
 
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
-                                                                  penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     kernels->type = PM_SUBTRACTION_KERNEL_GUNK;
     psStringPrepend(&kernels->description, "GUNK=");
@@ -893,5 +900,5 @@
 // RINGS --- just what it says
 pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder,
-                                                float penalty, pmSubtractionMode mode)
+                                                float penalty, psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_INT_POSITIVE(size, NULL);
@@ -924,5 +931,5 @@
 
     pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, size,
-                                                              spatialOrder, penalty, mode); // The kernels
+                                                              spatialOrder, penalty, bounds, mode); // Kernels
     kernels->inner = inner;
     psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d,%.2e)", size, inner, ringsOrder, spatialOrder,
@@ -1046,26 +1053,26 @@
 pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder,
                                                    const psVector *fwhms, const psVector *orders, int inner,
-                                                   int binning, int ringsOrder, float penalty,
+                                                   int binning, int ringsOrder, float penalty, psRegion bounds,
                                                    pmSubtractionMode mode)
 {
     switch (type) {
       case PM_SUBTRACTION_KERNEL_POIS:
-        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, mode);
+        return pmSubtractionKernelsPOIS(size, spatialOrder, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_ISIS:
-        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, mode);
+        return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_ISIS_RADIAL:
-        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, mode);
+        return pmSubtractionKernelsISIS_RADIAL(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_HERM:
-        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, mode);
+        return pmSubtractionKernelsHERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_DECONV_HERM:
-        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, mode);
+        return pmSubtractionKernelsDECONV_HERM(size, spatialOrder, fwhms, orders, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_SPAM:
-        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, mode);
+        return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_FRIES:
-        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, mode);
+        return pmSubtractionKernelsFRIES(size, spatialOrder, inner, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_GUNK:
-        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, mode);
+        return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, penalty, bounds, mode);
       case PM_SUBTRACTION_KERNEL_RINGS:
-        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, mode);
+        return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, penalty, bounds, mode);
       default:
         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type);
@@ -1103,5 +1110,5 @@
 
 pmSubtractionKernels *pmSubtractionKernelsFromDescription(const char *description, int bgOrder,
-                                                          pmSubtractionMode mode)
+                                                          psRegion bounds, pmSubtractionMode mode)
 {
     PS_ASSERT_STRING_NON_EMPTY(description, NULL);
@@ -1156,7 +1163,5 @@
         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
         penalty = parseStringFloat(ptr);
-
-        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
-
+        break;
       case PM_SUBTRACTION_KERNEL_RINGS:
         PARSE_STRING_NUMBER(size, ptr, ',', parseStringInt);
@@ -1165,9 +1170,11 @@
         PARSE_STRING_NUMBER(spatialOrder, ptr, ',', parseStringInt);
         PARSE_STRING_NUMBER(penalty, ptr, ')', parseStringInt);
-        return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning, ringsOrder, penalty, mode);
+        break;
       default:
         psAbort("Deciphering kernels other than ISIS, HERM, DECONV_HERM or RINGS is not currently supported.");
     }
-    return NULL;
+
+    return pmSubtractionKernelsGenerate(type, size, spatialOrder, fwhms, orders, inner, binning,
+                                        ringsOrder, penalty, bounds, mode);
 }
 
@@ -1245,6 +1252,8 @@
     out->bgOrder = in->bgOrder;
     out->mode = in->mode;
-    out->numCols = in->numCols;
-    out->numRows = in->numRows;
+    out->xMin = in->xMin;
+    out->xMax = in->xMax;
+    out->yMin = in->yMin;
+    out->yMax = in->yMax;
     out->solution1 = in->solution1 ? psVectorCopy(NULL, in->solution1, PS_TYPE_F64) : NULL;
     out->solution2 = in->solution2 ? psVectorCopy(NULL, in->solution2, PS_TYPE_F64) : NULL;
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionKernels.h	(revision 26739)
@@ -12,5 +12,5 @@
     PM_SUBTRACTION_KERNEL_ISIS_RADIAL,  ///< ISIS + higher-order radial Hermitians
     PM_SUBTRACTION_KERNEL_HERM,         ///< Hermitian polynomial kernels
-    PM_SUBTRACTION_KERNEL_DECONV_HERM,	///< Deconvolved Hermitian polynomial kernels
+    PM_SUBTRACTION_KERNEL_DECONV_HERM,  ///< Deconvolved Hermitian polynomial kernels
     PM_SUBTRACTION_KERNEL_SPAM,         ///< Summed Pixels for Advanced Matching --- summed delta functions
     PM_SUBTRACTION_KERNEL_FRIES,        ///< Fibonacci Radius Increases Excellence of Subtraction
@@ -32,4 +32,5 @@
     pmSubtractionKernelsType type;      ///< Type of kernels --- allowing the use of multiple kernels
     psString description;               ///< Description of the kernel parameters
+    int xMin, xMax, yMin, yMax;         ///< Bounds of image (for normalisation)
     long num;                           ///< Number of kernel components (not including the spatial ones)
     psVector *u, *v;                    ///< Offset (for POIS) or polynomial order (for ISIS, HERM or DECONV_HERM)
@@ -44,27 +45,26 @@
     int bgOrder;                        ///< The order for the background fitting
     pmSubtractionMode mode;             ///< Mode for subtraction
-    int numCols, numRows;               ///< Size of image (for normalisation), or zero to use image provided
     psVector *solution1, *solution2;    ///< Solution for the PSF matching
     // Quality information
     float mean, rms;                    ///< Mean and RMS of chi^2 from stamps
     int numStamps;                      ///< Number of good stamps
-    float fSigResMean;			///< mean fractional stdev of residuals
-    float fSigResStdev;			///< stdev of fractional stdev of residuals
-    float fMaxResMean;			///< mean fractional positive swing in residuals
-    float fMaxResStdev;			///< stdev of fractional positive swing in residuals
-    float fMinResMean;			///< mean fractional negative swing in residuals
-    float fMinResStdev;			///< stdev of fractional negative swing in residuals
-    psArray *sampleStamps;		///< array of brightest set of stamps for output visualizations
+    float fSigResMean;                  ///< mean fractional stdev of residuals
+    float fSigResStdev;                 ///< stdev of fractional stdev of residuals
+    float fMaxResMean;                  ///< mean fractional positive swing in residuals
+    float fMaxResStdev;                 ///< stdev of fractional positive swing in residuals
+    float fMinResMean;                  ///< mean fractional negative swing in residuals
+    float fMinResStdev;                 ///< stdev of fractional negative swing in residuals
+    psArray *sampleStamps;              ///< array of brightest set of stamps for output visualizations
 } pmSubtractionKernels;
 
 // pmSubtractionKernels->preCalc is an array of pmSubtractionKernelPreCalc structures
 typedef struct {
-    psVector *uCoords;			// used by RINGS
-    psVector *vCoords;			// used by RINGS
-    psVector *poly;			// used by RINGS
-
-    psVector *xKernel;			// used by ISIS, HERM, DECONV_HERM
-    psVector *yKernel;			// used by ISIS, HERM, DECONV_HERM
-    psKernel *kernel;			// used by ISIS, HERM, DECONV_HERM
+    psVector *uCoords;                  // used by RINGS
+    psVector *vCoords;                  // used by RINGS
+    psVector *poly;                     // used by RINGS
+
+    psVector *xKernel;                  // used by ISIS, HERM, DECONV_HERM
+    psVector *yKernel;                  // used by ISIS, HERM, DECONV_HERM
+    psKernel *kernel;                   // used by ISIS, HERM, DECONV_HERM
 } pmSubtractionKernelPreCalc;
 
@@ -162,4 +162,5 @@
                                                 int spatialOrder, ///< Order of spatial variations
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -168,8 +169,8 @@
 pmSubtractionKernelPreCalc *pmSubtractionKernelPreCalcAlloc(
     pmSubtractionKernelsType type, ///< type of kernel to allocate (not all can be pre-calculated)
-    int uOrder,			   ///< order in x-direction 
-    int vOrder,			   ///< order in x-direction 
-    int size,			   ///< Half-size of the kernel
-    float sigma			   ///< sigma of gaussian kernel
+    int uOrder,                    ///< order in x-direction
+    int vOrder,                    ///< order in x-direction
+    int size,                      ///< Half-size of the kernel
+    float sigma                    ///< sigma of gaussian kernel
     );
 
@@ -179,4 +180,5 @@
                                                int spatialOrder, ///< Order of spatial variations
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -188,4 +190,5 @@
                                                     const psVector *orders, ///< Polynomial order of gaussians
                                                     float penalty, ///< Penalty for wideness
+                                                    psRegion bounds,       ///< Bounds for validity
                                                     pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -197,4 +200,5 @@
                                                const psVector *orders, ///< Polynomial order of gaussians
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
                                                );
@@ -202,9 +206,10 @@
 /// Generate ISIS + RADIAL_HERM kernels
 pmSubtractionKernels *pmSubtractionKernelsISIS_RADIAL(int size, ///< Half-size of the kernel
-						      int spatialOrder, ///< Order of spatial variations
-						      const psVector *fwhms, ///< Gaussian FWHMs
-						      const psVector *orders, ///< Polynomial order of gaussians
-						      float penalty, ///< Penalty for wideness
-						      pmSubtractionMode mode ///< Mode for subtraction
+                                                      int spatialOrder, ///< Order of spatial variations
+                                                      const psVector *fwhms, ///< Gaussian FWHMs
+                                                      const psVector *orders, ///< Polynomial order of gaussians
+                                                      float penalty, ///< Penalty for wideness
+                                                      psRegion bounds,       ///< Bounds for validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
                                                );
 
@@ -215,4 +220,5 @@
                                                const psVector *orders, ///< order of hermitian polynomials
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
                                                );
@@ -220,9 +226,10 @@
 /// Generate DECONV_HERM kernels
 pmSubtractionKernels *pmSubtractionKernelsDECONV_HERM(int size, ///< Half-size of the kernel
-						     int spatialOrder, ///< Order of spatial variations
-						     const psVector *fwhms, ///< Gaussian FWHMs
-						     const psVector *orders, ///< order of hermitian polynomials
-						     float penalty, ///< Penalty for wideness
-						     pmSubtractionMode mode ///< Mode for subtraction
+                                                      int spatialOrder, ///< Order of spatial variations
+                                                      const psVector *fwhms, ///< Gaussian FWHMs
+                                                      const psVector *orders, ///< order of hermitian polynomials
+                                                      float penalty, ///< Penalty for wideness
+                                                      psRegion bounds,       ///< Bounds for validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
     );
 
@@ -233,4 +240,5 @@
                                                int binning, ///< Kernel binning factor
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -241,4 +249,5 @@
                                                 int inner, ///< Inner radius to preserve unbinned
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -251,4 +260,5 @@
                                                int inner, ///< Inner radius containing grid of delta functions
                                                float penalty, ///< Penalty for wideness
+                                               psRegion bounds,       ///< Bounds for validity
                                                pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -260,4 +270,5 @@
                                                 int ringsOrder, ///< Polynomial order
                                                 float penalty, ///< Penalty for wideness
+                                                psRegion bounds,       ///< Bounds for validity
                                                 pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -274,4 +285,5 @@
                                                    int ringsOrder, ///< Polynomial order for RINGS
                                                    float penalty, ///< Penalty for wideness
+                                                   psRegion bounds,       ///< Bounds for validity
                                                    pmSubtractionMode mode ///< Mode for subtraction
     );
@@ -281,4 +293,5 @@
     const char *description,            ///< Description of kernel
     int bgOrder,                        ///< Polynomial order for background fitting
+    psRegion bounds,                    ///< Bounds for validity
     pmSubtractionMode mode              ///< Mode for subtraction
     );
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.c	(revision 26739)
@@ -38,14 +38,15 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psImageMaskType maskVal,
-                           int size, int footprint, float badFrac, pmSubtractionMode mode)
-{
-    PS_ASSERT_IMAGE_NON_NULL(mask1, NULL);
-    PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_IMAGE_MASK, NULL);
-    if (mask2) {
-        PS_ASSERT_IMAGE_NON_NULL(mask2, NULL);
-        PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_IMAGE_MASK, NULL);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL);
-    }
+psImage *pmSubtractionMask(psRegion *bounds, const pmReadout *ro1, const pmReadout *ro2,
+                           psImageMaskType maskVal, int size, int footprint, float badFrac,
+                           pmSubtractionMode mode)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro1, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro1, NULL);
+    PM_ASSERT_READOUT_MASK(ro1, NULL);
+    PM_ASSERT_READOUT_NON_NULL(ro2, NULL);
+    PM_ASSERT_READOUT_IMAGE(ro2, NULL);
+    PM_ASSERT_READOUT_MASK(ro2, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL);
     PS_ASSERT_INT_NONNEGATIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
@@ -55,28 +56,38 @@
     }
 
-    int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images
+    int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images
 
     // Dereference inputs for convenience
-    psImageMaskType **data1 = mask1->data.PS_TYPE_IMAGE_MASK_DATA;
-    psImageMaskType **data2 = NULL;
-    if (mask2) {
-        data2 = mask2->data.PS_TYPE_IMAGE_MASK_DATA;
-    }
+    psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32;
+    psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA,
+        **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA;
 
     // First, a pass through to determine the fraction of bad pixels
-    if (isfinite(badFrac) && badFrac != 1.0) {
+    if (bounds || (isfinite(badFrac) && badFrac != 1.0)) {
+        int xMin = numCols, xMax = 0, yMin = numRows, yMax = 0; // Bounds of good pixels
         int numBad = 0;                 // Number of bad pixels
         for (int y = 0; y < numRows; y++) {
             for (int x = 0; x < numCols; x++) {
-                if (data1[y][x] & maskVal) {
+                if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) {
                     numBad++;
                     continue;
                 }
-                if (data2 && data2[y][x] & maskVal) {
+                if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) {
                     numBad++;
+                    continue;
                 }
-            }
-        }
-        if (numBad > badFrac * numCols * numRows) {
+                xMin = PS_MIN(xMin, x);
+                xMax = PS_MAX(xMax, x);
+                yMin = PS_MIN(yMin, y);
+                yMax = PS_MAX(yMax, y);
+            }
+        }
+        if (bounds) {
+            bounds->x0 = xMin;
+            bounds->x1 = xMax;
+            bounds->y0 = yMin;
+            bounds->y1 = yMax;
+        }
+        if (isfinite(badFrac) && badFrac != 1.0 && numBad > badFrac * numCols * numRows) {
             psError(PM_ERR_SMALL_AREA, true,
                     "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n",
@@ -117,8 +128,8 @@
     for (int y = 0; y < numRows; y++) {
         for (int x = 0; x < numCols; x++) {
-            if (data1[y][x] & maskVal) {
+            if (maskData1[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1;
             }
-            if (data2 && data2[y][x] & maskVal) {
+            if (maskData2[y][x] & maskVal) {
                 maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2;
             }
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMask.h	(revision 26739)
@@ -5,11 +5,13 @@
 
 /// Generate a mask for use in the subtraction process
-psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved)
-                           const psImage *inMask, ///< Mask for the input image, or NULL
-                           psImageMaskType maskVal, ///< Value to mask out
-                           int size, ///< Half-size of the kernel (pmSubtractionKernels.size)
-                           int footprint, ///< Half-size of the kernel footprint
-                           float badFrac, ///< Maximum fraction of bad input pixels to accept
-                           pmSubtractionMode mode  ///< Subtraction mode
+psImage *pmSubtractionMask(
+    psRegion *bounds,                   ///< Bounds of valid pixels (or NULL), returned
+    const pmReadout *ro1,               ///< Readout 1
+    const pmReadout *ro2,               ///< Readout 2
+    psImageMaskType maskVal,            ///< Value to mask out
+    int size,                           ///< Half-size of the kernel (pmSubtractionKernels.size)
+    int footprint,                      ///< Half-size of the kernel footprint
+    float badFrac,                      ///< Maximum fraction of bad input pixels to accept
+    pmSubtractionMode mode              ///< Subtraction mode
     );
 
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionMatch.c	(revision 26739)
@@ -70,5 +70,5 @@
                                  const psImage *subMask, // Mask for subtraction, or NULL
                                  psImage *variance,  // Variance map
-                                 const psRegion *region, // Region of interest, or NULL
+                                 const psRegion *region, // Region of interest
                                  float thresh1,  // Threshold for stamp finding on readout 1
                                  float thresh2,  // Threshold for stamp finding on readout 2
@@ -82,4 +82,11 @@
     )
 {
+    PS_ASSERT_PTR_NON_NULL(stamps, false);
+    PM_ASSERT_READOUT_NON_NULL(ro1, false);
+    PM_ASSERT_READOUT_NON_NULL(ro2, false);
+    PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+    PS_ASSERT_IMAGE_NON_NULL(variance, false);
+    PS_ASSERT_PTR_NON_NULL(region, false);
+
     psTrace("psModules.imcombine", 3, "Finding stamps...\n");
 
@@ -96,5 +103,5 @@
 
     psTrace("psModules.imcombine", 3, "Extracting stamps...\n");
-    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, variance, size)) {
+    if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2->image, variance, size, *region)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps.");
         return false;
@@ -319,5 +326,7 @@
     pmSubtractionMaskInvalid(ro2, maskVal);
 
-    psImage *subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, 0,
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+
+    psImage *subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, 0,
                                          badFrac, mode); // Subtraction mask
     if (!subMask) {
@@ -441,5 +450,5 @@
     // Putting important variable declarations here, since they are freed after a "goto" if there is an error.
     psImage *subMask = NULL;            // Mask for subtraction
-    psRegion *region = NULL;            // Iso-kernel region
+    psRegion *region = psRegionAlloc(NAN, NAN, NAN, NAN); // Iso-kernel region
     psString regionString = NULL;       // String for region
     pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF
@@ -457,6 +466,7 @@
     pmSubtractionMaskInvalid(ro2, maskVal);
 
-    subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskVal, size, footprint,
-                                badFrac, subMode);
+    psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels
+
+    subMask = pmSubtractionMask(&bounds, ro1, ro2, maskVal, size, footprint, badFrac, subMode);
     if (!subMask) {
         psError(psErrorCodeLast(), false, "Unable to generate subtraction mask.");
@@ -468,11 +478,13 @@
     // Get region of interest
     int xRegions = 1, yRegions = 1;     // Number of iso-kernel regions
-    float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions
+    float xRegionSize = NAN, yRegionSize = NAN; // Size of iso-kernel regions
     if (isfinite(regionSize) && regionSize != 0.0) {
-        xRegions = numCols / regionSize + 1;
-        yRegions = numRows / regionSize + 1;
-        xRegionSize = (float)numCols / (float)xRegions;
-        yRegionSize = (float)numRows / (float)yRegions;
-        region = psRegionAlloc(NAN, NAN, NAN, NAN);
+        xRegions = (bounds.x1 - bounds.x0) / regionSize + 1;
+        yRegions = (bounds.y1 - bounds.y0) / regionSize + 1;
+        xRegionSize = (float)(bounds.x1 - bounds.x0) / (float)xRegions;
+        yRegionSize = (float)(bounds.y1 - bounds.y0) / (float)yRegions;
+    } else {
+        xRegionSize = bounds.x1 - bounds.x0;
+        yRegionSize = bounds.y1 - bounds.y0;
     }
 
@@ -480,5 +492,5 @@
     float stampThresh1 = NAN, stampThresh2 = NAN; // Stamp thresholds for images
     {
-        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for backgroun
+        psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background
         if (ro1) {
             psStatsInit(bg);
@@ -504,4 +516,45 @@
     }
 
+    // Just in case the iso-kernel region doesn't cover the entire image...
+    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv1->image) {
+            conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv1->image, NAN);
+        if (ro1->variance) {
+            if (!conv1->variance) {
+                conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv1->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv1->mask) {
+                conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv1->mask, maskBad);
+        }
+    }
+    if (subMode == PM_SUBTRACTION_MODE_1 || subMode == PM_SUBTRACTION_MODE_UNSURE ||
+        subMode == PM_SUBTRACTION_MODE_DUAL) {
+        if (!conv2->image) {
+            conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+        }
+        psImageInit(conv2->image, NAN);
+        if (ro2->variance) {
+            if (!conv2->variance) {
+                conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32);
+            }
+            psImageInit(conv2->variance, NAN);
+        }
+        if (subMask) {
+            if (!conv2->mask) {
+                conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK);
+            }
+            psImageInit(conv2->mask, maskBad);
+        }
+    }
+
+
     // Iterate over iso-kernel regions
     for (int j = 0; j < yRegions; j++) {
@@ -509,12 +562,12 @@
             psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n",
                     j * xRegions + i + 1, xRegions * yRegions);
-            if (region) {
-                *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize),
-                                      (int)(j * yRegionSize), (int)((j + 1) * yRegionSize));
-                psFree(regionString);
-                regionString = psRegionToString(*region);
-                psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n",
-                        regionString, numCols, numRows);
-            }
+            *region = psRegionSet(bounds.x0 + (int)(i * xRegionSize),
+                                  bounds.x0 + (int)((i + 1) * xRegionSize),
+                                  bounds.y0 + (int)(j * yRegionSize),
+                                  bounds.y0 + (int)((j + 1) * yRegionSize));
+            psFree(regionString);
+            regionString = psRegionToString(*region);
+            psLogMsg("psModules.imcombine", PS_LOG_DETAIL, "Iso-kernel region: %s out of %d,%d\n",
+                    regionString, numCols, numRows);
 
             if (stampsName && strlen(stampsName) > 0) {
@@ -530,5 +583,5 @@
             // We get the stamps here; we will also attempt to get stamps at the first iteration, but it
             // doesn't matter.
-            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, NULL, stampThresh1, stampThresh2,
+            if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region, stampThresh1, stampThresh2,
                                       stampSpacing, normFrac, sysError, skyError, size, footprint, subMode)) {
                 goto MATCH_ERROR;
@@ -544,6 +597,7 @@
             // Define kernel basis functions
             if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) {
-                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder,
-                                                          stamps, footprint, optThreshold, penalty, subMode);
+                kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder,
+                                                          optFWHMs, optOrder, stamps, footprint,
+                                                          optThreshold, penalty, bounds, subMode);
                 if (!kernels) {
                     psErrorClear();
@@ -554,5 +608,5 @@
                 // Not an ISIS/GUNK kernel, or the optimum kernel search failed
                 kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders,
-                                                       inner, binning, ringsOrder, penalty, subMode);
+                                                       inner, binning, ringsOrder, penalty, bounds, subMode);
                 // pmSubtractionVisualShowKernels(kernels);
             }
@@ -603,5 +657,4 @@
                 psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k);
 
-#if 0
                 if (!subtractionGetStamps(&stamps, ro1, ro2, subMask, variance, region,
                                           stampThresh1, stampThresh2, stampSpacing, normFrac,
@@ -609,5 +662,4 @@
                     goto MATCH_ERROR;
                 }
-#endif
 
                 // generate the window function from the set of stamps
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.c	(revision 26739)
@@ -204,5 +204,6 @@
                                                       int spatialOrder, const psVector *fwhms, int maxOrder,
                                                       const pmSubtractionStampList *stamps, int footprint,
-                                                      float tolerance, float penalty, pmSubtractionMode mode)
+                                                      float tolerance, float penalty, psRegion bounds,
+                                                      pmSubtractionMode mode)
 {
     if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) {
@@ -232,5 +233,5 @@
     psVectorInit(orders, maxOrder);
     pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, fwhms, orders,
-                                                                  penalty, mode); // Kernels
+                                                                  penalty, bounds, mode); // Kernels
     psFree(orders);
     psFree(kernels->description);
@@ -483,10 +484,10 @@
     if (type == PM_SUBTRACTION_KERNEL_ISIS) {
 
-	// XXX in r26035, this code was just wrong.  we had:
-
-	// psKernel *subtract = kernels->preCalc->data[0]
-
-	// but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
-	// an array of pmSubtractionKernelPreCalc.
+        // XXX in r26035, this code was just wrong.  we had:
+
+        // psKernel *subtract = kernels->preCalc->data[0]
+
+        // but, kernels->preCalc was an array of psArray, not an array of kernels.  It is now
+        // an array of pmSubtractionKernelPreCalc.
 
         pmSubtractionKernelPreCalc *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionParams.h	(revision 26739)
@@ -17,5 +17,6 @@
                                                       float tolerance, ///< Maximum difference in chi^2
                                                       float penalty, ///< Penalty for wideness
-                                                      pmSubtractionMode mode // Mode for subtraction
+                                                      psRegion bounds,       ///< Bounds of validity
+                                                      pmSubtractionMode mode ///< Mode for subtraction
     );
 
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.c	(revision 26739)
@@ -709,5 +709,4 @@
         }
     }
-    fprintf(stderr, "Normalization window radius: %d\n", stamps->normWindow);
 
     psTrace("psModules.imcombine", 3, "Normalisation window radius set to %d\n", stamps->normWindow);
@@ -741,5 +740,5 @@
 
 bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2,
-                                psImage *variance, int kernelSize)
+                                psImage *variance, int kernelSize, psRegion bounds)
 {
     PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false);
@@ -758,5 +757,4 @@
     }
 
-    int numCols = image1->numCols, numRows = image1->numRows; // Size of images
     int size = kernelSize + stamps->footprint; // Size of postage stamps
 
@@ -767,14 +765,10 @@
         }
 
-        if (isnan(stamp->xNorm)) {
-            stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols;
-        }
-        if (isnan(stamp->yNorm)) {
-            stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows;
-        }
+        p_pmSubtractionPolynomialNormCoords(&stamp->xNorm, &stamp->yNorm, stamp->x, stamp->y,
+                                            bounds.x0, bounds.x1, bounds.y0, bounds.y1);
 
         int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates
-        if (x < size || x > numCols - size || y < size || y > numRows - size) {
-            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y);
+        if (x < bounds.x0 + size || x > bounds.x1 - size || y < bounds.y0 + size || y > bounds.y1 - size) {
+            psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the region border.\n", i, x, y);
             return false;
         }
Index: /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h
===================================================================
--- /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h	(revision 26738)
+++ /branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionStamps.h	(revision 26739)
@@ -164,5 +164,6 @@
                                 psImage *image2, ///< Input image (or NULL)
                                 psImage *variance, ///< Variance map
-                                int kernelSize ///< Kernel half-size
+                                int kernelSize, ///< Kernel half-size
+                                psRegion bounds ///< Bounds of validity
     );