Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 19286) +++ trunk/psphot/src/psphotSourceSize.c (revision 19915) @@ -2,5 +2,5 @@ # include -float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro); +float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro); // we need to call this function after sources have been fitted to the PSF model and @@ -51,5 +51,5 @@ // check for extendedness: measure the delta flux significance at the 1 sigma contour - source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, source->modelPSF, 1.0); + source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, maskVal, source->modelPSF, 1.0); // XXX prevent a source from being both CR and EXT? @@ -78,19 +78,35 @@ if (!keep) continue; + if (mask[yPeak][xPeak] & maskVal) continue; + // XXX need to deal with edge peaks... and mask - float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; - float cY = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak+0] - resid[yPeak+1][xPeak+0]; - float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; - float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; - - float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1]; - float dcY = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak+0] + weight[yPeak+1][xPeak+0]; - float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1]; - float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1]; - - float nX = cX / sqrt(dcX); - float nY = cY / sqrt(dcY); - float nL = cL / sqrt(dcL); - float nR = cR / sqrt(dcR); + float nX = 0; + float nY = 0; + float nL = 0; + float nR = 0; + + if (!(mask[yPeak+0][xPeak-1] & maskVal) && !(mask[yPeak+0][xPeak+1] & maskVal)) { + float cX = 2*resid[yPeak][xPeak] - resid[yPeak+0][xPeak-1] - resid[yPeak+0][xPeak+1]; + float dcX = 4*weight[yPeak][xPeak] + weight[yPeak+0][xPeak-1] + weight[yPeak+0][xPeak+1]; + nX = cX / sqrt(dcX); + } + + if (!(mask[yPeak-1][xPeak+0] & maskVal) && !(mask[yPeak+1][xPeak+0] & maskVal)) { + float cY = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak+0] - resid[yPeak+1][xPeak+0]; + float dcY = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak+0] + weight[yPeak+1][xPeak+0]; + nY = cY / sqrt(dcY); + } + + if (!(mask[yPeak-1][xPeak-1] & maskVal) && !(mask[yPeak+1][xPeak+1] & maskVal)) { + float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1]; + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + nL = cL / sqrt(dcL); + } + + if (!(mask[yPeak+1][xPeak-1] & maskVal) && !(mask[yPeak-1][xPeak+1] & maskVal)) { + float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1]; + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + nR = cR / sqrt(dcR); + } // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) @@ -100,46 +116,29 @@ // 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) + // 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 (nX) + PS_SQR (nR); float fCR = 0.0; - float fEXT = 0.0; int nCR = 0; - int nEXT = 0; if (nX > 0.0) { fCR += nX; nCR ++; - } else { - fEXT += nX; - nEXT ++; } if (nY > 0.0) { fCR += nY; nCR ++; - } else { - fEXT += nY; - nEXT ++; } if (nL > 0.0) { fCR += nL; nCR ++; - } else { - fEXT += nL; - nEXT ++; } if (nR > 0.0) { fCR += nR; nCR ++; - } else { - fEXT += nR; - nEXT ++; } source->crNsigma = (nCR > 0) ? fCR / nCR : 0.0; - // NOTE: abs needed to make the Nsigma value positive - - if (!isfinite(source->crNsigma)) { - fprintf (stderr, "."); - source->crNsigma = -6.0; - } + assert (isfinite(source->crNsigma)); // this source is thought to be a cosmic ray. flag the detection and mask the pixels @@ -212,4 +211,5 @@ } + // now that we have masked pixels associated with CRs, we can grow the mask if (grow > 0) { psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow); @@ -225,50 +225,14 @@ sources->n - first, psTimerMark ("psphot")); + psphotVisualPlotSourceSize (sources); + psphotVisualShowSourceSize (sources); + return true; } - - -// some possible classes: -// all 4 curvatures are highly negative : CR -// all 4 curvatures are highly positive : EXT - -// at least 2 are significantly negative, none are significantly positive : CR -// at least 2 are significantly positive, none are significantly negative : EXT - -// any are significantly negative, some may be significantly positive : CR -// any are significantly positive, none may be significantly positive : EXT - -/* Nn Np No - 4 0 0 CR_1 - 3 1 0 CR_1 - 3 0 1 CR_1 - 2 2 0 CR_2 - 2 1 1 CR_2 - 2 0 2 CR_2 - 1 3 0 CR_3 - 1 2 1 CR_3 - 1 1 2 CR_3 - 1 0 3 CR_3 - 0 4 0 EXT - 0 3 1 EXT - 0 2 2 EXT - 0 1 3 PSF - 0 0 4 PSF -*/ - -/* Alternatively, write a f(CR) = Sum(nX,etc if >0) */ - - -/* I can write the formal probability that the 4 measurements are consistent with a PSF - * based on how large the error values are (Nsigma -> erf -> P(Nsigma)). - * I should examine this value as a function of flux and also examine the distribution of - * the probabilities for the PSF sources. - */ - // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total // deviation in sigmas. This is measure on the residual image - should we ignore negative // deviations? -float psphotModelContour (psImage *image, psImage *weight, psImage *mask, pmModel *model, float Ro) { + float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro) { psF32 *PAR = model->params->data.F32; @@ -309,5 +273,5 @@ if ((yPixM >= 0) && (yPixM < image->numRows)) { - if (!mask || !mask->data.U8[yPixM][xPix]) { + if (!mask || !(mask->data.U8[yPixM][xPix] & maskVal)) { float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]); nSigma += dSigma; @@ -319,5 +283,5 @@ if ((yPixP >= 0) && (yPixP < image->numRows)) { - if (!mask || !mask->data.U8[yPixP][xPix]) { + if (!mask || !(mask->data.U8[yPixP][xPix] & maskVal)) { float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]); nSigma += dSigma;