Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 37766) +++ /trunk/psphot/src/psphot.h (revision 37767) @@ -172,4 +172,5 @@ // in psphotChoosePSF.c: bool psphotPSFstats (pmReadout *readout, pmPSF *psf); +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf); bool psphotMomentsStats (pmReadout *readout, psArray *sources); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 37766) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 37767) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf); // generate a PSF model for inputs without PSF models already loaded @@ -388,5 +389,7 @@ psFree (models); - if (!psphotPSFstats (readout, psf)) { + // if (!psphotPSFstats (readout, psf)) { + + if (!psphotPSFstatsSources (readout, sources, psf)) { psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); psFree(options); @@ -466,13 +469,13 @@ if (modelPSF->params->n > 7) { - psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); + psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); } if (modelPSF->params->n > 8) { - psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); + psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); } psFree (modelPSF); - // float fwhmtest = pmPSFtoFWHM(psf, xc, yc); - // fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); } } @@ -522,19 +525,172 @@ if (psfExtra1->n) { - if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + if (psfExtra2->n) { + if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + 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); + + // psLogMsg ("psphot", PS_LOG_DETAIL, "PSF angle: %f, nstars: %d", axes.theta, psf->nPSFstars); + + char *psfModelName = pmModelClassGetName(psf->type); + psMetadataAddStr(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "PSF Model Name", psfModelName); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSF_OK", PS_META_REPLACE, "Valid PSF Model?", true); + + int nParams = pmModelClassParameterCount(psf->type); + psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "PSF_NPAR", PS_META_REPLACE, "Number of PSF parameters", nParams); + + psFree (fwhmMajor); + psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); + psFree (stats); + return true; + +escape: + psFree (fwhmMajor); + psFree (fwhmMinor); + psFree (psfExtra1); + psFree (psfExtra2); + psFree (stats); + return false; +} + +// measure average parameters of the PSF model +bool psphotPSFstatsSources (pmReadout *readout, psArray *sources, pmPSF *psf) { + + psEllipseShape shape; + psEllipseAxes axes; + + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(psf, false); + + psImage *image = readout->image; + PS_ASSERT_PTR_NON_NULL(image, false); + + // XXX a minor hack: measure the values for a grid of points across the image, determine mean, UQ, LQ, etc: + psVector *fwhmMajor = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *fwhmMinor = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra1 = psVectorAllocEmpty (100, PS_DATA_F32); + psVector *psfExtra2 = psVectorAllocEmpty (100, PS_DATA_F32); + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + + if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue; + + float xc = source->peak->xf; + float yc = source->peak->yf; + + // create modelPSF from this model + pmModel *modelPSF = pmModelFromPSFforXY (psf, xc, yc, 1.0); + if (!modelPSF) { + fprintf (stderr, "?"); + continue; + } + + // get the model full-width at half-max + float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5); + + // XXX make sure this is consistent with the re-definition of PM_PAR_SXX + shape.sx = modelPSF->params->data.F32[PM_PAR_SXX]; + shape.sy = modelPSF->params->data.F32[PM_PAR_SYY]; + shape.sxy = modelPSF->params->data.F32[PM_PAR_SXY]; + axes = psEllipseShapeToAxes (shape, 20.0); + + float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); + if (!isfinite(FWHM_MAJOR) || !isfinite(FWHM_MINOR)) { + fprintf (stderr, "!"); + psFree (modelPSF); + continue; + } + psVectorAppend (fwhmMajor, FWHM_MAJOR); + psVectorAppend (fwhmMinor, FWHM_MINOR); + + if (modelPSF->params->n > 7) { + psVectorAppend (psfExtra1, modelPSF->params->data.F32[7]); + } + if (modelPSF->params->n > 8) { + psVectorAppend (psfExtra2, modelPSF->params->data.F32[8]); + } + psFree (modelPSF); + + float fwhmtest = pmPSFtoFWHM(psf, xc, yc); + fprintf (stderr, "fwhm: %f, %f : %f\n", FWHM_MAJOR, FWHM_MINOR, fwhmtest); + } + + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE); + if (!psVectorStats (stats, fwhmMajor, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MAJOR"); goto escape; - } - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); - psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); + } + + 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)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for FWHM MINOR"); + goto escape; + } + 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) { + + // we now have 2 possible measurements of the seeing : the PSF model based version and + // the Moments based value we need to define a definitive "CHIP.SEEING" value. This + // value is used by the PSF-matching programs (ppSub and ppStack). The moments-based + // value is probably more representative of the value needed by those tools (it is also + // more stable?) + bool status = false; + float fwhmMajorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + float fwhmMinorMoments = psMetadataLookupF32(&status, readout->analysis, "IQ_FW1"); + + 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 * (fwhmMajorMoments + fwhmMinorMoments); + + psLogMsg ("psphot", PS_LOG_DETAIL, "fwhm (psf): %f,%f (moments): %f,%f", fwhmMaj, fwhmMin, fwhmMajorMoments, fwhmMinorMoments); + } + + if (psfExtra1->n) { + if (!psVectorStats (stats, psfExtra1, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 1"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT1", PS_META_REPLACE, "PSF extra param 1", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 1: %f +/- %f", stats->sampleMean, stats->sampleStdev); } if (psfExtra2->n) { - if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); - goto escape; - } - psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); - psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); + if (!psVectorStats (stats, psfExtra2, NULL, NULL, 0)) { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats for PSF EXTRA 2"); + goto escape; + } + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PSF_EXT2", PS_META_REPLACE, "PSF extra param 2", stats->sampleMean); + psLogMsg ("psphot", PS_LOG_DETAIL, "PSF extra parameter 2: %f +/- %f", stats->sampleMean, stats->sampleStdev); } Index: /trunk/psphot/src/psphotReadoutCleanup.c =================================================================== --- /trunk/psphot/src/psphotReadoutCleanup.c (revision 37766) +++ /trunk/psphot/src/psphotReadoutCleanup.c (revision 37767) @@ -57,5 +57,5 @@ // use the psf-model to measure FWHM stats if (psf) { - if (!psphotPSFstats (readout, psf)) { + if (!psphotPSFstatsSources (readout, sources, psf)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters"); return false;