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

Timestamp:

Jun 19, 2007, 4:40:44 PM (19 years ago)

Author:

Paul Price

Message:

Extensive changes to APIs to allow use of a nominated value to mask
against (the maskVal). Previously, the mask values were either
hard-coded (e.g., PM_MASK_SAT) or taken as anything non-zero. The
code is tested with psModules (which has similar changes) and does not
crash, but neither is it successful in marking all bad pixels. For
this reason, I have left the "gutter" pixels (cell gaps) set to 0
instead of NAN in pmFPAMosaic.

File:

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

trunk/psphot/src/psphotFindPeaks.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-// In this function, we smooth the image, then search for the peaks 
+// In this function, we smooth the image, then search for the peaks
 psArray *psphotFindPeaks (pmReadout *readout, psMetadata *recipe,
-                          bool returnFootprints,
-                          const int pass) {
+                          bool returnFootprints,
+                          const int pass, psMaskType maskVal) {
 
     float SIGMA_SMTH, NSIGMA_SMTH, NSIGMA_PEAK;
@@ -13 +13 @@
 
     if (pass == 1) {
-        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
-        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        SIGMA_SMTH  = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_SIGMA");
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
     } else {
-        bool status_x, status_y;
-        float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
-        float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
-        if (!status_x | !status_y) {
-            psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
-            return false;
-        }
-        SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
-        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
+        bool status_x, status_y;
+        float FWHM_X = psMetadataLookupF32 (&status_x, recipe, "FWHM_X");
+        float FWHM_Y = psMetadataLookupF32 (&status_y, recipe, "FWHM_Y");
+        if (!status_x | !status_y) {
+            psError(PSPHOT_ERR_CONFIG, false, "FWHM_X or FWHM_Y not defined");
+            return false;
+        }
+        SIGMA_SMTH  = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrt(2.0*log(2.0)));
+        NSIGMA_SMTH = psMetadataLookupF32 (&status, recipe, "PEAKS_SMOOTH_NSIGMA");
     }
 
     // smooth the image, applying the mask as we go
     psImage *smooth_im = psImageCopy (NULL, readout->image, PS_TYPE_F32);
-    psImageSmoothMaskF32 (smooth_im, readout->mask, 0xff, SIGMA_SMTH, NSIGMA_SMTH);
+    psImageSmoothMaskF32 (smooth_im, readout->mask, maskVal, SIGMA_SMTH, NSIGMA_SMTH);
     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);
-    psImageSmoothMaskF32 (smooth_wt, readout->mask, 0xff, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
+    psImageSmoothMaskF32 (smooth_wt, readout->mask, maskVal, SIGMA_SMTH/sqrt(2), NSIGMA_SMTH);
     psLogMsg ("psphot", PS_LOG_MINUTIA, "smooth weight: %f sec\n", psTimerMark ("psphot"));
 
     psImage *mask = readout->mask;
 
-    // optionally save example images under trace 
+    // optionally save example images under trace
     if (psTraceGetLevel("psphot") > 5) {
-        char name[64];
-        sprintf (name, "imsmooth.v%d.fits", pass);
-        psphotSaveImage (NULL, smooth_im, name);
-        sprintf (name, "wtsmooth.v%d.fits", pass);
-        psphotSaveImage (NULL, smooth_wt, name);
+        char name[64];
+        sprintf (name, "imsmooth.v%d.fits", pass);
+        psphotSaveImage (NULL, smooth_im, name);
+        sprintf (name, "wtsmooth.v%d.fits", pass);
+        psphotSaveImage (NULL, smooth_wt, name);
     }
 
     // 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++) {
-            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;
-            } else {
-                smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
-            }
-        }
+        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.U8[j][i] & maskVal)) {
+                smooth_im->data.F32[j][i] = 0.0;
+            } else {
+                smooth_im->data.F32[j][i] = PS_SQR(value) / smooth_wt->data.F32[j][i];
+            }
+        }
     }
     psLogMsg ("psphot", PS_LOG_INFO, "built smoothed signficance image: %f sec\n", psTimerMark ("psphot"));
 
-    // optionally save example images under trace 
+    // optionally save example images under trace
     if (psTraceGetLevel("psphot") > 5) {
-        char name[64];
-        sprintf (name, "snsmooth.v%d.fits", pass);
-        psphotSaveImage (NULL, smooth_im, name);
+        char name[64];
+        sprintf (name, "snsmooth.v%d.fits", pass);
+        psphotSaveImage (NULL, smooth_im, name);
     }
 
@@ -73 +73 @@
     // signal/noise limit for the detected peaks
     if (pass == 1) {
-        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
+        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT");
     } else {
-        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
-    }   
-    
+        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2");
+    }
+
     // we need to define the threshold based on the value of NSIGMA_PEAK and the applied smoothing
     // gaussian SIGMA.  a peak in the significance image has an effective S/N for faint sources
@@ -87 +87 @@
     // 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: 
+    // thus, the threshold is:
     float effArea = 4.0*M_PI*PS_SQR(SIGMA_SMTH);
     if (effArea < 1) {
-        effArea = 1;                    // never less than a pixel
+        effArea = 1;                    // never less than a pixel
     }
     float threshold = PS_SQR(NSIGMA_PEAK) / effArea;
@@ -97 +97 @@
     psArray *peaks = pmFindImagePeaks (smooth_im, threshold);
     if (peaks == NULL) {
-        // XXX this may also be due to a programming or config error
-        // XXX do we need to set something in the readout->analysis to indicate that 
-        // we tried and failed to find peaks (something in the header data)
-        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
-        return false;
+        // XXX this may also be due to a programming or config error
+        // XXX do we need to set something in the readout->analysis to indicate that
+        // we tried and failed to find peaks (something in the header data)
+        psError(PSPHOT_ERR_DATA, false, "no peaks found in this image");
+        return false;
     }
-    
+
     // correct the peak values to S/N = sqrt(value*effArea)
     // get the peak flux from the unsmoothed image
@@ -110 +110 @@
     int col0 = readout->image->col0;
     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-row0][peak->x-col0];
+        pmPeak *peak = peaks->data[i];
+        peak->SN = sqrt(peak->value*effArea);
+        peak->flux = readout->image->data.F32[peak->y-row0][peak->x-col0];
     }
 
     // limit the total number of returned peak as specified
     if (pass == 1) {
-        psArraySort (peaks, pmPeakSortBySN);
-        int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX");
-        if (NMAX && (peaks->n > NMAX)) {
-            psArray *tmpPeaks = psArrayAllocEmpty (NMAX);
-            for (int i = 0; i < NMAX; i++) {
-                psArrayAdd (tmpPeaks, 100, peaks->data[i]);
-            }
-            psFree (peaks);
-            peaks = tmpPeaks;
-        }
-    }    
+        psArraySort (peaks, pmPeakSortBySN);
+        int NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX");
+        if (NMAX && (peaks->n > NMAX)) {
+            psArray *tmpPeaks = psArrayAllocEmpty (NMAX);
+            for (int i = 0; i < NMAX; i++) {
+                psArrayAdd (tmpPeaks, 100, peaks->data[i]);
+            }
+            psFree (peaks);
+            peaks = tmpPeaks;
+        }
+    }
 
     // optional dump of all peak data
     char *output = psMetadataLookupStr (&status, recipe, "PEAKS_OUTPUT_FILE");
     if (status && (output != NULL) && (output[0])) {
-        pmPeaksWriteText (peaks, output);
+        pmPeaksWriteText (peaks, output);
     }
     psLogMsg ("psphot", PS_LOG_INFO, "%ld peaks: %f sec\n", peaks->n, psTimerMark ("psphot"));
@@ -142 +142 @@
     // want to do that
     //
-    if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
-        int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
-        if (!status) {
-            npixMin = 1;
-        }
-        float FOOTPRINT_NSIGMA_LIMIT =
-            psMetadataLookupS32(&status, recipe,
-                                (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
-        if (!status) {
-            FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
-        }
-        threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;
+    if (returnFootprints) {     // We want an array of pmFootprint, not pmPeak
+        int npixMin = psMetadataLookupS32(&status, recipe, "FOOTPRINT_NPIXMIN");
+        if (!status) {
+            npixMin = 1;
+        }
+        float FOOTPRINT_NSIGMA_LIMIT =
+            psMetadataLookupS32(&status, recipe,
+                                (pass == 1) ? "FOOTPRINT_NSIGMA_LIMIT" : "FOOTPRINT_NSIGMA_LIMIT_2");
+        if (!status) {
+            FOOTPRINT_NSIGMA_LIMIT = NSIGMA_PEAK;
+        }
+        threshold = PS_SQR(FOOTPRINT_NSIGMA_LIMIT)/effArea;
 
-        psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
-        pmPeaksAssignToFootprints(footprints, peaks);
+        psArray *footprints = pmFindFootprints(smooth_im, threshold, npixMin);
+        pmPeaksAssignToFootprints(footprints, peaks);
 
-        psFree(peaks);
-        peaks = footprints;             // well, you know what I mean
+        psFree(peaks);
+        peaks = footprints;             // well, you know what I mean
     }
 
@@ -174 +174 @@
  */
 psErrorCode
-psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
-                const psImage *weight,  // corresponding variance image
-                const psMetadata *recipe,
-                psArray *footprints) {  // array of pmFootprints
+psphotCullPeaks(const psImage *image,   // the image wherein lives the footprint
+                const psImage *weight,  // corresponding variance image
+                const psMetadata *recipe,
+                psArray *footprints) {  // array of pmFootprints
     bool status = false;
     float nsigma_delta = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_DELTA");
     if (!status) {
-        nsigma_delta = 0; // min. 
+        nsigma_delta = 0; // min.
     }
     float nsigma_min = psMetadataLookupF32(&status, recipe, "FOOTPRINT_CULL_NSIGMA_MIN");
     if (!status) {
-        nsigma_min = 0;
+        nsigma_min = 0;
     }
     const float skyStdev = psMetadataLookupF32(NULL, recipe, "SKY_STDEV");
 
     return pmFootprintArrayCullPeaks(image, weight, footprints,
-                                     nsigma_delta, nsigma_min*skyStdev);
+                                     nsigma_delta, nsigma_min*skyStdev);
 }