Changeset 26646

Timestamp:

Jan 20, 2010, 2:49:40 PM (16 years ago)

Author:

eugene

Message:

clean up of CZW code assuming that footprints are not corrupted, but keep it commented out until the size checks can be improved

File:

• branches/eam_branches/20091201/psphot/src/psphotSourceSize.c (modified) (13 diffs)

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

@@ -15 +15 @@
 } psphotSourceSizeOptions;
 
+// local functions:
 bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf);
 bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
@@ -21 +22 @@
 bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
+bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
+
+bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        pmPeak *peak = source->peak;
+        pmFootprint *footprint = peak->footprint;
+        if (!footprint) continue;
+        for (int j = 0; j < footprint->spans->n; j++) {
+            pmSpan *sp = footprint->spans->data[j];
+            psAssert (sp, "missing span");
+        }
+    }
+    return true;
+}
 
 // we need to call this function after sources have been fitted to the PSF model and
@@ -70 +87 @@
 
     // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness
-    // of and object.  We need to model this distribution for the PSF stars before we can test
+    // of an object.  We need to model this distribution for the PSF stars before we can test
     // the significance for a specific object
     // XXX move this to the code that generates the PSF?
@@ -90 +107 @@
 }
 
-bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
+// This attempt to mask the cosmic rays used the isophotal boundary
+bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
 
     // replace the source flux
@@ -103 +121 @@
     pmFootprint *footprint = peak->footprint;
     if (!footprint) {
+      psTrace("psphot.czw",2,"Using isophot CR mask code.");
+      
         // if we have not footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
@@ -109 +129 @@
 
     if (!footprint->spans) {
+      psTrace("psphot.czw",2,"Using isophot CR mask code.");
+      
         // if we have no footprint, use the old code to mask by isophot
         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
         return true;
     }
-
+    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     // 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.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
-        }
-    }
+    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.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; 
+        } 
+    } 
     return true;
 }
@@ -232 +254 @@
         float dMag = source->psfMag - apMag;
 
-        psVectorAppend (Ap, dMag);
-        psVectorAppend (ApErr, source->errMag);
+        psVectorAppend (Ap, 100, dMag);
+        psVectorAppend (ApErr, 100, source->errMag);
     }
 
@@ -247 +269 @@
     options->ApResid = stats->robustMedian;
     options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
-    // XXX this is quite arbitrary...
-    if (!isfinite(options->ApSysErr)) options->ApSysErr = 0.01;
     psLogMsg ("psphot", PS_LOG_DETAIL, "psf - Sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
 
@@ -375 +395 @@
         // Anything within this region is a probably PSF-like object. Saturated stars may land
         // in this region, but are detected elsewhere on the basis of their peak value.
-        bool isPSF = (fabs(nSigma) < options->nSigmaApResid) && (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) && (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY);
+        bool isPSF = ((fabs(nSigma) < options->nSigmaApResid) &&
+                      (fabs(Mxx - psfClump->X) < options->nSigmaMoments*psfClump->dX) &&
+                      (fabs(Myy - psfClump->Y) < options->nSigmaMoments*psfClump->dY));
         if (isPSF) {
-            Npsf ++;
-            continue;
+          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g PSF\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
+                  options->nSigmaApResid,options->nSigmaMoments);
+          Npsf ++;
+          continue;
         }
 
@@ -385 +410 @@
         // XXX this rule is not great
         if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) {
+
+          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g CR\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
+                  options->nSigmaApResid,options->nSigmaMoments);
             source->mode |= PM_SOURCE_MODE_DEFECT;
             Ncr ++;
@@ -393 +422 @@
         // just large saturated regions.
         if (source->mode & PM_SOURCE_MODE_SATSTAR) {
+          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g SAT\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
+                  options->nSigmaApResid,options->nSigmaMoments);
+          
             Nsat ++;
             continue;
@@ -400 +433 @@
         bool isEXT = (nSigma > options->nSigmaApResid) || ((Mxx > psfClump->X) && (Myy > psfClump->Y));
         if (isEXT) {
+          psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Ext\t%g %g\n",
+                  source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
+                  options->nSigmaApResid,options->nSigmaMoments);
+          
             source->mode |= PM_SOURCE_MODE_EXT_LIMIT;
             Next ++;
             continue;
         }
-
+        psTrace("psphot.czw",4,"CLASS: %g %g\t%g %g  %g %g  %g %g\t%g %g\t%g Unk\t%g %g\n",
+                source->peak->xf,source->peak->yf,Mxx,Myy,psfClump->X,psfClump->Y,psfClump->dX,psfClump->dY,apMag,dMag,nSigma,
+                options->nSigmaApResid,options->nSigmaMoments);
+        
         psWarning ("sourse size was missed for %f,%f : %f %f -- %f\n", source->peak->xf, source->peak->yf, Mxx, Myy, nSigma);
         Nmiss ++;
@@ -488 +528 @@
 }
 
+// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
+// if enough of the object is detected as 'sharp', consider the object a cosmic ray
 bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
+
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+
+        // Integer position of peak
+        int xPeak = source->peak->xf - source->pixels->col0 + 0.5;
+        int yPeak = source->peak->yf - source->pixels->row0 + 0.5;
+
+        // Skip sources which are too close to a boundary.  These are mostly caught as DEFECT
+        if (xPeak < 1 || xPeak > source->pixels->numCols - 2 ||
+            yPeak < 1 || yPeak > source->pixels->numRows - 2) {
+            psTrace("psphot", 7, "Not calculating crNsigma due to edge\n");
+            //      psTrace("psphot.czw", 2, "Not calculating crNsigma due to edge\n");
+            continue;
+        }
+        
+        // this source is thought to be a cosmic ray.  flag the detection and mask the pixels
+        if (source->mode & PM_SOURCE_MODE_DEFECT) {
+            // XXX this is running slowly and is too agressive, but it more-or-less works
+            // psphotMaskCosmicRayCZW(readout, source, options->crMask);
+            
+            // XXX these are the old versions which we are not using
+            // psphotMaskCosmicRay (readout->mask, source, maskVal, crMask);
+            // psphotMaskCosmicRay_Old (source, maskVal, crMask);
+        }
+    }
+    
+    // now that we have masked pixels associated with CRs, we can grow the mask
+    if (options->grow > 0) {
+        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask
+        psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow);
+        psImageConvolveSetThreads(oldThreads);
+        if (!newMask) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask");
+            return false;
+        }
+        psFree(readout->mask);
+        readout->mask = newMask;
+    }
+
+    // XXX test : save the mask image
+    if (0) {
+        psphotSaveImage (NULL, readout->mask,   "mask.fits");
+    }
+    return true;
+}
+
+// Comments by CZW 20091209 : Mechanics of how to identify CR pixels taken from "Cosmic-Ray
+// Rejection by Laplacian Edge Detection" by Pieter van Dokkum, arXiv:astro-ph/0108003.  This
+// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
+// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
+bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
+
+    // Get the actual images and information about the peak.
+    psImage *mask = readout->mask;
+    pmPeak *peak = source->peak;
+    pmFootprint *footprint = peak->footprint;
+
+    int xm = footprint->bbox.x0;
+    int xM = footprint->bbox.x1;
+    int ym = footprint->bbox.y0;
+    int yM = footprint->bbox.y1;
+
+    if (xm < 0) { xm = 0; }
+    if (ym < 0) { ym = 0; }
+    if (xM > mask->numCols) { xM = mask->numCols; }
+    if (yM > mask->numRows) { yM = mask->numRows; }
+    int dx = xM - xm;
+    int dy = yM - ym;
+
+    // Bounding boxes are inclusive of final pixel
+    // XXX ensure dx,dy do not become too large here
+    dx++;
+    dy++;
+
+    psImage *image= readout->image;
+    psImage *variance = readout->variance;
+      
+    int binning = 1;
+    float sigma_thresh = 2.0;
+    int iteration = 0;
+    int max_iter = 5;
+    float noise_factor = 5.0 / 4.0;  // Intrinsic to the Laplacian making noise spikes spikier.
+
+    // Temporary images.
+    psImage *mypix  = psImageAlloc(dx,dy,image->type.type);
+    psImage *myvar  = psImageAlloc(dx,dy,image->type.type);
+    psImage *binned = psImageAlloc(dx * binning,dy * binning,image->type.type);
+    psImage *conved = psImageAlloc(dx * binning,dy * binning,image->type.type);
+    psImage *edges  = psImageAlloc(dx,dy,image->type.type);
+    psImage *mymask = psImageAlloc(dx,dy,PS_TYPE_IMAGE_MASK);
+      
+    int x,y;
+    // Load my copy of things.
+    for (x = 0; x < dx; x++) {
+        for (y = 0; y < dy; y++) {
+            psImageSet(mypix,x,y,psImageGet(image,x+xm,y+ym));
+            psImageSet(myvar,x,y,psImageGet(variance,x+xm,y+ym));
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0x00;
+        }
+    }
+    // Mask so I can see on the output image where the footprint is.
+    for (int i = 0; i < footprint->spans->n; i++) {
+        pmSpan *sp = footprint->spans->data[i];
+        for (int j = sp->x0; j <= sp->x1; j++) {
+            y = sp->y - ym;
+            x = j - xm;
+            mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x01;
+        }
+    }
+
+    int CRpix_count = 0;      
+    do {
+        CRpix_count = 0;
+        // Zero out my temp images.
+        for (x = 0; x < binning * dx; x++) {
+            for (y = 0; y < binning * dy; y++) {
+                psImageSet(binned,x,y,0.0);
+                psImageSet(conved,x,y,0.0);
+                psImageSet(edges,x/binning,y/ binning,0.0);
+            }
+        }      
+        // Make subsampled image. Maybe this should be called "unbinned" or something
+        for (x = 0; x < binning * dx; x++) {
+            for (y = 0; y < binning * dy; y++) {
+                psImageSet(binned,x,y,psImageGet(mypix,x / binning,y / binning));
+            }
+        }
+        // Apply Laplace transform (kernel = [[0 -0.25 0][-0.25 1 -0.25][0 -0.25 0]]), clipping at zero
+        for (x = 1; x < dx - 1; x++) {
+            for (y = 1; y < dy - 1; y++) {
+                psImageSet(conved,x,y,psImageGet(binned,x,y) - 0.25 *
+                           (psImageGet(binned,x-1,y) + psImageGet(binned,x+1,y) +
+                            psImageGet(binned,x,y-1) + psImageGet(binned,x,y+1)));
+                if (psImageGet(conved,x,y) < 0.0) {
+                    psImageSet(conved,x,y,0.0);
+                }
+            }
+        }
+        // Create an edge map by rebinning
+        for (x = 0; x < binning * dx; x++) {
+            for (y = 0; y < binning * dy; y++) {
+                psImageSet(edges,x / binning, y / binning,
+                           psImageGet(edges, x / binning, y / binning) +
+                           psImageGet(conved,x,y));
+            }
+        }
+        // Modify my mask if we're above the significance threshold
+        for (x = 0; x < dx; x++) {
+            for (y = 0; y < dy; y++) {
+                if ( psImageGet(edges,x,y) / (binning * sqrt(noise_factor * psImageGet(myvar,x,y))) > sigma_thresh ) {
+                    if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) {
+                        mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= 0x40;
+                        CRpix_count++;
+                    }
+                }
+            }
+        }
+
+        // "Repair" Masked pixels for the next round.
+        for (x = 1; x < dx - 1; x++) {
+            for (y = 1; y < dy - 1; y++) {
+                if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
+                    psImageSet(mypix,x,y,0.25 *
+                               (psImageGet(mypix,x-1,y) + psImageGet(mypix,x+1,y) +
+                                psImageGet(mypix,x,y-1) + psImageGet(mypix,x,y+1)));
+                }
+            }
+        }
+        
+# if 0
+        if ((psTraceGetLevel("psphot.czw") >= 2)&&(abs(xm - 2770) < 5)&&(abs(ym - 2581) < 5)&&(iteration == 0)) {
+          psTrace("psphot.czw",2,"TRACEMOTRON %d %d %d %d %d\n",xm,ym,dx,dy,iteration);
+          psphotSaveImage (NULL, mypix,   "czw.pix.fits");
+          psphotSaveImage (NULL, myvar,   "czw.var.fits");
+          psphotSaveImage (NULL, binned,  "czw.binn.fits");
+          psphotSaveImage (NULL, conved,  "czw.conv.fits");
+          psphotSaveImage (NULL, edges,   "czw.edge.fits");
+          psphotSaveImage (NULL, mymask,  "czw.mask.fits");
+        }
+# endif
+
+        psTrace("psphot.czw",2,"Iter: %d Count: %d",iteration,CRpix_count);
+        iteration++;
+    } while ((iteration < max_iter)&&(CRpix_count > 0));
+
+    // A solitary masked pixel is likely a lie. Remove those.
+    for (x = 0; x < dx; x++) {
+        for (y = 0; y < dy; y++) {
+            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
+                if ((x-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x-1] & 0x40)) {
+                    continue;
+                }
+                if ((y-1 >= 0)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y-1][x] & 0x40)) {
+                    continue;
+                }
+                if ((x+1 < dx)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x+1] & 0x40)) {
+                    continue;
+                }
+                if ((y+1 < dy)&&(mymask->data.PS_TYPE_IMAGE_MASK_DATA[y+1][x] & 0x40)) {
+                    continue;
+                }
+                mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] ^= 0x40;
+            }
+        }
+    }
+
+    // transfer temporary mask to real mask
+    for (x = 0; x < dx; x++) {
+        for (y = 0; y < dy; y++) {
+            if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x40) {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= maskVal;
+            }
+        }
+    }
+      
+    // XXX if we decide this REALLY is a cosmic ray, set the CR_LIMIT bit
+    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
+
+    psFree(mypix);
+    psFree(myvar);
+    psFree(binned);
+    psFree(conved);
+    psFree(edges);
+    psFree(mymask);
+      
+    return true;
+}
+
+// this was an old attempt to identify cosmic rays based on the peak curvature
+bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
 
     // classify the sources based on the CR test (place this in a function?)