Index: trunk/psModules/src/imcombine/pmStackReject.c
===================================================================
--- trunk/psModules/src/imcombine/pmStackReject.c	(revision 16476)
+++ trunk/psModules/src/imcombine/pmStackReject.c	(revision 16479)
@@ -11,54 +11,41 @@
 #define PIXEL_LIST_BUFFER 100           // Number of pixels to add to list at a time
 
-psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold,
-                        const psArray *subRegions, const psArray *kernels)
+psPixels *pmStackReject(const psPixels *in, float threshold, const psArray *regions,
+                        const psArray *solutions, const pmSubtractionKernels *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(subRegions, NULL);
-    PS_ASSERT_ARRAY_NON_NULL(kernels, NULL);
-    PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(regions, NULL);
+    PS_ASSERT_ARRAY_NON_NULL(solutions, NULL);
+    PS_ASSERT_ARRAYS_SIZE_EQUAL(regions, solutions, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernels, NULL);
 
     // Get the original image size
-    int numRegions = subRegions->n;        // Number of regions
+    int numRegions = regions->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 *subRegion = subRegions->data[i]; // Region of interest
-        if (subRegion->x0 < minCols) {
-            minCols = subRegion->x0;
+        psRegion *region = regions->data[i]; // Region of interest
+        if (region->x0 < minCols) {
+            minCols = region->x0;
         }
-        if (subRegion->y0 < minRows) {
-            minRows = subRegion->y0;
+        if (region->y0 < minRows) {
+            minRows = region->y0;
         }
-        if (subRegion->x1 > numCols) {
-            numCols = subRegion->x1;
+        if (region->x1 > numCols) {
+            numCols = region->x1;
         }
-        if (subRegion->y1 > numRows) {
-            numRows = subRegion->y1;
-        }
-
-        pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest
-        if (size == 0) {
-            size = kernel->size;
-        } else if (kernel->size != size) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d",
-                    size, kernel->size);
-            return NULL;
+        if (region->y1 > numRows) {
+            numRows = region->y1;
         }
     }
-
-    // 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);
+    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;
     }
 
-    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1),
-                                   0x01); // Mask
+    psImage *mask = psPixelsToMask(NULL, in, psRegionSet(0, numCols, 0, numRows), 0x01); // Mask image
     psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
     psFree(mask);
@@ -68,15 +55,7 @@
     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 = subRegions->data[i]; // Region of interest
-        if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 ||
-                      region->y0 > valid->y1 || region->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)) {
+        psRegion *region = regions->data[i]; // Region of interest
+        if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernels, true)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
             psFree(convRO);
@@ -89,10 +68,8 @@
 
         // Image of the kernel at the centre of the region
-        float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) /
-            (float)kernel->numCols;
-        float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) /
-            (float)kernel->numRows;
-        psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false);
-        if (!image) {
+        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;
+        psImage *kernel = pmSubtractionKernelImage(kernels, xNorm, yNorm, false);
+        if (!kernel) {
             psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image.");
             psFree(convRO);
@@ -101,23 +78,13 @@
         }
         float sum = 0.0;
-        for (int y = 0; y < image->numRows; y++) {
-            for (int x = 0; x < image->numCols; x++) {
-                sum += image->data.F32[y][x];
+        for (int y = 0; y < kernel->numRows; y++) {
+            for (int x = 0; x < kernel->numCols; x++) {
+                sum += kernel->data.F32[y][x];
             }
         }
-        psFree(image);
+        psFree(kernel);
 
-        // Range for normalisation
-        int yMin = PS_MAX(minRows, region->y0) - inRO->row0;
-        int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0;
-        int xMin = PS_MAX(minCols, region->x0) - inRO->col0;
-        int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0;
-        psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n",
-                sum, xMin, xMax, yMin, yMax);
-        for (int y = yMin; y <= yMax; y++) {
-            for (int x = xMin; x <= xMax; x++) {
-                convRO->image->data.F32[y][x] /= sum;
-            }
-        }
+        psImage *subConv = psImageSubset(convRO->image, *region); // Sub-image of convolved image
+        psBinaryOp(subConv, subConv, "*", psScalarAlloc(1.0 / sum, PS_TYPE_F32));
     }
     psFree(inRO);
@@ -127,6 +94,6 @@
     // Threshold the convolved image
     psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked
-    for (int y = 0; y < convolved->numRows; y++) {
-        for (int x = 0; x < convolved->numCols; x++) {
+    for (int y = 0; y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
             if (convolved->data.F32[y][x] > threshold) {
                 bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y);
@@ -136,7 +103,7 @@
     psFree(convolved);
 
-    // 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);
+    // Now, we want to convolve the original pixels properly
+    mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols, 0, numRows), 0xff);
+    int size = kernels->size;           // Size of kernels
     for (int i = 0; i < bad->n; i++) {
         int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest
@@ -148,5 +115,4 @@
         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;
             }