Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 17869) +++ /trunk/psphot/src/psphot.h (revision 17870) @@ -51,5 +51,5 @@ // used by psphotFindDetections psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); -psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); +psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax); bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal); psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints); Index: /trunk/psphot/src/psphotFindDetections.c =================================================================== --- /trunk/psphot/src/psphotFindDetections.c (revision 17869) +++ /trunk/psphot/src/psphotFindDetections.c (revision 17870) @@ -6,4 +6,6 @@ bool status; int pass; + float NSIGMA_PEAK = 25.0; + int NMAX = 0; psTimerStart ("psphot"); @@ -20,7 +22,13 @@ detections = pmDetectionsAlloc(); pass = 1; + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); PS_ASSERT (status, NULL); + NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); PS_ASSERT (status, NULL); } else { pass = 2; + NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL); + NMAX = 0; // unlimited number of peaks in final pass: allow a limit (PEAKS_NMAX_2) ? } + + float threshold = PS_SQR(NSIGMA_PEAK); // move the old peaks array (if it exists) to oldPeaks @@ -34,5 +42,6 @@ // detect the peaks in the significance image - detections->peaks = psphotFindPeaks (significance, readout, recipe, pass, maskVal); + detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); // optionally merge peaks into footprints @@ -47,5 +56,2 @@ // if we use the footprints, the output peaks list contains both old and new peaks, // otherwise it only contains the new peaks. - -// if psf is defined, we should treat this as pass "2", not "1". re-define this boolean to -// be "have PSF" Index: /trunk/psphot/src/psphotFindFootprints.c =================================================================== --- /trunk/psphot/src/psphotFindFootprints.c (revision 17869) +++ /trunk/psphot/src/psphotFindFootprints.c (revision 17870) @@ -18,8 +18,8 @@ PS_ASSERT (status, false); - float effArea = psMetadataLookupF32 (&status, recipe, "EFFECTIVE_AREA"); - psAssert (status, "EFFECTIVE_AREA missing: call psphotFindPeaks first"); + // XXX do we need to use the same threshold here as for peaks? does it make sense for + // these to be different? - float threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea; + float threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT); int growRadius = 0; @@ -31,7 +31,4 @@ PS_ASSERT (status, false); - // XXX do we need to use the same threshold here as for peaks? does it make sense for - // these to be different? - // find the raw footprints & assign the peaks to those footprints psArray *footprints = pmFootprintsFind (significance, threshold, npixMin); @@ -39,6 +36,5 @@ // XXX handle the error conditions here - // footprints now owns the peaks; after culling (below), we will rebuild - // the peaks array + // footprints now owns the peaks; after culling (below), we will rebuild the peaks array psFree (detections->peaks); Index: /trunk/psphot/src/psphotFindPeaks.c =================================================================== --- /trunk/psphot/src/psphotFindPeaks.c (revision 17869) +++ /trunk/psphot/src/psphotFindPeaks.c (revision 17870) @@ -1,50 +1,12 @@ # include "psphotInternal.h" -// 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 +// Find peaks in the significance image above a threshold significance level. The significance +// image must be constructed to represent (S/N)^2. If nMax is non-zero, only return a maximum +// of nMax peaks +psArray *psphotFindPeaks (psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int nMax) { -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"); - - // 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"); - } else { - 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 - // gaussian SIGMA. a peak in the significance image has an effective S/N for faint sources - // of (S = Io*2pi*sigma_eff^2) / sqrt(Var), where sigma_eff^2 = sigma_obs^2 + sigma_sm^2. - // the smoothing of the varience map does not affect the faint-source S/N since the noise - // is constant under a faint source. - - // if sigma_sm = sigma_obs, then sigma_eff^2 = 2 sigma_sm^2. in this case, the threshold in - // 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: - float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH); - if (effArea < 1) { - effArea = 1; // never less than a pixel - } - if (pass == 1) { - effArea = 4.0*M_PI; // XXX better than nothing but still fairly poor - } - - 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 @@ -59,5 +21,5 @@ } - // correct the peak values to S/N = sqrt(value*effArea) + // correct the peak values to S/N = sqrt(significance) // get the peak flux from the unsmoothed image // the peak pixel coords are guaranteed to be on the image @@ -66,26 +28,22 @@ for (int i = 0; i < peaks->n; i++) { pmPeak *peak = peaks->data[i]; - peak->SN = sqrt(peak->value*effArea); + peak->SN = sqrt(peak->value); 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; - } + // limit the total number of returned peaks as specified + psArraySort (peaks, pmPeakSortBySN); + 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 - // XXX need a better check for this option char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE"); - if (status && (output != NULL) && (output[0])) { + if (output && strcasecmp (output, "NONE")) { pmPeaksWriteText (peaks, output); } @@ -95,3 +53,2 @@ } - Index: /trunk/psphot/src/psphotOutput.c =================================================================== --- /trunk/psphot/src/psphotOutput.c (revision 17869) +++ /trunk/psphot/src/psphotOutput.c (revision 17870) @@ -44,7 +44,7 @@ // optional dump of all rough source data char *output = psMetadataLookupStr (&status, recipe, "MOMENTS_OUTPUT_FILE"); - if (!status) return false; - if (output == NULL) return false; + if (!output) return false; if (output[0] == 0) return false; + if (!strcasecmp (output, "NONE")) return false; pmMomentsWriteText (sources, output); Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 17869) +++ /trunk/psphot/src/psphotReadout.c (revision 17870) @@ -117,4 +117,24 @@ psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); + if (0) { + FILE *out = fopen ("out.dat", "w"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + pmPeak *peak = source->peak; + if (!peak) continue; + pmModel *model = source->modelPSF; + if (!model) continue; + if (!model->params) continue; + + fprintf (out, "%d %f %f %f %f %f %f\n", + i, peak->xf, peak->yf, peak->flux, peak->SN, + model->params->data.F32[PM_PAR_I0], + model->dparams->data.F32[PM_PAR_I0]); + } + fclose (out); + } + // identify CRs and extended sources psphotSourceSize (readout, sources, recipe, 0); Index: /trunk/psphot/src/psphotSignificanceImage.c =================================================================== --- /trunk/psphot/src/psphotSignificanceImage.c (revision 17869) +++ /trunk/psphot/src/psphotSignificanceImage.c (revision 17870) @@ -1,6 +1,7 @@ # include "psphotInternal.h" -// In this function, we smooth the image, then search for the peaks -// if FWMH_X,Y have been recorded, use them; otherwise use PEAKS_SMOOTH_SIGMA +// In this function, we smooth the image and weight, then generate the significance image : +// (S/N)^2. If FWMH_X,Y have been recorded, use them, otherwise use PEAKS_SMOOTH_SIGMA for the +// smoothing kernel. psImage *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, const int pass, psMaskType maskVal) { @@ -34,5 +35,10 @@ psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("smooth")); - // smooth the weight, applying the mask as we go + // Smooth the weight, applying the mask as we go. The variance *should* be smoothed by the + // PSF^2, which does not have unity normalization (variance decreases as we smooth). + // Instead, we are smoothing with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity + // normalization. The resulting variance is a factor of 4*pi*SIGMA_SMTH^2 too large. We + // correct for this effect, and the effective area, in the calculation of the (S/N)^2 image + // below. psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32); psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH); @@ -51,5 +57,29 @@ } - // build the significance image on top of smooth_im + // We have an input image with PSF size sigma_r. We have smoothed it with a kernel of size + // sigma_s. The result is an image with PSF size sigma_o: sigma_o^2 = sigma_r^2 + + // sigma_s^2. Ideally, we are choosing sigma_s = sigma_r, in which case sigma_o^2 = 2 + // sigma_s^2. If we do not know the input image PSF size (initial detection stage), then + // we are assuming that sigma_r = sigma_s. + + // Build the significance image on top of smooth_im. We need to correct the ratio im/wt by + // two factors: 1) the error in the variance normalization above and 2) a factor to account + // for the relationship the peak value and the integrated flux, and the relationship + // between the per-pixel variance (var_i) and the total variance included in the flux + // measurement (effective area). These latter correction comes from: flux = Io * + // 2\pi\sigma_o^2 and total variance = var_i * 4\pi\sigma_o^2, thus (S/N)^2 = flux^2 / var + // = var_i \pi sigma_o^2 + + // thus: + // (S/N)^2 = (im^2 / wt) * (\pi \sigma_o^2 * 4 \sigma_s^2) + // (S/N)^2 = (im^2 / wt) * (\pi 2 \sigma_s^2 * 4 \sigma_s^2) + // (S/N)^2 = (im^2 / wt) * (\pi 8 \sigma_s^4) + + float factor = 8.0 * PS_SQR(M_PI) * pow(SIGMA_SMTH, 4.0); + // record the effective area and significance scaling factor + float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); + for (int j = 0; j < smooth_im->numRows; j++) { for (int i = 0; i < smooth_im->numCols; i++) { @@ -58,5 +88,5 @@ smooth_im->data.F32[j][i] = 0.0; } else { - smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i]; + smooth_im->data.F32[j][i] = factor * PS_SQR(value) / smooth_wt->data.F32[j][i]; } } @@ -75,2 +105,3 @@ return smooth_im; } +