Changeset 18632

Timestamp:

Jul 20, 2008, 1:06:21 PM (18 years ago)

Author:

eugene

Message:

mask the cosmic rays

File:

• branches/eam_branch_20080719/psphot/src/psphotSourceSize.c (modified) (2 diffs)

branches/eam_branch_20080719/psphot/src/psphotSourceSize.c

@@ -10 +10 @@
 // deviation from the psf model at the r = FWHM/2 position
 
-bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, long first) {
+bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) {
 
     bool status;
 
     psTimerStart ("psphot");
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    psMaskType maskVal = psMetadataLookupU8(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // bit to mask the cosmic-ray pixels
+    psMaskType crMask  = pmConfigMaskGet("CR", config); // Mask value for cosmic rays
 
     float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT");
@@ -129 +136 @@
         }
 
+        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
         if (source->crNsigma > CR_NSIGMA_LIMIT) {
-          source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+            source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+            pmPeak *peak = source->peak;
+            psAssert (peak, "NULL peak");
+
+            // replace the source flux
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+            source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
+
+            psImage *mask   = source->maskView;
+            psImage *pixels = source->pixels;
+            psImage *weight = source->weight;
+
+            # define SN_LIMIT 5.0
+
+            int xo = peak->x - pixels->col0;
+            int yo = peak->y - pixels->row0;
+
+            // mark the pixels in this row to the left, then the right
+            for (int ix = xo; ix >= 0; ix--) {
+                float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
+                if (SN > SN_LIMIT) {
+                    mask->data.U8[yo][ix] |= crMask;
+                }
+            }
+            for (int ix = xo + 1; ix < pixels->numCols; ix++) {
+                float SN = pixels->data.F32[yo][ix] / sqrt(weight->data.F32[yo][ix]);
+                if (SN > SN_LIMIT) {
+                    mask->data.U8[yo][ix] |= crMask;
+                }
+            }
+
+            // for each of the neighboring rows, mark the high pixels if they have a marked neighbor 
+            // first go up:
+            for (int iy = yo; iy >= 0; iy--) {
+                // mark the pixels in this row to the left, then the right
+                for (int ix = 0; ix < pixels->numCols; ix++) {
+                    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+                    if (SN < SN_LIMIT) continue;
+                    
+                    bool valid = false;
+                    valid |= (mask->data.U8[iy+1][ix] & crMask);
+                    valid |= (ix > 0) ? (mask->data.U8[iy+1][ix-1] & crMask) : 0;
+                    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy+1][ix+1] & crMask) : 0;
+
+                    if (!valid) continue;
+                    mask->data.U8[iy][ix] |= crMask;
+                }
+            }
+            // next go down:
+            for (int iy = yo+1; iy < pixels->numRows; iy++) {
+                // mark the pixels in this row to the left, then the right
+                for (int ix = 0; ix < pixels->numCols; ix++) {
+                    float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+                    if (SN < SN_LIMIT) continue;
+                    
+                    bool valid = false;
+                    valid |= (mask->data.U8[iy-1][ix] & crMask);
+                    valid |= (ix > 0) ? (mask->data.U8[iy-1][ix-1] & crMask) : 0;
+                    valid |= (ix <= mask->numCols) ? (mask->data.U8[iy-1][ix+1] & crMask) : 0;
+
+                    if (!valid) continue;
+                    mask->data.U8[iy][ix] |= crMask;
+                }
+            }
         }
     }