Index: trunk/psphot/src/pmFootprint.c
===================================================================
--- trunk/psphot/src/pmFootprint.c	(revision 13372)
+++ trunk/psphot/src/pmFootprint.c	(revision 13441)
@@ -448,7 +448,7 @@
     void *data;				// saved data
     void *data_home;			// where the data was saved from
-} Startspan;
-
-static void startspanFree(Startspan *sspan) {
+} Old_Spartspan;
+
+static void Old_startspanFree(Old_Spartspan *sspan) {
     if (sspan->nbyte > 0) {
 	memcpy(sspan->data_home, sspan->data, sspan->nbyte); // restore image's pixels to their pristine state
@@ -459,11 +459,11 @@
 }
 
-static Startspan *
-StartspanAlloc(const pmSpan *span, // The span in question
+static Old_Spartspan *
+Old_SpartspanAlloc(const pmSpan *span, // The span in question
 	       psImage *img,// the data wherein lives the span
 	       const PM_SSPAN_DIR dir	// Should we continue searching towards the top of the image?
     ) {
-    Startspan *sspan = psAlloc(sizeof(Startspan));
-    psMemSetDeallocator(sspan, (psFreeFunc)startspanFree);
+    Old_Spartspan *sspan = psAlloc(sizeof(Old_Spartspan));
+    psMemSetDeallocator(sspan, (psFreeFunc)Old_startspanFree);
 
     sspan->span = psMemIncrRefCounter((void *)span);
@@ -507,15 +507,15 @@
 
 //
-// Add a new Startspan to an array of Startspans
+// Add a new Old_Spartspan to an array of Old_Spartspans
 //
-static void add_startspan(psArray *startspans, // the saved Startspans
+static void old_add_startspan(psArray *startspans, // the saved Old_Spartspans
 			  const pmSpan *sp, // the span in question
 			  const psImage *img, // image to save/restore (or NULL)
 			  const PM_SSPAN_DIR dir) { // the desired direction to search
     if (dir == PM_SSPAN_RESTART) {
-	add_startspan(startspans, sp, img,  PM_SSPAN_UP);
-	add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN);
+	old_add_startspan(startspans, sp, img,  PM_SSPAN_UP);
+	old_add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN);
     } else {
-	Startspan *sspan = StartspanAlloc(sp, (psImage *)img, dir);
+	Old_Spartspan *sspan = Old_SpartspanAlloc(sp, (psImage *)img, dir);
 	psArrayAdd(startspans, 1, sspan);
 	psFree(sspan);			// as it's now owned by startspans
@@ -525,12 +525,12 @@
 /************************************************************************************************************/
 /*
- * Search the image for pixels above threshold, starting at a single Startspan.
+ * Search the image for pixels above threshold, starting at a single Old_Spartspan.
  * We search the array looking for one to process; it'd be better to move the
  * ones that we're done with to the end, but it probably isn't worth it for
  * the anticipated uses of this routine.
  *
- * This is the guts of pmFindFootprintAtPoint
- */
-static bool do_startspan(pmFootprint *fp, // the footprint that we're building
+ * This is the guts of pmOldFindFootprintAtPoint
+ */
+static bool old_do_startspan(pmFootprint *fp, // the footprint that we're building
 			 const psImage *img, // the psImage we're working on
 			 const float threshold,	// Threshold
@@ -556,5 +556,5 @@
     /********************************************************************************************************/
     
-    Startspan *sspan = NULL;
+    Old_Spartspan *sspan = NULL;
     for (int i = 0; i < startspans->n; i++) {
 	sspan = startspans->data[i];
@@ -563,5 +563,5 @@
 	}
     }
-    if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process
+    if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Old_Spartspans to process
 	return false;
     }
@@ -615,5 +615,5 @@
 	    const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
 	    
-	    add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
+	    old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
 	}
 	//
@@ -646,15 +646,15 @@
 		    if (sx1 <= x1) {
 			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
-			add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);
+			old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);
 		    } else {		// overhangs to right
 			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
-			add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);
+			old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_DONE);
 			sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
-			add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
+			old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
 		    }
 		    first = false;
 		} else {
 		    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
-		    add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
+		    old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
 		}
 	    }
@@ -689,5 +689,5 @@
  */
 pmFootprint *
-pmFindFootprintAtPoint(const psImage *img,	// image to search
+pmOldFindFootprintAtPoint(const psImage *img,	// image to search
 		       const float threshold,	// Threshold
 		       int row, int col) { // starting position (in img's parent's coordinate system)
@@ -753,13 +753,395 @@
        const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1);
 
-       add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
+       old_add_startspan(startspans, sp, (psImage *)img, PM_SSPAN_RESTART);
    }
    /*
-    * Now workout from those Startspans, searching for pixels above threshold
+    * Now workout from those Old_Spartspans, searching for pixels above threshold
     */
-   while (do_startspan(fp, img, threshold, startspans)) continue;
+   while (old_do_startspan(fp, img, threshold, startspans)) continue;
    /*
     * Cleanup
     */
+   psFree(startspans);			// restores the image pixel
+
+   return fp;				// pmFootprint really
+}
+
+/************************************************************************************************************/
+/*
+ * A data structure to hold the starting point for a search for pixels above threshold,
+ * used by pmFindFootprintAtPoint
+ *
+ * We don't want to find this span again --- it's already part of the footprint ---
+ * so we set the pixels in the image to as small a value as we can.  N.b. We may want
+ * to reconsider this when dealing with difference images
+ */
+#if 0
+typedef enum {PM_SSPAN_DOWN = 0,	// scan down from this span
+	      PM_SSPAN_UP,		// scan up from this span
+	      PM_SSPAN_RESTART,		// restart scanning from this span
+	      PM_SSPAN_DONE		// this span is processed
+} PM_SSPAN_DIR;				// How to continue searching
+#endif
+
+//
+// An enum for mask's pixel values.  We're looking for pixels that are above threshold, and
+// we keep extra book-keeping information in the PM_SSPAN_STOP plane.  It's simpler to be
+// able to check for 
+enum {
+    PM_SSPAN_INITIAL = 0x0,		// initial state of pixels.
+    PM_SSPAN_DETECTED = 0x1,		// we've seen this pixel
+    PM_SSPAN_STOP = 0x2			// you may stop searching when you see this pixel
+};
+
+typedef struct {
+    const pmSpan *span;			// save the pixel range
+    PM_SSPAN_DIR direction;		// How to continue searching
+    bool stop;				// should we stop searching?
+} Spartspan;
+
+static void startspanFree(Spartspan *sspan) {
+    psFree((void *)sspan->span);
+}
+
+static Spartspan *
+SpartspanAlloc(const pmSpan *span,	// The span in question
+	       psImage *mask,		// Pixels that we've already detected
+	       const PM_SSPAN_DIR dir	// Should we continue searching towards the top of the image?
+    ) {
+    Spartspan *sspan = psAlloc(sizeof(Spartspan));
+    psMemSetDeallocator(sspan, (psFreeFunc)startspanFree);
+
+    sspan->span = psMemIncrRefCounter((void *)span);
+    sspan->direction = dir;
+    sspan->stop = false;
+    
+    if (mask != NULL) {			// remember that we've detected these pixels
+	psMaskType *mpix = &mask->data.PS_TYPE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];
+
+	for (int i = 0; i <= span->x1 - span->x0; i++) {
+	    mpix[i] |= PM_SSPAN_DETECTED;
+	    if (mpix[i] & PM_SSPAN_STOP) {
+		sspan->stop = true;
+	    }
+	}
+    }
+    
+    return sspan;
+}
+
+//
+// Add a new Spartspan to an array of Spartspans.  Iff we see a stop bit, return true
+//
+static bool add_startspan(psArray *startspans, // the saved Spartspans
+			  const pmSpan *sp, // the span in question
+			  psImage *mask, // mask of detected/stop pixels
+			  const PM_SSPAN_DIR dir) { // the desired direction to search
+    if (dir == PM_SSPAN_RESTART) {
+	if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
+	    add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
+	    return true;
+	}
+    } else {
+	Spartspan *sspan = SpartspanAlloc(sp, mask, dir);
+	if (sspan->stop) {		// we detected a stop bit
+	    psFree(sspan);		// don't allocate new span
+
+	    return true;
+	} else {
+	    psArrayAdd(startspans, 1, sspan);
+	    psFree(sspan);		// as it's now owned by startspans
+	}
+    }
+
+    return false;
+}
+
+/************************************************************************************************************/
+/*
+ * Search the image for pixels above threshold, starting at a single Spartspan.
+ * We search the array looking for one to process; it'd be better to move the
+ * ones that we're done with to the end, but it probably isn't worth it for
+ * the anticipated uses of this routine.
+ *
+ * This is the guts of pmFindFootprintAtPoint
+ */
+static bool do_startspan(pmFootprint *fp, // the footprint that we're building
+			 const psImage *img, // the psImage we're working on
+			 psImage *mask, // the associated masks
+			 const float threshold,	// Threshold
+			 psArray *startspans) {	// specify which span to process next
+    bool F32 = false;			// is this an F32 image?
+    if (img->type.type == PS_TYPE_F32) {
+	F32 = true;
+    } else if (img->type.type == PS_TYPE_S32) {
+	F32 = false;
+    } else {				// N.b. You can't trivially add more cases here; F32 is just a bool
+	psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
+	return NULL;
+    }
+
+    psF32 *imgRowF32 = NULL;		// row pointer if F32
+    psS32 *imgRowS32 = NULL;		//  "   "   "  "  !F32
+    psMaskType *maskRow = NULL;		//  masks's row pointer
+    
+    const int row0 = img->row0;
+    const int col0 = img->col0;
+    const int numRows = img->numRows;
+    const int numCols = img->numCols;
+    
+    /********************************************************************************************************/
+    
+    Spartspan *sspan = NULL;
+    for (int i = 0; i < startspans->n; i++) {
+	sspan = startspans->data[i];
+	if (sspan->direction != PM_SSPAN_DONE) {
+	    break;
+	}
+	if (sspan->stop) {
+	    break;
+	}
+    }
+    if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Spartspans to process
+	return false;
+    }
+    if (sspan->stop) {			// they don't want any more spans processed
+	return false;
+    }
+    /*
+     * Work
+     */
+    const PM_SSPAN_DIR dir = sspan->direction;
+    /*
+     * Set initial span to the startspan
+     */
+    int x0 = sspan->span->x0 - col0, x1 = sspan->span->x1 - col0;
+    /*
+     * Go through image identifying objects
+     */
+    int nx0, nx1 = -1;			// new values of x0, x1
+    const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row
+    bool stop = false;			// should I stop searching for spans?
+
+    for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) {
+	imgRowF32 = img->data.F32[i];	// only one of
+	imgRowS32 = img->data.S32[i];	//      these is valid!
+	maskRow = mask->data.PS_TYPE_MASK_DATA[i];
+	//
+	// Search left from the pixel diagonally to the left of (i - di, x0). If there's
+	// a connected span there it may need to grow up and/or down, so push it onto
+	// the stack for later consideration
+	//
+	nx0 = -1;
+	for (int j = x0 - 1; j >= -1; j--) {
+	    double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+	    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+		if (j < x0 - 1) {	// we found some pixels above threshold
+		    nx0 = j + 1;
+		}
+		break;
+	    }
+	}
+
+	if (nx0 < 0) {			// no span to the left
+	    nx1 = x0 - 1;		// we're going to resume searching at nx1 + 1
+	} else {
+	    //
+	    // Search right in leftmost span
+	    //
+	    //nx1 = 0;			// make gcc happy
+	    for (int j = nx0 + 1; j <= numCols; j++) {
+		double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+		if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+		    nx1 = j - 1;
+		    break;
+		}
+	    }
+	    
+	    const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
+	    
+	    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+		stop = true;
+		break;
+	    }
+	}
+	//
+	// Now look for spans connected to the old span.  The first of these we'll
+	// simply process, but others will have to be deferred for later consideration.
+	//
+	// In fact, if the span overhangs to the right we'll have to defer the overhang
+	// until later too, as it too can grow in both directions
+	//
+	// Note that column numCols exists virtually, and always ends the last span; this
+	// is why we claim below that sx1 is always set
+	//
+	bool first = false;		// is this the first new span detected?
+	for (int j = nx1 + 1; j <= x1 + 1; j++) {
+	    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+	    if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
+		int sx0 = j++;		// span that we're working on is sx0:sx1
+		int sx1 = -1;		// We know that if we got here, we'll also set sx1
+		for (; j <= numCols; j++) {
+		    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
+		    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
+			sx1 = j;
+			break;
+		    }
+		}
+		assert (sx1 >= 0);
+
+		const pmSpan *sp;
+		if (first) {
+		    if (sx1 <= x1) {
+			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
+			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
+			    stop = true;
+			    break;
+			}
+		    } else {		// overhangs to right
+			sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
+			if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
+			    stop = true;
+			    break;
+			}
+			sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
+			if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+			    stop = true;
+			    break;
+			}
+		    }
+		    first = false;
+		} else {
+		    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
+		    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
+			stop = true;
+			break;
+		    }
+		}
+	    }
+	}
+
+	if (stop || first == false) {	// we're done
+	    break;
+	}
+
+	x0 = nx0; x1 = nx1;
+    }
+    /*
+     * Cleanup
+     */
+
+    sspan->direction = PM_SSPAN_DONE;
+    return stop ? false : true;
+}
+
+/*
+ * Go through an image, starting at (row, col) and assembling all the pixels
+ * that are connected to that point (in a chess kings-move sort of way) into
+ * a pmFootprint.
+ *
+ * This is much slower than pmFindFootprints if you want to find lots of
+ * footprints, but if you only want a small region about a given point it
+ * can be much faster
+ *
+ * N.b. The returned pmFootprint is not in "normal form"; that is the pmSpans
+ * are not sorted by increasing y, x0, x1.  If this matters to you, call
+ * pmFootprintNormalize()
+ */
+pmFootprint *
+pmFindFootprintAtPoint(const psImage *img,	// image to search
+		       const float threshold,	// Threshold
+		       const psArray *peaks, // array of peaks; finding one terminates search for footprint
+		       int row, int col) { // starting position (in img's parent's coordinate system)
+   assert(img != NULL);
+
+   bool F32 = false;			// is this an F32 image?
+   if (img->type.type == PS_TYPE_F32) {
+       F32 = true;
+   } else if (img->type.type == PS_TYPE_S32) {
+       F32 = false;
+   } else {				// N.b. You can't trivially add more cases here; F32 is just a bool
+       psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
+       return NULL;
+   }
+   psF32 *imgRowF32 = NULL;		// row pointer if F32
+   psS32 *imgRowS32 = NULL;		//  "   "   "  "  !F32
+   
+   const int row0 = img->row0;
+   const int col0 = img->col0;
+   const int numRows = img->numRows;
+   const int numCols = img->numCols;
+/*
+ * Is point in image, and above threshold?
+ */
+   row -= row0; col -= col0;
+   if (row < 0 || row >= numRows ||
+       col < 0 || col >= numCols) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",
+		row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
+       return NULL;
+   }
+
+   double pixVal = F32 ? img->data.F32[row][col] : img->data.S32[row][col];
+   if (pixVal < threshold) {
+       return pmFootprintAlloc(0, img);
+   }
+   
+   pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img);
+/*
+ * We need a mask for two purposes; to indicate which pixels are already detected,
+ * and to store the "stop" pixels --- those that, once reached, should stop us
+ * looking for the rest of the pmFootprint.  These are generally set from peaks.
+ */
+   psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK);
+   P_PSIMAGE_SET_ROW0(mask, row0);
+   P_PSIMAGE_SET_COL0(mask, col0);
+   psImageInit(mask, PM_SSPAN_INITIAL);
+   //
+   // Set stop bits from peaks list
+   //
+   assert (peaks == NULL || peaks->n == 0 || pmIsPeak(peaks->data[0]));
+   if (peaks != NULL) {
+       for (int i = 0; i < peaks->n; i++) {
+	   pmPeak *peak = peaks->data[i];
+	   mask->data.PS_TYPE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
+       }
+   }
+/*
+ * Find starting span passing through (row, col)
+ */
+   psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search
+   
+   imgRowF32 = img->data.F32[row];	// only one of
+   imgRowS32 = img->data.S32[row];	//      these is valid!
+   psMaskType *maskRow = mask->data.PS_TYPE_MASK_DATA[row];
+   {
+       int i;
+       for (i = col; i >= 0; i--) {
+	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
+	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+	       break;
+	   }
+       }
+       int i0 = i;
+       for (i = col; i < numCols; i++) {
+	   pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
+	   if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
+	       break;
+	   }
+       }
+       int i1 = i;
+       const pmSpan *sp = pmFootprintAddSpan(fp, row + row0, i0 + col0 + 1, i1 + col0 - 1);
+
+       (void)add_startspan(startspans, sp, mask, PM_SSPAN_RESTART);
+   }
+   /*
+    * Now workout from those Spartspans, searching for pixels above threshold
+    */
+   while (do_startspan(fp, img, mask, threshold, startspans)) continue;
+   /*
+    * Cleanup
+    */
+   psFree(mask);
    psFree(startspans);			// restores the image pixel
 
@@ -1058,5 +1440,5 @@
   * starting point, discard the peak.
   */
-psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
+psErrorCode pmOldFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
 				 const psImage *weight,	// corresponding variance image
 				 pmFootprint *fp, // Footprint containing mortal peaks
@@ -1117,5 +1499,5 @@
 #else
 	const int peak_id = 1; // the ID for the peak of interest
-	pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, peak->y, peak->x);
+	pmFootprint *peakFootprint = pmOldFindFootprintAtPoint(subImg, threshold, peak->y, peak->x);
 	psImage *idImg = pmSetFootprintID(peakFootprint, peak_id);
 	psFree(peakFootprint);
@@ -1143,4 +1525,103 @@
     }
 
+    psFree((psImage *)subImg);
+    psFree((psImage *)subWt);
+
+    return PS_ERR_NONE;
+}
+
+/************************************************************************************************************/
+ /*
+  * 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.
+  */
+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;
+    const psImage *subImg = psImageSubset((psImage *)img, subRegion);
+    const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
+    assert (subImg != NULL && subWt != NULL);
+    //
+    // We need a psArray of peaks brighter than the current peak.  We'll fake this
+    // for efficiency --- shhh
+    //
+    psArray *brightPeaks = psArrayAlloc(0);
+    psFree(brightPeaks->data);
+    brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
+    //
+    // 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) {
+#if 1					// min_threshold is assumed to be below the detection threshold,
+					// so all the peaks are pmFootprint, and this isn't the brightest
+	    (void)psArrayRemoveIndex(fp->peaks, i);
+	    i--;			// we moved everything down one
+	    continue;
+#else
+#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
+	    threshold = min_threshold;
+#endif
+	}
+	if (threshold > subImg->data.F32[y][x]) {
+	    threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
+	}
+
+	const int peak_id = 1; // the ID for the peak of interest
+	brightPeaks->n = i;		// only stop at a peak brighter than we are
+	pmFootprint *peakFootprint = pmFindFootprintAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
+	brightPeaks->n = 0;		// don't double free
+	psImage *idImg = pmSetFootprintID(peakFootprint, peak_id);
+	psFree(peakFootprint);
+
+	int j;
+	for (j = 0; j < i; j++) {
+	    const pmPeak *peak2 = fp->peaks->data[j];
+	    int x2 = peak2->x - subImg->col0;
+	    int y2 = peak2->y - subImg->row0;
+	    const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
+
+	    if (peak2_id == peak_id) {	// There's a brighter peak within the footprint above
+		;			// threshold; so cull our initial peak
+		(void)psArrayRemoveIndex(fp->peaks, i);
+		i--;			// we moved everything down one
+		break;
+	    }
+	}
+	if (j == i) {
+	    j++;
+	}
+
+	psFree(idImg);
+    }
+
+    brightPeaks->n = 0; psFree(brightPeaks);
     psFree((psImage *)subImg);
     psFree((psImage *)subWt);