Changeset 17870 for trunk/psphot/src/psphotSignificanceImage.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/psphotSignificanceImage.c (modified) (5 diffs)

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;
 }
+