Index: branches/tap_branches/psphot/src/psphotChoosePSF.c =================================================================== --- branches/tap_branches/psphot/src/psphotChoosePSF.c (revision 25900) +++ branches/tap_branches/psphot/src/psphotChoosePSF.c (revision 27838) @@ -1,18 +1,69 @@ # include "psphotInternal.h" +// 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 @@ -61,7 +112,10 @@ // assert (status); - // We have calculated a Gaussian window function, use that for both the PSF fit radius and - // the aperture radius (scaling SIGMA) - float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA"); + // 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"); @@ -70,10 +124,6 @@ // XXX use the same radii for standard analysis as for the PSF creation - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius); - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius); - - // 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); + 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 @@ -116,4 +166,11 @@ 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 psArray *modelNames = NULL; @@ -155,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"); @@ -230,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"); @@ -287,5 +344,5 @@ psFree (models); - if (!psphotPSFstats (readout, recipe, psf)) { + if (!psphotPSFstats (readout, psf)) { psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); psFree(options); @@ -299,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; @@ -309,5 +377,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(psf, false); @@ -329,7 +396,7 @@ pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0); if (!modelPSF) { - fprintf (stderr, "?"); - continue; - } + fprintf (stderr, "?"); + continue; + } // get the model full-width at half-max @@ -345,7 +412,7 @@ float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { - fprintf (stderr, "!"); - continue; - } + fprintf (stderr, "!"); + continue; + } psVectorAppend (fwhmMajor, FWHM_MAJOR); psVectorAppend (fwhmMinor, FWHM_MINOR); @@ -359,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)) { @@ -368,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); @@ -384,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 @@ -398,5 +475,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(sources, false); @@ -426,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;