Changeset 10022 for trunk/psModules/src/detrend/pmMaskBadPixels.c

Timestamp:

Nov 16, 2006, 2:40:10 PM (20 years ago)

Author:

Paul Price

Message:

Adding additional option for identifying bad pixels.

File:

• trunk/psModules/src/detrend/pmMaskBadPixels.c (modified) (2 diffs)

trunk/psModules/src/detrend/pmMaskBadPixels.c

@@ -138 +138 @@
     PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, NULL);
 
-    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics
-    if (!psImageStats(stats, suspects, NULL, 0) || !isfinite(stats->sampleMean) ||
-            !isfinite(stats->sampleStdev)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics.\n");
-        psFree(stats);
+    float limit = NAN;                  // Limit for masking
+    if (thresh > 0) {
+        psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV); // Statistics
+        stats->clipSigma = 5.0;
+        stats->clipIter = 1;
+        if (!psImageStats(stats, suspects, NULL, 0) || !isfinite(stats->clippedStdev)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics.\n");
+            psFree(stats);
+            return NULL;
+        }
+        limit = thresh * stats->clippedStdev;
+        psFree(stats);
+    } else if (thresh < 0) {
+        psStats *stats = psStatsAlloc(PS_STAT_MAX); // Statistics
+        if (!psImageStats(stats, suspects, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics.\n");
+            psFree(stats);
+            return NULL;
+        }
+        psHistogram *histo = psHistogramAlloc(-0.5, stats->max + 0.5, stats->max + 1);
+        psFree(stats);
+        if (!psImageHistogram(histo, suspects, NULL, 0)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram.\n");
+            psFree(histo);
+            return NULL;
+        }
+
+        float target = histo->nums->data.F32[0] * powf(10.0, thresh); // Level to which histogram must fall
+        int index;
+        for (index = 0; histo->nums->data.F32[index] > target && index < histo->nums->n; index++)
+            ;
+        limit = histo->bounds->data.F32[index];
+    } else {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "Threshold (%f) must be a positive or negative real number.\n", thresh);
         return NULL;
     }
-
-    float mean = stats->robustMedian;     // Mean value
-    float stdev = stats->robustStdev;   // Standard deviation
-    psFree(stats);
 
     psImage *badpix = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); // Bad pixel mask
@@ -164 +190 @@
         psFree(stats);
         psFree(histo);
-        printf("Threshold: %f\n", mean + thresh * stdev);
+        printf("Threshold: %f\n", limit);
     }
 
     for (int y = 0; y < suspects->numRows; y++) {
         for (int x = 0; x < suspects->numCols; x++) {
-            if (suspects->data.S32[y][x] - mean >= thresh * stdev) {
+            if (suspects->data.S32[y][x] >= limit) {
                 badpix->data.PS_TYPE_MASK_DATA[y][x] = maskVal;
             }