Changeset 17516

Timestamp:

May 4, 2008, 2:10:40 PM (18 years ago)

Author:

eugene

Message:

some efficiency improvements and cleanups to the footprint code

Location:

trunk/psphot/src

Files:

• pmFootprint.h (modified) (1 diff)
• pmFootprintArrayGrow.c (modified) (4 diffs)
• pmFootprintIDs.c (modified) (2 diffs)
• pmFootprintsAssignPeaks.c (modified) (2 diffs)
• psphotCullPeaks.c (modified) (3 diffs)
• psphotFindFootprints.c (modified) (1 diff)
• psphotFitSourcesLinear.c (modified) (9 diffs)

trunk/psphot/src/pmFootprint.h

@@ -47 +47 @@
 
 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

trunk/psphot/src/pmFootprintArrayGrow.c

@@ -8 +8 @@
     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 +18 @@
      */
     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 +36 @@
     }
 
+# 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 +62 @@
     psFree((psArray *)peaks);
 
+    psLogMsg ("psphot", PS_LOG_DETAIL, "finished grow: %f sec\n", psTimerMark ("grow"));
+
     return grown;
     

trunk/psphot/src/pmFootprintIDs.c

@@ -74 +74 @@
  * 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 +84 @@
    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);

trunk/psphot/src/pmFootprintsAssignPeaks.c

@@ -8 +8 @@
 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 +54 @@
     
     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;

trunk/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;
+}
+

trunk/psphot/src/psphotFindFootprints.c

@@ -1 +1 @@
 # 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) {

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -158 +158 @@
     // 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 +175 @@
         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 +194 @@
         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 +202 @@
         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 +216 @@
 }
 
+// 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 +228 @@
     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 +239 @@
         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 +276 @@
         }
     }
+    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 +297 @@
         psSparseBorderElementT (border, 2, 2, y2);
     }
+    psLogMsg ("psphot.ensemble", PS_LOG_INFO, "step 5: %f sec\n", psTimerMark ("psphot"));
+
     return true;
 }