Index: /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25358) +++ /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25359) @@ -2,6 +2,13 @@ # include -static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro); +bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal); + +bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr); + +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSyserr); +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr); + +float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, + psImageMaskType maskVal, const pmModel *model, float Ro); bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); @@ -14,5 +21,5 @@ // 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, pmPSF *psf, long first) { bool status; @@ -30,6 +37,6 @@ assert (status); - float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); - assert (status); + // float EXT_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); + // assert (status); int grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs @@ -45,144 +52,154 @@ } - // loop over all source + // 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 + // the significance for a specific object + // XXX move this to the code that generates the PSF? + // XXX store the results on pmPSF? + float ApResid, ApSysErr; + psphotSourceSizePSF (&ApResid, &ApSysErr, sources, psf, maskVal); + + // classify the sources based on ApResid and Moments (extended sources) + psphotSourceClass(readout, sources, recipe, psf, maskVal, ApResid, ApSysErr); + + // classify the sources based on the CR test (place this in a function?) for (int i = first; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip source if it was already measured - if (isfinite(source->crNsigma)) { - psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n"); - continue; - } - - // source must have been subtracted - if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) { + pmSource *source = sources->data[i]; + + // skip source if it was already measured + if (isfinite(source->crNsigma)) { + psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already 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 extNsigma,crNsigma 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; - - // check for extendedness: measure the delta flux significance at the 1 sigma contour - source->extNsigma = psphotModelContour(source->pixels, source->variance, 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; - } - - // Integer position of peak - int xPeak = source->peak->xf - source->pixels->col0 + 0.5; - int yPeak = source->peak->yf - source->pixels->row0 + 0.5; - - // XXX for now, skip sources which are too close to a boundary - // XXX raise a flag? - if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || - yPeak < 1 || yPeak > source->pixels->numRows - 2) { + psTrace("psphot", 7, "Not calculating extNsigma,crNsigma 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 + int xPeak = source->peak->xf - source->pixels->col0 + 0.5; + int yPeak = source->peak->yf - source->pixels->row0 + 0.5; + + // XXX for now, skip sources which are too close to a boundary + // XXX raise a flag? + if (xPeak < 1 || xPeak > source->pixels->numCols - 2 || + yPeak < 1 || yPeak > source->pixels->numRows - 2) { source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); - continue; - } - - // XXX for now, just skip any sources with masked pixels - // XXX raise a flag? - 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] & maskVal) { - keep = false; - } - } - } - if (!keep) { - psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); + psTrace("psphot", 7, "Not calculating crNsigma due to edge\n"); + continue; + } + + // XXX for now, just skip any sources with masked pixels + // XXX raise a flag? + 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] & maskVal) { + keep = false; + } + } + } + if (!keep) { + psTrace("psphot", 7, "Not calculating crNsigma due to masked pixels\n"); source->mode |= PM_SOURCE_MODE_SIZE_SKIPPED; - continue; - } - - // Compare the central pixel with those on either side, for the four possible lines through it. - - // Soften variances (add systematic error) - float softening = 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; - if (!isfinite(source->crNsigma)) { - continue; - } - - // this source is thought to be a cosmic ray. flag the detection and mask the pixels - if (source->crNsigma > CR_NSIGMA_LIMIT) { - // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); - // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); - } - } - + continue; + } + + // Compare the central pixel with those on either side, for the four possible lines through it. + + // Soften variances (add systematic error) + float softening = 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; + if (!isfinite(source->crNsigma)) { + continue; + } + + // this source is thought to be a cosmic ray. flag the detection and mask the pixels + if (source->crNsigma > CR_NSIGMA_LIMIT) { + // XXX still testing... : psphotMaskCosmicRay_New (readout->mask, source, maskVal, crMask); + // XXX acting strange... psphotMaskCosmicRay_Old (source, maskVal, crMask); + } + } + + // pause and wait for user input: + // continue, save (provide name), ?? + char key[10]; + fprintf (stdout, "[c]ontinue? "); + if (!fgets(key, 8, stdin)) { + psWarning("Unable to read option"); + } + // now that we have masked pixels associated with CRs, we can grow the mask if (grow > 0) { - bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask - psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow); - psImageConvolveSetThreads(oldThreads); - if (!newMask) { - psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); - return false; - } - psFree(readout->mask); - readout->mask = newMask; + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask + psImage *newMask = psImageConvolveMask(NULL, readout->mask, crMask, crMask, -grow, grow, -grow, grow); + psImageConvolveSetThreads(oldThreads); + if (!newMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); + return false; + } + psFree(readout->mask); + readout->mask = newMask; } psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", - sources->n - first, psTimerMark ("psphot.size")); + sources->n - first, psTimerMark ("psphot.size")); psphotVisualPlotSourceSize (sources); @@ -195,6 +212,6 @@ // deviation in sigmas. This is measured on the residual image - should we ignore negative // deviations? -static float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro) +float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, + psImageMaskType maskVal, const pmModel *model, float Ro) { psF32 *PAR = model->params->data.F32; // Model parameters @@ -211,6 +228,6 @@ float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); if (Q < 0.0) { - // ellipse is imaginary - return NAN; + // ellipse is imaginary + return NAN; } @@ -220,44 +237,44 @@ for (int x = -radius; x <= radius; x++) { - // Polynomial coefficients - // XXX Should we be using the centre of the pixel as x or x+0.5? - 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; - } - - // 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 || xPix >= image->numCols) { - continue; - } - - if (yPixM >= 0 && yPixM < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); - nSigma += dSigma; - nPts++; - } - - if (yPixM == yPixP) { - continue; - } - - if (yPixP >= 0 && yPixP < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); - nSigma += dSigma; - nPts++; - } + // Polynomial coefficients + // XXX Should we be using the centre of the pixel as x or x+0.5? + 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; + } + + // 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 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; + } } nSigma /= nPts; @@ -279,26 +296,26 @@ pmFootprint *footprint = peak->footprint; if (!footprint) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; } if (!footprint->spans) { - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRay_Old (source, maskVal, crMask); - return true; + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRay_Old (source, maskVal, crMask); + return true; } // 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; - } + 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; @@ -323,14 +340,14 @@ // mark the pixels in this row to the left, then the right for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } } @@ -338,34 +355,211 @@ // first go up: for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } } // next go down: for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } - return true; -} + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +// is this source extended or not? +bool psphotSourceSizeExtended (FILE *output, pmSource *source, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) { + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + pmSourceMagnitudes (source, psf, photMode, maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + + float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr); + + fprintf (output, "%f %f : %f %f : %f %f : %f %f\n", source->peak->xf, source->peak->yf, nSigma, dMag, source->psfMag, source->errMag, source->moments->Mxx, source->moments->Myy); + + return true; +} + +// model the apmifit distribution for the psf stars: +bool psphotSourceSizePSF (float *ApResid, float *ApSysErr, psArray *sources, pmPSF *psf, psImageMaskType maskVal) { + + // select the psf stars + psArray *psfstars = psArrayAllocEmpty (100); + + psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32); + + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + psArrayAdd (psfstars, 100, source); + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + + psVectorAppend (Ap, 100, dMag); + psVectorAppend (ApErr, 100, source->errMag); + } + + // model the distribution as a mean or median value and a systematic error from that value: + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + psVectorStats (stats, Ap, NULL, NULL, 0); + + psVector *dAp = psVectorAlloc (Ap->n, PS_TYPE_F32); + for (int i = 0; i < Ap->n; i++) { + dAp->data.F32[i] = Ap->data.F32[i] - stats->robustMedian; + } + + *ApResid = stats->robustMedian; + *ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); + fprintf (stderr, "psf-sum: %f +/- %f\n", *ApResid, *ApSysErr); + + return true; +} + +// classify sources based on the combination of psf-mag, Mxx, Myy +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) { + + bool status; + pmPSFClump psfClump; + char regionName[64]; + + int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS"); + for (int i = 0; i < nRegions; i ++) { + snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", i); + psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName); + psAssert (regionMD, "regions must be defined by earlier call to psphotRoughClassRegion"); + + psRegion *region = psMetadataLookupPtr (&status, regionMD, "REGION"); + psAssert (region, "regions must be defined by earlier call to psphotRoughClassRegion"); + + // pull FWHM_X,Y from the recipe, use to define psfClump.X,Y + psfClump.X = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X"); psAssert (status, "missing PSF.CLUMP.X"); + psfClump.Y = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y"); psAssert (status, "missing PSF.CLUMP.Y"); + psfClump.dX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX"); psAssert (status, "missing PSF.CLUMP.DX"); + psfClump.dY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY"); psAssert (status, "missing PSF.CLUMP.DY"); + + if ((psfClump.X < 0) || (psfClump.Y < 0) || !psfClump.X || !psfClump.Y || isnan(psfClump.X) || isnan(psfClump.Y)) { + psLogMsg ("psphot", 4, "Failed to find a valid PSF clump for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + + if (!psphotSourceClassRegion (region, &psfClump, sources, recipe, psf, maskVal, ApResid, ApSysErr)) { + psLogMsg ("psphot", 4, "Failed to determine source classification for region %f,%f - %f,%f\n", region->x0, region->y0, region->x1, region->y1); + continue; + } + } + + return true; +} + +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal, float ApResid, float ApSysErr) { + + PS_ASSERT_PTR_NON_NULL(sources, false); + PS_ASSERT_PTR_NON_NULL(recipe, false); + + int Nsat = 0; + int Next = 0; + int Npsf = 0; + int Ncr = 0; + int Nmiss = 0; + int Nskip = 0; + + pmSourceMode noMoments = PM_SOURCE_MODE_MOMENTS_FAILURE | PM_SOURCE_MODE_SKYVAR_FAILURE | PM_SOURCE_MODE_SKY_FAILURE | PM_SOURCE_MODE_BELOW_MOMENTS_SN; + pmSourcePhotometryMode photMode = PM_SOURCE_PHOT_WEIGHT; + + for (psS32 i = 0 ; i < sources->n ; i++) { + + pmSource *source = (pmSource *) sources->data[i]; + + // psfClumps are found for image subregions: + // skip sources not in this region + if (source->peak->x < region->x0) continue; + if (source->peak->x >= region->x1) continue; + if (source->peak->y < region->y0) continue; + if (source->peak->y >= region->y1) continue; + + // we are basically classifying by moments; use the default if not found + psAssert (source->moments, "why is this source missing moments?"); + if (source->mode & noMoments) { + Nskip ++; + continue; + } + + psF32 Mxx = source->moments->Mxx; + psF32 Myy = source->moments->Myy; + + // saturated star (determined in PSF fit) + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + Nsat ++; + continue; + } + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + float nSigma = (dMag - ApResid) / hypot(source->errMag, ApSysErr); + + source->extNsigma = nSigma; + + // XXX these sigma cuts should be user-configured parameters + bool isPSF = (fabs(nSigma) < 3.0) && (fabs(Mxx - psfClump->X) < 2.0*psfClump->dX) && (fabs(Myy - psfClump->Y) < 2.0*psfClump->dY); + if (isPSF) { + Npsf ++; + continue; + } + + // XXX get the sign on nSigma right here and above + bool isEXT = (nSigma < 3.0) && (Mxx > psfClump->X) && (Myy > psfClump->Y); + if (isEXT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + Next ++; + continue; + } + + + // XXX possible or likely? need to compare my standard CR test with this + // Mark in both cases? + if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { + Ncr ++; + continue; + } + + Nmiss ++; + } + + psLogMsg("psModules.objects", PS_LOG_INFO, "Rough classifications: %d %d %d %d %d %d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); + + return true; +}