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

Timestamp:

Feb 14, 2006, 6:33:17 AM (20 years ago)

Author:

eugene

Message:

mods to change from threshold above sky to S/N thresholds

File:

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

trunk/psphot/src/psphotFindPeaks.c

@@ -2 +2 @@
 
 // 2006.02.02 : no leaks
-psArray *psphotFindPeaks (pmReadout *readout, psMetadata *config, psStats *sky) {
+psArray *psphotFindPeaks (pmReadout *readout, psMetadata *config) {
 
     bool  status = false;
@@ -9 +9 @@
     float threshold;
 
-    // smooth the image 
+    // smooth the image and weight map
+
+    psphotSaveImage (NULL, readout->image, "image.fits");
 
     psTimerStart ("psphot");
@@ -16 +18 @@
     NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_SMOOTH_NSIGMA");
 
-    psImage *smooth = psImageCopy (NULL, readout->image, PS_TYPE_F32);
-    psphotSaveImage (NULL, smooth, "presmooth.fits");
+    psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
+    psImageSmooth (smooth_im, SIGMA, NSIGMA);
+    psLogMsg ("psphot", 4, "smooth image: %f sec\n", psTimerMark ("psphot"));
 
-    psImageSmooth (smooth, SIGMA, NSIGMA);
-    psLogMsg ("psphot", 4, "smooth: %f sec\n", psTimerMark ("psphot"));
+    psphotSaveImage (NULL, smooth_im, "smooth_im.fits");
 
-    psphotSaveImage (NULL, smooth, "smooth.fits");
+    psImage *smooth_wt = psImageCopy (NULL, readout->weight, PS_TYPE_F32);
+    psImageSmooth (smooth_wt, SIGMA, NSIGMA);
+    psLogMsg ("psphot", 4, "smooth weight: %f sec\n", psTimerMark ("psphot"));
+
+    psphotSaveImage (NULL, smooth_wt, "smooth_wt.fits");
+
+    psTimerStart ("psphot");
+    for (int j = 0; j < smooth_im->numRows; j++) {
+        for (int i = 0; i < smooth_im->numCols; i++) {
+            if (smooth_wt->data.F32[j][i] <= 0) continue;
+            smooth_im->data.F32[j][i] = PS_SQR(smooth_im->data.F32[j][i]) / smooth_wt->data.F32[j][i];
+        }
+    }
+    psLogMsg ("psphot", 4, "built S/N map: %f sec\n", psTimerMark ("psphot"));
 
     psTimerStart ("psphot");
 
     // set peak threshold
-    NSIGMA    = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
-    
-    threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
-    // threshold = NSIGMA*sky->sampleStdev;
+    NSIGMA = psMetadataLookupF32 (&status, config, "PEAKS_NSIGMA_LIMIT");
+    // threshold = NSIGMA*sky->sampleStdev + sky->sampleMean; 
+    threshold = PS_SQR(NSIGMA);
     psLogMsg ("psphot", 4, "threshold: %f DN\n", threshold);
 
     // find the peaks in the smoothed image 
-    psArray *peaks = pmFindImagePeaks (smooth, threshold);
+    psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
     if (peaks == NULL) psAbort ("find peaks", "no peaks found");
-    psFree (smooth);
+    psFree (smooth_im);
+    psFree (smooth_wt);
 
     // optional dump of all peak data