IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

updates from trunk

Location:
branches/simmosaic_branches
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • branches/simmosaic_branches

  • branches/simmosaic_branches/psphot

  • branches/simmosaic_branches/psphot/src

    • Property svn:ignore

      •  

        old new  
        1818psphotVersionDefinitions.h
        1919psphotMomentsStudy
         20psphotPetrosianStudy
         21psphotForced
         22psphotMakePSF
         23psphotStack

  • branches/simmosaic_branches/psphot/src/psphotCullPeaks.c

    r17516 r27839  
    55 */
    66psErrorCode
    7 psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
    8                 const psImage *weight,  // corresponding variance image
     7psphotCullPeaks(const pmReadout *readout,   // the image wherein lives the footprint
    98                const psMetadata *recipe,
    109                psArray *footprints) {  // array of pmFootprints
     10
    1111    bool status = false;
     12
    1213    float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA");
    1314    if (!status) {
    … …  
    1819        nsigma_min = 0;
    1920    }
    20     const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");
     21    float fPadding = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_PAD");
     22    if (!status) {
     23        fPadding = 0;
     24    }
     25    const float skyStdev = psMetadataLookupF32(&status, readout->analysis, "SKY_STDEV");
     26    psAssert (status, "cannot find SKY_STDEV in readout->analysis");
    2127
    22     return pmFootprintArrayCullPeaks(image, weight, footprints,
    23                                      nsigma_delta, nsigma_min*skyStdev);
    24 }
    25 
    26 
    27 /*
    28  * Cull an entire psArray of pmFootprints
    29  * XXX drop this intermediate level function?
    30  */
    31 psErrorCode
    32 pmFootprintArrayCullPeaks(const psImage *img, // the image wherein lives the footprint
    33                           const psImage *weight,        // corresponding variance image
    34                           psArray *footprints, // array of pmFootprints
    35                           const float nsigma_delta, // how many sigma above local background a peak
    36                                         // needs to be to survive
    37                           const float min_threshold) { // minimum permitted coll height
     28    const float min_threshold = nsigma_min*skyStdev;
     29   
    3830    for (int i = 0; i < footprints->n; i++) {
     31        // if (i % 50 == 0) fprintf (stderr, "cull %d\n", i);
    3932        pmFootprint *fp = footprints->data[i];
    40         if (pmFootprintCullPeaks(img, weight, fp, nsigma_delta, min_threshold) != PS_ERR_NONE) {
     33        if (fp->npix > 30000) {
     34            fprintf (stderr, "big footprint: %f %f to %f %f (%d pix)\n", fp->bbox.x0, fp->bbox.y0, fp->bbox.x1, fp->bbox.y1, fp->npix);
     35        }
     36        psTimerStart ("psphot.cull.footprints");
     37        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, min_threshold) != PS_ERR_NONE) {
    4138            return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
    4239        }
     40        float dtime = psTimerMark ("psphot.cull.footprints");
     41        if (dtime > 1.0) {
     42            fprintf (stderr, "slow cull for %d (%f sec)\n", i, dtime);
     43        }
    4344    }
    44    
    4545    return PS_ERR_NONE;
    4646}
    4747
    48  /*
    49   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    50   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    51   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    52   * starting point, discard the peak.
    53   */
    54 psErrorCode pmFootprintCullPeaks_OLD(const psImage *img, // the image wherein lives the footprint
    55                                  const psImage *weight, // corresponding variance image
    56                                  pmFootprint *fp, // Footprint containing mortal peaks
    57                                  const float nsigma_delta, // how many sigma above local background a peak
    58                                         // needs to be to survive
    59                                  const float min_threshold) { // minimum permitted coll height
    60     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    61     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    62     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    63     assert (fp != NULL);
    64 
    65     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    66         return PS_ERR_NONE;
    67     }
    68 
    69     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    70     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    71     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    72     const psImage *subImg = psImageSubset((psImage *)img, subRegion);
    73     const psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    74     assert (subImg != NULL && subWt != NULL);
    75     //
    76     // We need a psArray of peaks brighter than the current peak.  We'll fake this
    77     // by reusing the fp->peaks but lying about n.
    78     //
    79     // We do this for efficiency (otherwise I'd need two peaks lists), and we are
    80     // rather too chummy with psArray in consequence.  But it works.
    81     //
    82     psArray *brightPeaks = psArrayAlloc(0);
    83     psFree(brightPeaks->data);
    84     brightPeaks->data = psMemIncrRefCounter(fp->peaks->data);// use the data from fp->peaks
    85     //
    86     // The brightest peak is always safe; go through other peaks trying to cull them
    87     //
    88     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    89         const pmPeak *peak = fp->peaks->data[i];
    90         int x = peak->x - subImg->col0;
    91         int y = peak->y - subImg->row0;
    92         //
    93         // Find the level nsigma below the peak that must separate the peak
    94         // from any of its friends
    95         //
    96         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
    97         const float stdev = sqrt(subWt->data.F32[y][x]);
    98         float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    99         if (isnan(threshold) || threshold < min_threshold) {
    100 #if 1     // min_threshold is assumed to be below the detection threshold,
    101           // so all the peaks are pmFootprint, and this isn't the brightest
    102             // XXX mark peak to be dropped
    103             (void)psArrayRemoveIndex(fp->peaks, i);
    104             i--;                        // we moved everything down one
    105             continue;
    106 #else
    107 #error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
    108             threshold = min_threshold;
    109 #endif
    110         }
    111 
    112         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    113         if (threshold > subImg->data.F32[y][x]) {
    114             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    115         }
    116 
    117         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    118         // perhaps this should alloc a single ID image above and pass it in to be set.
    119 
    120         const int peak_id = 1;          // the ID for the peak of interest
    121         brightPeaks->n = i;             // only stop at a peak brighter than we are
    122 
    123         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    124 
    125         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    126         brightPeaks->n = 0;             // don't double free
    127         psImage *idImg = pmSetFootprintID(NULL, peakFootprint, peak_id);
    128         psFree(peakFootprint);
    129 
    130         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    131         // If so, the current peak is bogus; drop it.
    132         int j;
    133         for (j = 0; j < i; j++) {
    134             const pmPeak *peak2 = fp->peaks->data[j];
    135             int x2 = peak2->x - subImg->col0;
    136             int y2 = peak2->y - subImg->row0;
    137             const int peak2_id = idImg->data.S32[y2][x2]; // the ID for some other peak
    138 
    139             if (peak2_id == peak_id) {  // There's a brighter peak within the footprint above
    140                 ;                       // threshold; so cull our initial peak
    141                 (void)psArrayRemoveIndex(fp->peaks, i);
    142                 i--;                    // we moved everything down one
    143                 break;
    144             }
    145         }
    146         if (j == i) {
    147             j++;
    148         }
    149 
    150         psFree(idImg);
    151     }
    152 
    153     brightPeaks->n = 0; psFree(brightPeaks);
    154     psFree((psImage *)subImg);
    155     psFree((psImage *)subWt);
    156 
    157     return PS_ERR_NONE;
    158 }
    159 
    160  /*
    161   * Examine the peaks in a pmFootprint, and throw away the ones that are not sufficiently
    162   * isolated.  More precisely, for each peak find the highest coll that you'd have to traverse
    163   * to reach a still higher peak --- and if that coll's more than nsigma DN below your
    164   * starting point, discard the peak.
    165   */
    166 
    167 # define IN_PEAK 1
    168 psErrorCode pmFootprintCullPeaks(const psImage *img, // the image wherein lives the footprint
    169                                  const psImage *weight, // corresponding variance image
    170                                  pmFootprint *fp, // Footprint containing mortal peaks
    171                                  const float nsigma_delta, // how many sigma above local background a peak
    172                                  // needs to be to survive
    173                                  const float min_threshold) { // minimum permitted coll height
    174     assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    175     assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    176     assert (img->row0 == weight->row0 && img->col0 == weight->col0);
    177     assert (fp != NULL);
    178 
    179     if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    180         return PS_ERR_NONE;
    181     }
    182 
    183     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    184     subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    185     subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    186 
    187     psImage *subImg = psImageSubset((psImage *)img, subRegion);
    188     psImage *subWt = psImageSubset((psImage *)weight, subRegion);
    189     assert (subImg != NULL && subWt != NULL);
    190 
    191     psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
    192 
    193     // We need a psArray of peaks brighter than the current peak. 
    194     // We reject peaks which either:
    195     // 1) are below the local threshold
    196     // 2) have a brighter peak within their threshold
    197 
    198     // allocate the full-sized array.  if the final array is much smaller, we can realloc
    199     // at that point.
    200     psArray *brightPeaks = psArrayAllocEmpty(fp->peaks->n);
    201     psArrayAdd (brightPeaks, 128, fp->peaks->data[0]);
    202 
    203     // The brightest peak is always safe; go through other peaks trying to cull them
    204     for (int i = 1; i < fp->peaks->n; i++) { // n.b. fp->peaks->n can change within the loop
    205         const pmPeak *peak = fp->peaks->data[i];
    206         int x = peak->x - subImg->col0;
    207         int y = peak->y - subImg->row0;
    208         //
    209         // Find the level nsigma below the peak that must separate the peak
    210         // from any of its friends
    211         //
    212         assert (x >= 0 && x < subImg->numCols && y >= 0 && y < subImg->numRows);
    213         const float stdev = sqrt(subWt->data.F32[y][x]);
    214         float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    215         if (isnan(threshold) || threshold < min_threshold) {
    216             // min_threshold is assumed to be below the detection threshold,
    217             // so all the peaks are pmFootprint, and this isn't the brightest
    218             continue;
    219         }
    220 
    221         // XXX EAM : if stdev >= 0, i'm not sure how this can ever be true?
    222         if (threshold > subImg->data.F32[y][x]) {
    223             threshold = subImg->data.F32[y][x] - 10*FLT_EPSILON;
    224         }
    225 
    226         // XXX this is a bit expensive: psImageAlloc for every peak contained in this footprint
    227         // perhaps this should alloc a single ID image above and pass it in to be set.
    228 
    229         // XXX optionally use the faster pmFootprintsFind if the subimage size is large (eg, M31)
    230 
    231         // at this point brightPeaks only has the peaks brighter than the current
    232         pmFootprint *peakFootprint = pmFootprintsFindAtPoint(subImg, threshold, brightPeaks, peak->y, peak->x);
    233 
    234         // XXX need to supply the image here
    235         // we set the IDs to either 1 (in peak) or 0 (not in peak)
    236         pmSetFootprintID (idImg, peakFootprint, IN_PEAK);
    237         psFree(peakFootprint);
    238 
    239         // Check if any of the previous (brighter) peaks are within the footprint of this peak
    240         // If so, the current peak is bogus; drop it.
    241         bool keep = true;
    242         for (int j = 0; keep && (j < brightPeaks->n); j++) {
    243             const pmPeak *peak2 = fp->peaks->data[j];
    244             int x2 = peak2->x - subImg->col0;
    245             int y2 = peak2->y - subImg->row0;
    246             if (idImg->data.S32[y2][x2] == IN_PEAK)
    247                 // There's a brighter peak within the footprint above threshold; so cull our initial peak
    248                 keep = false;
    249         }
    250         if (!keep) continue;
    251 
    252         psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
    253     }
    254 
    255     psFree (fp->peaks);
    256     fp->peaks = brightPeaks;
    257 
    258     psFree(idImg);
    259     psFree(subImg);
    260     psFree(subWt);
    261 
    262     return PS_ERR_NONE;
    263 }
    264 
     48// psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n,
Note: See TracChangeset for help on using the changeset viewer.