Changeset 13441

Timestamp:

May 21, 2007, 7:49:18 AM (19 years ago)

Author:

rhl

Message:

Rewrote pmFindFootprintAtPoint for speed using a mask and stop-pixels

Location:

trunk/psphot/src

Files:

• pmFootprint.c (modified) (13 diffs)
• pmFootprint.h (modified) (1 diff)

trunk/psphot/src/pmFootprint.c

@@ -448 +448 @@
     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 +459 @@
 }
 
-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 +507 @@
 
 //
-// 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 +525 @@
 /************************************************************************************************************/
 /*
- * 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 +556 @@
     /********************************************************************************************************/
     
-    Startspan *sspan = NULL;
+    Old_Spartspan *sspan = NULL;
     for (int i = 0; i < startspans->n; i++) {
         sspan = startspans->data[i];
@@ -563 +563 @@
         }
     }
-    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 +615 @@
             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 +646 @@
                     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 +689 @@
  */
 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 +753 @@
        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 +1440 @@
   * 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 +1499 @@
 #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 +1525 @@
     }
 
+    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);

trunk/psphot/src/pmFootprint.h

@@ -31 +31 @@
 void pmFootprintSetBBox(pmFootprint *fp);
 psArray *pmFindFootprints(const psImage *img, const float threshold, const int npixMin);
+pmFootprint *pmOldFindFootprintAtPoint(const psImage *img,
+                                    const float threshold,
+                                    int row, int col);
 pmFootprint *pmFindFootprintAtPoint(const psImage *img,
                                     const float threshold,
+                                    const psArray *peaks,
                                     int row, int col);