Index: branches/tap_branches/ppSub/src/ppSubThreshold.c
===================================================================
--- branches/tap_branches/ppSub/src/ppSubThreshold.c	(revision 25900)
+++ branches/tap_branches/ppSub/src/ppSubThreshold.c	(revision 27838)
@@ -27,4 +27,5 @@
     psImageMaskType maskIgnore,         // Ignore pixels with this mask
     psImageMaskType maskThresh,         // Give pixels this mask if below threshold
+    psRegion *region,                   // Region of interest
     const char *description             // Description of image
     )
@@ -37,10 +38,12 @@
     psImage *image = ro->image;         // Image
     psImage *mask = ro->mask;           // Mask
-    int numCols = ro->image->numCols, numRows = ro->image->numRows; // Size of image
+
+    psImage *subImage = psImageSubset(image, *region); // Image with region of interest
+    psImage *subMask = psImageSubset(mask, *region);  // Maks with region of interest
 
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);                               // Random number generator
-    if (!psImageBackground(stats, NULL, image, mask, maskIgnore, rng)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold.");
+    if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) {
+        psError(PPSUB_ERR_DATA, false, "Unable to determine threshold.");
         psFree(rng);
         psFree(stats);
@@ -49,10 +52,13 @@
     psFree(rng);
 
+    psFree(subImage);
+    psFree(subMask);
+
     float threshold = stats->robustMedian - thresh * stats->robustStdev; // Threshold below which to clip
     psFree(stats);
     psLogMsg("ppSub", PS_LOG_INFO, "Masking pixels below %f in %s", threshold, description);
 
-    for (int y = 0; y < numRows; y++) {
-        for (int x = 0; x < numCols; x++) {
+    for (int y = region->y0; y < region->y1; y++) {
+        for (int x = region->x0; x < region->x1; x++) {
             if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskIgnore) {
                 continue;
@@ -91,9 +97,5 @@
     pmReadout *in = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT.CONV"); // Input image
     if (!in) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout.");
-        return false;
-    }
-    if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+        psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout.");
         return false;
     }
@@ -101,11 +103,24 @@
     pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image
     if (!ref) {
-        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find readout.");
+        psError(PPSUB_ERR_UNKNOWN, false, "Unable to find readout.");
         return false;
     }
-    if (!lowThreshold(ref, thresh, maskVal, maskThresh, "reference convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
-        return false;
+
+    psMetadataIterator *regIter = psMetadataIteratorAlloc(in->analysis, PS_LIST_HEAD,
+                                                          "^" PM_SUBTRACTION_ANALYSIS_REGION "$");
+    psMetadataItem *regItem;        // Item with region
+    while ((regItem = psMetadataGetAndIncrement(regIter))) {
+        psAssert(regItem->type == PS_DATA_REGION && regItem->data.V, "Expect region type");
+        psRegion *region = regItem->data.V; // Region of interest
+        if (!lowThreshold(in, thresh, maskVal, maskThresh, region, "input convolved image")) {
+            psError(psErrorCodeLast(), false, "Unable to threshold input image.");
+            return false;
+        }
+        if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) {
+            psError(psErrorCodeLast(), false, "Unable to threshold input image.");
+            return false;
+        }
     }
+    psFree(regIter);
 
     psFree(view);