Changeset 9981

Timestamp:

Nov 14, 2006, 2:27:52 PM (19 years ago)

Author:

Paul Price

Message:

Adding functions to create bad pixel masks from flattened data.

Location:

trunk/psModules/src/detrend

Files:

• pmMaskBadPixels.c (modified) (2 diffs)
• pmMaskBadPixels.h (modified) (2 diffs)

trunk/psModules/src/detrend/pmMaskBadPixels.c

@@ -4 +4 @@
 
 #include <stdio.h>
+#include <math.h>
 #include <pslib.h>
 
@@ -67 +68 @@
     return true;
 }
+
+
+psImage *pmMaskFlagSuspectPixels(psImage *out, const pmReadout *readout, float rej,
+                                 psMaskType maskVal, float frac, psRandom *rng)
+{
+    PS_ASSERT_PTR_NON_NULL(readout, NULL);
+    PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, NULL);
+    PS_ASSERT_FLOAT_WITHIN_RANGE(frac, 0.0, 1.0, NULL);
+    PS_ASSERT_IMAGE_NON_NULL(readout->image, NULL);
+    PS_ASSERT_IMAGE_NON_EMPTY(readout->image, NULL);
+    PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, NULL);
+    if (readout->mask) {
+        PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, NULL);
+        PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, NULL);
+    }
+    if (out) {
+        PS_ASSERT_IMAGE_NON_EMPTY(out, NULL);
+        PS_ASSERT_IMAGE_TYPE(out, PS_TYPE_S32, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, out, NULL);
+    }
+
+    psImage *image = readout->image;    // Image of interest
+    psImage *mask = readout->mask;      // Corresponding mask
+
+    if (rng) {
+        psMemIncrRefCounter(rng);
+    } else {
+        rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    }
+
+    psStats *stats = psImageBackground(image, mask, maskVal, 0.25, 0.75,
+                                       frac * image->numCols * image->numRows, rng); // Image statistics
+    if (!stats || !isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) ||
+            !isfinite(stats->robustLQ)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to measure image statistics.\n");
+        psFree(stats);
+        psFree(rng);
+        return NULL;
+    }
+    psFree(rng);
+
+    float median = stats->robustMedian; // Median value
+    float stdev = 0.74*(stats->robustUQ - stats->robustLQ); // Estimate of the standard deviation
+    psFree(stats);
+
+    if (!out) {
+        out = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S32);
+        psImageInit(out, 0);
+    }
+
+    for (int y = 0; y < image->numRows; y++) {
+        for (int x = 0; x < image->numCols; x++) {
+            if (fabs((image->data.F32[y][x] - median) / stdev) >= rej &&
+                    (!mask || !(mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) {
+                out->data.S32[y][x]++;
+            }
+        }
+    }
+
+    return out;
+}
+
+
+psImage *pmMaskIdentifyBadPixels(const psImage *suspects, float thresh, psMaskType maskVal)
+{
+    PS_ASSERT_IMAGE_NON_NULL(suspects, NULL);
+    PS_ASSERT_IMAGE_NON_EMPTY(suspects, NULL);
+    PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, NULL);
+
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_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);
+        return NULL;
+    }
+
+    float mean = stats->sampleMean;     // Mean value
+    float stdev = stats->sampleStdev;   // Standard deviation
+
+    psImage *badpix = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); // Bad pixel mask
+    psImageInit(badpix, 0);
+
+    if (psTraceGetLevel("psModules.detrend") > 9) {
+        psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics
+        psImageStats(stats, suspects, NULL, 0);
+        psHistogram *histo = psHistogramAlloc(-0.5, stats->max + 0.5, stats->max + 1);
+        psImageHistogram(histo, suspects, NULL, 0);
+        for (int i = 0; i < histo->nums->n; i++) {
+            printf("%f --> %f : %f\n", histo->bounds->data.F32[i], histo->bounds->data.F32[i + 1],
+                   histo->nums->data.F32[i]);
+        }
+        psFree(stats);
+        psFree(histo);
+        printf("Threshold: %f\n", mean + thresh * stdev);
+    }
+
+    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) {
+                badpix->data.PS_TYPE_MASK_DATA[y][x] = maskVal;
+            }
+        }
+    }
+
+    return badpix;
+}

trunk/psModules/src/detrend/pmMaskBadPixels.h

@@ -8 +8 @@
 /// @author Eugene Magnier, IfA
 ///
-/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
-/// @date $Date: 2006-10-17 22:17:38 $
+/// @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+/// @date $Date: 2006-11-15 00:27:52 $
 ///
 /// Copyright 2004 Institute for Astronomy, University of Hawaii
@@ -33 +33 @@
                     );
 
+/// Find pixels outlying from the background, flagging suspect pixels
+///
+/// Pixels more than "rej" standard deviations from the background level (in flat-fielded images) have the
+/// corresponding pixel in the "suspect pixels" image incremented.  After accumulating over a suitable sample
+/// of images, bad pixels should have a high value in the suspect pixels image, allowing them to be
+/// identified.  The suspect pixels image is of type S32.
+psImage *pmMaskFlagSuspectPixels(psImage *out, ///< Suspected bad pixels image, or NULL
+                                 const pmReadout *readout, ///< Readout to inspect
+                                 float rej, ///< Rejection threshold (standard deviations)
+                                 psMaskType maskVal, ///< Mask value for statistics
+                                 float frac, ///< Fraction of pixels to consider
+                                 psRandom *rng ///< Random number generator
+                                );
+
+/// Identify bad pixels from the suspect pixels image
+///
+/// Bad pixels are identified from the suspect pixels image (accumulated over a large number of images).
+/// Pixels marked as suspect in more than "thresh" standard deviations from the mean are identified as bad
+/// pixels (output image).
+psImage *pmMaskIdentifyBadPixels(const psImage *suspects, ///< Accumulated suspect pixels image
+                                 float thresh, ///< Threshold for bad pixel (standard deviations)
+                                 psMaskType maskVal ///< Value to set for bad pixels
+                                );
+
+
+
+
 #endif