Changeset 42088

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/psphot.h (modified) (1 diff)
• src/psphotFindDetections.c (modified) (1 diff)
• src/psphotFindPeaks.c (modified) (2 diffs)
• src/psphotFitSourcesLinear.c (modified) (7 diffs)
• src/psphotSignificanceImage.c (modified) (2 diffs)
• src/psphotSourceFits.c (modified) (2 diffs)

trunk/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);

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

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

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -217 +217 @@
     int Nsat = 0;
 
+    // track number kept at several stages of the analysis
+    int Nstep_0 = 0;
+    int Nstep_1 = 0;
+    int Nstep_2 = 0;
+    int Nstep_3 = 0;
+    int Nstep_4 = 0;
+    int Nstep_5 = 0;
+
     // select the sources which will be used for the fitting analysis
     for (int i = 0; i < sources->n; i++) {
@@ -246 +254 @@
         if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
 
+        Nstep_0++;
+
         // XXX count saturated stars
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
@@ -255 +265 @@
             if (!pmSourceCacheModel (source, maskVal)) continue;
         }
+        Nstep_1++;
 
         // save the original coords
@@ -265 +276 @@
         if (x > AnalysisRegion.x1) continue;
         if (y > AnalysisRegion.y1) continue;
+        Nstep_2++;
 
         // check the integral of the model : is it large enough?
@@ -283 +295 @@
             continue;
         }
+        Nstep_3++;
+
         bool isPSF = false;
         pmModel *model = pmSourceGetModel (&isPSF, source);
@@ -306 +320 @@
         } 
         if (modelSum  < cutModelSum * normFlux) continue;
+        Nstep_4++;
+
         if (maskedSum < cutMaskedSum * normFlux) continue;
+        Nstep_5++;
 
         // clear the 'mark' pixels and remask on the fit aperture 
@@ -322 +339 @@
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "Number of sources kept at each stage: %d %d %d %d %d %d\n", Nstep_0, Nstep_1, Nstep_2, Nstep_3, Nstep_4, Nstep_5);
 
     if (fitSources->n == 0) {

trunk/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);

trunk/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