Changeset 20938

Timestamp:

Dec 7, 2008, 4:52:04 PM (17 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

Location:

trunk/psphot/src

Files:

• psphot.h (modified) (1 diff)
• psphotFitSourcesLinear.c (modified) (1 diff)
• psphotSourceSize.c (modified) (4 diffs)
• psphotVisual.c (modified) (3 diffs)

trunk/psphot/src/psphot.h

@@ -165 +165 @@
 bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
 bool psphotVisualShowFlags (psArray *sources);
-bool psphotVisualShowSourceSize (psArray *sources);
+bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
 bool psphotVisualShowResidualImage (pmReadout *readout);
 bool psphotVisualPlotApResid (psArray *sources);

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -68 +68 @@
 
         // if (source->type == PM_SOURCE_TYPE_STAR &&
-        // source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
 
         if (final) {

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;
+}

trunk/psphot/src/psphotVisual.c

@@ -26 +26 @@
 
     isVisual = mode;
+    return true;
+}
+
+bool psphotVisualShowMask (int kapaFD, psImage *inImage, const char *name, int channel) {
+
+    KiiImage image;
+    KapaImageData data;
+    Coords coords;
+
+    strcpy (coords.ctype, "RA---TAN");
+
+    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0);
+    if (!psImageBackground(stats, NULL, inImage, NULL, 0, rng)) {
+        fprintf (stderr, "failed to get background values\n");
+        return false;
+    }
+
+    image.Nx = inImage->numCols;
+    image.Ny = inImage->numRows;
+
+    ALLOCATE (image.data2d, float *, image.Ny);
+    for (int iy = 0; iy < image.Ny; iy++) {
+        ALLOCATE (image.data2d[iy], float, image.Nx);
+        for (int ix = 0; ix < image.Nx; ix++) {
+            image.data2d[iy][ix] = inImage->data.U8[iy][ix];
+        }
+    }
+
+    strcpy (data.name, name);
+    strcpy (data.file, name);
+    data.zero = -1;
+    data.range = 32;
+    data.logflux = 0;
+  
+    KiiSetChannel (kapaFD, channel);
+    KiiNewPicture2D (kapaFD, &image, &data, &coords);
+
+    for (int iy = 0; iy < image.Ny; iy++) {
+        free (image.data2d[iy]);
+    }
+    free (image.data2d);
+
+    psFree (stats);
+    psFree (rng);
+
     return true;
 }
@@ -1336 +1382 @@
 }
 
-bool psphotVisualShowSourceSize (psArray *sources) {
+bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources) {
 
     int Noverlay, NOVERLAY;
@@ -1398 +1444 @@
     KiiLoadOverlay (kapa, overlay, Noverlay, "blue");
     FREE (overlay);
+
+    psphotVisualShowMask (kapa, readout->mask, "mask", 2);
 
     // pause and wait for user input: