Index: /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c (revision 26645) +++ /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c (revision 26646) @@ -15,4 +15,5 @@ } psphotSourceSizeOptions; +// local functions: bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf); bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); @@ -21,4 +22,20 @@ 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,5 +87,5 @@ // 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,5 +107,6 @@ } -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,4 +121,6 @@ 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,21 +129,23 @@ 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,6 +254,6 @@ float dMag = source->psfMag - apMag; - psVectorAppend (Ap, dMag); - psVectorAppend (ApErr, source->errMag); + psVectorAppend (Ap, 100, dMag); + psVectorAppend (ApErr, 100, source->errMag); } @@ -247,6 +269,4 @@ 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,8 +395,13 @@ // 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,4 +410,8 @@ // 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,4 +422,8 @@ // 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,9 +433,16 @@ 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,5 +528,238 @@ } +// 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?)