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

Timestamp:

Sep 27, 2021, 11:11:33 AM (5 years ago)

Author:

eugene

Message:

added binned vs unbinned version of pattern fitting; stricter cuts on thresholding

File:

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

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

@@ -137 +137 @@
 
 // Comparison and swap functions for sorting values directly
-#define SORT_COMPARE(A,B) (value[A] < value[B])
-#define SORT_SWAP(A,B) { \
+#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B])
+#define SORT_SWAP(TYPE,A,B) {                   \
     if (A != B) { \
-        float temp = value[A]; \
-        value[A] = value[B]; \
-        value[B] = temp; \
+        TYPE temp = sampleArray[A];                     \
+        sampleArray[A] = sampleArray[B]; \
+        sampleArray[B] = temp; \
     } \
 }
@@ -150 +150 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                          psStatsOptions clipMean, psStatsOptions clipStdev,
+                          psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
+                        psStatsOptions clipMean, psStatsOptions clipStdev,
+                        psImageMaskType maskVal, psImageMaskType maskBad);
+
+
+bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
+                  psStatsOptions clipMean, psStatsOptions clipStdev,
+                  psImageMaskType maskVal, psImageMaskType maskBad) {
+
+//  bool status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+    bool status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+    return status;
+}
+
 // USE_BACKGROUND_STDEV: if true, the analysis will use the measured robust stdev to clip the out-of-range pixles
 // if false, the stdev will be estimated based on the Poisson statistics of the median (sky level).
 # define USE_BACKGROUND_STDEV 0
 
-bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh,
+bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh,
                   psStatsOptions clipMean, psStatsOptions clipStdev,
                   psImageMaskType maskVal, psImageMaskType maskBad)
@@ -233 +252 @@
 
     psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics
-    // XXX clip->clipIter = iter;
-    clip->clipIter = 1; // XXX skip iteration for a test
+    clip->clipIter = iter;
     clip->clipSigma = rej;
     psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping
@@ -440 +458 @@
 }
 
-// bin by Npix in the x-direction to reduce the number of calculations needed to measure
+# define NPIX 31    
+
+// bin by NPIX in the x-direction to reduce the number of calculations needed to measure
 // the pattern
 bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh,
@@ -474 +494 @@
     }
 
+    // if USE_BACKGROUND_STDEV is TRUE, the observed standard deviation is used to set the
+    // thresholds.  this is going to be an overestimate if there is any structure in the
+    // image.  If FALSE, the thresholds are set based on poisson stats for the background
+    // level.  We assume the gain is 1, so this is an overestimate if the gain is > 1
+
 # if (USE_BACKGROUND_STDEV) 
     float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data
@@ -497 +522 @@
     # 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 delta = PS_MIN (thresh * sigma, 40);
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
     float background = stats->robustMedian;
 # endif    
@@ -513 +539 @@
     psFree(rng);
 
-# define NPIX 15    
-
-    // Indices are distributed [-1:1] [-1 = 0, +1 = numCols = indices[nSamples]
+    // the vector 'indices' maps the x-coordinate to a range [-1:1].  the element number (i) of indices
+    // related to the x-coordinate (column number) by x = (i + 0.5) * NPIX
+
     int nSamples = numCols / NPIX;
-    psVector *indices = psVectorAlloc(nSamples, PS_TYPE_F32); // Indices for fitting
-    psVector *fitData = psVectorAlloc(nSAmples, PS_TYPE_F32); // Data to fit
-
-    // indices elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX
+
+    psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fit solutions
+    psVector *xFit = psVectorAlloc(nSamples, PS_TYPE_F32); // x-coordinate for fitting
+    psVector *yFit = psVectorAlloc(nSamples, PS_TYPE_F32); // flux values for fitting
+
+    // xFit elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX
 
     float norm = 2.0 / (float)numCols;  // Normalisation for indices
     for (int sample = 0; sample < nSamples; sample ++) {
         int x = (sample + 0.5)*NPIX;
-        indices->data.F32[sample] = x * norm - 1.0;
+        xFit->data.F32[sample] = x * norm - 1.0;
+    }
+    for (int x = 0; x < numCols; x ++) {
+        indices->data.F32[x] = x * norm - 1.0;
     }
 
     psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics
-    // XXX clip->clipIter = iter;
-    clip->clipIter = 1; // XXX skip iteration for a test
+    clip->clipIter = iter;
     clip->clipSigma = rej;
-    psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping
+    psVector *clipMask = psVectorAlloc(nSamples, 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
@@ -539 +569 @@
     psImageInit(corr, NAN);
 
-#ifdef PATTERN_ROW_BKG_FIX
     // CZW: 2011-11-30
     // Define the vectors to hold the "x" and "y" slope trends.
@@ -556 +585 @@
     psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term
     psVectorInit(yaxisMask, 0);
-#endif
-
-    // we really need more than order + 1 points (= 4).
+
+    // validNmin is the minimum number of samples needed to measure the trend.
     // this should be tunable, but let's try 5 - 10%
-    int validNmin = numCols * 0.1;
+    int validNmin = PS_MAX (nSamples * 0.1, order + 2);
 
     for (int y = 0; y < numRows; y++) {
@@ -572 +600 @@
         float validXmax = -1;
 
-        // XXX bookkeeping might be easier if this loop is over elements of 'indices'
-        // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY)
-        // (x % 3) ||
         for (int sample = 0; sample < nSamples; sample ++) {
 
-            // store valid samples in the array
+            // store valid samples in the array to be sorted
             float sampleArray[NPIX];
             int seq = 0;
             for (int j = 0; j < NPIX; j++) {
-                int pix = sample  * NPIX + j;
+                int pix = sample  * NPIX + j; // real pixel elements in x-dir
                 if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][pix] & maskVal)) continue;
                 if (data->data.F32[pix] < lower || data->data.F32[pix] > upper) continue;
-                sampleArray[seq] = data->data.F32[pix];
+                sampleArray[seq] = data->data.F32[pix]; // store the value to be sorted
                 seq ++;
             } 
             if (seq < 1) {
                 clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF;
-            } else {
-                clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0;
-                num++;
-                validXmin = PS_MIN(indices->data.F32[sample], validXmin);
-                validXmax = PS_MAX(indices->data.F32[sample], validXmax);
-            }
-
-            PSSORT (seq, sampleArray
-
-        }
-
-        // XXX how much time is spent in the fitting
+                yFit->data.F32[sample] = NAN;
+                continue;
+            }
+            // note that we are treating the x-coordinate as the center
+            // of the binned pixel group, even if some or most pixels have
+            // been masked.  compared to the amplitude of the slope, this
+            // error is small
+            clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0;
+            validXmin = PS_MIN(xFit->data.F32[sample], validXmin);
+            validXmax = PS_MAX(xFit->data.F32[sample], validXmax);
+            num++;
+
+            // PSSORT operates on sampleArray (see define of macro SORT_SWAP above)
+            PSSORT (seq, SORT_COMPARE, SORT_SWAP, float);
+
+            int midPt = 0.5 * seq;
+            float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]);
+            yFit->data.F32[sample] = medValue;
+        }
+
+        // XXX how much time is spent in the fitting? to test: make this always true
         if (num < validNmin) {
             // Not enough points to fit
@@ -609 +643 @@
         }
         // 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)) {
+        if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, yFit, NULL, xFit)) {
             psWarning("Unable to fit polynomial to row %d", y);
             psErrorClear();
@@ -628 +662 @@
             psErrorClear();
             patternMaskRow(ro, y, maskBad);
-#ifdef PATTERN_ROW_BKG_FIX
             yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
-#endif
             continue;
         }
@@ -641 +673 @@
     }
 
-#ifdef PATTERN_ROW_BKG_FIX
     // Put the global trends back that were removed by the PATTERN.ROW correction.
     // Set up the indices for the polynomial
@@ -665 +696 @@
         corr->data.F64[y][0]  -= background;
       }
-    }
-    else {
+    } else {
       psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices);
       if (!solution) {
@@ -679 +709 @@
           corr->data.F64[y][0]  -= background;
         }
-      }
-      else {
+      } else {
         for (int y = 0; y < numRows; y++) {
           for (int x = 0; x < numCols; x++) {
@@ -699 +728 @@
       psWarning("Unable to fit polynomial to x-axis trend");
       psErrorClear();
-    }
-    else {
+    } else {
       psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices);
       if (!solution) {
         psWarning("Unable to evaluate polynomial");
         psErrorClear();
-      }
-      else {
+      } else {
         for (int y = 0; y < numRows; y++) {
           for (int x = 0; x < numCols; x++) {
@@ -723 +750 @@
     psFree(yaxisData);
     psFree(xaxisData);
-    // End PATTERN_ROW_BKG_FIX global trend replacement
-#endif 
     
     psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE,
@@ -731 +756 @@
 
     psFree(indices);
+    psFree(xFit);
+    psFree(yFit);
     psFree(clip);
     psFree(clipMask);