Changeset 41806 for branches/eam_branches/ipp-dev-20210817/psModules

Timestamp:

Sep 8, 2021, 10:03:37 AM (5 years ago)

Author:

eugene

Message:

use poisson stats of sky background (median) to clip data values; require more valid pixels to fit the trend

File:

• branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c (modified) (4 diffs)

branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c

@@ -171 +171 @@
         return true;
     }
+
+# if (0) 
     float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
     float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data
     float background = stats->robustMedian;
+# else
+
+    // the signal we are looking for is a small variation on top of the background.  if
+    // the background is uniform with only read noise + sky noise, then the pixel-to-pixel
+    // stdev should only be due to known noise sources and predictable.  If the
+    // pixel-to-pixel variations are from other features, then those variations will
+    // probably dominate the row-by-row bias variations.
+
+    // instead of using the image pixel statistics to measure the stdev, lets assume only
+    // dark noise plus poisson sky noise.  we are not carrying in the read noise, but it is
+    // fairly modest for GPC1 (~10 DN)
+
+    // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have
+    // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17
+
+    // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2)
+    // so smaller by a factor of 1.4 than what we predict, which is not very large
+
+    # define READNOISE 10
+    float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    float lower = stats->robustMedian - thresh * sigma; // Lower bound for data
+    float upper = stats->robustMedian + thresh * sigma; // Upper bound for data
+    float background = stats->robustMedian;
+# endif    
+    
+    // the expected amplitude of the 2nd order term is 
+    // if the predicted sky poisson signal is larger than 
+
+    // I want to add a constraint to the fit so the amplitude or coeffs are |value| < x
+
     psFree(stats);
     psFree(rng);
@@ -212 +244 @@
     psVectorInit(yaxisMask, 0);
 #endif
+
+    // we really need more than order + 1 points (= 4).
+    // this should be tunable, but let's try 5 - 10%
+    int validNmin = numCols * 0.1;
+
     for (int y = 0; y < numRows; y++) {
         psVectorInit(clipMask, 0);
         data = psImageRow(data, image, y);
         int num = 0;                    // Number of good pixels
+
+        // if the unmasked pixels only span a small range in x then we cannot fit the
+        // 2nd order polynomial variations very well.  Require a minimum fractional range
+        float validXmin = +1;
+        float validXmax = -1;
+
         for (int x = 0; x < numCols; x++) {
             if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) ||
@@ -223 +266 @@
                 clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0;
                 num++;
+                validXmin = PS_MIN(indices->data.F32[x], validXmin);
+                validXmax = PS_MAX(indices->data.F32[x], validXmax);
             }
         }
-        if (num < order + 1) {
+
+        if (num < validNmin) {
             // Not enough points to fit
             patternMaskRow(ro, y, maskBad);
@@ -234 +280 @@
             continue;
         }
+        // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise?
         if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) {
             psWarning("Unable to fit polynomial to row %d", y);