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

Timestamp:

Oct 23, 2021, 1:48:02 PM (5 years ago)

Author:

eugene

Message:

adding PATTERN.ROW.AMP I/O functions, clean up pmPatternRow code

File:

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

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

@@ -160 +160 @@
 
 
+// XXX allow user choice of binned vs unbinned analysis?
 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;
+  bool status = false;
+  if (true) {
+    status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad);
+  } else {
+    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).
+// 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
 
@@ -184 +189 @@
     PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
 
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
     psImage *image = ro->image;         // Image to correct
     psImage *mask = ro->mask;           // Mask for image
@@ -191 +203 @@
     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; masking entire readout.");
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
         psErrorClear();
         psFree(stats);
         psFree(rng);
-        psImageInit(image, NAN);
-        if (mask) {
-            psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK));
-        }
-        if (ro->variance) {
-            psImageInit(image, NAN);
-        }
         return true;
     }
@@ -226 +231 @@
     // so smaller by a factor of 1.4 than what we predict, which is not very large
 
-    # define READNOISE 10
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
+    // XXX retrieve noise and gain from 'concepts' if possible
+# define READNOISE 10 /* arbitrary number */
     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, 2*nominalAmplitude); 
+    float lower = stats->robustMedian - delta; // Lower bound for data
+    float upper = stats->robustMedian + delta; // Upper bound for data
     float background = stats->robustMedian;
+
+    // if the noise from the background is too large 
+    float significance = nominalAmplitude / sigma;
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      return true; 
+    }
+    fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
 # endif    
     
-    pmCell *cell = ro->parent;
-    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
-    fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName,lower, background, upper);
-
-
-    // XXX add code to skip fit if background stdev is large compared to expected signal
-    // XXX add a constraint to the fit so the |amplitude| < x (x ~ 20 - 25)
-
     psFree(stats);
     psFree(rng);
@@ -458 +473 @@
 }
 
-# define NPIX 31    
+# define NPIX 15
 
 // bin by NPIX in the x-direction to reduce the number of calculations needed to measure
@@ -473 +488 @@
     PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false);
 
+    bool mdok;                          // Status of MD lookup
+
+    pmCell *cell = ro->parent;
+    pmChip *chip = cell->parent;
+    const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip
+    const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell
+
     psImage *image = ro->image;         // Image to correct
     psImage *mask = ro->mask;           // Mask for image
     int numCols = image->numCols, numRows = image->numRows; // Size of image
 
-    // use the background and typical amplitude value to choose to correct or not
-    // pmCell *cell = ro->parent;
-
     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; skipping pattern correction.");
+        psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName);
         psErrorClear();
         psFree(stats);
@@ -525 +544 @@
     // so smaller by a factor of 1.4 than what we predict, which is not very large
 
+    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+    float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP");
+    if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number
+    // If we cannot determine the nominal amplitude, we fall-back on the worst case
+
     // retrieve noise and gain from 'concepts' if possible
 # define READNOISE 10 /* arbitrary number */
     float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
-    float delta = PS_MIN (thresh * sigma, 40); // XXX EAM : arbitrary number
+    float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); 
     float lower = stats->robustMedian - delta; // Lower bound for data
     float upper = stats->robustMedian + delta; // Upper bound for data
     float background = stats->robustMedian;
 
-    // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
-    float nominalAmplitude = 20; // XXX EAM : arbitrary number
-
     // if the noise from the background is too large 
     float significance = nominalAmplitude / sigma;
-    if (significance < 0.5) return true; // XXX EAM : arbitrary number
-    
-    pmCell *cell = ro->parent;
-    const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
-    fprintf (stderr, "pattern row background %s: %f - %f - %f\n", cellName, lower, background, upper);
+   
+    // XXX EAM : arbitrary number
+    if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      psFree(stats);
+      psFree(rng);
+      return true; 
+    }
+    psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
 # endif    
 
@@ -643 +669 @@
         }
 
-        // XXX how much time is spent in the fitting? to test: make this always true
+        // If not enough points are valid, skip
         if (num < validNmin) {
-            // Not enough points to fit
-            // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
-            // patternMaskRow(ro, y, maskBad);
             // Ignore this row in our subsequent fits, because the fit failed.
             yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
@@ -656 +679 @@
             psWarning("Unable to fit polynomial to row %d", y);
             psErrorClear();
-            // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
-            // patternMaskRow(ro, y, maskBad);
             // Ignore this row in our subsequent fits, because the fit failed.
             yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
@@ -672 +693 @@
             psWarning("Unable to evaluate polynomial for row %d", y);
             psErrorClear();
-            // EAM 20211011 : we used to mask rows which could not be fit, but this is excessive.  just skip.
-            // patternMaskRow(ro, y, maskBad);
             yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
             continue;