Changeset 26739 for branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c

Timestamp:

Jan 29, 2010, 6:02:38 PM (16 years ago)

Author:

Paul Price

Message:

Fairly extensive change to make subtractions more stable on images with limited lit area in the presence of spatial variation. Took the polynomials which were scaled over the entire image and scale them only over the lit area. Removed some fprintf statements. Fixed the problem causing the background term to always come out 1.0.

File:

• branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c (modified) (5 diffs)

branches/eam_branches/20091201/ppSub/src/ppSubThreshold.c

@@ -27 +27 @@
     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 +38 @@
     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)) {
+    if (!psImageBackground(stats, NULL, subImage, subMask, maskIgnore, rng)) {
         psError(PS_ERR_UNKNOWN, false, "Unable to determine threshold.");
         psFree(rng);
@@ -49 +52 @@
     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;
@@ -94 +100 @@
         return false;
     }
-    if (!lowThreshold(in, thresh, maskVal, maskThresh, "input convolved image")) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
-        return false;
-    }
 
     pmReadout *ref = pmFPAfileThisReadout(config->files, view, "PPSUB.REF.CONV"); // Reference image
@@ -104 +106 @@
         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(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+            return false;
+        }
+        if (!lowThreshold(ref, thresh, maskVal, maskThresh, region, "reference convolved image")) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to threshold input image.");
+            return false;
+        }
     }
+    psFree(regIter);
 
     psFree(view);