Index: /trunk/psModules/src/imcombine/pmStackReject.c
===================================================================
--- /trunk/psModules/src/imcombine/pmStackReject.c	(revision 16422)
+++ /trunk/psModules/src/imcombine/pmStackReject.c	(revision 16423)
@@ -11,31 +11,32 @@
 #define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
 
-psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions, const psArray *kernels)
+psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
+                        const psArray *subRegions, const psArray *kernels)
 {
     PS_ASSERT_PIXELS_NON_NULL(in, NULL);
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL);
     PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL);
-    PS_ASSERT_ARRAY_NON_NULL(regions, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL);
     PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, kernels, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
 
     // Get the original image size
-    int numRegions = regions->n;        // Number of regions
+    int numRegions = subRegions->n;        // Number of regions
     int numCols = 0, numRows = 0;       // Size of original image
     int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0
     int size = 0;                       // Size of kernel
     for (int i = 0; i < numRegions; i++) {
-        psRegion *region = regions->data[i]; // Region of interest
-        if (region->x0 < minCols) {
-            minCols = region->x0;
+        psRegion *subRegion = subRegions->data[i]; // Region of interest
+        if (subRegion->x0 < minCols) {
+            minCols = subRegion->x0;
         }
-        if (region->y0 < minRows) {
-            minRows = region->y0;
+        if (subRegion->y0 < minRows) {
+            minRows = subRegion->y0;
         }
-        if (region->x1 > numCols) {
-            numCols = region->x1;
+        if (subRegion->x1 > numCols) {
+            numCols = subRegion->x1;
         }
-        if (region->y1 > numRows) {
-            numRows = region->y1;
+        if (subRegion->y1 > numRows) {
+            numRows = subRegion->y1;
         }
 
@@ -49,11 +50,14 @@
         }
     }
-    if (minCols != 0 || minRows != 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                "Some error with image regions --- minimum coordinate is not 0,0");
-        return NULL;
+
+    // Adjust the size for the size of the subimage
+    if (valid) {
+        minCols = PS_MAX(valid->x0, minCols);
+        minRows = PS_MAX(valid->y0, minRows);
+        numCols = PS_MIN(valid->x1, numCols);
+        numRows = PS_MIN(valid->y1, numRows);
     }
 
-    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
+    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols, minRows, numRows), 0x01); // Mask
     psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
     psFree(mask);
@@ -63,8 +67,15 @@
     pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image
     inRO->image = image;
+    inRO->col0 = minCols;
+    inRO->row0 = minRows;
     for (int i = 0; i < numRegions; i++) {
-        psRegion *region = regions->data[i]; // Region of interest
+        psRegion *subRegion = subRegions->data[i]; // Region of interest
+        if (valid && (subRegion->x0 > valid->x1 || subRegion->x1 < valid->x0 ||
+                      subRegion->y0 > valid->y1 || subRegion->y1 < valid->y0)) {
+            // Region is outside of our sub-image
+            continue;
+        }
         pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
-        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) {
+        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, subRegion, kernel, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -77,6 +88,8 @@
 
         // Image of the kernel at the centre of the region
-        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - numCols/2.0) / (float)numCols;
-        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - numRows/2.0) / (float)numRows;
+        float xNorm = (subRegion->x0 + 0.5 * (subRegion->x1 - subRegion->x0) - kernel->numCols/2.0) /
+            (float)kernel->numCols;
+        float yNorm = (subRegion->y0 + 0.5 * (subRegion->y1 - subRegion->y0) - kernel->numRows/2.0) /
+            (float)kernel->numRows;
         psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
         if (!image) {
@@ -94,5 +107,5 @@
         psFree(image);
 
-        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
+        psImage *subConv = psImageSubset(convRO->image, *subRegion); // Sub-image of convolved image
         psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
@@ -114,4 +127,5 @@
     // Now, grow the mask to include everything that touches a bad pixel in the convolution
     mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff);
+    assert(mask->numCols == numCols && mask->numRows == numRows);
     for (int i = 0; i < bad->n; i++) {
         int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
@@ -123,4 +137,5 @@
         for (int y = yMin; y <= yMax; y++) {
             for (int x = xMin; x <= xMax; x++) {
+                assert(x < mask->numCols && y < mask->numRows);
                 mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff;
             }
Index: /trunk/psModules/src/imcombine/pmStackReject.h
===================================================================
--- /trunk/psModules/src/imcombine/pmStackReject.h	(revision 16422)
+++ /trunk/psModules/src/imcombine/pmStackReject.h	(revision 16423)
@@ -10,4 +10,5 @@
 /// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels
 psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image
+                        const psRegion *valid, ///< Valid region to consider
                         float threshold, ///< Threshold on convolved image, 0..1
                         const psArray *regions, ///< Array of image regions for image
Index: /trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 16422)
+++ /trunk/psModules/src/imcombine/pmSubtraction.c	(revision 16423)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.77 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-02-07 04:03:23 $
+ *  @version $Revision: 1.78 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-02-14 00:13:09 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -715,4 +715,5 @@
 
     int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions
+    int x0 = ro1->col0, y0 = ro1->row0; // Image offset
 
     psImage *convImage1 = out1->image;   // Convolved image
@@ -771,4 +772,9 @@
         yMax = PS_MIN(region->y1, yMax);
     }
+
+    // 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);
 
     psMaskType maskSource;              // Mask these pixels when convolving
@@ -798,8 +804,10 @@
     for (int j = yMin; j < yMax; j += fullSize) {
         int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest
-        float yNorm = 2.0 * (float)(j + size + 1 - numRows/2.0) / (float)numRows; // Normalised coordinate
+        float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) /
+            (float)yNormSize; // Normalised coordinate
         for (int i = xMin; i < xMax; i += fullSize) {
             int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest
-            float xNorm = 2.0 * (float)(i + size + 1 - numCols/2.0) / (float)numCols; // Normalised coordinate
+            float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) /
+                (float)xNormSize; // Normalised coordinate
 
             // Only generate polynomial values every kernel footprint, since we have already assumed
Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.c
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 16422)
+++ /trunk/psModules/src/imcombine/pmSubtractionKernels.c	(revision 16423)
@@ -180,4 +180,6 @@
     kernels->bgOrder = 0;
     kernels->mode = mode;
+    kernels->numCols = 0;
+    kernels->numRows = 0;
     kernels->solution1 = NULL;
     kernels->solution2 = NULL;
Index: /trunk/psModules/src/imcombine/pmSubtractionKernels.h
===================================================================
--- /trunk/psModules/src/imcombine/pmSubtractionKernels.h	(revision 16422)
+++ /trunk/psModules/src/imcombine/pmSubtractionKernels.h	(revision 16423)
@@ -39,4 +39,5 @@
     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
 } pmSubtractionKernels;