Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 26617) +++ trunk/psphot/src/psphotSourceSize.c (revision 26618) @@ -22,5 +22,4 @@ bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); -bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal,psphotSourceSizeOptions *options); // we need to call this function after sources have been fitted to the PSF model and // subtracted. To determine the CR-nature, this function examines the 9 pixels in the 3x3 @@ -71,5 +70,5 @@ // We are using the value psfMag - 2.5*log10(moment.Sum) as a measure of the extendedness - // of an object. We need to model this distribution for the PSF stars before we can test + // of and 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? @@ -104,6 +103,4 @@ 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); @@ -112,28 +109,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; */ -/* } */ -/* } */ + + // 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; + } + } return true; } - - - bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { @@ -240,6 +232,6 @@ float dMag = source->psfMag - apMag; - psVectorAppend (Ap, 100, dMag); - psVectorAppend (ApErr, 100, source->errMag); + psVectorAppend (Ap, dMag); + psVectorAppend (ApErr, source->errMag); } @@ -255,4 +247,6 @@ 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); @@ -381,13 +375,8 @@ // 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) { - 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; + Npsf ++; + continue; } @@ -396,8 +385,4 @@ // 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 ++; @@ -408,8 +393,4 @@ // 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; @@ -419,16 +400,9 @@ 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 ++; @@ -521,22 +495,20 @@ pmSource *source = sources->data[i]; -/* // skip source if it was already measured */ -/* if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { */ -/* psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); */ -/* psTrace("psphot.czw", 2, "Not calculating source size since it has already been measured\n"); */ -/* continue; */ -/* } */ - -/* // source must have been subtracted */ -/* if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { */ -/* source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; */ -/* psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); */ -/* psTrace("psphot.czw", 2, "Not calculating source size since source is not subtracted\n"); */ -/* continue; */ -/* } */ - -// psF32 **resid = source->pixels->data.F32; - // psF32 **variance = source->variance->data.F32; - // psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } + + // source must have been subtracted + if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + continue; + } + + psF32 **resid = source->pixels->data.F32; + psF32 **variance = source->variance->data.F32; + psImageMaskType **mask = source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA; // Integer position of peak @@ -548,111 +520,90 @@ 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; - } - -/* // Skip sources with masked pixels. These are mostly caught as DEFECT */ -/* bool keep = true; */ -/* for (int iy = -1; (iy <= +1) && keep; iy++) { */ -/* for (int ix = -1; (ix <= +1) && keep; ix++) { */ -/* if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { */ -/* keep = false; */ -/* } */ -/* } */ -/* } */ -/* if (!keep) { */ -/* psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); */ -/* // psTrace("psphot.czw", 2, "Not calculating crNsigma due to masked pixels\n"); */ -/* continue; */ -/* } */ - -/* // psTrace("psphot.czw", 2, "Actually doing something to mask a CR. \n"); */ -/* // Compare the central pixel with those on either side, for the four possible lines through it. */ - -/* // Soften variances (add systematic error) */ -/* float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances */ - -/* // Across the middle: y = 0 */ -/* float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; */ -/* float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; */ -/* float nX = cX / sqrtf(dcX + softening); */ - -/* // Up the centre: x = 0 */ -/* float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; */ -/* float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; */ -/* float nY = cY / sqrtf(dcY + softening); */ - -/* // Diagonal: x = y */ -/* float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; */ -/* float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; */ -/* float nL = cL / sqrtf(dcL + softening); */ - -/* // Diagonal: x = - y */ -/* float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; */ -/* float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; */ -/* float nR = cR / sqrtf(dcR + softening); */ - -/* // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) */ -/* // Ndof = 4 ? (four measurements, no free parameters) */ -/* // XXX this value is going to be biased low because of systematic errors. */ -/* // we need to calibrate it somehow */ -/* // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); */ - -/* // not strictly accurate: overcounts the chisq contribution from the center pixel (by */ -/* // factor of 4); also biases a bit low if any pixels are masked */ -/* // XXX I am not sure I want to keep this value... */ -/* source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); */ - -/* float fCR = 0.0; */ -/* int nCR = 0; */ -/* if (nX > 0.0) { */ -/* fCR += nX; */ -/* nCR ++; */ -/* } */ -/* if (nY > 0.0) { */ -/* fCR += nY; */ -/* nCR ++; */ -/* } */ -/* if (nL > 0.0) { */ -/* fCR += nL; */ -/* nCR ++; */ -/* } */ -/* if (nR > 0.0) { */ -/* fCR += nR; */ -/* nCR ++; */ -/* } */ -/* source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; */ -/* source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; */ - -/* if (!isfinite(source->crNsigma)) { */ -/* continue; */ -/* } */ - - //psImageMaskType crMask = 0x08; - + continue; + } + + // Skip sources with masked pixels. These are mostly caught as DEFECT + bool keep = true; + for (int iy = -1; (iy <= +1) && keep; iy++) { + for (int ix = -1; (ix <= +1) && keep; ix++) { + if (mask[yPeak+iy][xPeak+ix] & options->maskVal) { + keep = false; + } + } + } + if (!keep) { + psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + continue; + } + + // Compare the central pixel with those on either side, for the four possible lines through it. + + // Soften variances (add systematic error) + float softening = options->soft * PS_SQR(source->peak->flux); // Softening for variances + + // Across the middle: y = 0 + float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; + float dcX = 4*variance[yPeak][xPeak] + variance[yPeak+0][xPeak-1] + variance[yPeak+0][xPeak+1]; + float nX = cX / sqrtf(dcX + softening); + + // Up the centre: x = 0 + float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; + float dcY = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak+0] + variance[yPeak+1][xPeak+0]; + float nY = cY / sqrtf(dcY + softening); + + // Diagonal: x = y + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + float dcL = 4*variance[yPeak][xPeak] + variance[yPeak-1][xPeak-1] + variance[yPeak+1][xPeak+1]; + float nL = cL / sqrtf(dcL + softening); + + // Diagonal: x = - y + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + float dcR = 4*variance[yPeak][xPeak] + variance[yPeak+1][xPeak-1] + variance[yPeak-1][xPeak+1]; + float nR = cR / sqrtf(dcR + softening); + + // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) + // Ndof = 4 ? (four measurements, no free parameters) + // XXX this value is going to be biased low because of systematic errors. + // we need to calibrate it somehow + // source->psfProb = gsl_sf_gamma_inc_Q (2, 0.5*chisq); + + // not strictly accurate: overcounts the chisq contribution from the center pixel (by + // factor of 4); also biases a bit low if any pixels are masked + // XXX I am not sure I want to keep this value... + source->psfChisq = PS_SQR(nX) + PS_SQR(nY) + PS_SQR(nL) + PS_SQR(nR); + + float fCR = 0.0; + int nCR = 0; + if (nX > 0.0) { + fCR += nX; + nCR ++; + } + if (nY > 0.0) { + fCR += nY; + nCR ++; + } + if (nL > 0.0) { + fCR += nL; + nCR ++; + } + if (nR > 0.0) { + fCR += nR; + nCR ++; + } + source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + if (!isfinite(source->crNsigma)) { + continue; + } + // this source is thought to be a cosmic ray. flag the detection and mask the pixels - if (source->mode & PM_SOURCE_MODE_DEFECT) { - // if (source->crNsigma > options->nSigmaCR) { - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - pmFootprint *footprint = peak->footprint; - // if ((footprint) && (footprint->spans) && !(footprint->npix < 0)) { - psTrace("psphot.czw",2,"Footprint sync: %d %d\n",peak->x,peak->y); - psTrace("psphot.czw",2,"Footprint info: %p %p %d %d %d\n",(void *)footprint, (void *) footprint->spans,footprint->id,peak->x,peak->y); - if ((footprint) && (footprint->spans)) { - psTrace("psphot.czw",2," Footprint good?: %p %p %d %g %g\n",(void *)footprint, (void *) footprint->spans, - footprint->id,footprint->bbox.x0,footprint->bbox.y0); - } - psImageMaskType maskVal = 0xbf; // HACK - maskVal = 0x80; - psphotMaskCosmicRayCZW(readout, source, maskVal,options); - - // psphotMaskCosmicRay(readout->mask, source, maskVal, maskVal); - - // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); - // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); - } - } - + if (source->crNsigma > options->nSigmaCR) { + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + // XXX still testing... : psphotMaskCosmicRay (readout->mask, source, maskVal, crMask); + // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); + } + } + // now that we have masked pixels associated with CRs, we can grow the mask if (options->grow > 0) { @@ -667,208 +618,4 @@ readout->mask = newMask; } - - if (psTraceGetLevel("psphot.czw") >= 2) { - psphotSaveImage (NULL, readout->mask, "mask.fits"); - } return true; } - - -// 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,psphotSourceSizeOptions *options) { - // Get the actual images and information about the peak. - psImage *mask = readout->mask; - pmPeak *peak = source->peak; - pmFootprint *footprint = peak->footprint; - - // Try to skip sources with bad footprints. - if ((footprint) && (footprint->spans) && footprint->id > 0 && footprint->id < 100000) { - 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; - - // Further try to skip bad footprints. - if ((footprint->npix > 0) && ((dx)*(dy) > 0) && (dx < 4000) && (dy < 4000) && (dx > 0) && (dy > 0)) { - // Bounding boxes are inclusive of final pixel. - 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 teh 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 ((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"); */ - -/* } */ - - 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; - } - } - } - - - - 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; - } -/* // Look at the footprint */ -/* if (mymask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & 0x01) { */ -/* mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= 0x04; */ -/* } */ -/* // Look at the bounding box. */ -/* mask->data.PS_TYPE_IMAGE_MASK_DATA[y+ym+mask->row0][x+xm+mask->col0] |= 0x08; */ - } - } - - - psFree(mypix); - psFree(myvar); - psFree(binned); - psFree(conved); - psFree(edges); - psFree(mymask); - - } - else { - psTrace("psphot.czw", 2, "BBox does not match npix: %d %d",(xM - xm) * (yM - ym),footprint->npix); - } - } - else { - psTrace("psphot.czw", 2, "No valid footprint for source (%g %g)",peak->xf,peak->yf); - } - return(true); -} - - -