Index: trunk/psphot/src/psphotFindDetections.c =================================================================== --- trunk/psphot/src/psphotFindDetections.c (revision 17443) +++ trunk/psphot/src/psphotFindDetections.c (revision 17444) @@ -1,7 +1,5 @@ # include "psphotInternal.h" -// XXX let's make a better distinction between 'pass' (used for counting intput / output -// entries) and needing/having a PSF. In this function, we smooth the image, then search for -// the peaks +// smooth the image, search for peaks, optionally define footprints based on the peaks pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe) { @@ -27,48 +25,21 @@ // move the old peaks array (if it exists) to oldPeaks + // XXX generically, we should be able to call this function an arbitrary number of times assert (detections->oldPeaks == NULL); detections->oldPeaks = detections->peaks; detections->peaks = NULL; - // if footprints are not requested, find the peaks and return them - if (!useFootprints) { - detections->peaks = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal); - return detections; + // generate the smoothed significance image + psImage *significance = psphotSignificanceImage (readout, recipe, pass, maskVal); + + // detect the peaks in the significance image + detections->peaks = psphotFindPeaks (significance, readout, recipe, pass, maskVal); + + // optionally merge peaks into footprints + if (useFootprints) { + psphotFindFootprints (detections, significance, readout, recipe, pass, maskVal); } - // psphotFindPeaks returns a different set of objects if useFootprints is true! - // XXX fix this by splitting psphotFindPeaks into smooth / peaks / footprint stages? - psArray *footprints = psphotFindPeaks (readout, recipe, useFootprints, pass, maskVal); - - int growRadius = 0; - if (pass == 1) { - growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS"); - } else { - growRadius = psMetadataLookupS32(NULL, recipe, "FOOTPRINT_GROW_RADIUS_2"); - } - - if (growRadius > 0) { - psArray *tmp = pmGrowFootprintArray(footprints, growRadius); - psFree(footprints); - footprints = tmp; - } - - if (pass == 2) { - // merge in old peaks; - const int includePeaks = 0x0 | 0x2; // i.e. just from newFootprints - - psLogMsg ("psphot", PS_LOG_INFO, "merging %ld new footprints into %ld existing ones\n", footprints->n, detections->footprints->n); - psArray *mergedFootprints = pmMergeFootprintArrays(detections->footprints, footprints, includePeaks); - psFree(footprints); - psFree(detections->footprints); - detections->footprints = mergedFootprints; - } else { - detections->footprints = footprints; - } - - psphotCullPeaks(readout->image, readout->weight, recipe, detections->footprints); - detections->peaks = pmFootprintArrayToPeaks(detections->footprints); - psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot")); - + psFree (significance); return detections; } Index: trunk/psphot/src/psphotFindPeaks.c =================================================================== --- trunk/psphot/src/psphotFindPeaks.c (revision 17443) +++ trunk/psphot/src/psphotFindPeaks.c (revision 17444) @@ -3,75 +3,15 @@ // XXX let's make a better distinction between 'pass' and needing/having a PSF. // In this function, we smooth the image, then search for the peaks -psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe, - bool returnFootprints, - const int pass, psMaskType maskVal) { - float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK; - bool status = false; +psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) { + + float NSIGMA_PEAK; + bool status = false; // smooth the image and weight map psTimerStart ("peaks"); - // XXX if we have been supplied a PSF, we can use that to set the FWHM_X,FWHM_Y values - if (pass == 1) { - 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"); - } - - // smooth the image, applying the mask as we go - psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("peaks")); - - // smooth the weight, applying the mask as we go - psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); - psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); - psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("peaks")); - - psImage *mask = readout->mask; - - // 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); - } - - // 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] & 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 ("peaks")); - - // 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); - } - - psTimerStart ("peaks"); - // set peak threshold - // signal/noise limit for the detected peaks + // XXX this is a little lame: can we do something better than 'pass'? if (pass == 1) { NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); @@ -79,4 +19,8 @@ NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); } + PS_ASSERT (status, NULL); + + float SIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "SIGMA_SMOOTH"); + psAssert (status, "SIGMA_SMOOTH missing: call psphotSignificanceImage first"); // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing @@ -87,5 +31,5 @@ // if sigma_sm = sigma_obs, then sigma_eff^2 = 2 sigma_sm^2. in this case, the threshold in - // terms of smooth_im peak counts Io, for a desired S/N limit corresponds to + // terms of significance peak counts Io, for a desired S/N limit corresponds to // S/N = sqrt(Io)*4*pi*sigma_sm^2 // thus, the threshold is: @@ -100,7 +44,12 @@ float threshold = PS_SQR(NSIGMA_PEAK) / effArea; + // record the actual effective area and peak threshold + psMetadataAddF32 (recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); + // find the peaks in the smoothed image - psArray *peaks = pmPeaksInImage (smooth_im, threshold); + psArray *peaks = pmPeaksInImage (significance, threshold); if (peaks == NULL) { + // XXX should we be sending back an empty array instead of 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 @@ -135,5 +84,6 @@ } - // optional dump of all peak data + // optional dump of all peak data + // XXX need a better check for this option char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE"); if (status && (output != NULL) && (output[0])) { @@ -142,36 +92,6 @@ psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("peaks")); - // If they asked us to return a set of pmFootprints, not a raw pmPeak list, - // go ahead and find the footprints, and assign them their peaks. - // - // N.b. We're not culling this list; call pmFootprintCullPeaks if you - // 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; + return peaks; - psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin); - pmPeaksAssignToFootprints(footprints, peaks); - - psFree(peaks); - peaks = footprints; // well, you know what I mean - - psLogMsg ("psphot", PS_LOG_INFO, "%ld footprints: %f sec\n", peaks->n, psTimerMark ("peaks")); - } - - psFree (smooth_im); - psFree (smooth_wt); - - return (peaks); }