Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14581)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14602)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.43 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-08-21 18:16:54 $
+ *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-08-22 02:56:53 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -24,4 +24,7 @@
 
 #include "pmSubtraction.h"
+
+#define PIXEL_LIST_BUFFER 100           // Number of entries to add to pixel list at a time
+
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -880,5 +883,6 @@
     }
 
-    float background = solution->data.F64[solution->n-1]; // The difference in background
+    int numBackground = polyTerms(kernels->bgOrder); // Number of background terms
+    float background = solution->data.F64[solution->n - numBackground]; // The difference in background
     int numCols = inImage->numCols, numRows = inImage->numRows; // Image dimensions
 
@@ -1017,2 +1021,86 @@
     return true;
 }
+
+
+psPixels *pmSubtractionDeconvolveMask(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(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 = 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
+    for (int i = 0; i < numRegions; i++) {
+        psRegion *region = regions->data[i]; // Region of interest
+        if (region->x0 < minCols) {
+            minCols = region->x0;
+        }
+        if (region->y0 < minRows) {
+            minRows = region->y0;
+        }
+        if (region->x1 > numCols) {
+            numCols = region->x1;
+        }
+        if (region->y1 > numRows) {
+            numRows = region->y1;
+        }
+    }
+    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(0, numCols, 0, numRows), 0x01); // Mask image
+    psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve
+    psFree(mask);
+
+    // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding
+    psImage *convolved = NULL;          // Convolved mask image
+    psPixels *out = NULL;               // List of pixels that should be masked
+    for (int i = 0; i < numRegions; i++) {
+        psRegion *region = regions->data[i]; // Region of interest
+        psVector *solution = solutions->data[i]; // Solution of interest
+        if (!pmSubtractionConvolve(&convolved, NULL, NULL, image, NULL, NULL, 0, region, solution, kernels)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i);
+            psFree(convolved);
+            psFree(image);
+            return NULL;
+        }
+
+        // Need to adjust the thresholding level for the normalisation of the kernel --- the application of
+        // the kernel may scale the unit level that we've inserted.
+
+        // Image of the kernel at the centre
+        psImage *kernel = pmSubtractionKernelImage(solution, kernels,
+                                                   region->x0 + 0.5 * (region->x1 - region->x0),
+                                                   region->y0 + 0.5 * (region->y1 - region->y0));
+        float sum = 0.0;
+        for (int y = 0; y < kernel->numRows; y++) {
+            for (int x = 0; x < kernel->numCols; x++) {
+                sum += kernel->data.F32[y][x];
+            }
+        }
+        psFree(kernel);
+
+        // Threshold the convolved image
+        float adjustedThreshold = sum * threshold; // Threshold adjusted for the kernel normalisation
+        for (int y = region->y0; y < region->y1; y++) {
+            for (int x = region->x0; x < region->x1; x++) {
+                if (image->data.F32[y][x] > adjustedThreshold) {
+                    out = psPixelsAdd(out, PIXEL_LIST_BUFFER, x, y);
+                }
+            }
+        }
+    }
+    psFree(image);
+
+    return out;
+}