Changeset 20978

Timestamp:

Dec 14, 2008, 9:18:03 AM (17 years ago)

Author:

eugene

Message:

working on the CR masking : trying both footprints and re-find the outline by hand

File:

• branches/eam_branch_20081214/psphot/src/psphotSourceSize.c (modified) (7 diffs)

branches/eam_branch_20081214/psphot/src/psphotSourceSize.c

@@ -2 +2 @@
 # include <gsl/gsl_sf_gamma.h>
 
+// XXX This should be a recipe variable
+# define SN_LIMIT 5.0
+
 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_Old (pmReadout *readout, pmSource *source, psMaskType maskVal, psMaskType crMask);
 bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psMaskType maskVal, psMaskType crMask);
 
@@ -44 +47 @@
         soft = 0.0;
     }
+
+    psphotVisualShowSNLimit (readout, SN_LIMIT);
 
     // loop over all source
@@ -161 +166 @@
         if (source->crNsigma > CR_NSIGMA_LIMIT) {
             // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
-            psphotMaskCosmicRay_Old (source, maskVal, crMask);
+            // psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask);
+            psphotMaskCosmicRay_Old (readout, source, maskVal, crMask);
         }
     }
@@ -275 +281 @@
     if (!footprint) {
         // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        // XXX this is broken, but we probably will not use it
+        psphotMaskCosmicRay_Old (NULL, source, maskVal, crMask);
         return true;
     }
@@ -281 +288 @@
     if (!footprint->spans) {
         // if we have not footprint, use the old code to mask by isophot
-        psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        // XXX this is broken, but we probably will not use it
+        psphotMaskCosmicRay_Old (NULL, source, maskVal, crMask);
         return true;
     }
@@ -300 +308 @@
 }
 
-bool psphotMaskCosmicRay_Old (pmSource *source, psMaskType maskVal, psMaskType crMask) {
+bool psphotMaskCosmicRay_Old (pmReadout *readout, pmSource *source, psMaskType maskVal, psMaskType crMask) {
 
     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
@@ -306 +314 @@
     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
+    // XXX should we be doing this on the smoothed image??
+    psImage *pixels = readout->image;
+    psImage *mask   = readout->mask;
+    psImage *weight = readout->weight;
 
     int xo = peak->x - pixels->col0;
     int yo = peak->y - pixels->row0;
-
-    // mark the pixels in this row to the left, then the right
+    bool rowDone;
+
+    psVector *xf = psVectorAllocEmpty (64, PS_TYPE_F32);
+    psVector *yf = psVectorAllocEmpty (64, PS_TYPE_F32);
+
+    // mark the pixels in this row to the left, then the right.  record the bounds of the valid pixel (SN >= SN_LIMIT)
+    int xs = 0;
+    int xe = pixels->numCols;
     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;
-        }
+        if (SN < SN_LIMIT) {
+            xs = ix + 1;
+            break;
+        }
+        mask->data.U8[yo][ix] |= crMask;
+        psVectorAppend (xf, (float) ix);
+        psVectorAppend (yf, (float) yo);
     }
     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;
-        }
-    }
+        if (SN < SN_LIMIT) {
+            xe = ix;
+            break;
+        } 
+        mask->data.U8[yo][ix] |= crMask;
+        psVectorAppend (xf, (float) ix);
+        psVectorAppend (yf, (float) yo);
+    }
+    psphotVisualShowCRMask (xf, yf);
+
+    int startXs = xs;
+    int startXe = xe;
+    rowDone = false;
 
     // 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;
-        }
-    }
+    // first go down:
+    for (int iy = yo - 1; !rowDone && (iy >= 0); iy--) {
+        rowDone = true; // we are not done if any valid pixels are found a row
+
+        // given the range xs to xe (valid pixel range from previous row), find the valid
+        // pixels within this range, and extend if the ends are valid
+
+        // newXs,newXe will define the valid pixel range for this row
+        int newXs = xe;
+        int newXe = xs;
+
+        // first, try all of the pixels in the valid range (xs < x < xe guaranteed to be real pixels)
+        for (int ix = xs; ix < xe; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                continue;
+            }
+            newXs = PS_MIN (newXs, ix);
+            newXe = PS_MAX (newXe, ix + 1);
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        // next, try the pixels to the left of the valid range for the previous row
+        for (int ix = xs - 1; ix >= 0; ix--) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                newXs = ix + 1;
+                break;
+            }
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        // finally, try the pixels to the right of the valid range for the previous row
+        for (int ix = xe; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                newXe = ix;
+                break;
+            }
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        xs = newXs;
+        xe = newXe;
+    }
+    psphotVisualShowCRMask (xf, yf);
+
+    xs = startXs;
+    xe = startXe;
+    rowDone = false;
+
+    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
+    // first go down:
+    for (int iy = yo + 1; !rowDone && (iy < pixels->numRows); iy++) {
+        rowDone = true; // we are not done if any valid pixels are found a row
+
+        // given the range xs to xe (valid pixel range from previous row), find the valid
+        // pixels within this range, and extend if the ends are valid
+
+        // newXs,newXe will define the valid pixel range for this row
+        int newXs = xe;
+        int newXe = xs;
+
+        // first, try all of the pixels in the valid range (xs < x < xe guaranteed to be real pixels)
+        for (int ix = xs; ix < xe; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                continue;
+            }
+            newXs = PS_MIN (newXs, ix);
+            newXe = PS_MAX (newXe, ix + 1);
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        // next, try the pixels to the left of the valid range for the previous row
+        for (int ix = xs - 1; ix >= 0; ix--) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                newXs = ix + 1;
+                break;
+            }
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        // finally, try the pixels to the right of the valid range for the previous row
+        for (int ix = xe; ix < pixels->numCols; ix++) {
+            float SN = pixels->data.F32[iy][ix] / sqrt(weight->data.F32[iy][ix]);
+            if (SN < SN_LIMIT) {
+                newXe = ix;
+                break;
+            }
+            rowDone = false;
+            mask->data.U8[iy][ix] |= crMask;
+            psVectorAppend (xf, (float) ix);
+            psVectorAppend (yf, (float) iy);
+        }
+
+        xs = newXs;
+        xe = newXe;
+    }
+
+    psphotVisualShowCRMask (xf, yf);
+    psFree (xf);
+    psFree (yf);
+
     return true;
 }