Changeset 20938 for trunk/psphot/src/psphotSourceSize.c

Timestamp:

Dec 7, 2008, 4:52:04 PM (18 years ago)

Author:

eugene

Message:

do not fit known CRs in FitSourcesLinear; defined a new CR masking function based on footprints, but it is not yet ready

File:

• trunk/psphot/src/psphotSourceSize.c (modified) (4 diffs)

trunk/psphot/src/psphotSourceSize.c

@@ -4 +4 @@
 static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask,
                                 psMaskType maskVal, const pmModel *model, float Ro);
+
+bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask);
+bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -157 +160 @@
         // 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;
-            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;
-
-            // XXX This should be a recipe variable
-            # 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 = PS_MIN(yo, mask->numRows-2); 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 = PS_MIN(yo+1, mask->numRows-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;
-                }
-            }
+            // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
+            psphotMaskCosmicRay_Old (source, maskVal, crMask);
         }
     }
@@ -242 +182 @@
 
     psphotVisualPlotSourceSize (sources);
-    psphotVisualShowSourceSize (sources);
+    psphotVisualShowSourceSize (readout, sources);
 
     return true;
@@ -319 +259 @@
     return nSigma;
 }
+
+bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask) {
+
+    // replace the source flux
+    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+    source->mode &= ~PM_SOURCE_MODE_SUBTRACTED;
+
+    // flag this as a CR
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    // grab the matching footprint
+    pmFootprint *footprint = peak->footprint;
+    if (!footprint) {
+        // if we have not footprint, use the old code to mask by isophot
+        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        return true;
+    }
+
+    if (!footprint->spans) {
+        // if we have not footprint, use the old code to mask by isophot
+        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        return true;
+    }
+
+    // mask all of the pixels covered by the spans of the footprint
+    for (int j = 1; j < footprint->spans->n; j++) {
+        pmSpan *span1 = footprint->spans->data[j];
+
+        int iy = span1->y;
+        int xs = span1->x0;
+        int xe = span1->x1;
+        
+        for (int ix = xs; ix < xe; ix++) {
+            mask->data.U8[iy][ix] |= crMask;
+        }
+    }
+    return true;
+}
+
+bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask) {
+
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+    pmPeak *peak = source->peak;
+    psAssert (peak, "NULL peak");
+
+    psImage *mask   = source->maskView;
+    psImage *pixels = source->pixels;
+    psImage *weight = source->weight;
+
+    // XXX This should be a recipe variable
+# 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 = PS_MIN(yo, mask->numRows-2); 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 = PS_MIN(yo+1, mask->numRows-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;
+        }
+    }
+    return true;
+}