Changeset 41968 for branches/eam_branches/ipp-20211108/psphot/src/psphotFindPeaks.c

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

File:

• branches/eam_branches/ipp-20211108/psphot/src/psphotFindPeaks.c (modified) (2 diffs)

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