Changeset 17870 for trunk/psphot/src/psphotFindPeaks.c

Timestamp:

May 30, 2008, 4:09:46 PM (18 years ago)

Author:

eugene

Message:

fix errors in calculation of sig image; API change: sig image is now actual S/N²

File:

• trunk/psphot/src/psphotFindPeaks.c (modified) (4 diffs)

trunk/psphot/src/psphotFindPeaks.c

@@ -1 +1 @@
 # 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 +21 @@
     }
 
-    // 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 +28 @@
     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 +53 @@
 
 }
-