IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 10, 2010, 7:34:39 PM (16 years ago)
Author:
eugene
Message:

updates from eam_branches/20091201 (substantially changes to the kernel matching code; updates to pmDetection container, added pmPattern, etc)

Location:
trunk/psModules/src/objects
Files:
2 deleted
16 edited
2 copied

Legend:

Unmodified
Added
Removed

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

    r23487 r26893  
    2626  psFree (detections->peaks);
    2727  psFree (detections->oldPeaks);
     28  psFree (detections->oldFootprints);
     29
     30  psFree (detections->newSources);
     31  psFree (detections->allSources);
    2832  return;
    2933}
    … …  
    3539    psMemSetDeallocator(detections, (psFreeFunc) pmDetectionsFree);
    3640
    37     detections->footprints = NULL;
    38     detections->peaks      = NULL;
    39     detections->oldPeaks   = NULL;
    40     detections->last       = 0;
     41    detections->footprints    = NULL;
     42    detections->peaks         = NULL;
     43    detections->oldPeaks      = NULL;
     44    detections->oldFootprints = NULL;
     45    detections->newSources    = NULL;
     46    detections->allSources    = NULL;
     47    detections->last          = 0;
    4148
    4249    return (detections);

  • trunk/psModules/src/objects/pmDetections.h

    r23487 r26893  
    2121typedef struct {
    2222  psArray *footprints;        // collection of footprints in the image
     23  psArray *oldFootprints;     // collection of footprints previously found
    2324  psArray *peaks;             // collection of all peaks contained by the footprints
    2425  psArray *oldPeaks;          // collection of all peaks previously found
     26  psArray *newSources;        // collection of sources
     27  psArray *allSources;        // collection of sources
    2528  int last;
    2629} pmDetections;

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

    r24888 r26893  
    2929  */
    3030
    31 # define IN_PEAK 1 
     31# define IN_PEAK 1
    3232psErrorCode 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
     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
    3838    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    3939    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    … …  
    4242
    4343    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    44         return PS_ERR_NONE;
    45     }
    46 
    47     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
     44        return PS_ERR_NONE;
     45    }
     46
     47    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    4848    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    4949    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    … …  
    5555    psImage *idImg = psImageAlloc(subImg->numCols, subImg->numRows, PS_TYPE_S32);
    5656
    57     // We need a psArray of peaks brighter than the current peak. 
     57    // We need a psArray of peaks brighter than the current peak.
    5858    // We reject peaks which either:
    5959    // 1) are below the local threshold
    … …  
    6767    // The brightest peak is always safe; go through other peaks trying to cull them
    6868    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);
    77 
    78         // const float stdev = sqrt(subWt->data.F32[y][x]);
    79         // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
    80 
    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
    86 
    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
     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);
     77
     78        // const float stdev = sqrt(subWt->data.F32[y][x]);
     79        // float threshold = subImg->data.F32[y][x] - nsigma_delta*stdev;
     80
     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
     86
     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
     92            continue;
     93        }
     94
     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        // If this peak has not already been assigned to a source, then we can look for any
     120        // brighter peaks within its footprint. Check if any of the previous (brighter) peaks
     121        // are within the footprint of this peak If so, the current peak is bogus; drop it.
     122        bool keep = true;
     123        for (int j = 0; keep && !peak->assigned && (j < brightPeaks->n); j++) {
     124            const pmPeak *peak2 = fp->peaks->data[j];
     125            int x2 = peak2->x - subImg->col0;
     126            int y2 = peak2->y - subImg->row0;
     127            if (idImg->data.S32[y2][x2] == IN_PEAK) {
     128                // There's a brighter peak within the footprint above threshold; so cull our initial peak
     129                keep = false;
     130            }
     131        }
     132        if (!keep) {
     133            psAssert (!peak->assigned, "logic error: trying to drop a previously-assigned peak");  // we should not drop any already assigned peaks.
    92134            continue;
    93135        }
    94136
    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]);
     137        psArrayAdd (brightPeaks, 128, fp->peaks->data[i]);
    133138    }
    134139
    … …  
    151156  */
    152157psErrorCode 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                                 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
    158163    assert (img != NULL); assert (img->type.type == PS_TYPE_F32);
    159164    assert (weight != NULL); assert (weight->type.type == PS_TYPE_F32);
    … …  
    162167
    163168    if (fp->peaks == NULL || fp->peaks->n == 0) { // nothing to do
    164         return PS_ERR_NONE;
    165     }
    166 
    167     psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
     169        return PS_ERR_NONE;
     170    }
     171
     172    psRegion subRegion;                 // desired subregion; 1 larger than bounding box (grr)
    168173    subRegion.x0 = fp->bbox.x0; subRegion.x1 = fp->bbox.x1 + 1;
    169174    subRegion.y0 = fp->bbox.y0; subRegion.y1 = fp->bbox.y1 + 1;
    … …  
    185190    //
    186191    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);
    195 
    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;
     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);
     200
     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;
    212217#else
    213218#error n.b. We will be running LOTS of checks at this threshold, so only find the footprint once
    214             threshold = min_threshold;
     219            threshold = min_threshold;
    215220#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);
     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);
    257262    }
    258263
    259264    brightPeaks->n = 0; psFree(brightPeaks);
    260     psFree((psImage *)subImg);
    261     psFree((psImage *)subWt);
     265    psFree(subImg);
     266    psFree(subWt);
    262267
    263268    return PS_ERR_NONE;

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

    r21183 r26893  
    2929 * so we set appropriate mask bits
    3030 *
    31  * EAM : these function were confusingly using "startspan" and "spartspan" 
     31 * EAM : these function were confusingly using "startspan" and "spartspan"
    3232 * I've rationalized them all to 'startspan'
    3333 */
    … …  
    3636// An enum for what we should do with a Startspan
    3737//
    38 typedef enum {PM_SSPAN_DOWN = 0,        // scan down from this span
    39               PM_SSPAN_UP,              // scan up from this span
    40               PM_SSPAN_RESTART,         // restart scanning from this span
    41               PM_SSPAN_DONE             // this span is processed
    42 } PM_SSPAN_DIR;                         // How to continue searching
     38typedef enum {PM_SSPAN_DOWN = 0,        // scan down from this span
     39              PM_SSPAN_UP,              // scan up from this span
     40              PM_SSPAN_RESTART,         // restart scanning from this span
     41              PM_SSPAN_DONE             // this span is processed
     42} PM_SSPAN_DIR;                         // How to continue searching
    4343//
    4444// An enum for mask's pixel values.  We're looking for pixels that are above threshold, and
    4545// we keep extra book-keeping information in the PM_SSPAN_STOP plane.  It's simpler to be
    46 // able to check for 
     46// able to check for
    4747//
    4848enum {
    49     PM_SSPAN_INITIAL = 0x0,             // initial state of pixels.
    50     PM_SSPAN_DETECTED = 0x1,            // we've seen this pixel
    51     PM_SSPAN_STOP = 0x2                 // you may stop searching when you see this pixel
     49    PM_SSPAN_INITIAL = 0x0,             // initial state of pixels.
     50    PM_SSPAN_DETECTED = 0x1,            // we've seen this pixel
     51    PM_SSPAN_STOP = 0x2                 // you may stop searching when you see this pixel
    5252};
    5353//
    … …  
    5555//
    5656typedef struct {
    57     const pmSpan *span;                 // save the pixel range
    58     PM_SSPAN_DIR direction;             // How to continue searching
    59     bool stop;                          // should we stop searching?
     57    const pmSpan *span;                 // save the pixel range
     58    PM_SSPAN_DIR direction;             // How to continue searching
     59    bool stop;                          // should we stop searching?
    6060} Startspan;
    6161
    6262static void startspanFree(Startspan *sspan) {
    63     psFree((void *)sspan->span);
     63    psFree(sspan->span);
    6464}
    6565
    6666static Startspan *
    67 StartspanAlloc(const pmSpan *span,      // The span in question
    68                psImage *mask,           // Pixels that we've already detected
    69                const PM_SSPAN_DIR dir   // Should we continue searching towards the top of the image?
     67StartspanAlloc(const pmSpan *span,      // The span in question
     68               psImage *mask,           // Pixels that we've already detected
     69               const PM_SSPAN_DIR dir   // Should we continue searching towards the top of the image?
    7070    ) {
    7171    Startspan *sspan = psAlloc(sizeof(Startspan));
    … …  
    7575    sspan->direction = dir;
    7676    sspan->stop = false;
    77    
    78     if (mask != NULL) {                 // remember that we've detected these pixels
    79         psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];
    80 
    81         for (int i = 0; i <= span->x1 - span->x0; i++) {
    82             mpix[i] |= PM_SSPAN_DETECTED;
    83             if (mpix[i] & PM_SSPAN_STOP) {
    84                 sspan->stop = true;
    85             }
    86         }
    87     }
    88    
     77
     78    if (mask != NULL) {                 // remember that we've detected these pixels
     79        psImageMaskType *mpix = &mask->data.PS_TYPE_IMAGE_MASK_DATA[span->y - mask->row0][span->x0 - mask->col0];
     80
     81        for (int i = 0; i <= span->x1 - span->x0; i++) {
     82            mpix[i] |= PM_SSPAN_DETECTED;
     83            if (mpix[i] & PM_SSPAN_STOP) {
     84                sspan->stop = true;
     85            }
     86        }
     87    }
     88
    8989    return sspan;
    9090}
    … …  
    9494//
    9595static bool add_startspan(psArray *startspans, // the saved Startspans
    96                           const pmSpan *sp, // the span in question
    97                           psImage *mask, // mask of detected/stop pixels
    98                           const PM_SSPAN_DIR dir) { // the desired direction to search
     96                          const pmSpan *sp, // the span in question
     97                          psImage *mask, // mask of detected/stop pixels
     98                          const PM_SSPAN_DIR dir) { // the desired direction to search
    9999    if (dir == PM_SSPAN_RESTART) {
    100         if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
    101             add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
    102             return true;
    103         }
     100        if (add_startspan(startspans, sp, mask,  PM_SSPAN_UP) ||
     101            add_startspan(startspans, sp, NULL, PM_SSPAN_DOWN)) {
     102            return true;
     103        }
    104104    } else {
    105         Startspan *sspan = StartspanAlloc(sp, mask, dir);
    106         if (sspan->stop) {              // we detected a stop bit
    107             psFree(sspan);              // don't allocate new span
    108 
    109             return true;
    110         } else {
    111             psArrayAdd(startspans, 1, sspan);
    112             psFree(sspan);              // as it's now owned by startspans
    113         }
     105        Startspan *sspan = StartspanAlloc(sp, mask, dir);
     106        if (sspan->stop) {              // we detected a stop bit
     107            psFree(sspan);              // don't allocate new span
     108
     109            return true;
     110        } else {
     111            psArrayAdd(startspans, 1, sspan);
     112            psFree(sspan);              // as it's now owned by startspans
     113        }
    114114    }
    115115
    … …  
    127127 */
    128128static bool do_startspan(pmFootprint *fp, // the footprint that we're building
    129                         const psImage *img, // the psImage we're working on
    130                         psImage *mask, // the associated masks
    131                          const float threshold, // Threshold
    132                          psArray *startspans) { // specify which span to process next
    133     bool F32 = false;                   // is this an F32 image?
     129                        const psImage *img, // the psImage we're working on
     130                        psImage *mask, // the associated masks
     131                         const float threshold, // Threshold
     132                         psArray *startspans) { // specify which span to process next
     133    bool F32 = false;                   // is this an F32 image?
    134134    if (img->type.type == PS_TYPE_F32) {
    135         F32 = true;
     135        F32 = true;
    136136    } else if (img->type.type == PS_TYPE_S32) {
    137         F32 = false;
    138     } else {                            // N.b. You can't trivially add more cases here; F32 is just a bool
    139         psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
    140         return NULL;
    141     }
    142 
    143     psF32 *imgRowF32 = NULL;            // row pointer if F32
    144     psS32 *imgRowS32 = NULL;            //  "   "   "  "  !F32
    145     psImageMaskType *maskRow = NULL;            //  masks's row pointer
    146    
     137        F32 = false;
     138    } else {                            // N.b. You can't trivially add more cases here; F32 is just a bool
     139        psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
     140        return NULL;
     141    }
     142
     143    psF32 *imgRowF32 = NULL;            // row pointer if F32
     144    psS32 *imgRowS32 = NULL;            //  "   "   "  "  !F32
     145    psImageMaskType *maskRow = NULL;            //  masks's row pointer
     146
    147147    const int row0 = img->row0;
    148148    const int col0 = img->col0;
    149149    const int numRows = img->numRows;
    150150    const int numCols = img->numCols;
    151    
     151
    152152    /********************************************************************************************************/
    153    
     153
    154154    Startspan *sspan = NULL;
    155155    for (int i = 0; i < startspans->n; i++) {
    156         sspan = startspans->data[i];
    157         if (sspan->direction != PM_SSPAN_DONE) {
    158             break;
    159         }
    160         if (sspan->stop) {
    161             break;
    162         }
     156        sspan = startspans->data[i];
     157        if (sspan->direction != PM_SSPAN_DONE) {
     158            break;
     159        }
     160        if (sspan->stop) {
     161            break;
     162        }
    163163    }
    164164    if (sspan == NULL || sspan->direction == PM_SSPAN_DONE) { // no more Startspans to process
    165         return false;
    166     }
    167     if (sspan->stop) {                  // they don't want any more spans processed
    168         return false;
     165        return false;
     166    }
     167    if (sspan->stop) {                  // they don't want any more spans processed
     168        return false;
    169169    }
    170170    /*
    … …  
    179179     * Go through image identifying objects
    180180     */
    181     int nx0, nx1 = -1;                  // new values of x0, x1
     181    int nx0, nx1 = -1;                  // new values of x0, x1
    182182    const int di = (dir == PM_SSPAN_UP) ? 1 : -1; // how much i changes to get to the next row
    183     bool stop = false;                  // should I stop searching for spans?
     183    bool stop = false;                  // should I stop searching for spans?
    184184
    185185    for (int i = sspan->span->y -row0 + di; i < numRows && i >= 0; i += di) {
    186         imgRowF32 = img->data.F32[i];   // only one of
    187         imgRowS32 = img->data.S32[i];   //      these is valid!
    188         maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i];
    189         //
    190         // Search left from the pixel diagonally to the left of (i - di, x0). If there's
    191         // a connected span there it may need to grow up and/or down, so push it onto
    192         // the stack for later consideration
    193         //
    194         nx0 = -1;
    195         for (int j = x0 - 1; j >= -1; j--) {
    196             double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
    197             if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
    198                 if (j < x0 - 1) {       // we found some pixels above threshold
    199                     nx0 = j + 1;
    200                 }
    201                 break;
    202             }
    203         }
    204 
    205         if (nx0 < 0) {                  // no span to the left
    206             nx1 = x0 - 1;               // we're going to resume searching at nx1 + 1
    207         } else {
    208             //
    209             // Search right in leftmost span
    210             //
    211             //nx1 = 0;                  // make gcc happy
    212             for (int j = nx0 + 1; j <= numCols; j++) {
    213                 double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
    214                 if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
    215                     nx1 = j - 1;
    216                     break;
    217                 }
    218             }
    219            
    220             const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
    221            
    222             if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
    223                 stop = true;
    224                 break;
    225             }
    226         }
    227         //
    228         // Now look for spans connected to the old span.  The first of these we'll
    229         // simply process, but others will have to be deferred for later consideration.
    230         //
    231         // In fact, if the span overhangs to the right we'll have to defer the overhang
    232         // until later too, as it too can grow in both directions
    233         //
    234         // Note that column numCols exists virtually, and always ends the last span; this
    235         // is why we claim below that sx1 is always set
    236         //
    237         bool first = false;             // is this the first new span detected?
    238         for (int j = nx1 + 1; j <= x1 + 1; j++) {
    239             double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
    240             if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
    241                 int sx0 = j++;          // span that we're working on is sx0:sx1
    242                 int sx1 = -1;           // We know that if we got here, we'll also set sx1
    243                 for (; j <= numCols; j++) {
    244                     double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
    245                     if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
    246                         sx1 = j;
    247                         break;
    248                     }
    249                 }
    250                 assert (sx1 >= 0);
    251 
    252                 const pmSpan *sp;
    253                 if (first) {
    254                     if (sx1 <= x1) {
    255                         sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
    256                         if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
    257                             stop = true;
    258                             break;
    259                         }
    260                     } else {            // overhangs to right
    261                         sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
    262                         if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
    263                             stop = true;
    264                             break;
    265                         }
    266                         sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
    267                         if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
    268                             stop = true;
    269                             break;
    270                         }
    271                     }
    272                     first = false;
    273                 } else {
    274                     sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
    275                     if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
    276                         stop = true;
    277                         break;
    278                     }
    279                 }
    280             }
    281         }
    282 
    283         if (stop || first == false) {   // we're done
    284             break;
    285         }
    286 
    287         x0 = nx0; x1 = nx1;
     186        imgRowF32 = img->data.F32[i];   // only one of
     187        imgRowS32 = img->data.S32[i];   //      these is valid!
     188        maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[i];
     189        //
     190        // Search left from the pixel diagonally to the left of (i - di, x0). If there's
     191        // a connected span there it may need to grow up and/or down, so push it onto
     192        // the stack for later consideration
     193        //
     194        nx0 = -1;
     195        for (int j = x0 - 1; j >= -1; j--) {
     196            double pixVal = (j < 0) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
     197            if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
     198                if (j < x0 - 1) {       // we found some pixels above threshold
     199                    nx0 = j + 1;
     200                }
     201                break;
     202            }
     203        }
     204
     205        if (nx0 < 0) {                  // no span to the left
     206            nx1 = x0 - 1;               // we're going to resume searching at nx1 + 1
     207        } else {
     208            //
     209            // Search right in leftmost span
     210            //
     211            //nx1 = 0;                  // make gcc happy
     212            for (int j = nx0 + 1; j <= numCols; j++) {
     213                double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
     214                if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) {
     215                    nx1 = j - 1;
     216                    break;
     217                }
     218            }
     219
     220            const pmSpan *sp = pmFootprintAddSpan(fp, i + row0, nx0 + col0, nx1 + col0);
     221
     222            if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
     223                stop = true;
     224                break;
     225            }
     226        }
     227        //
     228        // Now look for spans connected to the old span.  The first of these we'll
     229        // simply process, but others will have to be deferred for later consideration.
     230        //
     231        // In fact, if the span overhangs to the right we'll have to defer the overhang
     232        // until later too, as it too can grow in both directions
     233        //
     234        // Note that column numCols exists virtually, and always ends the last span; this
     235        // is why we claim below that sx1 is always set
     236        //
     237        bool first = false;             // is this the first new span detected?
     238        for (int j = nx1 + 1; j <= x1 + 1; j++) {
     239            double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
     240            if (!(maskRow[j] & PM_SSPAN_DETECTED) && pixVal >= threshold) {
     241                int sx0 = j++;          // span that we're working on is sx0:sx1
     242                int sx1 = -1;           // We know that if we got here, we'll also set sx1
     243                for (; j <= numCols; j++) {
     244                    double pixVal = (j >= numCols) ? threshold - 100 : (F32 ? imgRowF32[j] : imgRowS32[j]);
     245                    if ((maskRow[j] & PM_SSPAN_DETECTED) || pixVal < threshold) { // end of span
     246                        sx1 = j;
     247                        break;
     248                    }
     249                }
     250                assert (sx1 >= 0);
     251
     252                const pmSpan *sp;
     253                if (first) {
     254                    if (sx1 <= x1) {
     255                        sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
     256                        if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
     257                            stop = true;
     258                            break;
     259                        }
     260                    } else {            // overhangs to right
     261                        sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, x1 + col0);
     262                        if (add_startspan(startspans, sp, mask, PM_SSPAN_DONE)) {
     263                            stop = true;
     264                            break;
     265                        }
     266                        sp = pmFootprintAddSpan(fp, i + row0, x1 + 1 + col0, sx1 + col0 - 1);
     267                        if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
     268                            stop = true;
     269                            break;
     270                        }
     271                    }
     272                    first = false;
     273                } else {
     274                    sp = pmFootprintAddSpan(fp, i + row0, sx0 + col0, sx1 + col0 - 1);
     275                    if (add_startspan(startspans, sp, mask, PM_SSPAN_RESTART)) {
     276                        stop = true;
     277                        break;
     278                    }
     279                }
     280            }
     281        }
     282
     283        if (stop || first == false) {   // we're done
     284            break;
     285        }
     286
     287        x0 = nx0; x1 = nx1;
    288288    }
    289289    /*
    … …  
    309309 */
    310310pmFootprint *
    311 pmFootprintsFindAtPoint(const psImage *img,     // image to search
    312                        const float threshold,   // Threshold
    313                        const psArray *peaks, // array of peaks; finding one terminates search for footprint
    314                        int row, int col) { // starting position (in img's parent's coordinate system)
     311pmFootprintsFindAtPoint(const psImage *img,     // image to search
     312                       const float threshold,   // Threshold
     313                       const psArray *peaks, // array of peaks; finding one terminates search for footprint
     314                       int row, int col) { // starting position (in img's parent's coordinate system)
    315315   assert(img != NULL);
    316316
    317    bool F32 = false;                    // is this an F32 image?
     317   bool F32 = false;                    // is this an F32 image?
    318318   if (img->type.type == PS_TYPE_F32) {
    319319       F32 = true;
    320320   } else if (img->type.type == PS_TYPE_S32) {
    321321       F32 = false;
    322    } else {                             // N.b. You can't trivially add more cases here; F32 is just a bool
     322   } else {                             // N.b. You can't trivially add more cases here; F32 is just a bool
    323323       psError(PS_ERR_UNKNOWN, true, "Unsupported psImage type: %d", img->type.type);
    324324       return NULL;
    325325   }
    326    psF32 *imgRowF32 = NULL;             // row pointer if F32
    327    psS32 *imgRowS32 = NULL;             //  "   "   "  "  !F32
    328    
     326   psF32 *imgRowF32 = NULL;             // row pointer if F32
     327   psS32 *imgRowS32 = NULL;             //  "   "   "  "  !F32
     328
    329329   const int row0 = img->row0;
    330330   const int col0 = img->col0;
    … …  
    339339        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    340340                "row/col == (%d, %d) are out of bounds [%d--%d, %d--%d]",
    341                 row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
     341                row + row0, col + col0, row0, row0 + numRows - 1, col0, col0 + numCols - 1);
    342342       return NULL;
    343343   }
    … …  
    347347       return pmFootprintAlloc(0, img);
    348348   }
    349    
     349
    350350   pmFootprint *fp = pmFootprintAlloc(1 + img->numRows/10, img);
    351351/*
    … …  
    364364   if (peaks != NULL) {
    365365       for (int i = 0; i < peaks->n; i++) {
    366            pmPeak *peak = peaks->data[i];
    367            mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
     366           pmPeak *peak = peaks->data[i];
     367           mask->data.PS_TYPE_IMAGE_MASK_DATA[peak->y - mask->row0][peak->x - mask->col0] |= PM_SSPAN_STOP;
    368368       }
    369369   }
    … …  
    372372 */
    373373   psArray *startspans = psArrayAllocEmpty(1); // spans where we have to restart the search
    374    
    375    imgRowF32 = img->data.F32[row];      // only one of
    376    imgRowS32 = img->data.S32[row];      //      these is valid!
     374
     375   imgRowF32 = img->data.F32[row];      // only one of
     376   imgRowS32 = img->data.S32[row];      //      these is valid!
    377377   psImageMaskType *maskRow = mask->data.PS_TYPE_IMAGE_MASK_DATA[row];
    378378   {
    379379       int i;
    380380       for (i = col; i >= 0; i--) {
    381            pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
    382            if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
    383                break;
    384            }
     381           pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
     382           if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
     383               break;
     384           }
    385385       }
    386386       int i0 = i;
    387387       for (i = col; i < numCols; i++) {
    388            pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
    389            if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
    390                break;
    391            }
     388           pixVal = F32 ? imgRowF32[i] : imgRowS32[i];
     389           if ((maskRow[i] & PM_SSPAN_DETECTED) || pixVal < threshold) {
     390               break;
     391           }
    392392       }
    393393       int i1 = i;
    … …  
    404404    */
    405405   psFree(mask);
    406    psFree(startspans);                  // restores the image pixel
    407 
    408    return fp;                           // pmFootprint really
     406   psFree(startspans);                  // restores the image pixel
     407
     408   return fp;                           // pmFootprint really
    409409}

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

    r24206 r26893  
    4242#include "pmErrorCodes.h"
    4343
     44
     45#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
     46
    4447/*****************************************************************************/
    4548/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    … …  
    405408    return psf;
    406409}
     410
     411
     412float pmPSFtoFWHM(const pmPSF *psf, float x, float y)
     413{
     414    PS_ASSERT_PTR_NON_NULL(psf, NAN);
     415
     416    pmModel *model = pmModelFromPSFforXY(psf, x, y, 1.0); // Model of source
     417    if (!model) {
     418        psError(PS_ERR_UNKNOWN, false, "Unable to determine PSF model at %f,%f\n", x, y);
     419        return NAN;
     420    }
     421    psF32 *params = model->params->data.F32; // Model parameters
     422    psEllipseAxes axes = pmPSF_ModelToAxes(params, MAX_AXIS_RATIO); // Ellipse axes
     423
     424    // Curiously, the minor axis can be larger than the major axis, so need to check.
     425    float fwhm = 2.355 * PS_MAX(axes.minor, axes.major); // FWHM, converted from sigma
     426
     427    psFree(model);
     428
     429    return fwhm;
     430}

  • trunk/psModules/src/objects/pmPSF.h

    r25754 r26893  
    111111psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR);
    112112
     113/// Calculate FWHM value from a PSF
     114float pmPSFtoFWHM(
     115    const pmPSF *psf,                   // PSF of interest
     116    float x, float y                    // Position of interest
     117    );
     118
     119
    113120/// @}
    114121# endif

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

    r26391 r26893  
    416416
    417417        } else {
    418             psError(PS_ERR_UNKNOWN, true, "peak specified valid column range.");
     418            psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
    419419        }
    420420    }
    … …  
    501501                }
    502502            } else {
    503                 psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
     503                psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
    504504            }
    505505
    … …  
    545545            }
    546546        } else {
    547             psError(PS_ERR_UNKNOWN, true, "peak specified outside valid column range.");
     547            psLogMsg ("psModules.objects", 5, "peak specified outside valid column range.");
    548548        }
    549549    }

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

    r26395 r26893  
    33 *  Functions to define and manipulate sources on images
    44 *
    5  *  @author GLG, MHPCC
    6  *  @author EAM, IfA: significant modifications.
     5 *  @author EAM, IfA
     6 *  @author GLG, MHPCC (initial code base)
    77 *
    88 *  @version $Revision: 1.70 $ $Name: not supported by cvs2svn $
    99 *  @date $Date: 2009-02-16 22:29:59 $
    10  *
    11  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    12  *
     10 *  Copyright 2009 Institute for Astronomy, University of Hawaii
    1311 */
    1412
    … …  
    277275// psphot-specific function which applies the recipe values
    278276// only apply selection to sources within specified region
    279 pmPSFClump pmSourcePSFClump(psRegion *region, psArray *sources, psMetadata *recipe)
     277pmPSFClump pmSourcePSFClump(psImage **savedImage, psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_GRID_SCALE, psF32 SX_MAX, psF32 SY_MAX, psF32 AR_MAX)
    280278{
    281279    psTrace("psModules.objects", 10, "---- begin ----\n");
    … …  
    287285
    288286    PS_ASSERT_PTR_NON_NULL(sources, errorClump);
    289     PS_ASSERT_PTR_NON_NULL(recipe, errorClump);
    290 
    291     bool status = true;                 // Status of MD lookup
    292     float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
    293     if (!status) {
    294         PSF_SN_LIM = 0;
    295     }
    296     float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
    297     if (!status) {
    298         PSF_CLUMP_GRID_SCALE = 0.1;
    299     }
    300287
    301288    // find the sigmaX, sigmaY clump
    302289    {
    303         psF32 SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
    304         if (!status) {
    305             psWarning("MOMENTS_SX_MAX not set in recipe");
    306             SX_MAX = 10.0;
    307         }
    308         psF32 SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
    309         if (!status) {
    310             psWarning("MOMENTS_SY_MAX not set in recipe");
    311             SY_MAX = 10.0;
    312         }
    313         psF32 AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX");
    314         if (!status) {
    315             psWarning("MOMENTS_AR_MAX not set in recipe");
    316             AR_MAX =  3.0;
    317         }
    318290        psF32 AR_MIN = 1.0 / AR_MAX;
    319291
    … …  
    401373        psfClump.nSigma = stats->sampleStdev;
    402374
    403         const bool keep_psf_clump = psMetadataLookupBool(NULL, recipe, "KEEP_PSF_CLUMP");
    404         if (keep_psf_clump)
    405         {
    406             psMetadataAdd(recipe, PS_LIST_TAIL,
    407                           "PSF_CLUMP", PS_DATA_IMAGE, "Image of PSF coefficients", splane);
     375        if (savedImage) {
     376            *savedImage = psMemIncrRefCounter(splane);
    408377        }
    409378        psFree (splane);
    … …  
    411380
    412381        // if we failed to find a valid peak, return the empty clump (failure signal)
    413         if (!peaks || !peaks->n) {
     382        if (peaks == NULL) {
     383            psError(PS_ERR_UNKNOWN, false, "failure in peak analysis for PSF clump.\n");
     384            psFree (peaks);
     385            return emptyClump;
     386        }
     387
     388        if (peaks->n == 0)
     389        {
    414390            psLogMsg ("psphot", 3, "failed to find a peak in the PSF clump image\n");
    415391            if (nValid == 0) {
    … …  
    418394                psLogMsg ("psphot", 3, "no significant peak\n");
    419395            }
     396            psFree (peaks);
    420397            return (emptyClump);
    421398        }
    … …  
    524501*****************************************************************************/
    525502
    526 bool pmSourceRoughClass(psRegion *region, psArray *sources, psMetadata *recipe, pmPSFClump clump, psImageMaskType maskSat)
     503bool pmSourceRoughClass(psRegion *region, psArray *sources, float PSF_SN_LIM, float PSF_CLUMP_NSIGMA, pmPSFClump clump, psImageMaskType maskSat)
    527504{
    528505    psTrace("psModules.objects", 10, "---- begin ----");
    529506
    530507    PS_ASSERT_PTR_NON_NULL(sources, false);
    531     PS_ASSERT_PTR_NON_NULL(recipe, false);
    532508
    533509    int Nsat     = 0;
    … …  
    542518    psVector *starsn_peaks = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    543519    psVector *starsn_moments = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    544 
    545     // get basic parameters, or set defaults
    546     bool status;
    547     float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM");
    548     if (!status) PSF_SN_LIM = 20.0;
    549     float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA");
    550     if (!status) PSF_CLUMP_NSIGMA = 1.5;
    551 
    552     // float INNER_RADIUS = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    553520
    554521    pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN;

  • trunk/psModules/src/objects/pmSource.h

    r25754 r26893  
    176176 *
    177177 * The return value indicates the success (TRUE) of the operation.
    178  *
    179  * XXX: Limit the S/N of the candidate sources (part of Metadata)? (TBD).
    180  * XXX: Save the clump parameters on the Metadata (TBD)
    181  *
    182  */
     178 */
     179
    183180pmPSFClump pmSourcePSFClump(
     181    psImage **savedImage,
    184182    psRegion *region,                   ///< restrict measurement to specified region
    185183    psArray *source,                    ///< The input pmSource
    186     psMetadata *metadata                ///< Contains classification parameters
     184    float PSF_SN_LIM,
     185    float PSF_CLUMP_GRID_SCALE,
     186    psF32 SX_MAX,
     187    psF32 SY_MAX,
     188    psF32 AR_MAX
    187189);
    188190
    … …  
    200202    psRegion *region,                   ///< restrict measurement to specified region
    201203    psArray *sources,                    ///< The input pmSources
    202     psMetadata *metadata,               ///< Contains classification parameters
     204    float PSF_SN_LIM,                    ///< min S/N for source to be used for PSF model
     205    float PSF_CLUMP_NSIGMA,              ///< size of region around peak of clump for PSF stars
    203206    pmPSFClump clump,                   ///< Statistics about the PSF clump
    204207    psImageMaskType maskSat             ///< Mask value for saturated pixels
    … …  
    220223    float radius,     ///< Use a circle of pixels around the peak
    221224    float sigma,      ///< size of Gaussian window function (<= 0.0 -> skip window)
    222     float minSN       ///< minimum pixel significance
     225    float minSN,              ///< minimum pixel significance
     226    psImageMaskType maskVal
    223227);
    224228

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

    r26393 r26893  
    4040#include "pmPSF.h"
    4141#include "pmModel.h"
     42#include "pmDetections.h"
    4243#include "pmSource.h"
    4344#include "pmModelClass.h"
    … …  
    344345
    345346    // if sources is NULL, write out an empty table
    346     // input / output sources are stored on the readout->analysis as "PSPHOT.SOURCES" -- a better name might be something like PM_SOURCE_DATA
    347     psArray *sources = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
     347    // input / output sources are stored on the readout->analysis as "PSPHOT.DETECTIONS"
     348
     349    psArray *sources = NULL;
     350    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     351    if (detections) {
     352        sources = detections->allSources;
     353    }
    348354    if (!sources) {
     355        detections = pmDetectionsAlloc();
    349356        sources = psArrayAlloc(0);
    350         psMetadataAddArray(readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "Blank array of sources", sources);
    351         psFree(sources); // Held onto by the metadata, so we can continue to use
     357        detections->allSources = sources;
     358        psMetadataAddPtr(readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_UNKNOWN | PS_META_REPLACE, "Blank array of sources", detections);
     359        psFree(detections); // Held onto by the metadata, so we can continue to use
    352360    }
    353361
    … …  
    10311039    }
    10321040    readout->data_exists = true;
    1033     status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "input sources", sources);
    1034     psFree (sources);
     1041
     1042    pmDetections *detections = pmDetectionsAlloc();
     1043    detections->allSources = sources;
     1044    status = psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY, "input sources", detections);
     1045    psFree (detections);
    10351046    return true;
    10361047}
    … …  
    11241135    bool status;
    11251136
    1126     // select the psf of interest
    1127     pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.SOURCES");
    1128     if (!psf) return false;
    1129     return true;
    1130 }
    1131 
    1132 
     1137    // select the detections of interest
     1138    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     1139    if (!detections) return false;
     1140    if (!detections->allSources) return false;
     1141    return true;
     1142}
     1143
     1144

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

    r25754 r26893  
    304304        source->peak->dx   = dPAR[PM_PAR_XPOS];
    305305        source->peak->dy   = dPAR[PM_PAR_YPOS];
     306        source->peak->SN   = sqrt(source->peak->flux); // XXX a proxy: various functions sort by peak S/N
    306307
    307308        source->pixWeight = psMetadataLookupF32 (&status, row, "PSF_QF");

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

    r26392 r26893  
    163163
    164164    // try find the MATCHED_REFS extension.  if non-existent, note that we tried, and move on.
     165    // It is not an error to lack this entry -- psFitsMoveExtNameClean does not raise an error
    165166    if (!psFitsMoveExtNameClean (fits, "MATCHED_REFS")) {
    166167        psMetadataAddBool (fpa->analysis, PS_LIST_TAIL, "READ.REFMATCH", PS_META_REPLACE, "attempted to read MATCHED_REFS", true);

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

    r25979 r26893  
    5454# define VALID_RADIUS(X,Y,RAD2) (((RAD2) >= (PS_SQR(X) + PS_SQR(Y))) ? 1 : 0)
    5555
    56 bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN)
     56bool pmSourceMoments(pmSource *source, psF32 radius, psF32 sigma, psF32 minSN, psImageMaskType maskVal)
    5757{
    5858    PS_ASSERT_PTR_NON_NULL(source, false);
    … …  
    114114        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    115115            if (vMsk) {
    116                 if (*vMsk) {
     116                if (*vMsk & maskVal) {
    117117                    vMsk++;
    118118                    continue;
    … …  
    135135            // stars.
    136136            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    137             if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
     137            // if (pDiff < 0) continue; // XXX : MWV says I should include < 0.0 valued points...
    138138
    139139            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    … …  
    226226        for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
    227227            if (vMsk) {
    228                 if (*vMsk) {
     228                if (*vMsk & maskVal) {
    229229                    vMsk++;
    230230                    continue;
    … …  
    249249            // stars.
    250250            if (PS_SQR(pDiff) < minSN2*wDiff) continue;
    251             if (pDiff < 0) continue;
     251            // if (pDiff < 0) continue;
    252252
    253253            // Apply a Gaussian window function.  Be careful with the window function.  S/N
    … …  
    315315    source->moments->Myyyy = YYYY/Sum;
    316316
    317     if (source->moments->Mxx < 0) {
    318         fprintf (stderr, "error: neg second moment??\n");
    319     }
    320     if (source->moments->Myy < 0) {
    321         fprintf (stderr, "error: neg second moment??\n");
    322     }
     317    // if (source->moments->Mxx < 0) {
     318    // fprintf (stderr, "error: neg second moment??\n");
     319    // }
     320    // if (source->moments->Myy < 0) {
     321    // fprintf (stderr, "error: neg second moment??\n");
     322    // }
    323323
    324324    psTrace ("psModules.objects", 4, "Mxx: %f  Mxy: %f  Myy: %f  Mxxx: %f  Mxxy: %f  Mxyy: %f  Myyy: %f  Mxxxx: %f  Mxxxy: %f  Mxxyy: %f  Mxyyy: %f  Mxyyy: %f\n",

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

    r20937 r26893  
    3434#include "pmPSF.h"
    3535#include "pmModel.h"
     36#include "pmDetections.h"
    3637#include "pmSource.h"
    3738#include "pmSourcePlots.h"
    … …  
    5354    PS_ASSERT_PTR_NON_NULL(layout, false);
    5455
     56    bool status;
    5557    Graphdata graphdata;
    5658    KapaSection section;
    … …  
    6163    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6264
    63     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    64     if (sources == NULL)
    65         return false;
     65    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     66    if (detections == NULL) return false;
     67
     68    psArray *sources = detections->allSources;
     69    if (sources == NULL) return false;
    6670
    6771    int kapa = pmKapaOpen (false);

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

    r20937 r26893  
    3737#include "pmPSF.h"
    3838#include "pmModel.h"
     39#include "pmDetections.h"
    3940#include "pmSource.h"
    4041#include "pmSourcePlots.h"
    … …  
    5455    PS_ASSERT_PTR_NON_NULL(layout, false);
    5556
     57    bool status;
    5658    Graphdata graphdata;
    5759    KapaSection section;
    … …  
    6264    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6365
    64     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    65     if (sources == NULL)
    66         return false;
     66    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     67    if (detections == NULL) return false;
     68
     69    psArray *sources = detections->allSources;
     70    if (sources == NULL) return false;
    6771
    6872    int kapa = pmKapaOpen (false);

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

    r20937 r26893  
    3737#include "pmPSF.h"
    3838#include "pmModel.h"
     39#include "pmDetections.h"
    3940#include "pmSource.h"
    4041#include "pmSourcePlots.h"
    … …  
    5657    PS_ASSERT_PTR_NON_NULL(layout, false);
    5758
     59    bool status;
    5860    Graphdata graphdata;
    5961    KapaSection section;
    … …  
    6466    pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    6567
    66     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
    67     if (sources == NULL)
    68         return false;
     68    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     69    if (detections == NULL) return false;
     70
     71    psArray *sources = detections->allSources;
     72    if (sources == NULL) return false;
    6973
    7074    int kapa = pmKapaOpen (false);
Note: See TracChangeset for help on using the changeset viewer.