Index: trunk/psphot/src/psphotSourceSize.c =================================================================== --- trunk/psphot/src/psphotSourceSize.c (revision 19915) +++ trunk/psphot/src/psphotSourceSize.c (revision 20144) @@ -2,5 +2,6 @@ # include -float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro); +static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask, + psMaskType maskVal, const pmModel *model, float Ro); // we need to call this function after sources have been fitted to the PSF model and @@ -10,7 +11,6 @@ // deviation from the psf model at the r = FWHM/2 position -bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) +bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, long first) { - bool status; @@ -41,8 +41,14 @@ // skip source if it was already measured - if (isfinite(source->crNsigma)) continue; + if (isfinite(source->crNsigma)) { + psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n"); + continue; + } // source must have been subtracted - if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) continue; + if (!(source->mode & PM_SOURCE_MODE_SUBTRACTED)) { + psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since source is not subtracted\n"); + continue; + } psF32 **resid = source->pixels->data.F32; @@ -51,11 +57,13 @@ // check for extendedness: measure the delta flux significance at the 1 sigma contour - source->extNsigma = psphotModelContour (source->pixels, source->weight, source->maskObj, maskVal, 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? if (source->extNsigma > EXT_NSIGMA_LIMIT) { - source->mode |= PM_SOURCE_MODE_EXT_LIMIT; - } - + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + } + + // Integer position of peak int xPeak = source->peak->xf - source->pixels->col0 + 0.5; int yPeak = source->peak->yf - source->pixels->row0 + 0.5; @@ -63,8 +71,9 @@ // XXX for now, skip sources which are too close to a boundary // XXX raise a flag? - if (xPeak < 1) continue; - if (xPeak > source->pixels->numCols - 2) continue; - if (yPeak < 1) continue; - if (yPeak > source->pixels->numRows - 2) continue; + 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"); + continue; + } // XXX for now, just skip any sources with masked pixels @@ -73,40 +82,35 @@ 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; - - if (mask[yPeak][xPeak] & maskVal) continue; - - // XXX need to deal with edge peaks... and mask - 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); - } + if (mask[yPeak+iy][xPeak+ix] & 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. + + // Across the middle: y = 0 + 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]; + float nX = cX / sqrtf(dcX); + + // Up the centre: x = 0 + 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]; + float nY = cY / sqrtf(dcY); + + // Diagonal: x = y + float cL = 2*resid[yPeak][xPeak] - resid[yPeak-1][xPeak-1] - resid[yPeak+1][xPeak+1]; + float dcL = 4*weight[yPeak][xPeak] + weight[yPeak-1][xPeak-1] + weight[yPeak+1][xPeak+1]; + float nL = cL / sqrtf(dcL); + + // Diagonal: x = - y + float cR = 2*resid[yPeak][xPeak] - resid[yPeak+1][xPeak-1] - resid[yPeak-1][xPeak+1]; + float dcR = 4*weight[yPeak][xPeak] + weight[yPeak+1][xPeak-1] + weight[yPeak-1][xPeak+1]; + float nR = cR / sqrtf(dcR); // P(chisq > chisq_obs; Ndof) = gamma_Q (Ndof/2, chisq/2) @@ -118,6 +122,6 @@ // 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); + // 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; @@ -156,4 +160,5 @@ psImage *weight = source->weight; + // XXX This should be a recipe variable # define SN_LIMIT 5.0 @@ -232,60 +237,69 @@ // 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 +// deviation in sigmas. This is measured on the residual image - should we ignore negative // deviations? - float psphotModelContour (psImage *image, psImage *weight, psImage *mask, psMaskType maskVal, pmModel *model, float Ro) { - - psF32 *PAR = model->params->data.F32; - - // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY; +static float psphotModelContour(const psImage *image, const psImage *weight, const psImage *mask, + psMaskType maskVal, const pmModel *model, float Ro) +{ + psF32 *PAR = model->params->data.F32; // Model parameters + float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters + + // We treat the contour as an ellipse: + // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; - // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]; - // y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A]; - - // min/max value of x is where T -> 0 - // solve this for x2: - float Q = Ro * PS_SQR(PAR[PM_PAR_SXX]) / (1.0 - PS_SQR(PAR[PM_PAR_SXX]*PAR[PM_PAR_SYY]*PAR[PM_PAR_SXY]) / 4.0); - if (Q < 0.0) return NAN; // ellipse is imaginary - - int xMax = sqrt(Q); - int xMin = -1.0*xMax; - - int nPts = 0; - float nSigma = 0.0; - - for (int x = xMin; x <= xMax; x++) { - float A = PS_SQR (1.0 / PAR[PM_PAR_SYY]); - float B = x * PAR[PM_PAR_SXY]; - float C = PS_SQR (x / PAR[PM_PAR_SXX]) - Ro; - + // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro + // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so: + // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2] + + // min/max value of x is where B^2 - 4AC = 0; solve this for x + float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); + if (Q < 0.0) { + // ellipse is imaginary + return NAN; + } + + int radius = sqrtf(Q) + 0.5; // Radius of ellipse + int nPts = 0; // Number of points in ellipse + float nSigma = 0.0; // + + for (int x = -radius; x <= radius; x++) { + // Polynomial coefficients + float A = PS_SQR (1.0 / syy); + float B = x * sxy; + float C = PS_SQR (x / sxx) - Ro; float T = PS_SQR(B) - 4*A*C; - if (T < 0.0) continue; - + if (T < 0.0) { + continue; + } + + // y position in source frame float yP = (-B + sqrt (T)) / (2.0 * A); float yM = (-B - sqrt (T)) / (2.0 * A); + // Get the closest pixel positions (image frame) int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - if (xPix < 0) continue; - if (xPix >= image->numCols) continue; - - if ((yPixM >= 0) && (yPixM < image->numRows)) { - if (!mask || !(mask->data.U8[yPixM][xPix] & maskVal)) { - float dSigma = image->data.F32[yPixM][xPix] / sqrt (weight->data.F32[yPixM][xPix]); - nSigma += dSigma; - nPts ++; - } - } - - if (yPixM == yPixP) continue; - - if ((yPixP >= 0) && (yPixP < image->numRows)) { - if (!mask || !(mask->data.U8[yPixP][xPix] & maskVal)) { - float dSigma = image->data.F32[yPixP][xPix] / sqrt (weight->data.F32[yPixP][xPix]); - nSigma += dSigma; - nPts ++; - } + if (xPix < 0 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(weight->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(weight->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; } }