IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 13, 2010, 4:42:21 PM (16 years ago)
Author:
eugene
Message:

modify pmFootprintCullPeaks to reduce repeated allocs

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmFootprintCullPeaks.c

    r26893 r27672  
    2020#include "pmSpan.h"
    2121#include "pmFootprint.h"
     22#include "pmFootprintSpans.h"
    2223#include "pmPeaks.h"
     24
     25bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp);
    2326
    2427 /*
    … …  
    6568    psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
    6669
     70    // allocate the peakFootprint and peakFPSpans containers -- these are re-used by pmFootprintsFindAtPoint to minimize allocs in this function
     71    pmFootprint *peakFootprint = pmFootprintAlloc(fp->nspans, subImg);
     72    pmFootprintSpans *peakFPSpans = pmFootprintSpansAlloc(2*fp->nspans);
     73
     74    // allocate empty spans for the footprints
     75    pmFootprintAllocEmptySpans(peakFootprint, fp->nspans);
     76
     77    psImage *subMask = psImageCopy(NULL, subImg, PS_TYPE_IMAGE_MASK);
     78
    6779    // The brightest peak is always safe; go through other peaks trying to cull them
    6880    for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    … …  
    98110        }
    99111
    100         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    101         // perhaps this should alloc a single ID image above and pass it in to be set.
    102 
    103         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    104 
    105         // XXX this is not quite there yet:
    106         // pmFootprint *peakFootprint = NULL;
    107         // int area = subImg->numCols * subImg->numRows;
    108         // if (area > 30000) {
    109 
    110         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
     112        // init peakFootprint here?
     113        pmFootprintsFindAtPoint(peakFootprint, peakFPSpans, subImg, subMask, threshold, brightPeaks, peak->y, peak->x);
     114        if (peakFPSpans->nStartSpans > 2000) {
     115            // dumpfootprints(peakFootprint, peakFPSpans);
     116            // fprintf (stderr, "big footprint %d : %d\n", peakFootprint->nspans, peakFPSpans->nStartSpans);
     117            fprintf (stderr, "big test footprint: %f %f to %f %f (%d pix)\n", peakFootprint->bbox.x0, peakFootprint->bbox.y0, peakFootprint->bbox.x1, peakFootprint->bbox.y1, peakFootprint->npix);
     118        }
    111119
    112120        // at this point brightPeaks only has the peaks brighter than the current
    113121
    114         // XXX need to supply the image here
    115122        // we set the IDs to either 1 (in peak) or 0 (not in peak)
    116123        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
    117         psFree(peakFootprint);
    118124
    119125        // If this peak has not already been assigned to a source, then we can look for any
    … …  
    144150    psFree(subImg);
    145151    psFree(subWt);
     152    psFree(subMask);
     153    psFree(peakFootprint);
     154    psFree(peakFPSpans);
    146155
    147156    return PS_ERR_NONE;
    … …  
    149158
    150159
    151  /*
    152   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    153   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    154   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    155   * starting point, discard the peak.
    156   */
    157 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
    158                                  const psImage *weight, // corresponding variance image
    159                                  pmFootprint *fp, // Footprint containing mortal peaks
    160                                  const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
    161                                  const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
    162                                  const float min_threshold) { // minimum permitted coll height
    163     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    164     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    165     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    166     assert (fp != NULL);
     160bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp) {
    167161
    168     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    169         return PS_ERR_NONE;
     162    FILE *f1 = fopen ("fp.dat", "w");
     163    if (!f1) return false;
     164
     165    for (int i = 0; i < fp->nspans; i++) {
     166        pmSpan *sp = fp->spans->data[i];
     167        fprintf (f1, "%d - %d : %d\n", sp->x0, sp->x1, sp->y);
    170168    }
     169    fclose (f1);
    171170
    172     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    173     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    174     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    175     const psImage *subImg = psImageSubset((psImage *)img, subRegion);
    176     const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    177     assert (subImg != NULL && subWt != NULL);
    178     //
    179     // We need a psArray of peaks brighter than the current peak.  We'll fake this
    180     // by reusing the fp->peaks but lying about n.
    181     //
    182     // We do this for efficiency (otherwise I'd need two peaks lists), and we are
    183     // rather too chummy with psArray in consequence.  But it works.
    184     //
    185     psArray *brightPeaks = psArrayAlloc(0);
    186     psFree(brightPeaks->data);
    187     brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
    188     //
    189     // The brightest peak is always safe; go through other peaks trying to cull them
    190     //
    191     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    192         const pmPeak *peak = fp->peaks->data[i];
    193         int x = peak->x - subImg->col0;
    194         int y = peak->y - subImg->row0;
    195         //
    196         // Find the level nsigma below the peak that must separate the peak
    197         // from any of its friends
    198         //
    199         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
     171    FILE *f2 = fopen ("fpSp.dat", "w");
     172    if (!f2) return false;
    200173
    201         // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations)
    202         // add additional padding in quadrature:
    203 
    204         const float stdev = sqrt(subWt->data.F32[y][x]);
    205         const float flux = subImg->data.F32[y][x];
    206         const float fStdev = stdev/flux;
    207         const float stdev_pad = flux * hypot(fStdev, fPadding);
    208         float threshold = flux - nsigma_delta*stdev_pad;
    209 
    210         if (isnan(threshold) || threshold < min_threshold) {
    211 #if 1     // min_threshold is assumed to be below the detection threshold,
    212           // so all the peaks are pmFootprint, and this isn't the brightest
    213             // XXX mark peak to be dropped
    214             (void)psArrayRemoveIndex(fp->peaks, i);
    215             i--;                        // we moved everything down one
    216             continue;
    217 #else
    218 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
    219             threshold = min_threshold;
    220 #endif
    221         }
    222 
    223         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    224         if (threshold > subImg->data.F32[y][x]) {
    225             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    226         }
    227 
    228         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    229         // perhaps this should alloc a single ID image above and pass it in to be set.
    230 
    231         const int peak_id = 1;          // the ID for the peak of interest
    232         brightPeaks->n = i;             // only stop at a peak brighter than we are
    233 
    234         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    235 
    236         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    237         brightPeaks->n = 0;             // don't double free
    238         psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
    239         psFree(peakFootprint);
    240 
    241         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    242         // If so, the current peak is bogus; drop it.
    243         int j;
    244         for (j = 0; j < i; j++) {
    245             const pmPeak *peak2 = fp->peaks->data[j];
    246             int x2 = peak2->x - subImg->col0;
    247             int y2 = peak2->y - subImg->row0;
    248             const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
    249 
    250             if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
    251                 ;                       // threshold; so cull our initial peak
    252                 (void)psArrayRemoveIndex(fp->peaks, i);
    253                 i--;                    // we moved everything down one
    254                 break;
    255             }
    256         }
    257         if (j == i) {
    258             j++;
    259         }
    260 
    261         psFree(idImg);
     174    for (int i = 0; i < fpSp->nStartSpans; i++) {
     175        pmStartSpan *sp = fpSp->startspans->data[i];
     176        if (!sp->span) continue;
     177        fprintf (f2, "%d - %d : %d\n", sp->span->x0, sp->span->x1, sp->span->y);
    262178    }
    263 
    264     brightPeaks->n = 0; psFree(brightPeaks);
    265     psFree(subImg);
    266     psFree(subWt);
    267 
    268     return PS_ERR_NONE;
     179    fclose (f2);
     180    return true;
    269181}
    270 
Note: See TracChangeset for help on using the changeset viewer.