IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 3, 2010, 8:50:52 AM (16 years ago)
Author:
eugene
Message:

updates from trunk

Location:
branches/simtest_nebulous_branches
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • branches/simtest_nebulous_branches

  • branches/simtest_nebulous_branches/psModules

  • branches/simtest_nebulous_branches/psModules/src/objects/pmFootprintCullPeaks.c

    r24888 r27840  
    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 /*
    … …  
    2932  */
    3033
    31 # define IN_PEAK 1 
     34# define IN_PEAK 1
    3235psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
    33                                  const psImage *weight, // corresponding variance image
    34                                 pmFootprint *fp, // Footprint containing mortal peaks
    35                                 const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
    36                                 const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
    37                                 const float min_threshold) { // minimum permitted coll height
     36                                 const psImage *weight, // corresponding variance image
     37                                pmFootprint *fp, // Footprint containing mortal peaks
     38                                const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
     39                                const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
     40                                const float min_threshold) { // minimum permitted coll height
    3841    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    3942    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    … …  
    4245
    4346    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    44         return PS_ERR_NONE;
     47        return PS_ERR_NONE;
    4548    }
    4649
    47     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
     50    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    4851    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    4952    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    … …  
    5558    psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
    5659
    57     // We need a psArray of peaks brighter than the current peak. 
     60    // We need a psArray of peaks brighter than the current peak.
    5861    // We reject peaks which either:
    5962    // 1) are below the local threshold
    … …  
    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
    69         const pmPeak *peak = fp->peaks->data[i];
    70         int x = peak->x - subImg->col0;
    71         int y = peak->y - subImg->row0;
    72         //
    73         // Find the level nsigma below the peak that must separate the peak
    74         // from any of its friends
    75         //
    76         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
     81        const pmPeak *peak = fp->peaks->data[i];
     82        int x = peak->x - subImg->col0;
     83        int y = peak->y - subImg->row0;
     84        //
     85        // Find the level nsigma below the peak that must separate the peak
     86        // from any of its friends
     87        //
     88        assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
    7789
    78         // const float stdev = sqrt(subWt->data.F32[y][x]);
    79         // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
     90        // const float stdev = sqrt(subWt->data.F32[y][x]);
     91        // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    8092
    81         const float stdev = sqrt(subWt->data.F32[y][x]);
    82         const float flux = subImg->data.F32[y][x];
    83         const float fStdev = fabs(stdev/flux);
    84         const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
    85         // if flux is negative, careful with fStdev
     93        const float stdev = sqrt(subWt->data.F32[y][x]);
     94        const float flux = subImg->data.F32[y][x];
     95        const float fStdev = fabs(stdev/flux);
     96        const float stdev_pad = fabs(flux * hypot(fStdev, fPadding));
     97        // if flux is negative, careful with fStdev
    8698
    87         float threshold = flux - nsigma_delta*stdev_pad;
    88        
    89         if (isnan(threshold) || threshold < min_threshold) {
    90             // min_threshold is assumed to be below the detection threshold,
    91             // so all the peaks are pmFootprint, and this isn't the brightest
     99        float threshold = flux - nsigma_delta*stdev_pad;
     100
     101        if (isnan(threshold) || threshold < min_threshold) {
     102            // min_threshold is assumed to be below the detection threshold,
     103            // so all the peaks are pmFootprint, and this isn't the brightest
     104            continue;
     105        }
     106
     107        // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
     108        if (threshold > subImg->data.F32[y][x]) {
     109            threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
     110        }
     111
     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        }
     119
     120        // at this point brightPeaks only has the peaks brighter than the current
     121
     122        // we set the IDs to either 1 (in peak) or 0 (not in peak)
     123        pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
     124
     125        // If this peak has not already been assigned to a source, then we can look for any
     126        // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
     127        // are within the footprint of this peak If so, the current peak is bogus; drop it.
     128        bool keep = true;
     129        for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
     130            const pmPeak *peak2 = fp->peaks->data[j];
     131            int x2 = peak2->x - subImg->col0;
     132            int y2 = peak2->y - subImg->row0;
     133            if (idImg->data.S32[y2][x2] == IN_PEAK) {
     134                // There's a brighter peak within the footprint above threshold; so cull our initial peak
     135                keep = false;
     136            }
     137        }
     138        if (!keep) {
     139            psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
    92140            continue;
    93141        }
    94142
    95         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    96         if (threshold > subImg->data.F32[y][x]) {
    97             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    98         }
    99 
    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);
    111 
    112         // at this point brightPeaks only has the peaks brighter than the current
    113 
    114         // XXX need to supply the image here
    115         // we set the IDs to either 1 (in peak) or 0 (not in peak)
    116         pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
    117         psFree(peakFootprint);
    118 
    119         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    120         // If so, the current peak is bogus; drop it.
    121         bool keep = true;
    122         for (int j = 0; keep && (j < brightPeaks->n); j++) {
    123             const pmPeak *peak2 = fp->peaks->data[j];
    124             int x2 = peak2->x - subImg->col0;
    125             int y2 = peak2->y - subImg->row0;
    126             if (idImg->data.S32[y2][x2] == IN_PEAK)
    127                 // There's a brighter peak within the footprint above threshold; so cull our initial peak
    128                 keep = false;
    129         }
    130         if (!keep) continue;
    131 
    132         psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
     143        psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
    133144    }
    134145
    … …  
    139150    psFree(subImg);
    140151    psFree(subWt);
     152    psFree(subMask);
     153    psFree(peakFootprint);
     154    psFree(peakFPSpans);
    141155
    142156    return PS_ERR_NONE;
    … …  
    144158
    145159
    146  /*
    147   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    148   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    149   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    150   * starting point, discard the peak.
    151   */
    152 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
    153                                  const psImage *weight, // corresponding variance image
    154                                  pmFootprint *fp, // Footprint containing mortal peaks
    155                                  const float nsigma_delta, // how many sigma above local background a peak needs to be to survive
    156                                  const float fPadding, // fractional padding added to stdev since bright peaks have unreasonably high significance
    157                                  const float min_threshold) { // minimum permitted coll height
    158     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    159     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    160     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    161     assert (fp != NULL);
     160bool dumpfootprints (pmFootprint *fp, pmFootprintSpans *fpSp) {
    162161
    163     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    164         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);
    165168    }
     169    fclose (f1);
    166170
    167     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    168     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    169     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    170     const psImage *subImg = psImageSubset((psImage *)img, subRegion);
    171     const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    172     assert (subImg != NULL && subWt != NULL);
    173     //
    174     // We need a psArray of peaks brighter than the current peak.  We'll fake this
    175     // by reusing the fp->peaks but lying about n.
    176     //
    177     // We do this for efficiency (otherwise I'd need two peaks lists), and we are
    178     // rather too chummy with psArray in consequence.  But it works.
    179     //
    180     psArray *brightPeaks = psArrayAlloc(0);
    181     psFree(brightPeaks->data);
    182     brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
    183     //
    184     // The brightest peak is always safe; go through other peaks trying to cull them
    185     //
    186     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    187         const pmPeak *peak = fp->peaks->data[i];
    188         int x = peak->x - subImg->col0;
    189         int y = peak->y - subImg->row0;
    190         //
    191         // Find the level nsigma below the peak that must separate the peak
    192         // from any of its friends
    193         //
    194         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
     171    FILE *f2 = fopen ("fpSp.dat", "w");
     172    if (!f2) return false;
    195173
    196         // stdev ~ sqrt(flux) : for bright regions (f ~ 1e6, sqrt(f) = 1e3 -- unrealistically small deviations)
    197         // add additional padding in quadrature:
    198 
    199         const float stdev = sqrt(subWt->data.F32[y][x]);
    200         const float flux = subImg->data.F32[y][x];
    201         const float fStdev = stdev/flux;
    202         const float stdev_pad = flux * hypot(fStdev, fPadding);
    203         float threshold = flux - nsigma_delta*stdev_pad;
    204 
    205         if (isnan(threshold) || threshold < min_threshold) {
    206 #if 1     // min_threshold is assumed to be below the detection threshold,
    207           // so all the peaks are pmFootprint, and this isn't the brightest
    208             // XXX mark peak to be dropped
    209             (void)psArrayRemoveIndex(fp->peaks, i);
    210             i--;                        // we moved everything down one
    211             continue;
    212 #else
    213 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
    214             threshold = min_threshold;
    215 #endif
    216         }
    217 
    218         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    219         if (threshold > subImg->data.F32[y][x]) {
    220             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    221         }
    222 
    223         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    224         // perhaps this should alloc a single ID image above and pass it in to be set.
    225 
    226         const int peak_id = 1;          // the ID for the peak of interest
    227         brightPeaks->n = i;             // only stop at a peak brighter than we are
    228 
    229         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    230 
    231         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    232         brightPeaks->n = 0;             // don't double free
    233         psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
    234         psFree(peakFootprint);
    235 
    236         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    237         // If so, the current peak is bogus; drop it.
    238         int j;
    239         for (j = 0; j < i; j++) {
    240             const pmPeak *peak2 = fp->peaks->data[j];
    241             int x2 = peak2->x - subImg->col0;
    242             int y2 = peak2->y - subImg->row0;
    243             const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
    244 
    245             if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
    246                 ;                       // threshold; so cull our initial peak
    247                 (void)psArrayRemoveIndex(fp->peaks, i);
    248                 i--;                    // we moved everything down one
    249                 break;
    250             }
    251         }
    252         if (j == i) {
    253             j++;
    254         }
    255 
    256         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);
    257178    }
    258 
    259     brightPeaks->n = 0; psFree(brightPeaks);
    260     psFree((psImage *)subImg);
    261     psFree((psImage *)subWt);
    262 
    263     return PS_ERR_NONE;
     179    fclose (f2);
     180    return true;
    264181}
    265 
Note: See TracChangeset for help on using the changeset viewer.