Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 13430) +++ trunk/psphot/src/psphotFindPeaks.c (revision 13900) @@ -1,8 +1,8 @@ # include "psphotInternal.h" -// In this function, we smooth the image, then search for the peaks +// In this function, we smooth the image, then search for the peaks psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, - bool returnFootprints, - const int pass) { + bool returnFootprints, + const int pass, psMaskType maskVal) { float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK; @@ -13,57 +13,57 @@ if (pass == 1) { - SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA"); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); } else { - bool status_x, status_y; - float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); - float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); - if (!status_x | !status_y) { - psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined"); - return false; - } - SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); - NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); + bool status_x, status_y; + float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X"); + float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y"); + if (!status_x | !status_y) { + psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined"); + return false; + } + SIGMA_SMTH = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0))); + NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); } // smooth the image, applying the mask as we go psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_im, readout->mask, 0xff, SIGMA_SMTH, NSIGMA_SMTH); + psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot")); // smooth the weight, applying the mask as we go psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_wt, readout->mask, 0xff, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); + psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot")); psImage *mask = readout->mask; - // optionally save example images under trace + // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - char name[64]; - sprintf (name, "imsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); - sprintf (name, "wtsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_wt, name); + char name[64]; + sprintf (name, "imsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); + sprintf (name, "wtsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_wt, name); } // build the significance image on top of smooth_im for (int j = 0; j < smooth_im->numRows; j++) { - for (int i = 0; i < smooth_im->numCols; i++) { - float value = smooth_im->data.F32[j][i]; - if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || mask->data.U8[j][i]) { - smooth_im->data.F32[j][i] = 0.0; - } else { - smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; - } - } + for (int i = 0; i < smooth_im->numCols; i++) { + float value = smooth_im->data.F32[j][i]; + if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.U8[j][i] & maskVal)) { + smooth_im->data.F32[j][i] = 0.0; + } else { + smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; + } + } } psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot")); - // optionally save example images under trace + // optionally save example images under trace if (psTraceGetLevel("psphot") > 5) { - char name[64]; - sprintf (name, "snsmooth.v%d.fits", pass); - psphotSaveImage (NULL, smooth_im, name); + char name[64]; + sprintf (name, "snsmooth.v%d.fits", pass); + psphotSaveImage (NULL, smooth_im, name); } @@ -73,9 +73,9 @@ // signal/noise limit for the detected peaks if (pass == 1) { - NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); } else { - NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); - } - + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); + } + // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing // gaussian SIGMA. a peak in the significance image has an effective S/N for faint sources @@ -87,8 +87,8 @@ // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to // S/N = sqrt(Io)*4*pi*sigma_sm^2 - // thus, the threshold is: + // thus, the threshold is: float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); if (effArea < 1) { - effArea = 1; // never less than a pixel + effArea = 1; // never less than a pixel } float threshold = PS_SQR(NSIGMA_PEAK) / effArea; @@ -97,11 +97,11 @@ psArray *peaks = pmFindImagePeaks (smooth_im, threshold); if (peaks == NULL) { - // XXX this may also be due to a programming or config error - // XXX do we need to set something in the readout->analysis to indicate that - // we tried and failed to find peaks (something in the header data) - psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); - return false; + // XXX this may also be due to a programming or config error + // XXX do we need to set something in the readout->analysis to indicate that + // we tried and failed to find peaks (something in the header data) + psError(PSPHOT_ERR_DATA, false, "no peaks found in this image"); + return false; } - + // correct the peak values to S/N = sqrt(value*effArea) // get the peak flux from the unsmoothed image @@ -110,27 +110,27 @@ int col0 = readout->image->col0; for (int i = 0; i < peaks->n; i++) { - pmPeak *peak = peaks->data[i]; - peak->SN = sqrt(peak->value*effArea); - peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; + pmPeak *peak = peaks->data[i]; + peak->SN = sqrt(peak->value*effArea); + peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0]; } // limit the total number of returned peak as specified if (pass == 1) { - psArraySort (peaks, pmPeakSortBySN); - int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); - if (NMAX && (peaks->n > NMAX)) { - psArray *tmpPeaks = psArrayAllocEmpty (NMAX); - for (int i = 0; i < NMAX; i++) { - psArrayAdd (tmpPeaks, 100, peaks->data[i]); - } - psFree (peaks); - peaks = tmpPeaks; - } - } + psArraySort (peaks, pmPeakSortBySN); + int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); + if (NMAX && (peaks->n > NMAX)) { + psArray *tmpPeaks = psArrayAllocEmpty (NMAX); + for (int i = 0; i < NMAX; i++) { + psArrayAdd (tmpPeaks, 100, peaks->data[i]); + } + psFree (peaks); + peaks = tmpPeaks; + } + } // optional dump of all peak data char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE"); if (status && (output != NULL) && (output[0])) { - pmPeaksWriteText (peaks, output); + pmPeaksWriteText (peaks, output); } psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot")); @@ -142,22 +142,22 @@ // want to do that // - if (returnFootprints) { // We want an array of pmFootprint, not pmPeak - int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); - if (!status) { - npixMin = 1; - } - float FOOTPRINT_NSIGMA_LIMIT = - psMetadataLookupS32(&status, recipe, - (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2"); - if (!status) { - FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK; - } - threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; + if (returnFootprints) { // We want an array of pmFootprint, not pmPeak + int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN"); + if (!status) { + npixMin = 1; + } + float FOOTPRINT_NSIGMA_LIMIT = + psMetadataLookupS32(&status, recipe, + (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2"); + if (!status) { + FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK; + } + threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; - psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin); - pmPeaksAssignToFootprints(footprints, peaks); + psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin); + pmPeaksAssignToFootprints(footprints, peaks); - psFree(peaks); - peaks = footprints; // well, you know what I mean + psFree(peaks); + peaks = footprints; // well, you know what I mean } @@ -174,20 +174,20 @@ */ psErrorCode -psphotCullPeaks(const psImage *image, // the image wherein lives the footprint - const psImage *weight, // corresponding variance image - const psMetadata *recipe, - psArray *footprints) { // array of pmFootprints +psphotCullPeaks(const psImage *image, // the image wherein lives the footprint + const psImage *weight, // corresponding variance image + const psMetadata *recipe, + psArray *footprints) { // array of pmFootprints bool status = false; float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA"); if (!status) { - nsigma_delta = 0; // min. + nsigma_delta = 0; // min. } float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN"); if (!status) { - nsigma_min = 0; + nsigma_min = 0; } const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV"); return pmFootprintArrayCullPeaks(image, weight, footprints, - nsigma_delta, nsigma_min*skyStdev); + nsigma_delta, nsigma_min*skyStdev); }