Changeset 30903

Timestamp:

Mar 14, 2011, 3:07:22 PM (15 years ago)

Author:

eugene

Message:

psphotSignificanceImage now returns a pmReadout so we can keep the smoothed signal image as well as the S/N image; psphotCullPeaks (and pmFootprintCullPeaks) now uses the smooted image to cull the peaks -- this deals better with low-significance detections are large, low-surface brightness areas

Location:

branches/eam_branches/ipp-20110213/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotCullPeaks.c (modified) (4 diffs)
• psphotFindDetections.c (modified) (1 diff)
• psphotFindFootprints.c (modified) (3 diffs)
• psphotSignificanceImage.c (modified) (3 diffs)

branches/eam_branches/ipp-20110213/psphot/src/psphot.h

@@ -173 +173 @@
 
 // used by psphotFindDetections
-psImage        *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal);
+pmReadout      *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType 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 float threshold, const int pass, psImageMaskType maskVal);
-psErrorCode     psphotCullPeaks(const pmReadout *readout, const psMetadata *recipe, psArray *footprints);
+bool            psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal);
+psErrorCode     psphotCullPeaks(const pmReadout *readout, const pmReadout *signifRO, const psMetadata *recipe, psArray *footprints);
 
 // in psphotApResid.c:

branches/eam_branches/ipp-20110213/psphot/src/psphotCullPeaks.c

@@ -6 +6 @@
 psErrorCode
 psphotCullPeaks(const pmReadout *readout,   // the image wherein lives the footprint
+                const pmReadout *signifR, // smoothed image and significance
                 const psMetadata *recipe,
                 psArray *footprints) {  // array of pmFootprints
@@ -26 +27 @@
     psAssert (status, "cannot find SKY_STDEV in readout->analysis");
 
-    const float MIN_THRESHOLD = nsigma_min*skyStdev;
+    // the sky has been smoothed, so we need to scale down the raw sky stdev by this amount:
+    float scaleFactor = psMetadataLookupF32(&status, readout->analysis, "SIGNIFICANCE_SCALE_FACTOR");
+    if (!status) scaleFactor = 1.0;
+
+    const float MIN_THRESHOLD = nsigma_min*skyStdev / sqrt(scaleFactor);
     
     // for saturated stars, we should be somewhat more agressive about culling: instead of
@@ -36 +41 @@
     // fStdev = 0.005, stdev_pad = 0.011*40k = 400
     // threshold = 40k - 4*400 = 38400 
-
-    // in practice, we should make the threshold in such a case more like 0.5 * saturation...
+    // this gives too tight of a tolerance on the bright stars
+    // in practice, we should make the threshold much lower.  
+    // below I am using 0.05 * saturation (eg, 2000 DN above sky in a GPC1 image)
 
     float SATURATION = NAN;
@@ -72 +78 @@
         }
 
-        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold) != PS_ERR_NONE) {
+        // if we cull using the significance image, then the thresholds are different
+# if (0)
+        if (pmFootprintCullPeaks(signifR->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) {
             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
         }
+# endif
+# if (1)
+        if (pmFootprintCullPeaks(signifR->image, signifR->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, false) != PS_ERR_NONE) {
+            return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
+        }
+# endif
+# if (0)
+        if (pmFootprintCullPeaks(readout->image, readout->variance, fp, nsigma_delta, fPadding, MIN_THRESHOLD, max_threshold, true) != PS_ERR_NONE) {
+             return psError(PS_ERR_UNKNOWN, false, "Culling pmFootprint %d", fp->id);
+        }
+# endif
+
         float dtime = psTimerMark ("psphot.cull.footprints");
         if (dtime > 1.0) {

branches/eam_branches/ipp-20110213/psphot/src/psphotFindDetections.c

@@ -84 +84 @@
 
     // generate the smoothed significance image
-    psImage *significance = psphotSignificanceImage (readout, recipe, maskVal);
+    pmReadout *significance = psphotSignificanceImage (readout, recipe, maskVal);
 
     // display the log significance image
-    psphotVisualShowLogSignificance (significance, 0.0, 4.5);
+    psphotVisualShowLogSignificance (significance->variance, 0.0, 4.5);
 
     // display the significance image
-    psphotVisualShowSignificance (significance, 0.98*threshold, 1.02*threshold);
+    psphotVisualShowSignificance (significance->variance, 0.98*threshold, 1.02*threshold);
 
     // detect the peaks in the significance image
-    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX);
+    detections->peaks = psphotFindPeaks (significance->variance, readout, recipe, threshold, NMAX);
     psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
     if (!detections->peaks) {

branches/eam_branches/ipp-20110213/psphot/src/psphotFindFootprints.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-bool psphotFindFootprints (pmDetections *detections, psImage *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) {
+bool psphotFindFootprints (pmDetections *detections, pmReadout *significance, pmReadout *readout, psMetadata *recipe, const float threshold, const int pass, psImageMaskType maskVal) {
 
     bool status;
@@ -18 +18 @@
     PS_ASSERT (status, false);
 
-    // find the raw footprints & assign the peaks to those footprints
-    psArray *footprints = pmFootprintsFind (significance, threshold, npixMin);
+    // find the raw footprints in the smoothed significance image & assign the peaks to those footprints
+    psArray *footprints = pmFootprintsFind (significance->variance, threshold, npixMin);
 
     if (pmFootprintsAssignPeaks(footprints, detections->peaks) != PS_ERR_NONE) {
@@ -56 +56 @@
     }
 
-    psphotCullPeaks(readout, recipe, detections->footprints);
+    psphotCullPeaks(readout, significance, recipe, detections->footprints);
     detections->peaks = pmFootprintArrayToPeaks(detections->footprints);
     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks, %ld total footprints: %f sec\n", detections->peaks->n, detections->footprints->n, psTimerMark ("psphot.footprints"));

branches/eam_branches/ipp-20110213/psphot/src/psphotSignificanceImage.c

@@ -4 +4 @@
 // (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, psImageMaskType maskVal) {
+pmReadout *psphotSignificanceImage (pmReadout *readout, psMetadata *recipe, psImageMaskType maskVal) {
 
     float SIGMA_SMTH, NSIGMA_SMTH;
@@ -94 +94 @@
     psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
 
-    // XXX multithread this if needed
+    // we are going to return both the image and the weight here: the image contains the signal
+    // while the 'weight' will contain the significance (NOTE the deviation from the usual
+    // definition)
+
+    // save the smoothed significance image in the weight array
     for (int j = 0; j < smooth_im->numRows; j++) {
         for (int i = 0; i < smooth_im->numCols; i++) {
             float value = smooth_im->data.F32[j][i];
             if (value < 0 || smooth_wt->data.F32[j][i] <= 0 || (mask->data.PS_TYPE_IMAGE_MASK_DATA[j][i] & maskVal)) {
-                smooth_im->data.F32[j][i] = 0.0;
+                smooth_wt->data.F32[j][i] = 0.0;
             } else {
                 // XXX the value of 100 here (or 1000 before) must depend on the FWHM of the smoothing kernel, right??
                 float v2 = value + PS_SQR(value/100.0);
-                smooth_im->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
+                smooth_wt->data.F32[j][i] = factor * PS_SQR(v2) / smooth_wt->data.F32[j][i];
             }
         }
@@ -115 +119 @@
         static int pass = 0;
         sprintf (name, "snsmooth.v%d.fits", pass);
-        psphotSaveImage (NULL, smooth_im, name);
+        psphotSaveImage (NULL, smooth_wt, name);
         pass ++;
     }
-
-    psFree(smooth_wt);
-
     psImageConvolveSetThreads(oldThreads);
 
-    return smooth_im;
+    pmReadout *significanceRO = pmReadoutAlloc(NULL);
+    significanceRO->variance = smooth_wt;    
+    significanceRO->image = smooth_im;    
+
+    return significanceRO;
 }