Index: branches/simmosaic_branches/psphot/src/psphotChoosePSF.c =================================================================== --- branches/simmosaic_branches/psphot/src/psphotChoosePSF.c (revision 24860) +++ branches/simmosaic_branches/psphot/src/psphotChoosePSF.c (revision 27839) @@ -1,33 +1,69 @@ # include "psphotInternal.h" -void psphotCountPSFStars (psArray *sources) { - - int nPSF = 0; - - // select the candidate PSF stars (pointers to original sources) - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - if (source->mode & PM_SOURCE_MODE_PSFSTAR) { - nPSF ++; - } - } - - fprintf (stderr, "N PSF: %d\n", nPSF); +// generate a PSF model for inputs without PSF models already loaded +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; } // try PSF models and select best option -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; psTimerStart ("psphot.choose.psf"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // do not generate a PSF if we already were supplied one + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index); + return true; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model"); + return false; + } // bit-masks to test for good/bad pixels psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); + psAssert (maskVal, "missing mask value?"); // bit-mask to mark pixels not used in analysis psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); + psAssert (markVal, "missing mark value?"); // maskVal is used to test for rejected pixels, and must include markVal @@ -73,10 +109,21 @@ // get the fixed PSF fit radius // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS - options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); - assert (status); - - // XXX ROBUST seems to be too agressive given the small numbers. - // options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - options->stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + // options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS"); + // assert (status); + + // values on readout->analysis if we have calculated a Gaussian window function, use that + // for both the PSF fit radius and the aperture radius (scaling SIGMA), otherwise base the value on the recipe value for MOMENTS_GAUSS_SIGMA + float gaussSigma = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + gaussSigma = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE"); + float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE"); + options->fitRadius = (int)(fitScale*gaussSigma); + options->apRadius = (int)(apScale*gaussSigma); + + // XXX use the same radii for standard analysis as for the PSF creation + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius); // dimensions of the field for which the PSF is defined @@ -99,6 +146,4 @@ psArray *stars = psArrayAllocEmpty (sources->n); - - // psphotCountPSFStars (sources); // select the candidate PSF stars (pointers to original sources) @@ -115,13 +160,16 @@ } - // psphotCountPSFStars (sources); - // check that the identified psf stars sufficiently cover the region; if not, extend the // limits somewhat psphotCheckStarDistribution (stars, sources, options); - // psphotCountPSFStars (sources); - psLogMsg ("psphot.pspsf", PS_LOG_DETAIL, "selected candidate %ld PSF objects\n", stars->n); + + if (stars->n < 50) { + // ROBUST is too agressive if we only have a small number of PSF stars + options->stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + } else { + options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + } // get the list pointers for the PSF_MODEL entries @@ -164,5 +212,5 @@ bool status = true; status &= psphotMakeFluxScale (readout->image, recipe, psf); - status &= psphotPSFstats (readout, recipe, psf); + status &= psphotPSFstats (readout, psf); if (!status) { psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF"); @@ -239,5 +287,5 @@ bool status = true; status &= psphotMakeFluxScale (readout->image, recipe, psf); - status &= psphotPSFstats (readout, recipe, psf); + status &= psphotPSFstats (readout, psf); if (!status) { psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF"); @@ -289,14 +337,12 @@ // XXX test dump of psf star data and psf-subtracted image if (psTraceGetLevel("psphot.psfstars") > 5) { - psphotDumpPSFStars (readout, try, options->radius, maskVal, markVal); + psphotDumpPSFStars (readout, try, options->fitRadius, maskVal, markVal); } // save only the best model; - // XXX we are not saving the fitted sources - // XXX do we want to keep them so we may optionally write them out? pmPSF *psf = psMemIncrRefCounter(try->psf); psFree (models); - if (!psphotPSFstats (readout, recipe, psf)) { + if (!psphotPSFstats (readout, psf)) { psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); psFree(options); @@ -305,6 +351,4 @@ } - // psphotCountPSFStars (sources); - char *modelName = pmModelClassGetName (psf->type); psLogMsg ("psphot.pspsf", PS_LOG_INFO, "select psf model: %f sec\n", psTimerMark ("psphot.choose.psf")); @@ -312,9 +356,20 @@ psFree (options); - return (psf); + + // save PSF on readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + psFree (psf); // XXX double-check + + // move into psphotChoosePSF + psphotVisualShowPSFModel (readout, psf); + + return true; } // measure average parameters of the PSF model -bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) { +bool psphotPSFstats (pmReadout *readout, pmPSF *psf) { psEllipseShape shape; @@ -322,5 +377,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(psf, false); @@ -341,5 +395,8 @@ // create modelPSF from this model pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0); - if (!modelPSF) continue; + if (!modelPSF) { + fprintf (stderr, "?"); + continue; + } // get the model full-width at half-max @@ -354,5 +411,8 @@ float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); - if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) continue; + if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { + fprintf (stderr, "!"); + continue; + } psVectorAppend (fwhmMajor, FWHM_MAJOR); psVectorAppend (fwhmMinor, FWHM_MINOR); @@ -366,8 +426,10 @@ } - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ", PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG", PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ", PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ", PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ", PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG", PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ", PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ", PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ); + + float fwhmMaj = stats->sampleMean; // FWHM on major axis if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) { @@ -375,12 +437,20 @@ return false; } - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG", PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ); - - psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); - psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG", PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ); + + float fwhmMin = stats->sampleMean; // FWHM on minor axis + if (readout->parent) { + pmChip *chip = readout->parent->parent; // Parent chip + psAssert(chip, "Cell should be attached to a chip."); + psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip + item->data.F32 = 0.5 * (fwhmMaj + fwhmMin); + } + + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", axes.theta); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true); psFree (fwhmMajor); @@ -391,5 +461,5 @@ // determine approximate PSF shape parameters based on the moments -bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) { +bool psphotMomentsStats (pmReadout *readout, psArray *sources) { // without the PSF model, we can only rely on the source->moments @@ -405,5 +475,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(sources, false); @@ -433,17 +502,22 @@ FWHM_T /= FWHM_N; - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ", PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG", PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ", PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ", PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0); - - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG", PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0); - - psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", FWHM_T); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0); - psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false); + if (readout->parent) { + pmChip *chip = readout->parent->parent; // Parent chip + psAssert(chip, "Cell should be attached to a chip."); + psMetadataItem *item = psMetadataLookup(chip->concepts, "CHIP.SEEING"); // Item with chip + item->data.F32 = 0.5 * (FWHM_X + FWHM_Y); + } + + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ", PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG", PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ", PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ", PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN", PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG", PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ", PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ", PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE", PS_META_REPLACE, "PSF angle", FWHM_T); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false); return true;