Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14125)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14195)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.16 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-11 01:28:39 $
+ *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-07-13 20:28:28 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -410,4 +410,108 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+psImage *pmSubtractionMask(const psImage *inMask, const psImage *refMask, psMaskType maskVal,
+                           int size, int footprint)
+{
+    PS_ASSERT_IMAGE_NON_NULL(inMask, NULL);
+    PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(refMask, NULL);
+    PS_ASSERT_IMAGE_TYPE(refMask, PS_TYPE_MASK, NULL);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, refMask, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(size, NULL);
+    PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
+
+    // Size of the images
+    int numCols = inMask->numCols;
+    int numRows = inMask->numRows;
+
+    // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask
+    psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask
+    psImageInit(mask, 0);
+
+    // Dereference for convenience
+    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA;
+    psMaskType **inData = inMask->data.PS_TYPE_MASK_DATA;
+    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;
+
+    // Block out a border around the edge of the image
+
+    // Bottom stripe
+    for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) {
+        for (int x = 0; x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+    // Either side
+    for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) {
+        for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+        for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+    // Top stripe
+    for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) {
+        for (int x = 0; x < numCols; x++) {
+            maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER;
+        }
+    }
+
+    // Mask the bad pixels, and around them
+    for (int y = footprint + size; y < numRows - footprint - size; y++) {
+        for (int x = 0; x < numCols; x++) {
+            if (inData[y][x] & maskVal) {
+                maskData[y][x] |= PM_SUBTRACTION_MASK_INPUT;
+                // Block out the entire stamp footprint around this pixel
+                for (int v = y - footprint; v <= y + footprint; v++) {
+                    for (int u = x - footprint; u <= x + footprint; u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
+                    }
+                }
+            }
+            if (refData[y][x] & maskVal) {
+                maskData[y][x] |= PM_SUBTRACTION_MASK_REF;
+
+                // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with
+                // this bad pixel (within 'size').  Then we want to block out with the FOOTPRINT mask
+                // everything within a footprint's distance of those (within 'footprint').
+
+                // Bottom stripe
+                for (int v = PS_MAX(y - footprint - size, 0); v < y - size; v++) {
+                    for (int u = PS_MAX(x - footprint - size, 0);
+                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
+                    }
+                }
+                // Middle
+                for (int v = PS_MAX(y - size, 0); v <= PS_MIN(y + size, numRows - 1); v++) {
+                    // Left side
+                    for (int u = PS_MAX(x - footprint - size, 0); u < x - size; u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
+                    }
+                    // Centre
+                    for (int u = PS_MAX(x - size, 0); u <= PS_MIN(x + size, numCols - 1); u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT | PM_SUBTRACTION_MASK_CONVOLVE;
+                    }
+                    // Right side
+                    for (int u = x + size + 1; u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
+                    }
+                }
+                // Top stripe
+                for (int v = y + size + 1; v <= PS_MIN(y + footprint + size, numRows - 1); v++) {
+                    for (int u = PS_MAX(x - footprint - size, 0);
+                         u <= PS_MIN(x + footprint + size, numCols - 1); u++) {
+                        maskData[v][u] |= PM_SUBTRACTION_MASK_FOOTPRINT;
+                    }
+                }
+            }
+        }
+    }
+
+    return mask;
+}
+
+
 bool pmSubtractionCalculateEquation(psArray *stamps, const psImage *reference, const psImage *input,
                                     const psImage *weight, const pmSubtractionKernels *kernels, int footprint)
@@ -507,13 +611,30 @@
             psFree(polyValues);
 
-            // Fill in lower diagonal of symmetric matrix
+            // Fill in lower diagonal of symmetric matrix, while checking for bad values
+            bool bad = false;           // Are there bad values?
             for (int i = 0; i < numKernels; i++) {
                 for (int j = 0; j < i; j++) {
                     stampMatrix->data.F64[i][j] = stampMatrix->data.F64[j][i];
+                    if (!isfinite(stampMatrix->data.F64[j][i])) {
+                        bad = true;
+                    }
                 }
                 stampMatrix->data.F64[bgIndex][i] = stampMatrix->data.F64[i][bgIndex];
-            }
-
-            stamp->status = PM_SUBTRACTION_STAMP_USED;
+                if (!isfinite(stampMatrix->data.F64[i][bgIndex]) ||
+                    !isfinite(stampMatrix->data.F64[i][i]) ||
+                    !isfinite(stampVector->data.F64[i])) {
+                    bad = true;
+                }
+            }
+            if (!isfinite(stampVector->data.F64[bgIndex])) {
+                bad = true;
+            }
+
+            if (bad) {
+                stamp->status = PM_SUBTRACTION_STAMP_REJECTED;
+                psTrace("psModules.imcombine", 3, "Rejecting stamp %d because of bad equation\n", i);
+            } else {
+                stamp->status = PM_SUBTRACTION_STAMP_USED;
+            }
 
             if (psTraceGetLevel("psModules.imcombine.equation") >= 10) {
@@ -610,6 +731,6 @@
 
 int pmSubtractionRejectStamps(psArray *stamps, const psImage *refImage, psImage *inImage,
-                              psImage *mask, psMaskType badStampMaskVal, const psVector *solution,
-                              int footprint, float sigmaRej, const pmSubtractionKernels *kernels)
+                              psImage *subMask, const psVector *solution, int footprint, float sigmaRej,
+                              const pmSubtractionKernels *kernels)
 {
     PS_ASSERT_ARRAY_NON_NULL(stamps, -1);
@@ -619,7 +740,7 @@
     PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, -1);
     PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, inImage, -1);
-    PS_ASSERT_IMAGE_NON_EMPTY(mask, -1);
-    PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, -1);
-    PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, mask, -1);
+    PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1);
+    PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1);
+    PS_ASSERT_IMAGES_SIZE_EQUAL(refImage, subMask, -1);
     PS_ASSERT_VECTOR_NON_NULL(solution, -1);
     PS_ASSERT_VECTOR_TYPE(solution, PS_TYPE_F64, -1);
@@ -714,5 +835,5 @@
             for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) {
                 for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) {
-                    mask->data.PS_TYPE_MASK_DATA[y][x] |= badStampMaskVal;
+                    subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ;
                 }
             }
@@ -767,14 +888,13 @@
 
 bool pmSubtractionConvolve(psImage **outImage, psImage **outWeight, psImage **outMask,
-                           const psImage *inImage, const psImage *inWeight, const psImage *inMask,
-                           psMaskType maskVal, psMaskType blank,
-                           const psVector *solution, const pmSubtractionKernels *kernels)
+                           const psImage *inImage, const psImage *inWeight, const psImage *subMask,
+                           psMaskType blank, const psVector *solution, const pmSubtractionKernels *kernels)
 {
     PS_ASSERT_IMAGE_NON_NULL(inImage, false);
     PS_ASSERT_IMAGE_TYPE(inImage, PS_TYPE_F32, false);
-    if (inMask) {
-        PS_ASSERT_IMAGE_NON_NULL(inMask, false);
-        PS_ASSERT_IMAGE_TYPE(inMask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGES_SIZE_EQUAL(inMask, inImage, false);
+    if (subMask) {
+        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
+        PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, inImage, false);
     }
     if (inWeight) {
@@ -793,5 +913,5 @@
         PS_ASSERT_IMAGES_SIZE_EQUAL(*outMask, inImage, false);
         PS_ASSERT_IMAGE_TYPE(*outMask, PS_TYPE_MASK, false);
-        PS_ASSERT_IMAGE_NON_NULL(inMask, false);
+        PS_ASSERT_IMAGE_NON_NULL(subMask, false);
     }
     if (outWeight && *outWeight) {
@@ -823,5 +943,5 @@
     }
     psImage *convMask = NULL;           // Convolved mask image
-    if (outMask && inMask) {
+    if (outMask && subMask) {
         if (!*outMask) {
             *outMask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
@@ -872,17 +992,10 @@
 
                     // Check and propagate the kernel footprint, if required
-                    if (inMask) {
-                        for (int v = -size; v <= size; v++) {
-                            for (int u = -size; u <= size; u++) {
-                                convMask->data.PS_TYPE_MASK_DATA[y][x] |=
-                                    inMask->data.PS_TYPE_MASK_DATA[y + v][x + u];
-                            }
-                        }
-                        if (convMask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) {
-                            convImage->data.F32[y][x] = NAN;
-                            if (inWeight) {
-                                convWeight->data.F32[y][x] = NAN;
-                            }
-                            continue;
+                    if (subMask && (subMask->data.PS_TYPE_MASK_DATA[y][x] &
+                                    (PM_SUBTRACTION_MASK_INPUT | PM_SUBTRACTION_MASK_CONVOLVE))) {
+                        convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank;
+                        convImage->data.F32[y][x] = NAN;
+                        if (inWeight) {
+                            convWeight->data.F32[y][x] = NAN;
                         }
                     }