IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 17444


Ignore:
Timestamp:
Apr 12, 2008, 10:07:26 PM (18 years ago)
Author:
eugene
Message:

re-org of psphotFindDetections into significance, peaks, footprints; find and cull footprints in one function

Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotFindDetections.c

    r17396 r17444  
    11# include "psphotInternal.h"
    22
    3 // XXX let's make a better distinction between 'pass' (used for counting intput / output
    4 // entries) and needing/having a PSF.  In this function, we smooth the image, then search for
    5 // the peaks
     3// smooth the image, search for peaks, optionally define footprints based on the peaks
    64pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe) {
    75
    … …  
    2725
    2826    // move the old peaks array (if it exists) to oldPeaks
     27    // XXX generically, we should be able to call this function an arbitrary number of times
    2928    assert (detections->oldPeaks == NULL);
    3029    detections->oldPeaks = detections->peaks;
    3130    detections->peaks = NULL;
    3231
    33     // if footprints are not requested, find the peaks and return them
    34     if (!useFootprints) {
    35         detections->peaks = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal);
    36         return detections;
     32    // generate the smoothed significance image
     33    psImage *significance = psphotSignificanceImage (readout, recipe, pass, maskVal);
     34
     35    // detect the peaks in the significance image
     36    detections->peaks = psphotFindPeaks (significance, readout, recipe, pass, maskVal);
     37
     38    // optionally merge peaks into footprints
     39    if (useFootprints) {
     40        psphotFindFootprints (detections, significance, readout, recipe, pass, maskVal);
    3741    }
    3842
    39     // psphotFindPeaks returns a different set of objects if useFootprints is true!
    40     // XXX fix this by splitting psphotFindPeaks into smooth / peaks / footprint stages?
    41     psArray *footprints = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal);
    42 
    43     int growRadius = 0;
    44     if (pass == 1) {
    45         growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS");
    46     } else {
    47         growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2");
    48     }
    49 
    50     if (growRadius > 0) {
    51         psArray *tmp = pmGrowFootprintArray(footprints, growRadius);
    52         psFree(footprints);
    53         footprints = tmp;
    54     }
    55 
    56     if (pass == 2) {
    57         // merge in old peaks;
    58         const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints
    59        
    60         psLogMsg ("psphot", PS_LOG_INFO, "merging %ld new footprints into %ld existing ones\n", footprints->n, detections->footprints->n);
    61         psArray *mergedFootprints = pmMergeFootprintArrays(detections->footprints, footprints, includePeaks);
    62         psFree(footprints);
    63         psFree(detections->footprints);
    64         detections->footprints = mergedFootprints;
    65     } else {
    66         detections->footprints = footprints;
    67     }
    68 
    69     psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints);
    70     detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
    71     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot"));
    72 
     43    psFree (significance);
    7344    return detections;
    7445}

  • trunk/psphot/src/psphotFindPeaks.c

    r17396 r17444  
    33// XXX let's make a better distinction between 'pass' and needing/having a PSF.
    44// In this function, we smooth the image, then search for the peaks
    5 psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe,
    6                           bool returnFootprints,
    7                           const int pass, psMaskType maskVal) {
    85
    9     float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK;
    10     bool  status = false;
     6psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) {
     7
     8    float NSIGMA_PEAK;
     9    bool status = false;
    1110
    1211    // smooth the image and weight map
    1312    psTimerStart ("peaks");
    1413
    15     // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values
    16     if (pass == 1) {
    17         SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
    18         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    19     } else {
    20         bool status_x, status_y;
    21         float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
    22         float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
    23         if (!status_x | !status_y) {
    24             psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
    25             return false;
    26         }
    27         SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
    28         NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
    29     }
    30 
    31     // smooth the image, applying the mask as we go
    32     psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
    33     psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH);
    34     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks"));
    35 
    36     // smooth the weight, applying the mask as we go
    37     psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
    38     psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
    39     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks"));
    40 
    41     psImage *mask = readout->mask;
    42 
    43     // optionally save example images under trace
    44     if (psTraceGetLevel("psphot") > 5) {
    45         char name[64];
    46         sprintf (name, "imsmooth.v%d.fits", pass);
    47         psphotSaveImage (NULL, smooth_im, name);
    48         sprintf (name, "wtsmooth.v%d.fits", pass);
    49         psphotSaveImage (NULL, smooth_wt, name);
    50     }
    51 
    52     // build the significance image on top of smooth_im
    53     for (int j = 0; j < smooth_im->numRows; j++) {
    54         for (int i = 0; i < smooth_im->numCols; i++) {
    55             float value = smooth_im->data.F32[j][i];
    56             if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) {
    57                 smooth_im->data.F32[j][i] = 0.0;
    58             } else {
    59                 smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
    60             }
    61         }
    62     }
    63     psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("peaks"));
    64 
    65     // optionally save example images under trace
    66     if (psTraceGetLevel("psphot") > 5) {
    67         char name[64];
    68         sprintf (name, "snsmooth.v%d.fits", pass);
    69         psphotSaveImage (NULL, smooth_im, name);
    70     }
    71 
    72     psTimerStart ("peaks");
    73     // set peak threshold
    74 
    7514    // signal/noise limit for the detected peaks
     15    // XXX this is a little lame: can we do something better than 'pass'?
    7616    if (pass == 1) {
    7717        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
    … …  
    7919        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
    8020    }
     21    PS_ASSERT (status, NULL);
     22
     23    float SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "SIGMA_SMOOTH");
     24    psAssert (status, "SIGMA_SMOOTH missing: call psphotSignificanceImage first");
    8125
    8226    // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing
    … …  
    8731
    8832    // if sigma_sm = sigma_obs, then sigma_eff^2 = 2 sigma_sm^2. in this case, the threshold in
    89     // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to
     33    // terms of significance peak counts Io, for a desired S/N limit corresponds to
    9034    // S/N = sqrt(Io)*4*pi*sigma_sm^2
    9135    // thus, the threshold is:
    … …  
    10044    float threshold = PS_SQR(NSIGMA_PEAK) / effArea;
    10145
     46    // record the actual effective area and peak threshold
     47    psMetadataAddF32  (recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
     48    psMetadataAddF32  (recipe, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
     49
    10250    // find the peaks in the smoothed image
    103     psArray *peaks = pmPeaksInImage (smooth_im, threshold);
     51    psArray *peaks = pmPeaksInImage (significance, threshold);
    10452    if (peaks == NULL) {
     53        // XXX should we be sending back an empty array instead of NULL?
    10554        // XXX this may also be due to a programming or config error
    10655        // XXX do we need to set something in the readout->analysis to indicate that
    … …  
    13584    }
    13685
    137     // optional dump of all peak data
     86    // optional dump of all peak data
     87    // XXX need a better check for this option
    13888    char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE");
    13989    if (status && (output != NULL) && (output[0])) {
    … …  
    14292    psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks"));
    14393
    144     // If they asked us to return a set of pmFootprints, not a raw pmPeak list,
    145     // go ahead and find the footprints, and assign them their peaks.
    146     //
    147     // N.b. We're not culling this list; call pmFootprintCullPeaks if you
    148     // want to do that
    149     //
    150     if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
    151         int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
    152         if (!status) {
    153             npixMin = 1;
    154         }
    155         float FOOTPRINT_NSIGMA_LIMIT =
    156             psMetadataLookupS32(&status, recipe,
    157                                 (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
    158         if (!status) {
    159             FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
    160         }
    161         threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;
     94    return peaks;
    16295
    163         psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
    164         pmPeaksAssignToFootprints(footprints, peaks);
    165 
    166         psFree(peaks);
    167         peaks = footprints;             // well, you know what I mean
    168 
    169         psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks"));
    170     }
    171 
    172     psFree (smooth_im);
    173     psFree (smooth_wt);
    174 
    175     return (peaks);
    17696}
    17797
Note: See TracChangeset for help on using the changeset viewer.