Index: /branches/eam_branch_20080413/psphot/src/psphotCullPeaks.c
===================================================================
--- /branches/eam_branch_20080413/psphot/src/psphotCullPeaks.c	(revision 17461)
+++ /branches/eam_branch_20080413/psphot/src/psphotCullPeaks.c	(revision 17462)
@@ -52,5 +52,5 @@
   * starting point, discard the peak.
   */
-psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
+psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
 				 const psImage *weight,	// corresponding variance image
 				 pmFootprint *fp, // Footprint containing mortal peaks
@@ -125,5 +125,5 @@
 	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
 	brightPeaks->n = 0;		// don't double free
-	psImage *idImg = pmSetFootprintID(peakFootprint, peak_id);
+	psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
 	psFree(peakFootprint);
 
@@ -158,2 +158,107 @@
 }
 
+ /*
+  * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
+  * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
+  * to reach a still higher peak --- and if that coll's more than nsigma DN below your
+  * starting point, discard the peak.
+  */
+
+# define IN_PEAK 1 
+psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
+				 const psImage *weight,	// corresponding variance image
+				 pmFootprint *fp, // Footprint containing mortal peaks
+				 const float nsigma_delta, // how many sigma above local background a peak
+				 // needs to be to survive
+				 const float min_threshold) { // minimum permitted coll height
+    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
+    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
+    assert (img->row0 == weight->row0 && img->col0 == weight->col0);
+    assert (fp != NULL);
+
+    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
+	return PS_ERR_NONE;
+    }
+
+    psRegion subRegion;			// desired subregion; 1 larger than bounding box (grr)
+    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
+    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
+
+    psImage *subImg = psImageSubset((psImage *)img, subRegion);
+    psImage *subWt = psImageSubset((psImage *)weight, subRegion);
+    assert (subImg != NULL && subWt != NULL);
+
+    psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
+
+    // We need a psArray of peaks brighter than the current peak.  
+    // We reject peaks which either:
+    // 1) are below the local threshold
+    // 2) have a brighter peak within their threshold
+
+    // allocate the full-sized array.  if the final array is much smaller, we can realloc
+    // at that point.
+    psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n);
+    psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
+
+    // The brightest peak is always safe; go through other peaks trying to cull them
+    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
+	const pmPeak *peak = fp->peaks->data[i];
+	int x = peak->x - subImg->col0;
+	int y = peak->y - subImg->row0;
+	//
+	// Find the level nsigma below the peak that must separate the peak
+	// from any of its friends
+	//
+	assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
+	const float stdev = sqrt(subWt->data.F32[y][x]);
+	float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
+	if (isnan(threshold) || threshold < min_threshold) {
+	    // min_threshold is assumed to be below the detection threshold,
+	    // so all the peaks are pmFootprint, and this isn't the brightest
+	    continue;
+	}
+
+	// XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
+	if (threshold > subImg->data.F32[y][x]) {
+	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
+	}
+
+	// XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
+	// perhaps this should alloc a single ID image above and pass it in to be set.
+
+	// XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
+
+	// at this point brightPeaks only has the peaks brighter than the current
+	pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
+
+	// XXX need to supply the image here
+	// we set the IDs to either 1 (in peak) or 0 (not in peak)
+	pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
+	psFree(peakFootprint);
+
+	// Check if any of the previous (brighter) peaks are within the footprint of this peak
+	// If so, the current peak is bogus; drop it.
+	bool keep = true;
+	for (int j = 0; keep && (j < brightPeaks->n); j++) {
+	    const pmPeak *peak2 = fp->peaks->data[j];
+	    int x2 = peak2->x - subImg->col0;
+	    int y2 = peak2->y - subImg->row0;
+	    if (idImg->data.S32[y2][x2] == IN_PEAK) 
+		// There's a brighter peak within the footprint above threshold; so cull our initial peak
+		keep = false;
+	}
+	if (!keep) continue;
+
+	psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
+    }
+
+    psFree (fp->peaks);
+    fp->peaks = brightPeaks;
+
+    psFree(idImg);
+    psFree(subImg);
+    psFree(subWt);
+
+    return PS_ERR_NONE;
+}
+