Changeset 24903

Timestamp:

Jul 22, 2009, 5:13:36 PM (17 years ago)

Author:

Paul Price

Message:

Fixing pattern subtraction so that it ignores pixels away from the median, in order to reduce bad subtractions around bright stars. Changed API for pmPatternRow, added new recipe paramter.

Location:

trunk

Files:

• ippconfig/recipes/ppImage.config (modified) (1 diff)
• ppImage/src/ppImage.h (modified) (1 diff)
• ppImage/src/ppImageDetrendReadout.c (modified) (1 diff)
• ppImage/src/ppImageOptions.c (modified) (2 diffs)
• psModules/src/detrend/pmPattern.c (modified) (5 diffs)
• psModules/src/detrend/pmPattern.h (modified) (1 diff)

trunk/ippconfig/recipes/ppImage.config

@@ -85 +85 @@
 PATTERN.ITER       S32     4               # Rejection iterations
 PATTERN.REJ        F32     2.5             # Rejection threshold
+PATTERN.THRESH     F32     10.0            # Ignore threshold
 PATTERN.MEAN       STR     SAMPLE_MEAN     # Statistic for mean
 PATTERN.STDEV      STR     SAMPLE_STDEV    # Statistic for standard deviation

trunk/ppImage/src/ppImage.h

@@ -98 +98 @@
     int patternIter;                    // Clipping iterations
     float patternRej;                   // Clipping threshold
+    float patternThresh;                // Ignore threshold
     psStatsOptions patternMean;         // Statistic for mean
     psStatsOptions patternStdev;        // Statistic for stdev

trunk/ppImage/src/ppImageDetrendReadout.c

@@ -127 +127 @@
     if (options->doPattern) {
         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej,
-                          options->patternMean, options->patternStdev, options->darkMask)) {
+                          options->patternThresh, options->patternMean, options->patternStdev,
+                          options->maskValue, options->darkMask)) {
             psFree(detview);
             return false;

trunk/ppImage/src/ppImageOptions.c

@@ -93 +93 @@
     options->patternIter     = 0;       // Clipping iterations
     options->patternRej      = NAN;     // Clipping rejection threshold
+    options->patternThresh   = NAN;     // Threshold for ignoring pixels
     options->patternMean     = PS_STAT_NONE; // Statistic for mean
     options->patternStdev    = PS_STAT_NONE; // Statistic for standard deviation
@@ -312 +313 @@
     options->patternIter = psMetadataLookupS32(NULL, recipe, "PATTERN.ITER");
     options->patternRej = psMetadataLookupF32(NULL, recipe, "PATTERN.REJ");
+    options->patternThresh = psMetadataLookupF32(NULL, recipe, "PATTERN.THRESH");
     options->patternMean = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "PATTERN.MEAN"));
     options->patternStdev = psStatsOptionFromString(psMetadataLookupStr(NULL, recipe, "PATTERN.STDEV"));

trunk/psModules/src/detrend/pmPattern.c

@@ -30 +30 @@
 }
 
-bool pmPatternRow(pmReadout *ro, int order, int iter, float rej,
+bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
                   psStatsOptions clipMean, psStatsOptions clipStdev,
-                  psImageMaskType maskBad)
+                  psImageMaskType maskVal, psImageMaskType maskBad)
 {
     PM_ASSERT_READOUT_NON_NULL(ro, false);
@@ -39 +39 @@
     PS_ASSERT_INT_NONNEGATIVE(iter, false);
     PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false);
+    PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
 
     psImage *image = ro->image;         // Image to correct
+    psImage *mask = ro->mask;           // Mask for image
     int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+    if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) {
+        psWarning("Unable to calculate statistics on readout.");
+        psFree(stats);
+        psFree(rng);
+        return false;
+    }
+    float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
+    float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
+    psFree(stats);
+    psFree(rng);
 
     // Indices are distributed [-1:1)
@@ -53 +68 @@
     clip->clipIter = iter;
     clip->clipSigma = rej;
-    psVector *mask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for fitting
+    psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping
     psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit
     psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit
 
     for (int y = 0; y < numRows; y++) {
-        psVectorInit(mask, 0);
+        psVectorInit(clipMask, 0);
         data = psImageRow(data, image, y);
         int num = 0;                    // Number of good pixels
         for (int x = 0; x < numCols; x++) {
-            if (ro->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] > 0) {
-                mask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
+            if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
+                data->data.F32[x] < lower || data->data.F32[x] > upper) {
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF;
             } else {
-                mask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0;
+                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0;
                 num++;
             }
@@ -74 +90 @@
             continue;
         }
-        if (!psVectorClipFitPolynomial1D(poly, clip, mask, 0xFF, data, NULL, indices)) {
+        if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) {
             psWarning("Unable to fit polynomial to row %d", y);
             psErrorClear();
@@ -96 +112 @@
     psFree(indices);
     psFree(clip);
-    psFree(mask);
+    psFree(clipMask);
     psFree(poly);
     psFree(data);

trunk/psModules/src/detrend/pmPattern.h

@@ -24 +24 @@
     int iter,                           ///< Number of clipping iterations
     float rej,                          ///< Rejection threshold for clipping
+    float thresh,                       ///< Threshold for ignoring pixels
     psStatsOptions clipMean,            ///< Statistic to use for mean
     psStatsOptions clipStdev,           ///< Statistic to use for standard deviation
+    psImageMaskType maskVal,            ///< Mask value to use
     psImageMaskType maskBad             ///< Mask value to give bad pixels
     );