Changeset 31154 for trunk/psphot/src/psphotCullPeaks.c

Timestamp:

Apr 4, 2011, 1:12:39 PM (15 years ago)

Author:

eugene

Message:

use the smoothed image for footprint culling; use a more stringent culling for saturated stars; more tweaks to get sat stars to be fitted; various updates to psphotStack; unify psphotImageLoop varients; be a bit careful about number of stars in a PSF clump

Location:

trunk/psphot

Files:

• . (modified) (2 props)
• src/psphotCullPeaks.c (modified) (2 diffs)

trunk/psphot

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM

trunk/psphot/src/psphotCullPeaks.c

@@ -6 +6 @@
 psErrorCode
 psphotCullPeaks(const pmReadout *readout,   // the image wherein lives the footprint
+                const pmReadout *signifR, // smoothed image and significance
                 const psMetadata *recipe,
                 psArray *footprints) {  // array of pmFootprints
@@ -26 +27 @@
     psAssert (status, "cannot find SKY_STDEV in readout->analysis");
 
-    const float min_threshold = nsigma_min*skyStdev;
+    // the sky has been smoothed, so we need to scale down the raw sky stdev by this amount:
+    float scaleFactor = psMetadataLookupF32(&status, readout->analysis, "SIGNIFICANCE_SCALE_FACTOR");
+    if (!status) scaleFactor = 1.0;
+
+    const float MIN_THRESHOLD = nsigma_min*skyStdev / sqrt(scaleFactor);
     
+    // for saturated stars, we should be somewhat more agressive about culling: instead of
+    // letting the height of the intervening col (saddle point) be set by the S/N of the image
+    // pixel, it should be set to a fraction of the saturation value.  example for a
+    // near-saturation pixel:
+    // flux = 40k, sigma = 200
+    // nsigma_delta = 4.0, nsigma_min = 1.0, fPadding = 0.01, 
+    // fStdev = 0.005, stdev_pad = 0.011*40k = 400
+    // threshold = 40k - 4*400 = 38400 
+    // this gives too tight of a tolerance on the bright stars
+    // in practice, we should make the threshold much lower.  
+    // below I am using 0.05 * saturation (eg, 2000 DN above sky in a GPC1 image)
+
+    float SATURATION = NAN;
+
+    // do not completely trust the values in the header...
+    float CELL_SATURATION = psMetadataLookupF32 (&status, readout->parent->concepts, "CELL.SATURATION");
+    float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION");
+    if (!status || !isfinite(MIN_SATURATION)) {
+        MIN_SATURATION = 40000.0;
+    }
+    if (!isfinite(CELL_SATURATION)) {
+        SATURATION = MIN_SATURATION;
+    } else {
+        SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION);
+    }
+    float SAT_TEST_LEVEL = 0.50*SATURATION;
+    float SAT_THRESHOLD  = 0.05*SATURATION;
+
+# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
+    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the smoothed image");
+# else
+    psLogMsg ("psphot", PS_LOG_INFO, "Culling peaks from footprints using the raw (unsmoothed) image");
+# endif
+
     for (int i = 0; i < footprints->n; i++) {
-        // if (i % 50 == 0) fprintf (stderr, "cull %d\n", i);
         pmFootprint *fp = footprints->data[i];
+        if (fp->peaks == NULL) continue;
+        if (fp->peaks->n == 0) continue;
+
         if (fp->npix > 30000) {
             fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix);
         }
         psTimerStart ("psphot.cull.footprints");
-        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) {
+
+        pmPeak *brightPeak = fp->peaks->data[0];
+        float max_threshold = SAT_TEST_LEVEL;
+        if (brightPeak->rawFlux > SAT_TEST_LEVEL) {
+            max_threshold = SAT_THRESHOLD;
+            brightPeak->type = PM_PEAK_SUSPECT_SATURATION;
+        }
+
+# if (PM_PEAKS_CULL_WITH_SMOOTHED_IMAGE)
+        // New (post r30869) style of culling using the smoothed image and variance (S/N)
+        // if we cull using the significance image, then the definition of variance is different (thus the bool in arg 8)
+        if (pmFootprintCullPeaks(signifR->image, signifR->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) {
             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
         }
+# else 
+        // Old (pre r30869) style of culling using the raw image and variance
+        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, false) != PS_ERR_NONE) {
+             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
+        }
+# endif
+
         float dtime = psTimerMark ("psphot.cull.footprints");
         if (dtime > 1.0) {

@@ -19 +19 @@
 psphot-config
 Doxyfile
+a.out.dSYM