Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 13081) +++ trunk/psphot/src/psphotSourceSize.c (revision 13086) @@ -1,18 +1,12 @@ # include "psphotInternal.h" +# include -// we need to call this function after sources have been fitted to the PSF -// model and subtracted. this function examines the 9 pixels in the 3x3 square containing the peak -// and compares the observed flux to the model. For the 4 pixel lines which contain the peak, we -// calculate the value +// we need to call this function after sources have been fitted to the PSF model and +// subtracted. this function examines the 9 pixels in the 3x3 square containing the peak and +// compares the observed flux to the model. -/* - * - * - * - * - * - */ +bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe) { -bool psphotSourceSize (pmReadout *readout, psArray *sources, psMetadata *recipe, bool add) { + FILE *f = fopen ("srcsize.dat", "w"); // loop over all source @@ -23,12 +17,21 @@ if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; - psF32 *resid = source->pixels->data.F32; - psF32 *weight = source->weight->data.F32; - psU8 *mask = source->maskObj->data.U8; + psF32 **resid = source->pixels->data.F32; + psF32 **weight = source->weight->data.F32; + psU8 **mask = source->maskObj->data.U8; - int xPeak = source->peak->xf - source->pixels->col0; - int yPeak = source->peak->yf - source->pixels->row0; + int xPeak = source->peak->xf - source->pixels->col0 + 0.5; + int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - // XXX need to deal with edge peaks... + // XXX for now, just skip any sources with masked pixels + 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]) { keep &= false; } + } + } + if (!keep) 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]; @@ -36,43 +39,101 @@ float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; - float dcX = 2*weight[yPeak][xPeak] - weight[yPeak+0][xPeak-1] - weight[yPeak+0][xPeak+1]; - float dcY = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak+0] - weight[yPeak+1][xPeak+0]; - float dcL = 2*weight[yPeak][xPeak] - weight[yPeak-1][xPeak-1] - weight[yPeak+1][xPeak+1]; - float dcR = 2*weight[yPeak][xPeak] - weight[yPeak+1][xPeak-1] - weight[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 chisq = PS_SQR (nX) + PS_SQR (nY) + PS_SQR (nX) + PS_SQR (nR); + + // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) + // Ndof = 4 ? (four measurements, no free parameters) + float Ppsf = gsl_sf_gamma_inc_Q (2, 0.5*chisq); + + 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 ++; + } + fCR = (nCR > 0) ? fCR / nCR : 0.0; + fEXT = (nEXT > 0) ? fEXT / nEXT : 0.0; // XXX plot this as a function of flux / magnitude? // XXX overlay outliers on an image? - // some possible classes: - // all 4 curvatures are highly negative : CR - // all 4 curvatures are highly positive : EXT + fprintf (f, "%f %f %f %f %f %f %f\n", source->peak->xf, source->peak->yf, source->peak->flux, + chisq, Ppsf, fCR, fEXT); + } + fclose (f); + return true; +} - // 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 +// 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 - */ +/* 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 +*/ - /* 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. - */ +/* 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. + */ +