Index: trunk/psModules/src/imcombine/pmSubtraction.c
===================================================================
--- trunk/psModules/src/imcombine/pmSubtraction.c	(revision 14870)
+++ trunk/psModules/src/imcombine/pmSubtraction.c	(revision 15305)
@@ -4,6 +4,6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.62 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-17 21:40:22 $
+ *  @version $Revision: 1.63 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-10-12 23:00:36 $
  *
  *  Copyright 2004-2007 Institute for Astronomy, University of Hawaii
@@ -459,5 +459,5 @@
 
 psImage *pmSubtractionMask(const psImage *refMask, const psImage *inMask, psMaskType maskVal,
-                           int size, int footprint)
+                           int size, int footprint, float badFrac)
 {
     PS_ASSERT_IMAGE_NON_NULL(refMask, NULL);
@@ -470,4 +470,8 @@
     PS_ASSERT_INT_NONNEGATIVE(size, NULL);
     PS_ASSERT_INT_NONNEGATIVE(footprint, NULL);
+    if (isfinite(badFrac)) {
+        PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL);
+        PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL);
+    }
 
     // Size of the images
@@ -475,15 +479,36 @@
     int numRows = refMask->numRows;
 
+    // Dereference inputs for convenience
+    psMaskType **inData = NULL;
+    if (inMask) {
+        inData = inMask->data.PS_TYPE_MASK_DATA;
+    }
+    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;
+
+    // First, a pass through to determine the fraction of bad pixels
+    if (isfinite(badFrac) && badFrac != 1.0) {
+        int numBad = 0;                 // Number of bad pixels
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if (inData && inData[y][x] & maskVal) {
+                    numBad++;
+                    continue;
+                }
+                if (refData[y][x] & maskVal) {
+                    numBad++;
+                }
+            }
+        }
+        if (numBad > badFrac * numCols * numRows) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Fraction of bad pixels (%d/%d) exceeds limit (%f)\n",
+                    numBad, numCols * numRows, badFrac);
+            return NULL;
+        }
+    }
+
     // 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 = NULL;
-    if (inMask) {
-        inData = inMask->data.PS_TYPE_MASK_DATA;
-    }
-    psMaskType **refData = refMask->data.PS_TYPE_MASK_DATA;
+    psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience
 
     // Block out a border around the edge of the image
@@ -510,4 +535,9 @@
         }
     }
+
+    // XXX Could do something smarter here --- we will get images that are predominantly masked (where the
+    // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad
+    // pixel is wasting time.  As a first cut, I've put in a check on the fraction of bad pixels, but we could
+    // imagine looking for the edge of big regions and convolving just at the edge.
 
     // Mask the bad pixels, and around them