Index: /trunk/psphot/src/psphotFindPeaks.c
===================================================================
--- /trunk/psphot/src/psphotFindPeaks.c	(revision 11168)
+++ /trunk/psphot/src/psphotFindPeaks.c	(revision 11169)
@@ -1,29 +1,23 @@
 # include "psphot.h"
 
-// 2006.02.02 : no leaks
+// In this function, we smooth the image, then search for the peaks 
 psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe) {
 
     bool  status = false;
-    float NSIGMA;
-    float SIGMA;
-    float threshold;
-    float value;
 
     // smooth the image and weight map
     psTimerStart ("psphot");
 
-    SIGMA  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
-    NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+    float SIGMA  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
+    float NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
 
-    // XXX we are somewhat doing the wrong thing here with the mask.
-    // it would be better to smooth the image using the mask information
+    // smooth the image, applying the mask as we go
     psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
-    psImageSmooth (smooth_im, SIGMA, NSIGMA);
+    psImageSmoothMaskF32 (smooth_im, readout->mask, 0xff, SIGMA, NSIGMA);
     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth image: %f sec\n", psTimerMark ("psphot"));
 
+    // smooth the weight, applying the mask as we go
     psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
-    psImageSmooth (smooth_wt, SIGMA/sqrt(2), NSIGMA);
-    smooth_wt = (psImage*)psBinaryOp(smooth_wt, smooth_wt, "*",
-				     psScalarAlloc(1/(4*M_PI*PS_SQR(SIGMA)), PS_TYPE_F32));
+    psImageSmoothMaskF32 (smooth_wt, readout->mask, 0xff, SIGMA/sqrt(2), NSIGMA);
     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot"));
 
@@ -36,7 +30,8 @@
     }
 
+    // build the significance image on top of smooth_im
     for (int j = 0; j < smooth_im->numRows; j++) {
 	for (int i = 0; i < smooth_im->numCols; i++) {
-	    value = smooth_im->data.F32[j][i];
+	    float value = smooth_im->data.F32[j][i];
 	    if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || mask->data.U8[j][i]) {
 		smooth_im->data.F32[j][i] = 0.0;
@@ -46,16 +41,43 @@
 	}
     }
-    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed image: %f sec\n", psTimerMark ("psphot"));
+    psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot"));
+
+    // optionally save example images under trace 
+    if (psTraceGetLevel("psphot") > 5) {
+	psphotSaveImage (NULL, smooth_im, "snsmooth.fits");
+    }
 
     psTimerStart ("psphot");
     // set peak threshold
+
+    // signal/noise limit for the detected peaks
     NSIGMA = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
-    // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
-    threshold = PS_SQR(NSIGMA);
-    psLogMsg ("psphot", PS_LOG_DETAIL, "threshold: %f DN\n", threshold);
 
-    // find the peaks in the smoothed image 
+    // we need to define the threshold based on the value of NSIGMA 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 smooth_im 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);
+    float threshold = PS_SQR(NSIGMA) / effArea;
+
+    // find the peaks in the smoothed image
     psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
     if (peaks == NULL) psAbort ("find peaks", "no peaks found");
+    
+    // correct the peak values to S/N = sqrt(value*effArea)
+    // get the peak flux from the unsmoothed image
+    // the peak pixel coords are guaranteed to be on the image
+    for (int i = 0; i < peaks->n; i++) {
+	pmPeak *peak = peaks->data[i];
+	peak->SN = sqrt(peak->value*effArea);
+	peak->flux = readout->image->data.F32[peak->y][peak->x];
+    }
+
     psFree (smooth_im);
     psFree (smooth_wt);
@@ -71,5 +93,2 @@
 }
 
-// In this function, we smooth the image, then search for the peaks 
-// Should we also subtract a super-binned image? (as an option?)
-// XXX : We need to gracefully handle no source detections