Changeset 42088 for trunk/psphot/src/psphotFindPeaks.c

Timestamp:

Feb 28, 2022, 2:40:32 PM (4 years ago)

Author:

eugene

Message:

for bright stars in stacks, the peak in the S/N image can have a central divot if the variance of the central pixels is relatively large compared to the ring of pixels further from the center. the detected peaks are then not well centered. this can be avoided if we detect on the signal image itself (after first determining an appropriate flux value for the desired S/N limit); added an option to use the signal image for the first pass peak detetion step

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src/psphotFindPeaks.c (modified) (2 diffs)

trunk/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");
+            psFree (SNvalues);
+            return NULL;
+        }
+        psFree (SNvalues);
+
+        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");
+            psFree (stats);
+            goto use_significance;
+        }
+        
+        float alt_threshold = stats->sampleMedian;
+        psFree (stats);
+
+        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) {