Changeset 24903 for trunk/psModules/src/detrend/pmPattern.c

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.

File:

• trunk/psModules/src/detrend/pmPattern.c (modified) (5 diffs)

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);