Index: /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c
===================================================================
--- /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c	(revision 30868)
+++ /branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c	(revision 30869)
@@ -26,14 +26,51 @@
     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);
 	}