Changeset 18005

Timestamp:

Jun 8, 2008, 2:47:42 PM (18 years ago)

Author:

eugene

Message:

another attempt to get the significance image right

File:

• trunk/psphot/src/psphotSignificanceImage.c (modified) (5 diffs)

trunk/psphot/src/psphotSignificanceImage.c

@@ -8 +8 @@
     float SIGMA_SMTH, NSIGMA_SMTH;
     bool status = false;
+    bool guess = false;
 
     // smooth the image and weight map
@@ -19 +20 @@
       SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
       NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+      guess = false;
     } else {
       // if we do not know the FWHM, use the guess smoothing kernel supplied.
@@ -26 +28 @@
       NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA"); 
       PS_ASSERT (status, NULL);
+      guess = true;
     }
     // record the actual smoothing sigma
@@ -35 +38 @@
     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("smooth"));
 
-    // 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.
+    // Smooth the weight, applying the mask as we go.  The variance is smoothed by the PSF^2,
+    // renomalized to maintain the input level of the variance.  We achieve this by smoothing
+    // with a Gaussian with sigma = SIGMA_SMTH/sqrt(2) with unity normalization.  Note that
+    // this process yields a smoothed image with correlated errors.  The pixel-to-pixel
+    // variations in smooth_im will be decreased by a factor of 4*pi*SIGMA_SMTH^2, but for
+    // measurements based on apertures comparable to or larger than the smoothing kernel, the
+    // effective per-pixel variance is maintained.
     psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
     psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
@@ -64 +68 @@
 
     // 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 factors: 1) the relationship the peak value and the integrated flux, and 2) 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
+    // = (Io/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)
+    // in fact, because of the smoothing and the resulting correlated errors, the per-pixel
+    // varience is a bit smaller than the value predicted here.  The total variance should be
+    // \alpha 4\pi(\sigma_r_^2 + \sigma_s^2) where \alpha approaches 1 for \sigma_s / \sigma_r
+    // << 1, and \alpha approaches 0.5 for \sigma_s / \sigma_r >> 1, with a value of 0.7
+    // (actually 1/1.5) for \sigma_s = \sigma_r
 
-    float factor = 8.0 * PS_SQR(M_PI) * pow(SIGMA_SMTH, 4.0);
+    // if the smoothing kernel size is known to be the psf size (\sigma_s = \sigma_r), then we
+    // can use the scaling factor of 1/1.5; otherwise, we want to be conservative and use \alpha
+    // = 0.5.  Thus, for the guess smoothing kernel, (S/N)^2 = 3.0 (Io/var_i) \pi sigma_s^2 for
+    // the guess case and (S/N)^2 = 2.0 (Io/var_i) \pi sigma_s^2 for the known case
+
+    // XXX hmm, somehow the simulations are coming up with 4pi, not 3pi (as if \alpha is 1/2.0,
+    // not 1/1.5 for sigma_s = sigma_r...
+
+    float factor = guess ? 4.0 * M_PI * PS_SQR(SIGMA_SMTH) : 4.0 * M_PI * PS_SQR(SIGMA_SMTH);
     // record the effective area and significance scaling factor
     float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);