Changeset 30869

Timestamp:

Mar 10, 2011, 5:29:34 PM (15 years ago)

Author:

eugene

Message:

modify the footprint analysis for potentially saturated objects -- use a max threshold of a smallish fraction of saturation so that we group all detections within the saturated region

File:

• branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c (modified) (1 diff)

branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c

@@ -26 +26 @@
     psAssert (status, "cannot find SKY_STDEV in readout->analysis");
 
-    const float min_threshold = nsigma_min*skyStdev;
+    const float MIN_THRESHOLD = nsigma_min*skyStdev;
     
+    // 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 
+
+    // in practice, we should make the threshold in such a case more like 0.5 * saturation...
+
+    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;
+
     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->flux > SAT_TEST_LEVEL) {
+            max_threshold = SAT_THRESHOLD;
+            brightPeak->type = PM_PEAK_SUSPECT_SATURATION;
+        }
+
+        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold) != PS_ERR_NONE) {
             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
         }