Changeset 27440 for branches/eam_branches/20100225/psphot/src/psphotSourceSize.c

Timestamp:

Mar 24, 2010, 4:56:11 PM (16 years ago)

Author:

eugene

Message:

updates for source sizes and CR masking

File:

• branches/eam_branches/20100225/psphot/src/psphotSourceSize.c (modified) (11 diffs)

branches/eam_branches/20100225/psphot/src/psphotSourceSize.c

@@ -23 +23 @@
 bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
 bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
+int  psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh);
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -252 +253 @@
 
 # define SIZE_SN_LIM 10
-
 bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) {
 
@@ -484 +484 @@
 
     // XXX test : save the mask image
-    if (1) {
+    if (0) {
         psphotSaveImage (NULL, readout->mask,   "mask.fits");
     }
@@ -491 +491 @@
 }
 
+# define DUMPPICS 0
 # define LIMIT_XRANGE(X, IMAGE) { X = PS_MIN(PS_MAX(0, X), IMAGE->numCols); }
 # define LIMIT_YRANGE(Y, IMAGE) { Y = PS_MIN(PS_MAX(0, Y), IMAGE->numRows); }
@@ -524 +525 @@
     int binning = 2;
     float sigma_thresh = 3.0;
-    int max_iter = 5;
-    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
+    int max_iter = 1; // XXX with isophot masking, we only want to do a single pass
 
     // Temporary images.
@@ -591 +591 @@
         }
 
-        // Modify my mask if we're above the significance threshold
-        for (int y = 0; y < dy; y++) {
-            for (int x = 0; x < dx; x++) {
-                if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01)) continue;
-
-                float noise = binning * sqrt(noise_factor * myvar->data.F32[y][x]);
-                float value = edges->data.F32[y][x] / noise;
-
-                if (value > sigma_thresh ) {
-                    mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
-                    nCRpix ++;
-                }
-            }
-        }
-
-# if 1
+        // coordinate of peak in subimage pixels:
+        int xPeak = peak->x - xs;
+        int yPeak = peak->y - ys;
+
+        // Modify my mask if we're above the significance threshold, but only for connected pixels
+        nCRpix = psphotMaskCosmicRayConnected (xPeak, yPeak, mymask, myvar, edges, binning, sigma_thresh);
+        
+# if DUMPPICS
         psphotSaveImage (NULL, mypix,   "crmask.pix.fits");
 # endif
-
+        
+// XXX do not repair the pixels in isophot version
+# if 0
         // "Repair" Masked pixels for the next round.
         for (int y = 1; y < dy - 1; y++) {
@@ -629 +623 @@
             }
         }
-
-# if 1
+# endif
+
+# if DUMPPICS
         fprintf (stderr, "CRMASK %d %d %d %d %d\n", xs, ys, dx, dy, iteration);
         psphotSaveImage (NULL, mypix,   "crmask.fix.fits");
@@ -642 +637 @@
     }
 
+# if 0
     // A solitary masked pixel is likely a lie. Remove those
     // XXX can't we use nCRpix == 1 to test for these?
@@ -662 +658 @@
         }
     }
+# endif
 
     // transfer temporary mask to real mask & count masked pixels
@@ -677 +674 @@
     if (nCRpix > 1) {
         source->mode |= PM_SOURCE_MODE_CR_LIMIT;
-    }
+        source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED;
+    }
+    // fprintf (stderr, "CRMASK %d %d %d %d %d\n", peak->x, peak->y, dx, dy, nCRpix);
 
     psFree(mypix);
@@ -753 +752 @@
 }
 
+# define VERBOSE 0
+int psphotMaskCosmicRayConnected (int xPeak, int yPeak, psImage *mymask, psImage *myvar, psImage *edges, int binning, float sigma_thresh) {
+    
+    int xLo, xRo;
+    int nCRpix = 0;
+
+    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
+    
+    // mark the pixels in this row to the left, then the right. stay within footprint
+    int xL = xPeak; // find the range of valid pixels in this row
+    int xR = xPeak;
+    for (int ix = xPeak; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix--) {
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xL = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    for (int ix = xPeak; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] & 0x01); ix++) {
+        float noise = binning * sqrt(noise_factor * myvar->data.F32[yPeak][ix]);
+        float value = edges->data.F32[yPeak][ix] / noise;
+        if (value < sigma_thresh ) break;
+        mymask->data.PS_TYPE_IMAGE_MASK_DATA[yPeak][ix] |= 0x40;
+        xR = ix;
+        nCRpix ++;
+        if (VERBOSE) fprintf (stderr, "mark %d,%d (%d) : %d - %d\n", ix, yPeak, nCRpix, xL, xR);
+    }
+    // xL and xR mark the first and last valid pixel in the row
+
+    // for each of the neighboring rows, mark the high pixels if they touch the range xL to xR 
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // first go down:
+    for (int iy = yPeak - 1; iy >= 0; iy--) {
+
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    xLo = PS_MAX(xL - 1, 0);
+    xRo = PS_MIN(xR + 1, mymask->numCols);
+
+    // next go up:
+    for (int iy = yPeak + 1; iy < mymask->numRows; iy++) {
+
+        int xLn = -1;
+        int xRn = -1;
+        int newPix = 0;
+
+        // mark the pixels in the good range
+        for (int ix = xLo; ix < xRo; ix++) {
+            if (!(mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01)) continue; // only use pixels in the footprint
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) continue;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;                // last valid pixel in this row
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark C %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the left of the good range
+        for (int ix = xLo; (ix >= 0) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix--) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xRn == -1) xRn = ix; // last valid pixel in this row
+            xLn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark L %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+
+        // mark the pixels to the right of the good range
+        for (int ix = xRo; (ix < mymask->numCols) && (mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & 0x01); ix++) {
+            float noise = binning * sqrt(noise_factor * myvar->data.F32[iy][ix]);
+            float value = edges->data.F32[iy][ix] / noise;
+            if (value < sigma_thresh ) break;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= 0x40;
+            if (xLn == -1) xLn = ix; // first valid pixel in this row
+            xRn = ix;
+            nCRpix ++;
+            newPix ++;
+            if (VERBOSE) fprintf (stderr, "mark R %d,%d (%d) : %d - %d | %d - %d | %d - %d \n", ix, iy, nCRpix, xL, xR, xLo, xRo, xLn, xRn);
+        }
+        if (newPix == 0) break;
+        xLo = PS_MAX(xLn - 1, 0);
+        xRo = PS_MIN(xRn + 1, mymask->numCols);
+    }
+
+    return nCRpix;
+}
+
 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {