Changeset 41968 for branches/eam_branches/ipp-20211108

Timestamp:

Dec 25, 2021, 2:22:11 PM (5 years ago)

Author:

eugene

Message:

allow option to use signal image on first pass for detections

Location:

branches/eam_branches/ipp-20211108/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotFindDetections.c (modified) (1 diff)
• psphotFindPeaks.c (modified) (2 diffs)
• psphotSignificanceImage.c (modified) (2 diffs)
• psphotSourceFits.c (modified) (2 diffs)

branches/eam_branches/ipp-20211108/psphot/src/psphot.h

@@ -216 +216 @@
 // used by psphotFindDetections
 pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
-psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks);
+psArray        *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks, bool firstPass);
 bool            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
 psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);

branches/eam_branches/ipp-20211108/psphot/src/psphotFindDetections.c

@@ -112 +112 @@
     // detect the peaks in the significance image
     int totalPeaks = 0;
-    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks);
+    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX, &totalPeaks, firstPass);
     psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
     if (!detections->peaks) {

branches/eam_branches/ipp-20211108/psphot/src/psphotFindPeaks.c

@@ -4 +4 @@
 // image must be constructed to represent (S/N)^2.  If nMax is non-zero, only return a maximum
 // of nMax peaks
-psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks) {
+psArray *psphotFindPeaks (pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax, int *totalPeaks, bool firstPass) {
 
     bool status = false;
@@ -10 +10 @@
     psTimerStart ("psphot.peaks");
 
-    // find the peaks in the smoothed image
-    // NOTE : significance->variance actually carries the detection S/N image
-    psArray *peaks = pmPeaksInImage (significance->variance, threshold);
-    if (peaks == NULL) {
-        // we only get a NULL peaks array due to a programming or config error. 
-        // this will result in a failure.
-        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
-        return NULL;
+    psArray *peaks = NULL;
+
+    bool useSignalImage = psMetadataLookupF32 (&status, recipe, "PEAKS_USE_SIGNAL_IMAGE"); PS_ASSERT (status, NULL);
+
+    if (firstPass && useSignalImage) {
+        // find the approximate smoothed signal image that corresponds to the desired threshold
+
+        // find all pixels within 25% of the threshold:
+        float minThresh = threshold * 0.80;
+        float maxThresh = threshold * 1.25;
+
+        psImage *smooth_sn = significance->variance;
+        psImage *smooth_im = significance->image;
+
+        psVector *SNvalues = psVectorAllocEmpty (1000, PS_TYPE_F32);
+        for (int iy = 0; iy < smooth_im->numRows; iy++) {
+            for (int ix = 0; ix < smooth_im->numCols; ix++) {
+                // select all valid pixels within the S/N range
+                if (!isfinite(smooth_sn->data.F32[iy][ix])) continue;
+                if (!isfinite(smooth_im->data.F32[iy][ix])) continue;
+                if (smooth_sn->data.F32[iy][ix] < minThresh) continue;
+                if (smooth_sn->data.F32[iy][ix] > maxThresh) continue;
+                psVectorAppend (SNvalues, smooth_im->data.F32[iy][ix]);
+            }
+        }
+
+        // what is the median of the selected values?
+        psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+        if (!psVectorStats(stats, SNvalues, NULL, NULL, 0)) {
+            // psVectorStats will only fail due to a programming error (e.g., invalid vector type)
+            psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
+            return NULL;
+        }
+        if (!isfinite(stats->sampleMedian)) {
+            // could not determine relationship between SN threshold and image values
+            // XXX fall-back could be the standard analysis above.
+            psLogMsg ("psphot", PS_LOG_INFO, "failure to map SN to image values");
+            goto use_significance;
+        }
+        
+        float alt_threshold = stats->sampleMedian;
+        peaks = pmPeaksInImage (significance->image, alt_threshold);
+        if (peaks == NULL) {
+            // we only get a NULL peaks array due to a programming or config error. 
+            // this will result in a failure.
+            psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
+            return NULL;
+        }
     }
+
+use_significance:
+    if (!peaks) {
+        // peaks is NULL at this point if we did not use signal image (by choice or by failure).
+        // find the peaks in the smoothed image
+        // NOTE : significance->variance actually carries the detection S/N image
+        peaks = pmPeaksInImage (significance->variance, threshold);
+        if (peaks == NULL) {
+            // we only get a NULL peaks array due to a programming or config error. 
+            // this will result in a failure.
+            psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
+            return NULL;
+        }
+    }
+
     // return the total number of peaks found before the nMax limit is applied
     if (totalPeaks) {

branches/eam_branches/ipp-20211108/psphot/src/psphotSignificanceImage.c

@@ -105 +105 @@
               //              float v2 = value + PS_SQR(value/100.0);
               // CZW 2013-06-20: I don't think this hack was helping.
+              // EAM 2021-12-21: see note below about divots in the significance image
               float v2 = value;
               smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
@@ -122 +123 @@
     }
     psImageConvolveSetThreads(oldThreads);
+
+    // We now have the significance image and the signal image.  In some cases (e.g.,
+    // stacks), the variance on pixels in the cores of stars is elevated compared to pure
+    // poisson statistics.  In this case, especially at high signal levels, the ratio of
+    // signal / noise in the core of the star can be lower than the surrounding ring of
+    // pixels.  This results in a divot in the center of the star in the significance
+    // image.  the apparent peak of the significance is then not centered on the star and
+    // chaos ensues.  A possible fix is to use the signal image for signficance for
+    // the high S/N detection pass.
 
     pmReadout *significanceRO = pmReadoutAlloc(NULL);

branches/eam_branches/ipp-20211108/psphot/src/psphotSourceFits.c

@@ -499 +499 @@
 
 # if (PS_TRACE_ON) 
-    if (psTraceGetLevel ("psphot") >= 6) {
+    if ( (source->peak->xf > 5500) &&
+         (source->peak->yf < 500) &&
+         psTraceGetLevel ("psphot") >= 6) {
 
       // Moments-based shapes parameters
@@ -546 +548 @@
 
 # if (PS_TRACE_ON) 
-    if (psTraceGetLevel ("psphot") >= 5) {
-      if (psTraceGetLevel ("psphot") >= 6) {
-        fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
-      }
+    if ((source->peak->xf > 5500) &&
+        (source->peak->yf < 500) &&
+        psTraceGetLevel ("psphot") >= 6) {
+      fprintf (stderr, "chisq: %f, nIter: %d, radius: %f, npix: %d\n", model->chisqNorm, model->nIter, model->fitRadius, model->nPix);
       fprintf (stderr, "--- fitted values ---\n");
       for (int i = 0; i < model->params->n; i++) {
         fprintf (stderr, "par %d: %f\n", i, model->params->data.F32[i]);
       }
-      psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
-    }
+    }
+    psTraceSetLevel("psLib.math.psMinimizeLMChi2", 0);
 # endif