Index: trunk/psphot/src/pmFootprint.h
===================================================================
--- trunk/psphot/src/pmFootprint.h	(revision 17444)
+++ trunk/psphot/src/pmFootprint.h	(revision 17516)
@@ -47,5 +47,5 @@
 
 psImage *pmSetFootprintArrayIDs(const psArray *footprints, const bool relativeIDs);
-psImage *pmSetFootprintID(const pmFootprint *fp, const int id);
+psImage *pmSetFootprintID(psImage *idImage, const pmFootprint *fp, const int id);
 void pmSetFootprintArrayIDsForImage(psImage *idImage,
 				    const psArray *footprints, // the footprints to insert
Index: trunk/psphot/src/pmFootprintArrayGrow.c
===================================================================
--- trunk/psphot/src/pmFootprintArrayGrow.c	(revision 17444)
+++ trunk/psphot/src/pmFootprintArrayGrow.c	(revision 17516)
@@ -8,4 +8,6 @@
     assert (footprints->n == 0 || pmFootprintTest(footprints->data[0]));
 
+    psTimerStart ("grow");
+
     if (footprints->n == 0) {		// we don't know the size of the footprint's region
 	return psArrayAlloc(0);
@@ -16,4 +18,6 @@
      */
     psImage *idImage = pmSetFootprintArrayIDs(footprints, true);
+    psLogMsg ("psphot", PS_LOG_DETAIL, "set footprint array IDs: %f sec\n", psTimerMark ("grow"));
+
     if (r <= 0) {
 	r = 1;				// r == 1 => no grow
@@ -32,8 +36,20 @@
     }
 
+# if (1)    
+    psImage *f32ImageIn = psImageCopy (NULL, idImage, PS_TYPE_F32);
+    psImage *f32ImageOut = psImageConvolveFFT(NULL, f32ImageIn, NULL, 0, circle);
+    psImage *grownIdImage = psImageCopy (NULL, f32ImageOut, PS_TYPE_S32);
+    psFree (f32ImageIn);
+    psFree (f32ImageOut);
+# else
     psImage *grownIdImage = psImageConvolveDirect(NULL, idImage, circle); // Here's the actual grow step
+# endif
+
+    psLogMsg ("psphot", PS_LOG_DETAIL, "convolved with grow disc: %f sec\n", psTimerMark ("grow"));
     psFree(circle);	
 
     psArray *grown = pmFootprintsFind(grownIdImage, 0.5, 1); // and here we rebuild the grown footprints
+    psLogMsg ("psphot", PS_LOG_DETAIL, "found grown footprints: %f sec\n", psTimerMark ("grow"));
+
     assert (grown != NULL);
     psFree(idImage); psFree(grownIdImage);
@@ -46,4 +62,6 @@
     psFree((psArray *)peaks);
 
+    psLogMsg ("psphot", PS_LOG_DETAIL, "finished grow: %f sec\n", psTimerMark ("grow"));
+
     return grown;
     
Index: trunk/psphot/src/pmFootprintIDs.c
===================================================================
--- trunk/psphot/src/pmFootprintIDs.c	(revision 17444)
+++ trunk/psphot/src/pmFootprintIDs.c	(revision 17516)
@@ -74,5 +74,6 @@
  * Set an image to the value of footprint's ID whever they may fall
  */
-psImage *pmSetFootprintID(const pmFootprint *fp, // the footprint to insert
+psImage *pmSetFootprintID(psImage *idImage,
+			  const pmFootprint *fp, // the footprint to insert
 			  const int id) {	// the desired ID
    assert(fp != NULL && pmFootprintTest((const psPtr)fp));
@@ -83,5 +84,11 @@
    assert (numCols >= 0 && numRows >= 0);
    
-   psImage *idImage = psImageAlloc(numCols, numRows, PS_TYPE_S32);
+   if (idImage == NULL) {
+       idImage = psImageAlloc(numCols, numRows, PS_TYPE_S32);
+   } else {
+       assert (idImage->numCols == numCols);
+       assert (idImage->numRows == numRows);
+       // XXX assert on type (S32)
+   }
    P_PSIMAGE_SET_ROW0(idImage, row0);
    P_PSIMAGE_SET_COL0(idImage, col0);
Index: trunk/psphot/src/pmFootprintsAssignPeaks.c
===================================================================
--- trunk/psphot/src/pmFootprintsAssignPeaks.c	(revision 17444)
+++ trunk/psphot/src/pmFootprintsAssignPeaks.c	(revision 17516)
@@ -8,5 +8,5 @@
 psErrorCode
 pmFootprintsAssignPeaks(psArray *footprints,	// the pmFootprints
-			  const psArray *peaks) { // the pmPeaks
+			const psArray *peaks) { // the pmPeaks
     assert (footprints != NULL);
     assert (footprints->n == 0 || pmFootprintTest(footprints->data[0]));
@@ -54,18 +54,61 @@
     
     psFree(ids);
-    //
+
     // Make sure that peaks within each footprint are sorted and unique
-    //
     for (int i = 0; i < footprints->n; i++) {
+	
 	pmFootprint *fp = footprints->data[i];
+
+	// XXX are we allowed to have peak-less footprints??
+	if (!fp->peaks->n) continue;
+
         fp->peaks = psArraySort(fp->peaks, pmPeakSortBySN);
 
-	for (int j = 1; j < fp->peaks->n; j++) { // check for duplicates
-	    if (fp->peaks->data[j] == fp->peaks->data[j-1]) {
-		(void)psArrayRemoveIndex(fp->peaks, j);
-		j--;			// we moved everything down one
+	// XXX EAM : the algorithm below should be much faster than using psArrayRemove if
+	// the number of peaks in the footprint is large, or if there are no duplicates.
+	// if we have a lot of small-number peak arrays with duplicates, this may be
+	// slower.
+
+	// track the number of good peaks in the footprint
+	int nGood = 1;
+
+	// check for duplicates
+	// on first pass, we set the index to NULL if peak is a duplicate
+	// XXX EAM : this can leave behind duplicates of the same S/N
+	// (if sorted list has A, B, A, B ...)
+	for (int j = 1; j < fp->peaks->n; j++) { 
+	    if (fp->peaks->data[j] == fp->peaks->data[nGood]) {
+		// everything on the array has its own mem reference; free and drop this one
+		psFree (fp->peaks->data[j]);
+		fp->peaks->data[j] = NULL;
+	    } else {
+		nGood ++;
 	    }
 	}
+
+	// no deleted peaks, go to next footprint
+	if (nGood == fp->peaks->n) continue;
+
+	int nKeep = 0;
+
+	psArray *goodPeaks = psArrayAlloc (nGood);
+	// on second pass, save the good peaks
+	for (int j = 0; j < fp->peaks->n; j++) { // check for duplicates
+	    if (fp->peaks->data[j] == NULL) continue;
+	    // transfer the data (NULL to avoid double free)
+	    // this is only slightly sleazy
+	    goodPeaks->data[nKeep] = fp->peaks->data[j];
+	    fp->peaks->data[j] = NULL;
+	    nKeep ++;
+	}
+	psAssert (nGood == nKeep, "mis-counted nKeep or nGood");
+
+	// free the old (now NULL-filled) array
+	psFree (fp->peaks);
+	fp->peaks = goodPeaks;
     }
+
+    // (void)psArrayRemoveIndex(fp->peaks, j);
+
 
     return PS_ERR_NONE;
Index: trunk/psphot/src/psphotCullPeaks.c
===================================================================
--- trunk/psphot/src/psphotCullPeaks.c	(revision 17444)
+++ trunk/psphot/src/psphotCullPeaks.c	(revision 17516)
@@ -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;
+}
+
Index: trunk/psphot/src/psphotFindFootprints.c
===================================================================
--- trunk/psphot/src/psphotFindFootprints.c	(revision 17444)
+++ trunk/psphot/src/psphotFindFootprints.c	(revision 17516)
@@ -1,6 +1,3 @@
 # include "psphotInternal.h"
-
-// N.b. We're not culling this list; call pmFootprintCullPeaks if you
-// want to do that
 
 bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) {
Index: trunk/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- trunk/psphot/src/psphotFitSourcesLinear.c	(revision 17444)
+++ trunk/psphot/src/psphotFitSourcesLinear.c	(revision 17516)
@@ -158,4 +158,5 @@
     // set the sky, sky_x, sky_y components of border matrix
     SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER);
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "set border: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
     psSparseConstraint constraint;
@@ -174,4 +175,5 @@
         skyfit = NULL;
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
     // adjust I0 for fitSources and subtract
@@ -192,4 +194,5 @@
         source->mode |= PM_SOURCE_MODE_SUBTRACTED;
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
     // measure chisq for each source
@@ -199,4 +202,5 @@
         pmSourceChisq (model, source->pixels, source->maskObj, source->weight, maskVal);
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f (%d elements)\n", psTimerMark ("psphot"), sparse->Nelem);
 
     // psFree (index);
@@ -212,4 +216,6 @@
 }
 
+// XXX do we need this?
+
 // Calculate the weight terms for the sky fit component of the matrix.  This function operates
 // on the pixels which correspond to all of the sources of interest.  These elements fill in
@@ -222,4 +228,5 @@
     fullArray = psRegionForImage (readout->mask, fullArray);
     psImageMaskRegion (readout->mask, fullArray, "OR", PM_MASK_MARK);
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 1: %f sec\n", psTimerMark ("psphot"));
 
     // turn off MARK for all object pixels
@@ -232,4 +239,5 @@
         psImageMaskCircle (source->maskView, x, y, model->radiusFit, "AND", PS_NOT_U8(PM_MASK_MARK));
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 2: %f sec\n", psTimerMark ("psphot"));
 
     // accumulate the image statistics from the masked regions
@@ -268,7 +276,9 @@
         }
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 3: %f sec\n", psTimerMark ("psphot"));
 
     // turn off MARK for all image pixels
     psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_U8(PM_MASK_MARK));
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 4: %f sec\n", psTimerMark ("psphot"));
 
     // set the Border T elements
@@ -287,4 +297,6 @@
         psSparseBorderElementT (border, 2, 2, y2);
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 5: %f sec\n", psTimerMark ("psphot"));
+
     return true;
 }