Index: /branches/eam_branches/psphot.stack.20100120/doc/footprints.txt =================================================================== --- /branches/eam_branches/psphot.stack.20100120/doc/footprints.txt (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/doc/footprints.txt (revision 26748) @@ -1,2 +1,20 @@ + +2010.01.18 + + * we are having some problems with footprints being freed to early in the second pass + + * footprint memory history: + + * entering psphotFindFootprints, on the first pass, we have an + array of detected peak in detections->peaks + + * in pmFootprintsFind, footprints array owns allocated footprints + footprints each own their spans + + * from pmFootprintAssignPeaks: + * the footprints have an array of peaks with valid memory references + * the peaks have just the pointer to the corresponding footprint (not owned) + * all peaks get assigned to a footprint + * at end of pmFootprintAssignPeaks, only a single copy of each peak is/should be available from the collection of footprints pmFootprintCullPeaks is very expensive (15 - 30 msec per object) and Index: /branches/eam_branches/psphot.stack.20100120/doc/stack.txt =================================================================== --- /branches/eam_branches/psphot.stack.20100120/doc/stack.txt (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/doc/stack.txt (revision 26748) @@ -1,2 +1,11 @@ + +20100126: + + * watch out for psphotSetMomentsWindow & MOMENTS_SX_MAX,etc + * watch out for psphotSignificanceImage.c:,psphotEfficiency using the FWHM_MAJ from psphotChoosePSF + * ppSimDetections.c : SIGMA_SMOOTH +ppSmooth/src/ppSmoothReadout.c: psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); +ppSmooth/src/ppSmoothReadout.c: psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); + 20100120 : more stack processing mods: Index: /branches/eam_branches/psphot.stack.20100120/src/Makefile.am =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/Makefile.am (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/Makefile.am (revision 26748) @@ -25,6 +25,6 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot -# bin_PROGRAMS = psphotPetrosianStudy psphotForced psphotMakePSF psphotTest psphotMomentsStudy +bin_PROGRAMS = psphot psphotForced psphotMakePSF +# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -32,23 +32,23 @@ psphot_LDADD = libpsphot.la -# psphotForced_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -# psphotForced_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -# psphotForced_LDADD = libpsphot.la -# -# psphotMakePSF_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -# psphotMakePSF_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -# psphotMakePSF_LDADD = libpsphot.la -# +psphotForced_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotForced_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotForced_LDADD = libpsphot.la + +psphotMakePSF_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotMakePSF_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotMakePSF_LDADD = libpsphot.la + +# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +# psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +# psphotMomentsStudy_LDADD = libpsphot.la + +# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +# psphotPetrosianStudy_LDADD = libpsphot.la + # psphotTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) # psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) # psphotTest_LDADD = libpsphot.la -# -# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -# psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -# psphotMomentsStudy_LDADD = libpsphot.la - -# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) -# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -# psphotPetrosianStudy_LDADD = libpsphot.la # standard psphot for generic photometry @@ -61,25 +61,24 @@ psphotCleanup.c -# # forced photometry of specified positions given a specified psf -# psphotForced_SOURCES = \ -# psphotForced.c \ -# psphotForcedArguments.c \ -# psphotForcedImageLoop.c \ -# psphotForcedReadout.c \ -# psphotParseCamera.c \ -# psphotMosaicChip.c \ -# psphotCleanup.c -# -# # a psphot-variant that simply generates the PSF model -# psphotMakePSF_SOURCES = \ -# psphotMakePSF.c \ -# psphotMakePSFArguments.c \ -# psphotMakePSFImageLoop.c \ -# psphotMakePSFReadout.c \ -# psphotParseCamera.c \ -# psphotMosaicChip.c \ -# psphotCleanup.c -# -# +# forced photometry of specified positions given a specified psf +psphotForced_SOURCES = \ + psphotForced.c \ + psphotForcedArguments.c \ + psphotForcedImageLoop.c \ + psphotForcedReadout.c \ + psphotParseCamera.c \ + psphotMosaicChip.c \ + psphotCleanup.c + +# a psphot-variant that simply generates the PSF model +psphotMakePSF_SOURCES = \ + psphotMakePSF.c \ + psphotMakePSFArguments.c \ + psphotMakePSFImageLoop.c \ + psphotMakePSFReadout.c \ + psphotParseCamera.c \ + psphotMosaicChip.c \ + psphotCleanup.c + # # psphot analysis of the detectability of specified positions # psphotDetect_SOURCES = \ @@ -91,5 +90,5 @@ # psphotMosaicChip.c \ # psphotCleanup.c -# + # psphotTest_SOURCES = \ # psphotTest.c @@ -97,5 +96,5 @@ # psphotMomentsStudy_SOURCES = \ # psphotMomentsStudy.c - +# # psphotPetrosianStudy_SOURCES = \ # psphotPetrosianStudy.c @@ -112,4 +111,7 @@ psphotDefineFiles.c \ psphotReadout.c \ + psphotReadoutFindPSF.c \ + psphotReadoutKnownSources.c \ + psphotReadoutMinimal.c \ psphotModelBackground.c \ psphotSubtractBackground.c \ @@ -170,8 +172,5 @@ psphotEfficiency.c -# XXX need to fix these for the new apis -# psphotReadoutFindPSF.c -# psphotReadoutKnownSources.c -# psphotReadoutMinimal.c +# XXX need to fix this for the new apis # psphotModelTest.c Index: /branches/eam_branches/psphot.stack.20100120/src/psphot.h =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphot.h (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphot.h (revision 26748) @@ -63,6 +63,6 @@ bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); -bool psphotFindDetections (pmConfig *config, const pmFPAview *view); -bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass); +bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass); bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow); @@ -140,6 +140,6 @@ // in psphotChoosePSF.c: -bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); -bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); +bool psphotPSFstats (pmReadout *readout, pmPSF *psf); +bool psphotMomentsStats (pmReadout *readout, psArray *sources); // in psphotGuessModel.c @@ -147,9 +147,11 @@ // in psphotMergeSources.c: -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources); +bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view); psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view); +bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view); +bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view); // generate the detection structure for the supplied array of sources -pmDetections *psphotDetectionsFromSources (pmConfig *config, psArray *sources); +bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources); // generate the detection structure for the supplied array of sources @@ -159,5 +161,5 @@ // Create a background model for a readout, without saving the result as a pmFPAfile on config->files. Otherwise identical to psphotModelBackgroundFileIndex. psImage *psphotModelBackgroundReadoutNoFile(pmReadout *readout, const pmConfig *config); -psImageBinning *psphotBackgroundBinning(const psImage *, const pmConfig *); +psImageBinning *psphotBackgroundBinning(const psImage *image, const pmConfig *config); // in psphotReplaceUnfit.c: @@ -186,5 +188,6 @@ // functions to set the correct source pixels -bool psphotInitRadiusPSF (const psMetadata *recipe, const pmModelType type); +bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type); + bool psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); bool psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR); @@ -224,5 +227,5 @@ bool psphotLoadPSF (pmConfig *config, const pmFPAview *view); -bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotApResid.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotApResid.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotApResid.c (revision 26748) @@ -67,8 +67,8 @@ if (!measureAptrend) { // save nan values since these were not calculated - psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); - psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); return true; } @@ -322,8 +322,8 @@ // save results for later output - psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); - psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid); - psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", psf->ApResid); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss); psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid); @@ -342,8 +342,8 @@ escape: // save nan values since these were not calculated - psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); - psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); - psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT", PS_DATA_F32 | PS_META_REPLACE, "aperture residual", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS", PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN); psFree (xPos); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotBlendFit.c (revision 26748) @@ -79,5 +79,5 @@ psphotInitLimitsPSF (recipe, readout); psphotInitLimitsEXT (recipe); - psphotInitRadiusPSF (recipe, psf->type); + psphotInitRadiusPSF (recipe, readout->analysis, psf->type); // starts the timer, sets up the array of fitSets Index: /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotChoosePSF.c (revision 26748) @@ -107,7 +107,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"); @@ -116,6 +119,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); + 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 @@ -204,5 +207,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"); @@ -279,5 +282,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"); @@ -336,5 +339,5 @@ psFree (models); - if (!psphotPSFstats (readout, recipe, psf)) { + if (!psphotPSFstats (readout, psf)) { psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms"); psFree(options); @@ -363,5 +366,5 @@ // measure average parameters of the PSF model -bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) { +bool psphotPSFstats (pmReadout *readout, pmPSF *psf) { psEllipseShape shape; @@ -369,5 +372,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(psf, false); @@ -419,8 +421,8 @@ } - 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); if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) { @@ -428,12 +430,12 @@ 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); + + 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); @@ -444,5 +446,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 @@ -458,5 +460,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(sources, false); @@ -486,17 +487,15 @@ 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); + 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; Index: /branches/eam_branches/psphot.stack.20100120/src/psphotEfficiency.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotEfficiency.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotEfficiency.c (revision 26748) @@ -35,8 +35,6 @@ psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing - psKernel *newCovar = psImageCovarianceCalculate(kernel, ro->covariance); // Covariance matrix + *covarFactor = psImageCovarianceCalculateFactor(kernel, ro->covariance); psFree(kernel); - *covarFactor = psImageCovarianceFactor(newCovar); // Variance factor - psFree(newCovar); psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance Index: /branches/eam_branches/psphot.stack.20100120/src/psphotFake.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotFake.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotFake.c (revision 26748) @@ -28,8 +28,6 @@ float smoothSigma = 0.5*(FWHM_X + FWHM_Y) / (2.0*sqrtf(2.0*log(2.0))); psKernel *kernel = psImageSmoothKernel(smoothSigma, smoothNsigma); // Kernel used for smoothing - psKernel *newCovar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + float factor = psImageCovarianceCalculateFactor(kernel, readout->covariance); // Covariance matrix psFree(kernel); - float factor = psImageCovarianceFactor(newCovar); // Variance factor - psFree(newCovar); psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for variance @@ -213,10 +211,7 @@ // Putting results on recipe because that appears to be the psphot standard, but it's not a good idea - psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE, - "Efficiency fractions", frac); - psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE, - "Efficiency magnitudes", magOffsets); - psMetadataAddF32(recipe, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, - "Efficiency reference magnitude", magLim); + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE, "Efficiency fractions", frac); + psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE, "Efficiency magnitudes", magOffsets); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, "Efficiency reference magnitude", magLim); psFree(frac); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotFindDetections.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotFindDetections.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotFindDetections.c (revision 26748) @@ -4,5 +4,5 @@ // peaks and new footprints. any old peaks are saved on oldPeaks. the resulting footprint set // contains all footprints (old and new) -bool psphotFindDetections (pmConfig *config, const pmFPAview *view) +bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass) { bool status = true; @@ -17,5 +17,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, recipe, firstPass)) { psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); return false; @@ -26,5 +26,5 @@ // smooth the image, search for peaks, optionally define footprints based on the peaks -bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass) { bool status; @@ -54,14 +54,17 @@ // save detections on the readout->analysis - if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detectinos", detections)) { + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) { psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); return false; } + } else { + psMemIncrRefCounter(detections); // so we can free the detections below + } + if (firstPass) { pass = 1; NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); PS_ASSERT (status, NULL); NMAX = psMetadataLookupS32 (&status, recipe, "PEAKS_NMAX"); PS_ASSERT (status, NULL); } else { - psMemIncrRefCounter(detections); // so we can free the detections below pass = 2; NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT_2"); PS_ASSERT (status, NULL); @@ -96,5 +99,5 @@ // detect the peaks in the significance image detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold); if (!detections->peaks) { // we only get a NULL peaks array due to a programming or config error. @@ -123,3 +126,2 @@ // if we use the footprints, the output peaks list contains both old and new peaks, // otherwise it only contains the new peaks. - Index: /branches/eam_branches/psphot.stack.20100120/src/psphotForcedReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotForcedReadout.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotForcedReadout.c (revision 26748) @@ -25,8 +25,4 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); - // optional break-point for processing char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); @@ -34,42 +30,34 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view, recipe); + psphotSetMaskAndVariance (config, view); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup(config, view); } - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); // generate a background model (median, smoothed image) if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } - // display the backsub and backgnd images - psphotVisualShowBackground (config, view, readout); - - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe - pmPSF *psf = psphotLoadPSF (config, view, recipe); - if (!psf) { - psError(PSPHOT_ERR_CONFIG, false, "unable to load psf model"); - return false; + if (!psphotLoadPSF (config, view)) { + // this only happens if we had a programming error in psphotLoadPSF + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view); } // include externally-supplied sources - psArray *sources = psArrayAllocEmpty(100); - psphotLoadExtSources (config, view, sources); + psphotLoadExtSources (config, view); // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); + psphotFitSourcesLinear (config, view, false); // identify CRs and extended sources @@ -79,5 +67,5 @@ // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); + psphotMagnitudes(config, view); // XXX do I want to do this? @@ -91,7 +79,7 @@ // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (config, view); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, NULL, psf, sources); + return psphotReadoutCleanup(config, view); } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotGuessModels.c (revision 26748) @@ -75,5 +75,5 @@ // setup the PSF fit radius details - psphotInitRadiusPSF (recipe, psf->type); + psphotInitRadiusPSF (recipe, readout->analysis, psf->type); // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) Index: /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotImageQuality.c (revision 26748) @@ -132,6 +132,5 @@ } - psMetadataAddS32(recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, - "Number of stars used for IQ measurements", M2->n); + psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, "Number of stars used for IQ measurements", M2->n); // XXX make this a recipe option @@ -145,8 +144,6 @@ float fwhm_major = stats->sampleMean; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, - "FWHM of Major Axis from moments", stats->sampleMean); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, - "FWHM scatter (Major) from moments", stats->sampleStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, "FWHM of Major Axis from moments", stats->sampleMean); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->sampleStdev); if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) { @@ -155,8 +152,6 @@ } float fwhm_minor = stats->sampleMean; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, - "FWHM of Minor Axis from moments", stats->sampleMean); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, - "FWHM scatter (Minor) from moments", stats->sampleStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->sampleMean); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->sampleStdev); if (!psVectorStats(stats, M2, NULL, NULL, 0)) { @@ -166,12 +161,8 @@ float vM2 = stats->sampleMean; float dM2 = stats->sampleStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, - "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, - "Stdev of M_2", dM2); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, - "Lower Quartile of M_2", stats->sampleLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, - "Upper Quartile of M_2", stats->sampleUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->sampleLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->sampleUQ); if (!psVectorStats(stats, M2c, NULL, NULL, 0)) { @@ -179,12 +170,8 @@ goto FAIL; } - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, - "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, - "Stdev of M_2c", stats->sampleStdev); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, - "Lower Quartile of M_2c", stats->sampleLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, - "Upper Quartile of M_2c", stats->sampleUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->sampleStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->sampleLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->sampleUQ); if (!psVectorStats(stats, M2s, NULL, NULL, 0)) { @@ -192,12 +179,8 @@ goto FAIL; } - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, - "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, - "Stdev of M_2s", stats->sampleStdev); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, - "Lower Quartile of M_2s", stats->sampleLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, - "Upper Quartile of M_2s", stats->sampleUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->sampleStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->sampleLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->sampleUQ); if (!psVectorStats(stats, M3, NULL, NULL, 0)) { @@ -207,12 +190,8 @@ float vM3 = stats->sampleMean; float dM3 = stats->sampleStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, - "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, - "Stdev of M_3", dM3); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, - "Lower Quartile of M_3", stats->sampleLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, - "Upper Quartile of M_3", stats->sampleUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->sampleLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->sampleUQ); if (!psVectorStats(stats, M4, NULL, NULL, 0)) { @@ -222,12 +201,8 @@ float vM4 = stats->sampleMean; float dM4 = stats->sampleStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, - "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, - "Stdev of M_4", dM4); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, - "Lower Quartile of M_4", stats->sampleLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, - "Upper Quartile of M_4", stats->sampleUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->sampleLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->sampleUQ); #else @@ -239,8 +214,6 @@ } float fwhm_major = stats->robustMedian; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, - "FWHM of Major Axis from moments", stats->robustMedian); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, - "FWHM scatter (Major) from moments", stats->robustStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, "FWHM of Major Axis from moments", stats->robustMedian); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->robustStdev); if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) { @@ -249,8 +222,6 @@ } float fwhm_minor = stats->robustMedian; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, - "FWHM of Minor Axis from moments", stats->robustMedian); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, - "FWHM scatter (Minor) from moments", stats->robustStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->robustMedian); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->robustStdev); if (!psVectorStats(stats, M2, NULL, NULL, 0)) { @@ -260,12 +231,8 @@ float vM2 = stats->robustMedian; float dM2 = stats->robustStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, - "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, - "Stdev of M_2", dM2); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, - "Lower Quartile of M_2", stats->robustLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, - "Upper Quartile of M_2", stats->robustUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->robustLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->robustUQ); if (!psVectorStats(stats, M2c, NULL, NULL, 0)) { @@ -273,12 +240,8 @@ goto FAIL; } - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, - "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, - "Stdev of M_2c", stats->robustStdev); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, - "Lower Quartile of M_2c", stats->robustLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, - "Upper Quartile of M_2c", stats->robustUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->robustStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->robustLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->robustUQ); if (!psVectorStats(stats, M2s, NULL, NULL, 0)) { @@ -286,12 +249,8 @@ goto FAIL; } - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, - "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, - "Stdev of M_2s", stats->robustStdev); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, - "Lower Quartile of M_2s", stats->robustLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, - "Upper Quartile of M_2s", stats->robustUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->robustStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->robustLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->robustUQ); if (!psVectorStats(stats, M3, NULL, NULL, 0)) { @@ -301,12 +260,8 @@ float vM3 = stats->robustMedian; float dM3 = stats->robustStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, - "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, - "Stdev of M_3", dM3); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, - "Lower Quartile of M_3", stats->robustLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, - "Upper Quartile of M_3", stats->robustUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->robustLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->robustUQ); if (!psVectorStats(stats, M4, NULL, NULL, 0)) { @@ -316,12 +271,8 @@ float vM4 = stats->robustMedian; float dM4 = stats->robustStdev; - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, - "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, - "Stdev of M_4", dM4); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, - "Lower Quartile of M_4", stats->robustLQ); - psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, - "Upper Quartile of M_4", stats->robustUQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->robustLQ); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->robustUQ); #endif Index: /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotLoadPSF.c (revision 26748) @@ -1,6 +1,15 @@ # include "psphotInternal.h" +// NOTE : pmPSF_IO.c functions must load the psf model onto the chip->analysis metadata because +// the I/O operation likely occurs before the readout exists. This implementation assumes that +// a single psf model is valid for the entire set of readouts (not valid for a time series of readouts) + +// XXX for now (2010.01.27), the supporting programs do not define multiple PSPHOT.PSF.LOAD +// files to go with multiple PSPHOT.INPUT files. as a result, the implementation below is +// currently going to work for the case of a single input file, but will fail if we try with a +// stack of images. + // load an externally supplied psf model -bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { bool status; @@ -25,5 +34,5 @@ } - if (!psphotPSFstats (readout, recipe, psf)) { + if (!psphotPSFstats (readout, psf)) { psAbort("cannot measure PSF shape terms"); } @@ -44,9 +53,5 @@ bool status = false; - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - psAssert (recipe, "missing recipe?"); - - // XXX PSPHOT.PSF.LOAD vs PSPHOT.INPUT?? + // XXX PSPHOT.PSF.LOAD vs PSPHOT.INPUT -- see note at top int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -56,5 +61,5 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotLoadPSFReadout (config, view, "PSPHOT.PSF.LOAD", i, recipe)) { + if (!psphotLoadPSFReadout (config, view, "PSPHOT.PSF.LOAD", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i); return false; Index: /branches/eam_branches/psphot.stack.20100120/src/psphotMakePSFReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotMakePSFReadout.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotMakePSFReadout.c (revision 26748) @@ -24,8 +24,4 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); - // optional break-point for processing char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); @@ -33,96 +29,48 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view, recipe); + psphotSetMaskAndVariance (config, view); if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup(config, view); } - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); // generate a background model (median, smoothed image) if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } - // display the backsub and backgnd images - psphotVisualShowBackground (config, view, readout); + psphotLoadExtSources (config, view); - pmDetections *detections = NULL; - - // If sources have been supplied, then these should be used to measure the PSF - // include externally-supplied sources - - // XXX sources loaded from a text file have no valid SN values, but psphotChoosePSF - // selected the top PSF_MAX_NSTARS to generate the PSF, excluding an arbitrary subset. - psArray *sources = psArrayAllocEmpty(100); - psphotLoadExtSources (config, view, sources); - if (sources->n) { - // the user wants to make the psf from these stars; define them as psf stars: - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->mode |= PM_SOURCE_MODE_PSFSTAR; - } - // force psphotChoosePSF to use all loaded sources - psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF_MAX_NSTARS", PS_META_REPLACE, "fit radius", sources->n); - - // measure stats of externally specified sources - if (!psphotSourceStatsUpdate (sources, config, readout)) { - psError(PSPHOT_ERR_CONFIG, false, "failure to measure stats of existing sources"); - return false; - } - - } else { - // find the detections (by peak and/or footprint) in the image. - detections = psphotFindDetections (NULL, readout, recipe); - if (!detections) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources); - } - - // construct sources and measure basic stats - psFree (sources); - sources = psphotSourceStats (config, readout, detections, true); - if (!sources) return false; - - // find blended neighbors of very saturated stars - // XXX merge this with Basic Deblend? - psphotDeblendSatstars (readout, sources, recipe); - - // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (sources, recipe)) { - psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources); - } - - // classify sources based on moments, brightness - if (!psphotRoughClass (readout, sources, recipe, false)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources); - } + // If sources have been supplied, then these should be used to measure the PSF include + // externally-supplied sources; if not, we need to generate a set of possible PSF sources. + // This function updates the SN entries for the loaded sources or generates a set of + // detections from the image, if no external ones have been supplied. Sources loaded from + // a text file have no valid SN values, but psphotChoosePSF needs to select the top + // PSF_MAX_NSTARS to generate the PSF. + if (!psphotCheckExtSources (config, view)) { + psLogMsg ("psphot", 3, "failure to select possible PSF sources (external or internal)"); + return psphotReadoutCleanup (config, view); } - // use bright stellar objects to measure PSF - // XXX if we do not have enough stars to generate the PSF, build one - // from the SEEING guess and model class - pmPSF *psf = psphotChoosePSF (readout, sources, recipe); - if (psf == NULL) { + // Use bright stellar objects to measure PSF. If we do not have enough stars to generate + // the PSF, build one from the SEEING guess and model class + if (!psphotChoosePSF (config, view)) { psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } - psphotVisualShowPSFModel (readout, psf); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); -} // measure aperture photometry corrections - // XXX isn't this part of the PSF?? - // if (!psphotApResid (config, readout, sources, psf)) { - // psLogMsg ("psphot", 3, "failed on psphotApResid"); - // return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - // } +# if 0 + if (!psphotApResid (config, view)) { + psLogMsg ("psphot", 3, "failed on psphotApResid"); + return psphotReadoutCleanup (config, view); + } +# endif + + return psphotReadoutCleanup (config, view); +} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotMaskReadout.c (revision 26748) @@ -36,7 +36,5 @@ psAssert (readout, "missing readout?"); - // ** Interpret the mask values: - // XXX drop the write to recipe and move config into psphotRoughClass? - // XXX alternatively, define a function to set the psphot recipe masks + // save maskSat and maskBad on the psphot recipe (mostly for psphotRoughClass) psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.SAT", PS_META_REPLACE, "user-defined mask", maskSat); @@ -44,5 +42,5 @@ psImageMaskType maskBad = pmConfigMaskGet("LOW", config); // Mask value for low pixels if (!maskBad) { - // XXX: for backward compatability look up old name + // for backward compatability look up old name maskBad = pmConfigMaskGet("BAD", config); } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotMergeSources.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotMergeSources.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotMergeSources.c (revision 26748) @@ -41,4 +41,12 @@ psAssert (newSources, "missing sources?"); + // XXX TEST: + if (detections->allSources) { + psphotMaskCosmicRayFootprintCheck(detections->allSources); + } + if (detections->newSources) { + psphotMaskCosmicRayFootprintCheck(detections->newSources); + } + if (!detections->allSources) { detections->allSources = psArrayAllocEmpty(newSources->n); @@ -51,36 +59,57 @@ } - // XXX free detections->newSources? - return true; -} - -// merge the externally supplied sources with the existing sources. mark them as having -// mode PM_SOURCE_MODE_EXTERNAL -// XXX this function needs to be updated to work with the new context of pshot inputs -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources) { - - psArray *extSourcesCMF = NULL; + psFree (detections->newSources); + detections->newSources = NULL; + + return true; +} + +// Merge the externally supplied sources with the existing sources. Mark them as having mode +// PM_SOURCE_MODE_EXTERNAL. + +// XXX This function needs to be updated to loop over set of input files. At the moment, we +// only expect a single entry for PSPHOT.INPUT.CMF and PSPHOT.SOURCES.TEXT, so we can only +// associate input sources with a single entry for PSPHOT.INPUT +bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view) { + + bool status; + pmDetections *extCMF = NULL; psArray *extSourcesTXT = NULL; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + // XXX allSources of newSources? + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + // load data from input CMF file: - { - pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF"); - if (!readoutCMF) goto loadTXT; - - extSourcesCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.SOURCES"); - if (extSourcesCMF) { - for (int i = 0; i < extSourcesCMF->n; i++) { - pmSource *source = extSourcesCMF->data[i]; + { + pmReadout *readoutCMF = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT.CMF"); + if (!readoutCMF) goto loadTXT; + + extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS"); + if (extCMF) { + for (int i = 0; i < extCMF->allSources->n; i++) { + pmSource *source = extCMF->allSources->data[i]; source->mode |= PM_SOURCE_MODE_EXTERNAL; - // the supplied peak flux needs to be re-normalized - source->peak->flux = 1.0; - source->peak->value = 1.0; + // the supplied peak flux needs to be re-normalized + source->peak->flux = 1.0; + source->peak->value = 1.0; // drop the loaded source modelPSF psFree (source->modelPSF); source->modelPSF = NULL; + + psArrayAdd (detections->allSources, 100, source); } - // XXX psphotMergeSources (sources, extSourcesCMF); } } @@ -89,19 +118,20 @@ // load data from input TXT file: - { - pmChip *chipTXT = pmFPAfileThisChip (config->files, view, "PSPHOT.INPUT"); - if (!chipTXT) goto finish; - - extSourcesTXT = psMetadataLookupPtr (NULL, chipTXT->analysis, "PSPHOT.SOURCES.TEXT"); - if (extSourcesTXT) { - for (int i = 0; i < extSourcesTXT->n; i++) { - pmSource *source = extSourcesTXT->data[i]; - source->mode |= PM_SOURCE_MODE_EXTERNAL; + { + pmChip *chipTXT = pmFPAfileThisChip (config->files, view, "PSPHOT.INPUT"); + if (!chipTXT) goto finish; + + extSourcesTXT = psMetadataLookupPtr (NULL, chipTXT->analysis, "PSPHOT.SOURCES.TEXT"); + if (extSourcesTXT) { + for (int i = 0; i < extSourcesTXT->n; i++) { + pmSource *source = extSourcesTXT->data[i]; + source->mode |= PM_SOURCE_MODE_EXTERNAL; // drop the loaded source modelPSF psFree (source->modelPSF); source->modelPSF = NULL; + + psArrayAdd (detections->allSources, 100, source); } - // XXX psphotMergeSources (sources, extSourcesTXT); } } @@ -109,14 +139,14 @@ finish: - if (!extSourcesTXT && !extSourcesTXT) { + if (!extCMF && !extSourcesTXT) { psLogMsg ("psphot", 3, "no external sources for this readout"); return true; } - int nCMF = extSourcesCMF ? extSourcesCMF->n : 0; + int nCMF = extCMF ? extCMF->allSources->n : 0; int nTXT = extSourcesTXT ? extSourcesTXT->n : 0; - psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d text) merged to yield %ld total sources", - nCMF + nTXT, nCMF, nTXT, sources->n); + psLogMsg ("psphot", 3, "%d external sources (%d cmf, %d text) merged to yield %ld total sources", + nCMF + nTXT, nCMF, nTXT, sources->n); return true; } @@ -125,4 +155,6 @@ // XXX this function needs to be updated to work with the new context of pshot inputs psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view) { + + bool status; // find the currently selected readout @@ -133,5 +165,10 @@ } - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + psLogMsg ("psphot", 3, "no psf sources for this readout"); + } + + psArray *sources = detections->allSources; if (!sources) { psLogMsg ("psphot", 3, "no psf sources for this readout"); @@ -141,6 +178,42 @@ } +bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + psError(PSPHOT_ERR_CONFIG, false, "missing detections"); + return false; + } + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + // peak flux is wrong : set based on previous image + // use the peak measured in the moments analysis: + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + source->peak->flux = source->moments->Peak; + } + + return true; +} + // generate the detection structure for the supplied array of sources -pmDetections *psphotDetectionsFromSources (pmConfig *config, psArray *sources) { +// XXX this currently assumes there is a single input file +bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); pmDetections *detections = pmDetectionsAlloc(); @@ -155,5 +228,5 @@ float snMin = psMetadataLookupF32(NULL, recipe, "MOMENTS_SN_MIN"); if (!isfinite(snMin)) { - return NULL; + return false; } @@ -166,9 +239,10 @@ } - // use the existing peak information, otherwise generate a new peak - if (source->peak) { - psArrayAdd (detections->peaks, 100, source->peak); - continue; - } + // use the existing peak information, otherwise generate a new peak + if (source->peak) { + source->peak->assigned = false; // So the moments will be measured + psArrayAdd (detections->peaks, 100, source->peak); + continue; + } float flux = powf(10.0, -0.4 * source->psfMag); @@ -181,9 +255,9 @@ peak->flux = flux; // this are being set wrong, but does it matter? - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - peak->SN = 1.0 / source->errMag; - } else { - peak->SN = 0.0; - } + if (isfinite (source->errMag) && (source->errMag > 0.0)) { + peak->SN = 1.0 / source->errMag; + } else { + peak->SN = 0.0; + } psArrayAdd (detections->peaks, 100, peak); @@ -193,5 +267,11 @@ psLogMsg ("psphot", 3, "%ld PSF sources loaded", detections->peaks->n); - return detections; + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detectinos", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree (detections); + return true; } @@ -222,9 +302,9 @@ peak->flux = flux; // this are being set wrong, but does it matter? - if (isfinite (source->errMag) && (source->errMag > 0.0)) { - peak->SN = 1.0 / source->errMag; - } else { - peak->SN = 0.0; - } + if (isfinite (source->errMag) && (source->errMag > 0.0)) { + peak->SN = 1.0 / source->errMag; + } else { + peak->SN = 0.0; + } source->peak = peak; @@ -235,2 +315,67 @@ return true; } + +bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view) { + + bool status; + + psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + // XXX allSources of newSources? + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (sources->n) { + // the user wants to make the psf from these stars; define them as psf stars: + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + source->mode |= PM_SOURCE_MODE_PSFSTAR; + } + // force psphotChoosePSF to use all loaded sources + psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF_MAX_NSTARS", PS_META_REPLACE, "max number of sources for PSF model", sources->n); + + // measure stats of externally specified sources + if (!psphotSourceStatsUpdate (sources, config, readout)) { + psError(PSPHOT_ERR_CONFIG, false, "failure to measure stats of existing sources"); + return false; + } + } else { + + // find the detections (by peak and/or footprint) in the image. + if (!psphotFindDetections (config, view, true)) { + psError(PSPHOT_ERR_CONFIG, false, "unable to find detections in this image"); + return psphotReadoutCleanup (config, view); + } + + // construct sources and measure basic stats + psphotSourceStats (config, view, true); + + // find blended neighbors of very saturated stars + psphotDeblendSatstars (config, view); + + // mark blended peaks PS_SOURCE_BLEND + if (!psphotBasicDeblend (config, view)) { + psLogMsg ("psphot", 3, "failed on deblend analysis"); + return psphotReadoutCleanup (config, view); + } + + // classify sources based on moments, brightness + if (!psphotRoughClass (config, view)) { + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view); + } + } + + return true; +} Index: /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotModelBackground.c (revision 26748) @@ -32,4 +32,9 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model +// +// NOTE that the 'analysis' metedata pass in here is used to store the binning information. +// This may be the analysis for this readout, but it may be the analysis for the pmFPAfile +// corresponding to the model. Other information about the background model is saved on the +// readout->analysis static bool psphotModelBackgroundReadout(psImage *model, // Model image psImage *modelStdev, // Model stdev image @@ -140,6 +145,5 @@ // we save the binning structure for use in psphotMagnitudes - psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", - PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); + psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning); psF32 **modelData = model->data.F32; @@ -296,6 +300,6 @@ psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background")); - psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); - psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev); psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev); @@ -306,16 +310,10 @@ PS_STAT_MAX); psImageStats (statsBck, model, NULL, 0); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN", - PS_META_REPLACE, "sky model mean", statsBck->sampleMean); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG", - PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX", - PS_META_REPLACE, "sky model maximum value", statsBck->max); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN", - PS_META_REPLACE, "sky model minimum value", statsBck->min); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX", - PS_META_REPLACE, "sky model size (x)", model->numCols); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY", - PS_META_REPLACE, "sky model size (y)", model->numRows); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean", statsBck->sampleMean); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev", statsBck->sampleStdev); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max); + psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)", model->numCols); + psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)", model->numRows); psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotOutput.c (revision 26748) @@ -197,8 +197,7 @@ // these values are saved in an output header stub - they are added to either the PHU header -// (CMP) or the MEF table header (CMF) XXX these are currently carried by the recipe -- this -// will not work in a Stack context XXX they should be place in the readout->analysis of the -// given image or directly on the output header -psMetadata *psphotDefineHeader (psMetadata *recipe) { +// (CMP) or the MEF table header (CMF). before the header is created, each readout has these +// values stored on readout->analysis +psMetadata *psphotDefineHeader (psMetadata *analysis) { bool status = true; @@ -207,75 +206,75 @@ // write necessary information to output header - psMetadataItemSupplement (&status, header, recipe, "ZERO_PT"); - psMetadataItemSupplement (&status, header, recipe, "PHOTCODE"); - - psMetadataItemSupplement (&status, header, recipe, "APMIFIT"); - psMetadataItemSupplement (&status, header, recipe, "DAPMIFIT"); - psMetadataItemSupplement (&status, header, recipe, "NAPMIFIT"); + psMetadataItemSupplement (&status, header, analysis, "ZERO_PT"); + psMetadataItemSupplement (&status, header, analysis, "PHOTCODE"); + + psMetadataItemSupplement (&status, header, analysis, "APMIFIT"); + psMetadataItemSupplement (&status, header, analysis, "DAPMIFIT"); + psMetadataItemSupplement (&status, header, analysis, "NAPMIFIT"); // PSF model parameters (shape values for image center) - psMetadataItemSupplement (&status, header, recipe, "NPSFSTAR"); - psMetadataItemSupplement (&status, header, recipe, "APLOSS"); - - psMetadataItemSupplement (&status, header, recipe, "FWHM_MAJ"); - psMetadataItemSupplement (&status, header, recipe, "FW_MJ_SG"); - psMetadataItemSupplement (&status, header, recipe, "FW_MJ_LQ"); - psMetadataItemSupplement (&status, header, recipe, "FW_MJ_UQ"); - - psMetadataItemSupplement (&status, header, recipe, "FWHM_MIN"); - psMetadataItemSupplement (&status, header, recipe, "FW_MN_SG"); - psMetadataItemSupplement (&status, header, recipe, "FW_MN_LQ"); - psMetadataItemSupplement (&status, header, recipe, "FW_MN_UQ"); - - psMetadataItemSupplement (&status, header, recipe, "ANGLE"); + psMetadataItemSupplement (&status, header, analysis, "NPSFSTAR"); + psMetadataItemSupplement (&status, header, analysis, "APLOSS"); + + psMetadataItemSupplement (&status, header, analysis, "FWHM_MAJ"); + psMetadataItemSupplement (&status, header, analysis, "FW_MJ_SG"); + psMetadataItemSupplement (&status, header, analysis, "FW_MJ_LQ"); + psMetadataItemSupplement (&status, header, analysis, "FW_MJ_UQ"); + + psMetadataItemSupplement (&status, header, analysis, "FWHM_MIN"); + psMetadataItemSupplement (&status, header, analysis, "FW_MN_SG"); + psMetadataItemSupplement (&status, header, analysis, "FW_MN_LQ"); + psMetadataItemSupplement (&status, header, analysis, "FW_MN_UQ"); + + psMetadataItemSupplement (&status, header, analysis, "ANGLE"); // Image Quality measurements - psMetadataItemSupplement (&status, header, recipe, "IQ_NSTAR"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_FW1"); - psMetadataItemSupplement (&status, header, recipe, "IQ_FW1_E"); - psMetadataItemSupplement (&status, header, recipe, "IQ_FW2"); - psMetadataItemSupplement (&status, header, recipe, "IQ_FW2_E"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_M2"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2_ER"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2_LQ"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2_UQ"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_M2C"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_E"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_L"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_U"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_M2S"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_E"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_L"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_U"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_M3"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M3_ER"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M3_LQ"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M3_UQ"); - - psMetadataItemSupplement (&status, header, recipe, "IQ_M4"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M4_ER"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M4_LQ"); - psMetadataItemSupplement (&status, header, recipe, "IQ_M4_UQ"); + psMetadataItemSupplement (&status, header, analysis, "IQ_NSTAR"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_FW1"); + psMetadataItemSupplement (&status, header, analysis, "IQ_FW1_E"); + psMetadataItemSupplement (&status, header, analysis, "IQ_FW2"); + psMetadataItemSupplement (&status, header, analysis, "IQ_FW2_E"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_M2"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2_ER"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2_LQ"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2_UQ"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_M2C"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_E"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_L"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_U"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_M2S"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_E"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_L"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_U"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_M3"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M3_ER"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M3_LQ"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M3_UQ"); + + psMetadataItemSupplement (&status, header, analysis, "IQ_M4"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M4_ER"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M4_LQ"); + psMetadataItemSupplement (&status, header, analysis, "IQ_M4_UQ"); // XXX these need to be defined from elsewhere psMetadataAdd (header, PS_LIST_TAIL, "FSATUR", PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG", 0.0); psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT", PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG", 0.0); - psMetadataItemSupplement (&status, header, recipe, "NSTARS"); - - psMetadataItemSupplement (&status, header, recipe, "NDET_EXT"); - psMetadataItemSupplement (&status, header, recipe, "NDET_CR"); + psMetadataItemSupplement (&status, header, analysis, "NSTARS"); + + psMetadataItemSupplement (&status, header, analysis, "NDET_EXT"); + psMetadataItemSupplement (&status, header, analysis, "NDET_CR"); // sky background model statistics - psMetadataItemSupplement (&status, header, recipe, "MSKY_MN"); - psMetadataItemSupplement (&status, header, recipe, "MSKY_SIG"); - psMetadataItemSupplement (&status, header, recipe, "MSKY_MIN"); - psMetadataItemSupplement (&status, header, recipe, "MSKY_MAX"); - psMetadataItemSupplement (&status, header, recipe, "MSKY_NX"); - psMetadataItemSupplement (&status, header, recipe, "MSKY_NY"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_MN"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_SIG"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_MIN"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_MAX"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_NX"); + psMetadataItemSupplement (&status, header, analysis, "MSKY_NY"); psMetadataAddF32 (header, PS_LIST_TAIL, "DT_PHOT", PS_META_REPLACE, "elapsed psphot time", psTimerMark ("psphotReadout")); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotRadiusChecks.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotRadiusChecks.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotRadiusChecks.c (revision 26748) @@ -4,18 +4,24 @@ static float PSF_FIT_NSIGMA; static float PSF_FIT_PADDING; +static float PSF_APERTURE = 0; // radius to use in PSF aperture mags static float PSF_FIT_RADIUS = 0; // radius to use in fitting (ignored if <= 0, // and a per-object radius is calculated) -static float PSF_APERTURE = 0; // radius to use in PSF aperture mags - - -bool psphotInitRadiusPSF(const psMetadata *recipe, const pmModelType type) { +bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) { bool status = true; - PSF_FIT_NSIGMA = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA"); + PSF_FIT_NSIGMA = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA"); PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING"); - PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); - PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + + PSF_FIT_RADIUS = psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS"); + if (!status) { + PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS"); + } + + PSF_APERTURE = psMetadataLookupF32(&status, analysis, "PSF_APERTURE"); + if (!status) { + PSF_APERTURE = psMetadataLookupF32(&status, recipe, "PSF_APERTURE"); + } return true; Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadout.c (revision 26748) @@ -50,19 +50,13 @@ } - // run a single-model test if desired (exits from here if test is run) - // XXX fix this in the new multi-input context (or drop this and only keep a stand-alone program) - // psphotModelTest (config, view, recipe); - - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe - // XXX the recipe is probably the wrong place to store those results... - // XXX this is using the wrong file... fix in the new context - // if (!psphotLoadPSF (config, view, recipe)) { - // // this only happens if we had a programming error in psphotLoadPSF - // psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); - // return psphotReadoutCleanup (config, view); - // } + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on + // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT + if (!psphotLoadPSF (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view); + } // find the detections (by peak and/or footprint) in the image. - if (!psphotFindDetections (config, view)) { // pass 1 + if (!psphotFindDetections (config, view, true)) { // pass 1 // this only happens if we had an error in psphotFindDetections psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); @@ -119,5 +113,5 @@ // include externally-supplied sources // XXX fix this in the new multi-input context - // psphotLoadExtSources (config, view, sources); // pass 1 + psphotLoadExtSources (config, view); // pass 1 // construct an initial model for each object, set the radius to fitRadius, set circular @@ -159,5 +153,5 @@ // find fainter sources // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections - psphotFindDetections (config, view); // pass 2 (detections->peaks, detections->footprints) + psphotFindDetections (config, view, false); // pass 2 (detections->peaks, detections->footprints) // remove noise for subtracted objects (ie, return to normal noise level) Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutCleanup.c (revision 26748) @@ -58,5 +58,5 @@ // use the psf-model to measure FWHM stats if (psf) { - if (!psphotPSFstats (readout, recipe, psf)) { + if (!psphotPSFstats (readout, psf)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters"); return false; @@ -65,5 +65,5 @@ // otherwise, use the source moments to measure FWHM stats if (!psf && sources) { - if (!psphotMomentsStats (readout, recipe, sources)) { + if (!psphotMomentsStats (readout, sources)) { psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters"); return false; @@ -82,8 +82,6 @@ } - // create an output header with stats results XXX this needs to be reworked : - // psphotImageQuality and others need to write these values to a header. that needs to be - // stored on the readout->analysis - psMetadata *header = psphotDefineHeader (recipe); + // create an output header with stats results currently saved on readout->analysis + psMetadata *header = psphotDefineHeader (readout->analysis); // write NSTARS to the image header Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutFindPSF.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutFindPSF.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutFindPSF.c (revision 26748) @@ -7,11 +7,4 @@ psTimerStart ("psphotReadout"); - // select the current recipe - psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); - if (!recipe) { - psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); - return false; - } - // set the photcode for the PSPHOT.INPUT if (!psphotAddPhotcode(config, view)) { @@ -20,13 +13,6 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); - // Generate the mask and variance images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view, recipe); - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); + psphotSetMaskAndVariance (config, view); // Note that in this implementation, we do NOT model the background and we do not @@ -34,39 +20,37 @@ // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF) - pmDetections *detections = psphotDetectionsFromSources (config, inSources); - if (!detections || !detections->peaks) { + // XXX we assume a single set of input sources is supplied + if (!psphotDetectionsFromSources (config, view, inSources)) { psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars"); - return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL); + return psphotReadoutCleanup(config, view); } // construct sources and measure basic stats (moments, local sky) - psArray *sources = psphotSourceStats(config, readout, detections, true); - if (!sources) return false; + if (!psphotSourceStats(config, view, true)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return false; + } - // peak flux is wrong : set based on previous image - // use the peak measured in the moments analysis: - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->peak->flux = source->moments->Peak; + // peak flux is wrong : use the peak measured in the moments analysis: + if (!psphotRepairLoadedSources(config, view)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources"); + return false; } // classify sources based on moments, brightness (psf is not known) - if (!psphotRoughClass (readout, sources, recipe, false)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources); + if (!psphotRoughClass (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class"); + return psphotReadoutCleanup (config, view); } - if (!psphotImageQuality (recipe, sources)) { - psLogMsg("psphot", 3, "failed to measure image quality"); - return psphotReadoutCleanup(config, readout, recipe, detections, NULL, sources); + if (!psphotImageQuality (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); + return psphotReadoutCleanup(config, view); } - pmPSF *psf = psphotChoosePSF(readout, sources, recipe); - if (!psf) { + if (!psphotChoosePSF(config, view)) { psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); - psFree(sources); - return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL); + return psphotReadoutCleanup(config, view); } - psphotVisualShowPSFModel(readout, psf); # if 0 @@ -75,18 +59,17 @@ // fits from that analysis, or run the linear PSF fit for all objects currently in hand // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { + if (!psphotApResid (config, view)) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } # endif // drop the references to the image pixels held by each source - psphotSourceFreePixels(sources); - psFree(sources); + psphotSourceFreePixels(config, view); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, NULL); + return psphotReadoutCleanup(config, view); } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutKnownSources.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutKnownSources.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutKnownSources.c (revision 26748) @@ -1,16 +1,9 @@ # include "psphotInternal.h" -// in this psphotReadout-variant, we are only measuring the photometry for known sources, -// using a PSF generated from this observation those sources +// in this psphotReadout-variant, we are only measuring the photometry for known sources, using +// a PSF generated for this observation from those sources bool psphotReadoutKnownSources(pmConfig *config, const pmFPAview *view, psArray *inSources) { psTimerStart ("psphotReadout"); - - // select the current recipe - psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); - if (!recipe) { - psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); - return false; - } // set the photcode for this image @@ -20,13 +13,6 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); - // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view, recipe); - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); + psphotSetMaskAndVariance (config, view); // Note that in this implementation, we do NOT model the background and we do not @@ -34,58 +20,51 @@ // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF) - pmDetections *detections = psphotDetectionsFromSources (config, inSources); - if (!detections || !detections->peaks) { + if (!psphotDetectionsFromSources (config, view, inSources)) { psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars"); - return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL); + return psphotReadoutCleanup(config, view); } // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, true); - if (!sources) return false; + if (!psphotSourceStats (config, view, true)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return false; + } - // peak flux is wrong : set based on previous image - // use the peak measured in the moments analysis: - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - source->peak->flux = source->moments->Peak; + // peak flux is wrong : use the peak measured in the moments analysis: + if (!psphotRepairLoadedSources(config, view)) { + psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources"); + return false; } // classify sources based on moments, brightness (psf is not known) - if (!psphotRoughClass (readout, sources, recipe, false)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources); + if (!psphotRoughClass (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class"); + return psphotReadoutCleanup(config, view); } - pmPSF *psf = psphotChoosePSF (readout, sources, recipe); - if (!psf) { - psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + if (!psphotChoosePSF (config, view)) { + psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model"); + return psphotReadoutCleanup(config, view); } - psphotVisualShowPSFModel (readout, psf); // construct an initial model for each object - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); // linear PSF fit to source peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - - // We have to place these visualizations here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); - psphotVisualShowSatStars (recipe, psf, sources); + psphotFitSourcesLinear (config, view, false); // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { + if (!psphotApResid (config, view)) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); + psphotMagnitudes(config, view); // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (config, view); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutMinimal.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutMinimal.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutMinimal.c (revision 26748) @@ -17,11 +17,4 @@ pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); - // select the current recipe - psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); - if (!recipe) { - psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); - return false; - } - // set the photcode for this image if (!psphotAddPhotcode(config, view)) { @@ -30,78 +23,61 @@ } - // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - PS_ASSERT_PTR_NON_NULL (readout, false); + // Generate the mask and weight images, including the user-defined analysis region of interest + psphotSetMaskAndVariance (config, view); - // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view, recipe); - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); - - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe - pmPSF *psf = psphotLoadPSF (config, view, recipe); - if (!psf) { + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved on readout->analysis + if (!psphotLoadPSF (config, view)) { psError (PSPHOT_ERR_CONFIG, false, "missing psf model"); - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } - // find the detections (by peak and/or footprint) in the image. - pmDetections *detections = pmDetectionsAlloc(); // New detections; allocated to ensure pass=2 - detections = psphotFindDetections(detections, readout, recipe); - if (!detections) { + // find the detections (by peak and/or footprint) in the image. (final pass) + if (!psphotFindDetections(config, view, false)) { psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); - } - if (!detections->peaks->n) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); + return psphotReadoutCleanup (config, view); } - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, false); - if (!sources) return false; + // construct sources and measure basic stats (saved on detections->newSources) + if (!psphotSourceStats (config, view, false)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); + } // find blended neighbors of very saturated stars - // XXX merge this with Basic Deblend? - psphotDeblendSatstars (readout, sources, recipe); + psphotDeblendSatstars (config, view); // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (sources, recipe)) { + if (!psphotBasicDeblend (config, view)) { psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // classify sources based on moments, brightness (use supplied psf shape parameters) - if (!psphotRoughClass (readout, sources, recipe, true)) { + if (!psphotRoughClass (config, view)) { psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // construct an initial model for each object - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); // linear PSF fit to source peaks - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - - // We have to place these visualizations here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); - psphotVisualShowSatStars (recipe, psf, sources); + psphotFitSourcesLinear (config, view, false); // XXX eventually, add the extended source fits here # if (0) // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, 0); + psphotSourceSize (config, view); - psphotExtendedSourceAnalysis (readout, sources, recipe); + psphotExtendedSourceAnalysis (config, view); - psphotExtendedSourceFits (readout, sources, recipe); + psphotExtendedSourceFits (config, view); # endif // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); + psphotMagnitudes(config, view); - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + // XXX ensure this is measured if the analysis succeeds (even if quality is low) + if (!psphotEfficiency(config, view)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -109,7 +85,7 @@ // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (config, view); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } Index: anches/eam_branches/psphot.stack.20100120/src/psphotReadoutStack.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotReadoutStack.c (revision 26747) +++ (revision ) @@ -1,207 +1,0 @@ -# include "psphotInternal.h" - -bool psphotReadoutStack (pmConfig *config, const pmFPAview *view) { - - psTimerStart ("psphotReadout"); - - pmModelClassSetLimits(PM_MODEL_LIMITS_LAX); // allow models to have ugly fits (eg, central cusp) - - // select the current recipe - psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE); - if (!recipe) { - psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE); - return false; - } - - // set the photcode for this image - if (!psphotAddPhotcode (recipe, config, view, "PSPHOT.INPUT")) { - psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); - return false; - } - - // generate a background model (median, smoothed image) - if (!psphotStackSetMaskAndVariance (config, view, recipe)) { - return psphotReadoutCleanupStack (config, recipe, NULL, NULL, NULL); - } - if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL); - } - - // optional break-point for processing - char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); - PS_ASSERT_PTR_NON_NULL (breakPt, false); - - // generate a background model (median, smoothed image) - if (!psphotStackModelBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); - } - if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); - } - if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); - } - - // find the detections (by peak and/or footprint) in the image. - pmDetections *detections = psphotFindDetections (NULL, readout, recipe); - if (!detections) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); - } - - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, true); - if (!sources) return false; - if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - // find blended neighbors of very saturated stars - // XXX merge this with Basic Deblend? - psphotDeblendSatstars (readout, sources, recipe); - - // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (sources, recipe)) { - psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // classify sources based on moments, brightness - if (!psphotRoughClass (readout, sources, recipe, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - if (!havePSF && !psphotImageQuality (recipe, sources)) { - psLogMsg("psphot", 3, "failed to measure image quality"); - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); - } - - // if we were not supplied a PSF, choose one here - if (psf == NULL) { - // use bright stellar objects to measure PSF - // XXX if we do not have enough stars to generate the PSF, build one - // from the SEEING guess and model class - psf = psphotChoosePSF (readout, sources, recipe); - if (psf == NULL) { - psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - havePSF = true; - } - if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - psphotVisualShowPSFModel (readout, psf); - - // include externally-supplied sources - psphotLoadExtSources (config, view, sources); - - // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, sources, psf); - - // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - - // We have to place these visualizations here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); - - // identify CRs and extended sources - psphotSourceSize (config, readout, sources, recipe, psf, 0); - if (!strcasecmp (breakPt, "ENSEMBLE")) { - goto finish; - } - psphotVisualShowSatStars (recipe, psf, sources); - - // non-linear PSF and EXT fit to brighter sources - // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) - psphotBlendFit (config, readout, sources, psf); - - // replace all sources - psphotReplaceAllSources (sources, recipe); - - // linear fit to include all sources (subtract again) - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); - - // if we only do one pass, skip to extended source analysis - if (!strcasecmp (breakPt, "PASS1")) { - goto pass1finish; - } - // NOTE: possibly re-measure background model here with objects subtracted - - // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe); - - // find fainter sources (pass 2) - detections = psphotFindDetections (detections, readout, recipe); - - // remove noise for subtracted objects (ie, return to normal noise level) - psphotSubNoise (readout, sources, recipe); - - // define new sources based on only the new peaks - psArray *newSources = psphotSourceStats (config, readout, detections, false); - - // set source type - if (!psphotRoughClass (readout, newSources, recipe, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // create full input models, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, newSources, psf); - - // replace all sources so fit below applies to all at once - psphotReplaceAllSources (sources, recipe); - - // merge the newly selected sources into the existing list - psphotMergeSources (sources, newSources); - psFree (newSources); - - // linear fit to all sources - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); - -pass1finish: - - // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, psf, 0); - - psphotExtendedSourceAnalysis (readout, sources, recipe); - - psphotExtendedSourceFits (readout, sources, recipe); - -finish: - - // plot positive sources - // psphotSourcePlots (readout, sources, recipe); - - // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { - psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); - - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { - psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); - psErrorClear(); - } - - // replace failed sources? - // psphotReplaceUnfitSources (sources); - - // replace background in residual image - psphotSkyReplace (config, view); - - // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); - - // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); -} - Index: /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotRoughClass.c (revision 26748) @@ -118,5 +118,26 @@ // determine the PSF parameters from the source moment values // XXX why not save the psfClump as a PTR? - psfClump = pmSourcePSFClump (region, sources, recipe); + + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); + float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); + + float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE"); + if (!status) { + PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE"); + psAssert (status, "missing PSF_CLUMP_GRID_SCALE"); + } + float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX"); + if (!status) { + MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX"); + psAssert (status, "missing MOMENTS_SX_MAX"); + } + float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX"); + if (!status) { + MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX"); + psAssert (status, "missing MOMENTS_SY_MAX"); + } + + psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); + psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y", PS_META_REPLACE, "psf clump center", psfClump.Y); @@ -146,6 +167,10 @@ psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); + // get basic parameters, or set defaults + float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); + float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); psAssert (status, "missing PSF_CLUMP_NSIGMA"); + // group into STAR, COSMIC, EXTENDED, SATURATED, etc. - if (!pmSourceRoughClass (region, sources, recipe, psfClump, maskSat)) { + if (!pmSourceRoughClass (region, sources, PSF_SN_LIM, PSF_CLUMP_NSIGMA, psfClump, maskSat)) { psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass"); return false; Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSetThreads.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSetThreads.c (revision 26748) @@ -25,5 +25,5 @@ psFree(task); - task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 5); + task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10); task->function = &psphotSourceStats_Threaded; psThreadTaskAdd(task); Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSignificanceImage.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSignificanceImage.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSignificanceImage.c (revision 26748) @@ -22,7 +22,8 @@ } + // if we have already determined the PSF model, then we have a better idea how to smooth this image bool statusMajor, statusMinor; - float fwhmMajor = psMetadataLookupF32(&statusMajor, recipe, "FWHM_MAJ"); - float fwhmMinor = psMetadataLookupF32(&statusMinor, recipe, "FWHM_MIN"); + float fwhmMajor = psMetadataLookupF32(&statusMajor, readout->analysis, "FWHM_MAJ"); + float fwhmMinor = psMetadataLookupF32(&statusMinor, readout->analysis, "FWHM_MIN"); if (statusMajor && statusMinor) { // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?) @@ -43,6 +44,5 @@ } // record the actual smoothing sigma - psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", - SIGMA_SMTH); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH); // smooth the image, applying the mask as we go @@ -59,6 +59,5 @@ // effective per-pixel variance is maintained. psImage *smooth_wt = psImageCopy(NULL, readout->variance, PS_TYPE_F32); - psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, - NSIGMA_SMTH, minGauss); + psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss); psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth")); @@ -77,14 +76,11 @@ // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing psKernel *kernel = psImageSmoothKernel(SIGMA_SMTH, NSIGMA_SMTH); // Kernel used for smoothing - psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + float factor = 1.0 / psImageCovarianceCalculateFactor(kernel, readout->covariance); psFree(kernel); - float factor = 1.0 / psImageCovarianceFactor(covar); - psFree(covar); // record the effective area and significance scaling factor float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH); - psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); - psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, - "Signicance scale factor", factor); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea); + psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor); // XXX multithread this if needed Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSourceSize.c (revision 26748) @@ -179,6 +179,6 @@ float dMag = source->psfMag - apMag; - psVectorAppend (Ap, 100, dMag); - psVectorAppend (ApErr, 100, source->errMag); + psVectorAppend (Ap, dMag); + psVectorAppend (ApErr, source->errMag); } Index: /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c (revision 26747) +++ /branches/eam_branches/psphot.stack.20100120/src/psphotSourceStats.c (revision 26748) @@ -42,4 +42,12 @@ psAssert (!detections->newSources, "new sources already defined?"); + // XXX TEST: + if (detections->allSources) { + psphotMaskCosmicRayFootprintCheck(detections->allSources); + } + if (detections->newSources) { + psphotMaskCosmicRayFootprintCheck(detections->newSources); + } + // determine the number of allowed threads int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads @@ -57,4 +65,15 @@ char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); psAssert (status, "missing BREAK_POINT?"); + + float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); psAssert (status, "missing SKY_INNER_RADIUS"); + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN"); + + // bit-masks to test for good/bad pixels + psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); + psAssert (maskVal, "missing MASK.PSPHOT"); + + // bit-mask to mark pixels not used in analysis + psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); + psAssert (markVal, "missing MARK.PSPHOT"); psArray *peaks = detections->peaks; @@ -112,4 +131,14 @@ } + // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis + float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA"); + if (!status) { + SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); + } + float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS"); + if (!status) { + RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); + } + // threaded measurement of the source magnitudes int Nfail = 0; @@ -133,5 +162,13 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, recipe); + + PS_ARRAY_ADD_SCALAR(job->args, INNER, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, MIN_SN, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, RADIUS, PS_TYPE_F32); + PS_ARRAY_ADD_SCALAR(job->args, SIGMA, PS_TYPE_F32); + + PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK); + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail @@ -152,4 +189,137 @@ psFree(sources); return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } else { + psScalar *scalar = NULL; + scalar = job->args->data[7]; + Nmoments += scalar->data.S32; + scalar = job->args->data[8]; + Nfail += scalar->data.S32; + scalar = job->args->data[9]; + Nfaint += scalar->data.S32; + } + psFree(job); + } + } + + psFree (cellGroups); + + psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats")); + + psphotVisualShowMoments (sources); + + // save the new sources on the detection structure: + detections->newSources = sources; + + + if (detections->allSources) { + psphotMaskCosmicRayFootprintCheck(detections->allSources); + } + if (detections->newSources) { + psphotMaskCosmicRayFootprintCheck(detections->newSources); + } + + return true; +} + +// this function is currently only called by psphotCheckExtSources +bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) { + + bool status = false; + + psTimerStart ("psphot.stats"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + + // determine properties (sky, moments) of initial sources + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + if (!status) return false; + + OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius + + char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); + if (!status) return false; + + for (int i = 0; i < sources->n; i++) { + + pmSource *source = sources->data[i]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + // allocate space for moments + if (!source->moments) { + source->moments = pmMomentsAlloc(); + } + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + source->peak->assigned = true; + } + + if (!strcasecmp (breakPt, "PEAKS")) { + psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats")); + psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n"); + psphotVisualShowMoments (sources); + return sources; + } + + // XXX how else could we get the window size in? + if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); + return NULL; + } + + // threaded measurement of the source magnitudes + int Nfail = 0; + int Nmoments = 0; + int Nfaint = 0; + + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); + + // XXX: this must match the above + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, recipe); + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail + PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint + + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); + psFree (job); + return NULL; + } + psFree(job); + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false)) { + psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); + return NULL; } @@ -178,127 +348,4 @@ psphotVisualShowMoments (sources); - // save the new sources on the detection structure: - detections->newSources = sources; - return true; -} - -// XXX fix the return status of this function... -bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) { - - bool status = false; - - psTimerStart ("psphot.stats"); - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); - - // determine the number of allowed threads - int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads - if (!status) { - nThreads = 0; - } - - // determine properties (sky, moments) of initial sources - float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); - if (!status) return false; - - OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius - - char *breakPt = psMetadataLookupStr (&status, recipe, "BREAK_POINT"); - if (!status) return false; - - for (int i = 0; i < sources->n; i++) { - - pmSource *source = sources->data[i]; - if (!source->peak) continue; // XXX how can we have a peak-less source? - - // allocate space for moments - if (!source->moments) { - source->moments = pmMomentsAlloc(); - } - - // allocate image, weight, mask arrays for each peak (square of radius OUTER) - pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); - source->peak->assigned = true; - } - - if (!strcasecmp (breakPt, "PEAKS")) { - psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats")); - psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n"); - psphotVisualShowMoments (sources); - return sources; - } - - // XXX how else could we get the window size in? - if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); - return NULL; - } - - // threaded measurement of the source magnitudes - int Nfail = 0; - int Nmoments = 0; - int Nfaint = 0; - - // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) - int Cx = 1, Cy = 1; - psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); - - psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); - - for (int i = 0; i < cellGroups->n; i++) { - - psArray *cells = cellGroups->data[i]; - - for (int j = 0; j < cells->n; j++) { - - // allocate a job -- if threads are not defined, this just runs the job - psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS"); - - psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, recipe); - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail - PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint - - if (!psThreadJobAddPending(job)) { - psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); - psFree (job); - return NULL; - } - psFree(job); - } - - // wait for the threads to finish and manage results - if (!psThreadPoolWait (false)) { - psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); - return NULL; - } - - // we have only supplied one type of job, so we can assume the types here - psThreadJob *job = NULL; - while ((job = psThreadJobGetDone()) != NULL) { - if (job->args->n < 1) { - fprintf (stderr, "error with job\n"); - } else { - psScalar *scalar = NULL; - scalar = job->args->data[2]; - Nmoments += scalar->data.S32; - scalar = job->args->data[3]; - Nfail += scalar->data.S32; - scalar = job->args->data[4]; - Nfaint += scalar->data.S32; - } - psFree(job); - } - } - - psFree (cellGroups); - - psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats")); - - psphotVisualShowMoments (sources); - return (sources); } @@ -311,22 +358,12 @@ psArray *sources = job->args->data[0]; - psMetadata *recipe = job->args->data[1]; - - float INNER = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); - if (!status) return false; - float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); - if (!status) return false; - float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS"); - if (!status) return false; - float SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA"); - if (!status) return false; - - // bit-masks to test for good/bad pixels - psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); - - // bit-mask to mark pixels not used in analysis - psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); + + float INNER = PS_SCALAR_VALUE(job->args->data[1],F32); + float MIN_SN = PS_SCALAR_VALUE(job->args->data[2],F32); + float RADIUS = PS_SCALAR_VALUE(job->args->data[3],F32); + float SIGMA = PS_SCALAR_VALUE(job->args->data[4],F32); + + psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA); + psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA); // maskVal is used to test for rejected pixels, and must include markVal @@ -392,11 +429,11 @@ // change the value of a scalar on the array (wrap this and put it in psArray.h) - scalar = job->args->data[2]; + scalar = job->args->data[7]; scalar->data.S32 = Nmoments; - scalar = job->args->data[3]; + scalar = job->args->data[8]; scalar->data.S32 = Nfail; - scalar = job->args->data[4]; + scalar = job->args->data[9]; scalar->data.S32 = Nfaint; @@ -410,6 +447,7 @@ bool status; - float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); - if (!status) return false; + float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN"); + float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM"); + psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX"); // XXX move this to a config file? @@ -441,10 +479,10 @@ // choose a grid scale that is a fixed fraction of the psf sigma^2 - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i])); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i])); + float PSF_CLUMP_GRID_SCALE = 0.1*PS_SQR(sigma[i]); + float MOMENTS_SX_MAX = 2.0*PS_SQR(sigma[i]); + float MOMENTS_SY_MAX = 2.0*PS_SQR(sigma[i]); // determine the PSF parameters from the source moment values - pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe); + pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX); psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); @@ -485,9 +523,9 @@ // choose a grid scale that is a fixed fraction of the psf sigma^2 - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma)); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); - psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma); - psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); + psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma)); + psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); + psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma)); + psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma); + psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma); psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);