Changeset 17462

Timestamp:

Apr 21, 2008, 8:21:54 AM (18 years ago)

Author:

eugene

Message:

changed usage of psArrayRemoveIndex to an algorithm which saves the new peaks in a separate array (cleaner on array api as well)

File:

• branches/eam_branch_20080413/psphot/src/psphotCullPeaks.c (modified) (3 diffs)

branches/eam_branch_20080413/psphot/src/psphotCullPeaks.c

@@ -52 +52 @@
   * 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 +125 @@
         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 +158 @@
 }
 
+ /*
+  * 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;
+}
+