Changeset 17870

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²

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotFindDetections.c (modified) (4 diffs)
• psphotFindFootprints.c (modified) (3 diffs)
• psphotFindPeaks.c (modified) (4 diffs)
• psphotOutput.c (modified) (1 diff)
• psphotReadout.c (modified) (1 diff)
• psphotSignificanceImage.c (modified) (5 diffs)

trunk/psphot/src/psphot.h

@@ -51 +51 @@
 // 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);

trunk/psphot/src/psphotFindDetections.c

@@ -6 +6 @@
     bool status;
     int pass;
+    float NSIGMA_PEAK = 25.0;
+    int NMAX = 0;
 
     psTimerStart ("psphot");
@@ -20 +22 @@
         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 +42 @@
 
     // 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 +56 @@
 // 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"

trunk/psphot/src/psphotFindFootprints.c

@@ -18 +18 @@
     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 +31 @@
     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 +36 @@
     // 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); 
 

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 @@
 
 }
-

trunk/psphot/src/psphotOutput.c

@@ -44 +44 @@
     // 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);

trunk/psphot/src/psphotReadout.c

@@ -117 +117 @@
     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);

trunk/psphot/src/psphotSignificanceImage.c

@@ -1 +1 @@
 # 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 +35 @@
     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 +57 @@
     }
 
-    // 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 +88 @@
                 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 +105 @@
     return smooth_im;
 }
+