Index: /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25389) +++ /branches/eam_branches/20090715/psphot/src/psphotSourceSize.c (revision 25390) @@ -2,16 +2,22 @@ # include -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); -bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); +typedef struct { + psImageMaskType maskVal; + psImageMaskType crMask; + float ApResid; + float ApSysErr; + float nSigmaApResid; + float nSigmaMoments; + float nSigmaCR; + float soft; + int grow; +} psphotSourceSizeOptions; + +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf); +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); +bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); +bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); +bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); // we need to call this function after sources have been fitted to the PSF model and @@ -21,33 +27,40 @@ // deviation from the psf model at the r = FWHM/2 position +// XXX use an internal flag to mark sources which have already been measured bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first) { bool status; + psphotSourceSizeOptions options; psTimerStart ("psphot.size"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + options.maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels + assert (options.maskVal); // bit to mask the cosmic-ray pixels - psImageMaskType crMask = pmConfigMaskGet("CR", config); // Mask value for cosmic rays - - float CR_NSIGMA_LIMIT = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.NSIGMA.LIMIT"); + options.crMask = pmConfigMaskGet("CR", config); // Mask value for cosmic rays + + options.nSigmaCR = psMetadataLookupF32 (&status, recipe, "PSPHOT.CR.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 - if (!status || grow < 0) { + // XXX recipe name is not great + options.nSigmaApResid = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.LIMIT"); + assert (status); + + // XXX recipe name is not great + options.nSigmaMoments = psMetadataLookupF32 (&status, recipe, "PSPHOT.EXT.NSIGMA.MOMENTS"); + assert (status); + + options.grow = psMetadataLookupS32(&status, recipe, "PSPHOT.CR.GROW"); // Growth size for CRs + if (!status || options.grow < 0) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSPHOT.CR.GROW is not positive."); return false; } - float soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter - if (!status || !isfinite(soft) || soft < 0.0) { + options.soft = psMetadataLookupF32(&status, recipe, "PSPHOT.CR.NSIGMA.SOFTEN"); // Softening parameter + if (!status || !isfinite(options.soft) || options.soft < 0.0) { psWarning("PSPHOT.CR.NSIGMA.SOFTEN not set; defaulting to zero."); - soft = 0.0; + options.soft = 0.0; } @@ -57,17 +70,232 @@ // 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); + psphotSourceSizePSF (&options, sources, psf); // 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++) { + psphotSourceClass(readout, sources, recipe, psf, &options); + + psphotSourceSizeCR (readout, sources, &options); + + psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size")); + + psphotVisualPlotSourceSize (recipe, sources); + psphotVisualShowSourceSize (readout, sources); + + return true; +} + +bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + // replace the source flux + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; + + // flag this as a CR + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + // grab the matching footprint + pmFootprint *footprint = peak->footprint; + if (!footprint) { + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + + if (!footprint->spans) { + // if we have no footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (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; + } + } + return true; +} + +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + psImage *mask = source->maskView; + psImage *pixels = source->pixels; + psImage *variance = source->variance; + + // XXX This should be a recipe variable +# define SN_LIMIT 5.0 + + int xo = peak->x - pixels->col0; + int yo = peak->y - pixels->row0; + + // 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; + } + } + 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; + } + } + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // 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; + } + } + // 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; +} + +// model the apmifit distribution for the psf stars: +bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) { + + // 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, options->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; + } + + options->ApResid = stats->robustMedian; + options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05); + fprintf (stderr, "psf-sum: %f +/- %f\n", options->ApResid, options->ApSysErr); + + return true; +} + +// classify sources based on the combination of psf-mag, Mxx, Myy +bool psphotSourceClass (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options) { + + 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, options)) { + 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, psphotSourceSizeOptions *options) { + + 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; // skip source if it was already measured - if (isfinite(source->crNsigma)) { - psTrace("psphot", 7, "Not calculating extNsigma,crNsigma since already measured\n"); + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); continue; } @@ -76,133 +304,63 @@ 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; - - // 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"); - 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); - } - } - - // 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; - } - - psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", - sources->n - first, psTimerMark ("psphot.size")); - - psphotVisualPlotSourceSize (sources); - psphotVisualShowSourceSize (readout, sources); + psTrace("psphot", 7, "Not calculating source size since source is not subtracted\n"); + 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; + + // XXX can we test if psfMag is set and calculate only if needed? + pmSourceMagnitudes (source, psf, photMode, options->maskVal); + + float apMag = -2.5*log10(source->moments->Sum); + float dMag = source->psfMag - apMag; + float nSigma = (dMag - options->ApResid) / hypot(source->errMag, options->ApSysErr); + + source->extNsigma = nSigma; + source->tmpFlags |= PM_SOURCE_TMPF_SIZE_MEASURED; + + // 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); + if (isPSF) { + Npsf ++; + continue; + } + + // Defects may not always match CRs from peak curvature analysis + // Defects may also be marked as SATSTAR -- XXX deactivate this flag? + if ((Mxx < psfClump->X) || (Myy < psfClump->Y)) { + source->mode |= PM_SOURCE_MODE_DEFECT; + Ncr ++; + continue; + } + + // saturated star (determined in PSF fit). These may also be saturated galaxies, or + // just large saturated regions. + if (source->mode & PM_SOURCE_MODE_SATSTAR) { + Nsat ++; + continue; + } + + // XXX allow the Mxx, Myy to be less than psfClump->X,Y (by some nSigma)? + bool isEXT = (nSigma > options->nSigmaApResid) && (Mxx > psfClump->X) && (Myy > psfClump->Y); + if (isEXT) { + source->mode |= PM_SOURCE_MODE_EXT_LIMIT; + Next ++; + continue; + } + + Nmiss ++; + } + + psLogMsg("psModules.objects", PS_LOG_INFO, "Source Size classifications: %d %d %d %d %d %d", Npsf, Next, Nsat, Ncr, Nmiss, Nskip); return true; @@ -211,5 +369,6 @@ // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total // deviation in sigmas. This is measured on the residual image - should we ignore negative -// deviations? +// deviations? NOTE: This function was an early attempt to classify extended objects, and is +// no longer used by psphot. float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, psImageMaskType maskVal, const pmModel *model, float Ro) @@ -282,284 +441,134 @@ } -bool psphotMaskCosmicRay_New (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - // replace the source flux - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED; - - // flag this as a CR - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - // grab the matching footprint - 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 (!footprint->spans) { - // 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; - } - } - return true; -} - -bool psphotMaskCosmicRay_Old (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - psImage *mask = source->maskView; - psImage *pixels = source->pixels; - psImage *variance = source->variance; - - // XXX This should be a recipe variable -# define SN_LIMIT 5.0 - - int xo = peak->x - pixels->col0; - int yo = peak->y - pixels->row0; - - // 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; - } - } - 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; - } - } - - // for each of the neighboring rows, mark the high pixels if they have a marked neighbor - // 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; - } - } - // 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; -} - -// 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; - +bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { + + // classify the sources based on the CR test (place this in a function?) + // XXX use an internal flag to mark sources which have already been measured 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; -} + + // 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 + int xPeak = source->peak->xf - source->pixels->col0 + 0.5; + int yPeak = source->peak->yf - source->pixels->row0 + 0.5; + + // Skip sources which are too close to a boundary. These are mostly caught as DEFECT + 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; + } + + // 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->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) { + bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading for psImageConvolveMask + psImage *newMask = psImageConvolveMask(NULL, readout->mask, options->crMask, options->crMask, -options->grow, options->grow, -options->grow, options->grow); + psImageConvolveSetThreads(oldThreads); + if (!newMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to grow CR mask"); + return false; + } + psFree(readout->mask); + readout->mask = newMask; + } + return true; +}