Index: /branches/eam_branches/20091201/psphot/doc/notes.20100131.txt =================================================================== --- /branches/eam_branches/20091201/psphot/doc/notes.20100131.txt (revision 26788) +++ /branches/eam_branches/20091201/psphot/doc/notes.20100131.txt (revision 26788) @@ -0,0 +1,22 @@ + +20100131 + + We have several stand-alone psphot versions, and a few different + psphotReadout versions called by other programs. Here is a complete + list: + + Stand-alone psphot versions | associated psphotReadout version + + o psphot | psphotReadout + o psphotForced | psphotForcedReadout + o psphotMakePSF | psphotMakePSFReadout + x psphotDetect | psphotDetectReadout + x psphotStack | psphotStackReadout + + External programs | associated psphotReadout version + + o ppImage | psphotReadout + o pswarp | psphotReadoutFindPSF + o ppSub | psphotReadoutFindPSF, psphotReadoutMinimal + o ppStack | psphotReadoutKnownSources + Index: /branches/eam_branches/20091201/psphot/doc/stack.txt =================================================================== --- /branches/eam_branches/20091201/psphot/doc/stack.txt (revision 26787) +++ /branches/eam_branches/20091201/psphot/doc/stack.txt (revision 26788) @@ -1,2 +1,42 @@ + +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: + + there are a number of data collections generated during psphotReadout: + + * detections (currently a structure containing multiple arrays) + * sources (a psArray) + * psf (the psf model structure) + + * new sources vs old sources: + there is a sequence on the second pass in which we need to distinguish the sources + from the first pass from those on the second pass: + + - add noise (old sources only) + - find detections + - subtract noise (old sources only) + - generate sources (new detections only) + - source classifications (new sources only) + - guess models (new sources only) + - replace sources (old sources only) + - merge sources (new + old -> sources) + + currently we are distiguishing the old vs new based on different arrays. + can we use the processing flags to distinguish the these cases and carry around + only a single source list? + + * detections->peaks holds the most recently detected set of peaks + detections->oldPeaks holds the previous collection + + detections->footprints holds the full set of merged footprints, + including assigned peaks from the old and new set. 20100107 : updates for stack processing Index: /branches/eam_branches/20091201/psphot/src/Makefile.am =================================================================== --- /branches/eam_branches/20091201/psphot/src/Makefile.am (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/Makefile.am (revision 26788) @@ -25,7 +25,6 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotTest psphotMomentsStudy -# bin_PROGRAMS = psphotPetrosianStudy -# bin_PROGRAMS = psphot +bin_PROGRAMS = psphot psphotForced psphotMakePSF +# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -41,15 +40,15 @@ psphotMakePSF_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 +# 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 # standard psphot for generic photometry @@ -82,21 +81,20 @@ psphotCleanup.c - -# psphot analysis of the detectability of specified positions -psphotDetect_SOURCES = \ - psphotDetect.c \ - psphotDetectArguments.c \ - psphotDetectParseCamera.c \ - psphotDetectImageLoop.c \ - psphotDetectReadout.c \ - psphotMosaicChip.c \ - psphotCleanup.c - -psphotTest_SOURCES = \ - psphotTest.c - -psphotMomentsStudy_SOURCES = \ - psphotMomentsStudy.c - +# # psphot analysis of the detectability of specified positions +# psphotDetect_SOURCES = \ +# psphotDetect.c \ +# psphotDetectArguments.c \ +# psphotDetectParseCamera.c \ +# psphotDetectImageLoop.c \ +# psphotDetectReadout.c \ +# psphotMosaicChip.c \ +# psphotCleanup.c + +# psphotTest_SOURCES = \ +# psphotTest.c +# +# psphotMomentsStudy_SOURCES = \ +# psphotMomentsStudy.c +# # psphotPetrosianStudy_SOURCES = \ # psphotPetrosianStudy.c @@ -146,5 +144,4 @@ psphotKernelFromPSF.c \ psphotPSFConvModel.c \ - psphotModelTest.c \ psphotFitSet.c \ psphotSourceFreePixels.c \ @@ -175,4 +172,7 @@ psphotEfficiency.c +# XXX need to fix this for the new apis +# psphotModelTest.c + # re-instate these # psphotIsophotal.c \ Index: /branches/eam_branches/20091201/psphot/src/psphot.h =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphot.h (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphot.h (revision 26788) @@ -12,4 +12,12 @@ #define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs + +// pmPCMData : PSF Convolved Model data storage structure +typedef struct { + psImage *model; + psArray *dmodels; + psImage *modelConv; + psArray *dmodelsConv; +} pmPCMData; // top-level psphot functions @@ -29,5 +37,7 @@ bool psphotReadoutMinimal(pmConfig *config, const pmFPAview *view); -bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources); +bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view); +bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + bool psphotDefineFiles (pmConfig *config, pmFPAfile *input); void psphotFilesActivate(pmConfig *config, bool state); @@ -38,66 +48,122 @@ // XXX test functions psArray *psphotFakeSources (void); -bool psphotMaskCosmicRayFootprintCheck (psArray *sources); +bool psphotMaskCosmicRayFootprintCheck (psArray *sources); // psphotReadout functions +bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view); +bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); + +bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view); +bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + bool psphotModelBackground (pmConfig *config, const pmFPAview *view); bool psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index); -// 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); - -// worker function for above background model functions -bool psphotModelBackgroundReadout(psImage *model, // Model image - psImage *modelStdev, // Model stdev image - psMetadata *analysis, // Analysis metadata for outputs - pmReadout *readout, // Readout for which to generate a background model - psImageBinning *binning, // Binning parameters - const pmConfig *config // Configuration - ); - -psImageBinning *psphotBackgroundBinning(const psImage *, const pmConfig *); - - - bool psphotSubtractBackground (pmConfig *config, const pmFPAview *view); -bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); - - -pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe); - -bool psphotRoughClass (pmReadout *readout, psArray *sources, psMetadata *recipe, const bool findPsfClump); -bool psphotBasicDeblend (psArray *sources, psMetadata *recipe); -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf); -bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final); -bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); - -bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe); -bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe); - -bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf); +bool psphotSubtractBackgroundReadout (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); +bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow); + +bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view); +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + +bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view); +bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + +bool psphotRoughClass (pmConfig *config, const pmFPAview *view); +bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); +bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *analysis, psMetadata *recipe, const bool havePSF); + +bool psphotImageQuality (pmConfig *config, const pmFPAview *view); +bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view); +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotGuessModels (pmConfig *config, const pmFPAview *view); +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + +bool psphotMergeSources (pmConfig *config, const pmFPAview *view); +bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final); +bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final); + +bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize); +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize); + +bool psphotBlendFit (pmConfig *config, const pmFPAview *view); +bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); bool psphotBlendFit_Threaded (psThreadJob *job); -bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout); -psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow); -bool psphotSourceStats_Threaded (psThreadJob *job); - -bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf); -bool psphotGuessModel_Threaded (psThreadJob *job); - -bool psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf); +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view); +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotAddNoise (pmConfig *config, const pmFPAview *view); +bool psphotSubNoise (pmConfig *config, const pmFPAview *view); +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add); +bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add); + +bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view); +bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view); +bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotApResid (pmConfig *config, const pmFPAview *view); +bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); + +bool psphotMagnitudes (pmConfig *config, const pmFPAview *view); +bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); bool psphotMagnitudes_Threaded (psThreadJob *job); + +bool psphotEfficiency (pmConfig *config, const pmFPAview *view); +bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe); bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources); bool psphotPSFWeights_Threaded (psThreadJob *job); -bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf); - -bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); -bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources); +bool psphotSkyReplace (pmConfig *config, const pmFPAview *view); +bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + +bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view); +bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index); + +// in psphotSourceStats.c: +bool psphotSourceStats_Threaded (psThreadJob *job); +bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout); +bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources); + +// in psphotChoosePSF.c: +bool psphotPSFstats (pmReadout *readout, pmPSF *psf); +bool psphotMomentsStats (pmReadout *readout, psArray *sources); + +// in psphotGuessModel.c +bool psphotGuessModel_Threaded (psThreadJob *job); + +// in psphotMergeSources.c: +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 +bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources); + +// generate the detection structure for the supplied array of sources +bool psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf); + +// in psphotModelBackground.c: +// 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 *image, const pmConfig *config); + +// in psphotReplaceUnfit.c: +bool psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe); +bool psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal); // thread-related: @@ -113,8 +179,5 @@ psErrorCode psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints); -// generate the detection structure for the supplied array of sources -pmDetections *psphotDetectionsFromSources (pmConfig *config, psArray *sources); - -// used by ApResid +// in psphotApResid.c: pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag); bool psphotApResidMags_Threaded (psThreadJob *job); @@ -123,10 +186,8 @@ void psphotModelClassInit (void); bool psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal); -bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe); -bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe); -void psphotSourceFreePixels (psArray *sources); // 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); @@ -134,8 +195,4 @@ bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal); float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal); - -// output functions -bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index); -bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view); bool psphotDumpMoments (psMetadata *recipe, psArray *sources); @@ -169,23 +226,12 @@ bool psphotFitSummary (void); -bool psphotMergeSources (psArray *oldSources, psArray *newSources); -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources); -psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view); - -pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe); +bool psphotLoadPSF (pmConfig *config, const pmFPAview *view); +bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index); + bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources); -bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe); -bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe); -bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add); bool psphotRadialPlot (int *kapa, const char *filename, pmSource *source); bool psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe); bool psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY, bool normalize); -bool psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal); -bool psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal); -bool psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal); -bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe); -bool psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first); - bool psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal); @@ -195,35 +241,33 @@ // functions related to extended source analysis -bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); -bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax); -float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); -bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax); -bool psphotEllipticalProfile (pmSource *source); -bool psphotEllipticalContour (pmSource *source); +bool psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); +bool psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax); +float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax); +bool psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax); +bool psphotEllipticalProfile (pmSource *source); +bool psphotEllipticalContour (pmSource *source); // psphotVisual functions -bool psphotVisualShowImage (pmReadout *readout); -bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout); -bool psphotVisualShowSignificance (psImage *image, float min, float max); -bool psphotVisualShowPeaks (pmDetections *detections); -bool psphotVisualShowFootprints (pmDetections *detections); -bool psphotVisualShowMoments (psArray *sources); -bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources); -bool psphotVisualShowRoughClass (psArray *sources); -bool psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources); -bool psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources); -bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf); -bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal); -bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources); -bool psphotVisualShowFlags (psArray *sources); -bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); -bool psphotVisualShowResidualImage (pmReadout *readout); -bool psphotVisualPlotApResid (psArray *sources, float mean, float error); -bool psphotVisualPlotChisq (psArray *sources); -bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources); -bool psphotVisualShowPetrosians (psArray *sources); - -// bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe); -// bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise); +bool psphotVisualShowImage (pmReadout *readout); +bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout); +bool psphotVisualShowSignificance (psImage *image, float min, float max); +bool psphotVisualShowPeaks (pmDetections *detections); +bool psphotVisualShowFootprints (pmDetections *detections); +bool psphotVisualShowMoments (psArray *sources); +bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources); +bool psphotVisualShowRoughClass (psArray *sources); +bool psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources); +bool psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources); +bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf); +bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal); +bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources); +bool psphotVisualShowFlags (psArray *sources); +bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources); +bool psphotVisualShowResidualImage (pmReadout *readout); +bool psphotVisualPlotApResid (psArray *sources, float mean, float error); +bool psphotVisualPlotChisq (psArray *sources); +bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources); +bool psphotVisualShowPetrosians (psArray *sources); +bool psphotVisualEraseOverlays (int channel, char *overlay); bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal); @@ -231,5 +275,5 @@ bool psphotPetrosianStats (pmSource *source); -// XXX old versions, currently disabled +// currently disabled: // bool psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); // bool psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal); @@ -246,16 +290,5 @@ // float petFlux, float radiusForFlux); -bool psphotImageQuality (psMetadata *recipe, psArray *sources); - // structures & functions to support psf-convolved model fitting - -// pmPCMData : PSF Convolved Model data storage structure -typedef struct { - psImage *model; - psArray *dmodels; - psImage *modelConv; - psArray *dmodelsConv; -} pmPCMData; - // psf-convolved model fitting @@ -287,5 +320,4 @@ int psphotKapaOpen (void); -bool psphotVisualEraseOverlays (int channel, char *overlay); bool psphotKapaClose (void); bool psphotImageBackgroundCellHistogram (psVector *values, float mean, float sigma, int ix, int iy); @@ -301,5 +333,4 @@ int psphotSupplementStars (psArray *stars, psArray *sources, psImageBinning *binning, int ix, int iy); - pmConfig *psphotForcedArguments(int argc, char **argv); bool psphotForcedImageLoop (pmConfig *config); Index: /branches/eam_branches/20091201/psphot/src/psphotAddNoise.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotAddNoise.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotAddNoise.c (revision 26788) @@ -1,13 +1,35 @@ # include "psphotInternal.h" -bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { - return psphotAddOrSubNoise (readout, sources, recipe, true); +bool psphotAddNoise (pmConfig *config, const pmFPAview *view) { + return psphotAddOrSubNoise (config, view, true); } -bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) { - return psphotAddOrSubNoise (readout, sources, recipe, false); +bool psphotSubNoise (pmConfig *config, const pmFPAview *view) { + return psphotAddOrSubNoise (config, view, false); } -bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) { bool status = false; @@ -16,7 +38,18 @@ psEllipseAxes axes; - PS_ASSERT (readout, false); - PS_ASSERT (readout->parent, false); - PS_ASSERT (readout->parent->concepts, false); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + psAssert (readout->parent, "missing cell?"); + psAssert (readout->parent->concepts, "missing concepts?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); psTimerStart ("psphot.noise"); @@ -24,5 +57,5 @@ // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + psAssert (maskVal, "missing mask value?"); // increase variance by factor*(object noise): Index: /branches/eam_branches/20091201/psphot/src/psphotApResid.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotApResid.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotApResid.c (revision 26788) @@ -4,5 +4,27 @@ // measure the aperture residual statistics and 2D variations -bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) +// for now, let's store the detections on the readout->analysis for each readout +bool psphotApResid (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { int Nfail = 0; @@ -13,14 +35,26 @@ pmSource *source; - PS_ASSERT_PTR_NON_NULL(config, false); - PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(sources, false); - PS_ASSERT_PTR_NON_NULL(psf, false); - psTimerStart ("psphot.apresid"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // determine the number of allowed threads @@ -33,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; } @@ -288,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); @@ -308,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); @@ -318,5 +352,6 @@ psFree (mag); psFree (dMag); - return false; + return true; + // this is a quality error, not a programming error } Index: /branches/eam_branches/20091201/psphot/src/psphotBasicDeblend.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotBasicDeblend.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotBasicDeblend.c (revision 26788) @@ -1,5 +1,23 @@ # include "psphotInternal.h" -bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) { int N; @@ -8,7 +26,29 @@ pmSource *source, *testSource; + psTimerStart ("psphot.deblend.basic"); + int Nblend = 0; - psTimerStart ("psphot.deblend.basic"); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // 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?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping basic deblend"); + return true; + } + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); float FRACTION = psMetadataLookupF32 (&status, recipe, "DEBLEND_PEAK_FRACTION"); Index: /branches/eam_branches/20091201/psphot/src/psphotBlendFit.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotBlendFit.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotBlendFit.c (revision 26788) @@ -1,6 +1,28 @@ # include "psphotInternal.h" +// for now, let's store the detections on the readout->analysis for each readout +bool psphotBlendFit (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // XXX I don't like this name -bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { +bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { int Nfit = 0; @@ -12,7 +34,24 @@ psTimerStart ("psphot.fit.nonlinear"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // determine the number of allowed threads @@ -40,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 @@ -47,4 +86,8 @@ // source analysis is done in S/N order (brightest first) sources = psArraySort (sources, pmSourceSortBySN); + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend"); + return true; + } // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) @@ -272,109 +315,2 @@ } -# if (0) -bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) { - - bool status = false; - int Nfit = 0; - int Npsf = 0; - int Next = 0; - int Nfail = 0; - - // 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); - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // S/N limit to perform full non-linear fits - float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM"); - - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region)); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - - // skip non-astronomical objects (very likely defects) - if (source->mode & PM_SOURCE_MODE_BLEND) continue; - if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue; - if (source->type == PM_SOURCE_TYPE_DEFECT) continue; - if (source->type == PM_SOURCE_TYPE_SATURATED) continue; - - // skip DBL second sources (ie, added by psphotFitBlob) - if (source->mode & PM_SOURCE_MODE_PAIR) continue; - - // limit selection to some SN limit - if (source->peak->SN < FIT_SN_LIM) continue; - - // exclude sources outside optional analysis region - if (source->peak->xf < AnalysisRegion.x0) continue; - if (source->peak->yf < AnalysisRegion.y0) continue; - if (source->peak->xf > AnalysisRegion.x1) continue; - if (source->peak->yf > AnalysisRegion.y1) continue; - - // if model is NULL, we don't have a starting guess - if (source->modelPSF == NULL) continue; - - // skip sources which are insignificant flux? - // XXX this is somewhat ad-hoc - if (source->modelPSF->params->data.F32[1] < 0.1) { - psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n", - source->modelPSF->params->data.F32[1], - source->modelPSF->params->data.F32[2], - source->modelPSF->params->data.F32[3]); - continue; - } - - // replace object in image - if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } - Nfit ++; - - // try fitting PSFs or extended sources depending on source->mode - // these functions subtract the resulting fitted source - if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) { - if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_EXTENDED; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf); - Next ++; - continue; - } - } else { - if (psphotFitBlend (readout, source, psf, maskVal, markVal)) { - source->type = PM_SOURCE_TYPE_STAR; - psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf); - Npsf ++; - continue; - } - } - - psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf); - Nfail ++; - - // re-subtract the object, leave local sky - pmSourceCacheModel (source, maskVal); - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) - *nfit = Nfit; - *npsf = Npsf; - *next = Next; - *nfail = Nfail; - - // moments are modified by the fit; re-display - psphotVisualPlotMoments (recipe, sources); - psphotVisualShowResidualImage (readout); - - return true; -} -# endif Index: /branches/eam_branches/20091201/psphot/src/psphotChoosePSF.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotChoosePSF.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotChoosePSF.c (revision 26788) @@ -1,18 +1,64 @@ # include "psphotInternal.h" +// generate a PSF model for inputs without PSF models already loaded +bool psphotChoosePSF (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // try PSF models and select best option -pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; psTimerStart ("psphot.choose.psf"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // do not generate a PSF if we already were supplied one + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index); + return true; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model"); + return true; + } // bit-masks to test for good/bad pixels psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); - assert (maskVal); + psAssert (maskVal, "missing mask value?"); // bit-mask to mark pixels not used in analysis psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT"); - assert (markVal); + psAssert (markVal, "missing mark value?"); // maskVal is used to test for rejected pixels, and must include markVal @@ -61,7 +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"); @@ -70,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 @@ -158,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"); @@ -233,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"); @@ -290,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); @@ -302,9 +351,20 @@ psFree (options); - return (psf); + + // save PSF on readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + psFree (psf); // XXX double-check + + // move into psphotChoosePSF + psphotVisualShowPSFModel (readout, psf); + + return true; } // measure average parameters of the PSF model -bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) { +bool psphotPSFstats (pmReadout *readout, pmPSF *psf) { psEllipseShape shape; @@ -312,5 +372,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(psf, false); @@ -362,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)) { @@ -371,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); @@ -387,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 @@ -401,5 +460,4 @@ PS_ASSERT_PTR_NON_NULL(readout, false); - PS_ASSERT_PTR_NON_NULL(recipe, false); PS_ASSERT_PTR_NON_NULL(sources, false); @@ -429,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/20091201/psphot/src/psphotDeblendSatstars.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotDeblendSatstars.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotDeblendSatstars.c (revision 26788) @@ -1,8 +1,27 @@ # include "psphotInternal.h" -bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) { int N; pmSource *source; + bool status; psTimerStart ("psphot.deblend.sat"); @@ -11,5 +30,26 @@ float SAT_MIN_RADIUS = 5.0; - bool status; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // 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?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend"); + return true; + } + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + pmCell *cell = readout->parent; float SATURATION = 0.75*psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION"); Index: /branches/eam_branches/20091201/psphot/src/psphotEfficiency.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotEfficiency.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotEfficiency.c (revision 26788) @@ -6,4 +6,6 @@ //#define TESTING + +# if 0 // Calculate the limiting magnitude for an image @@ -146,16 +148,56 @@ } +# endif + +bool psphotEfficiency (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} // Determine detection efficiency -bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, - psMetadata *recipe, const psArray *realSources) +bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { +# if 0 + bool status = true; + + psTimerStart("psphot.fake"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *realSources = detections->allSources; + psAssert (realSources, "missing sources?"); + + // XXX do we need to skip this step if we do not have a psf? + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); + PM_ASSERT_READOUT_NON_NULL(readout, false); PM_ASSERT_READOUT_IMAGE(readout, false); PS_ASSERT_PTR_NON_NULL(psf, false); PS_ASSERT_METADATA_NON_NULL(recipe, false); - PS_ASSERT_ARRAY_NON_NULL(realSources, false); - - psTimerStart("psphot.fake"); // Collect recipe information @@ -200,5 +242,5 @@ // remove all sources, adding noise for subtracted sources psphotRemoveAllSources(realSources, recipe); - // psphotAddNoise(readout, realSources, recipe); + // psphotAddNoise(readout, realSources, recipe); @@ -474,4 +516,5 @@ psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake")); +# endif return true; Index: /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceAnalysis.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceAnalysis.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceAnalysis.c (revision 26788) @@ -1,6 +1,34 @@ # include "psphotInternal.h" +// for now, let's store the detections on the readout->analysis for each readout +bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full non-linear fits / extended source analysis? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // aperture-like measurements for extended sources -bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -11,13 +39,25 @@ int Nkron = 0; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels assert (maskVal); - - // perform full non-linear fits / extended source analysis? - if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) { - psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); - return true; - } // XXX require petrosian analysis for non-linear fits? Index: /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceFits.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceFits.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotExtendedSourceFits.c (revision 26788) @@ -1,6 +1,34 @@ # include "psphotInternal.h" +// for now, let's store the detections on the readout->analysis for each readout +bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // perform full extended source non-linear fits? + if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); + return true; + } + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotExtendedSourceFitsReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // non-linear model fitting for extended sources -bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe) { +bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -12,4 +40,22 @@ bool savePics = false; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels @@ -21,10 +67,4 @@ // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; - - // perform full extended source non-linear fits? - if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) { - psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n"); - return true; - } // select the collection of desired models Index: /branches/eam_branches/20091201/psphot/src/psphotFake.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotFake.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotFake.c (revision 26788) @@ -211,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/20091201/psphot/src/psphotFindDetections.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotFindDetections.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotFindDetections.c (revision 26788) @@ -1,6 +1,30 @@ # include "psphotInternal.h" +// we store the detections on the readout->analysis for each readout this function finds new +// 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 firstPass) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + 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; + } + } + return true; +} + // smooth the image, search for peaks, optionally define footprints based on the peaks -pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe) { +bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass) { bool status; @@ -9,14 +33,34 @@ int NMAX = 0; + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + psAssert (maskVal, "missing mask value?"); // Use the new pmFootprints approach? const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS"); - // on first pass, detections have not yet been allocated + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + // on the initial pass, detections have not yet been allocated or saved on readout->analysis if (!detections) { + // create the container detections = pmDetectionsAlloc(); + + // save detections on the readout->analysis + 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); @@ -30,6 +74,10 @@ float threshold = PS_SQR(NSIGMA_PEAK); - // move the old peaks array (if it exists) to oldPeaks - // XXX generically, we should be able to call this function an arbitrary number of times + // move the old peaks array (if it exists) to oldPeaks XXX generically, we should be able + // to call this function an arbitrary number of times the old peaks are saved so they can + // be freed later -- the have to be freed after psphotFindFootprints is called below, since + // they are also owned by the oldFootprints, which are in turn merged into the new + // footprints. (what about the source->peak entry?) + assert (detections->oldPeaks == NULL); detections->oldPeaks = detections->peaks; @@ -51,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. @@ -57,5 +105,5 @@ psFree (detections); psError (PSPHOT_ERR_CONFIG, false, "failed on peak search"); - return NULL; + return false; } @@ -71,27 +119,9 @@ psphotVisualShowFootprints (detections); - return detections; + psFree (detections); + + return true; } // if we use the footprints, the output peaks list contains both old and new peaks, // otherwise it only contains the new peaks. - -# if 0 -// XXX where do we place the N sets of detections? -bool psphotFindDetections (pmConfig *config, const pmFPAview *view) -{ - bool status = true; - - int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAbort (!status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} -# endif Index: /branches/eam_branches/20091201/psphot/src/psphotFitSourcesLinear.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotFitSourcesLinear.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotFitSourcesLinear.c (revision 26788) @@ -12,5 +12,23 @@ static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal); -bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) { bool status; @@ -19,4 +37,29 @@ float f; // float r; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + assert (recipe); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (sources, "missing psf?"); psTimerStart ("psphot.linear"); @@ -248,4 +291,8 @@ psphotVisualPlotChisq (sources); // psphotVisualShowFlags (sources); + + // We have to place this visualization here because the models are not realized until + // psphotGuessModels or fitted until psphotFitSourcesLinear. + psphotVisualShowPSFStars (recipe, psf, sources); return true; Index: /branches/eam_branches/20091201/psphot/src/psphotForcedReadout.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotForcedReadout.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotForcedReadout.c (revision 26788) @@ -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/20091201/psphot/src/psphotGuessModels.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotGuessModels.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotGuessModels.c (revision 26788) @@ -7,10 +7,49 @@ // 3) define the threaded function to work with sources for a given cell -// construct an initial PSF model for each object -bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotGuessModels (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +// construct an initial PSF model for each object (new sources only) +bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { bool status; psTimerStart ("psphot.models"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping model guess"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (sources, "missing psf?"); // select the appropriate recipe information @@ -36,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/20091201/psphot/src/psphotImageQuality.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotImageQuality.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotImageQuality.c (revision 26788) @@ -4,7 +4,54 @@ // XXX Lots of code duplication here +// for now, let's store the detections on the readout->analysis for each readout +bool psphotImageQuality (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4 -bool psphotImageQuality(psMetadata *recipe, psArray *sources) +bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { + bool status = true; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // if we have a PSF, skip this analysis + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + return true; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping image quality"); + return true; + } + psVector *FWHM_MAJOR = psVectorAllocEmpty(100, PS_TYPE_F32); psVector *FWHM_MINOR = psVectorAllocEmpty(100, PS_TYPE_F32); @@ -81,10 +128,9 @@ if (num == 0) { - psWarning("Unable to find sources from which to measure image quality"); - return false; - } - - psMetadataAddS32(recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, - "Number of stars used for IQ measurements", M2->n); + psLogMsg ("psphot", PS_LOG_INFO, "no valid sources for image quality, skipping"); + return true; + } + + 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 @@ -98,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)) { @@ -108,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)) { @@ -119,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)) { @@ -132,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)) { @@ -145,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)) { @@ -160,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)) { @@ -175,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 @@ -192,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)) { @@ -202,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)) { @@ -213,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)) { @@ -226,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)) { @@ -239,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)) { @@ -254,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)) { @@ -269,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/20091201/psphot/src/psphotLoadPSF.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotLoadPSF.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotLoadPSF.c (revision 26788) @@ -1,28 +1,69 @@ # 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 -pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { +bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); // find the currently selected chip - pmChip *chip = pmFPAfileThisChip (config->files, view, "PSPHOT.PSF.LOAD"); - if (!chip) return NULL; + pmChip *chip = pmFPAviewThisChip (view, file->fpa); + if (!chip) return false; // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.PSF.LOAD"); - if (!readout) return NULL; + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + if (!readout) return false; // check if a PSF model is supplied by the user - pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSPHOT.PSF"); + pmPSF *psf = psMetadataLookupPtr (&status, chip->analysis, "PSPHOT.PSF"); if (psf == NULL) { psLogMsg ("psphot", 3, "no psf supplied for this chip"); - return NULL; + return true; } - if (!psphotPSFstats (readout, recipe, psf)) psAbort("cannot measure PSF shape terms"); + if (!psphotPSFstats (readout, psf)) { + psAbort("cannot measure PSF shape terms"); + } psLogMsg ("psphot", 3, "using externally supplied PSF model for this readout"); - // we return a psf which can be free (and which may be removed from the metadata) - psMemIncrRefCounter (psf); - return psf; + // save PSF on readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) { + psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout"); + return false; + } + + return true; } + +bool psphotLoadPSF (pmConfig *config, const pmFPAview *view) { + + bool status = false; + + // 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"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // Generate the mask and weight images, including the user-defined analysis region of interest + 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; + } + } + return true; +} Index: /branches/eam_branches/20091201/psphot/src/psphotMagnitudes.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotMagnitudes.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotMagnitudes.c (revision 26788) @@ -1,5 +1,26 @@ # include "psphotInternal.h" -bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) { +bool psphotMagnitudes (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotMagnitudesReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status = false; @@ -8,7 +29,24 @@ psTimerStart ("psphot.mags"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // determine the number of allowed threads @@ -64,5 +102,5 @@ psArrayAdd(job->args, 1, cells->data[j]); // sources - psArrayAdd(job->args, 1, (pmPSF*)psf); // Casting away const + psArrayAdd(job->args, 1, psf); psArrayAdd(job->args, 1, binning); psArrayAdd(job->args, 1, backModel); @@ -179,43 +217,4 @@ } -# if (0) -bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) { - - bool status; - int Nap = 0; - - for (int i = 0; i < sources->n; i++) { - pmSource *source = (pmSource *) sources->data[i]; - status = pmSourceMagnitudes (source, psf, photMode, maskVal); - if (status && isfinite(source->apMag)) Nap ++; - - if (backModel) { - psAssert (binning, "if backModel is defined, so should binning be"); - source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning); - if (isnan(source->sky) && false) { - psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky"); - } - } else { - source->sky = NAN; - } - - if (backStdev) { - psAssert (binning, "if backStdev is defined, so should binning be"); - source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning); - if (isnan(source->skyErr) && false) { - psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr"); - } - } else { - source->skyErr = NAN; - } - } - - // change the value of a scalar on the array (wrap this and put it in psArray.h) - *nap = Nap; - - return true; -} -# endif - bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) { Index: /branches/eam_branches/20091201/psphot/src/psphotMakePSFReadout.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotMakePSFReadout.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotMakePSFReadout.c (revision 26788) @@ -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/20091201/psphot/src/psphotMaskReadout.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotMaskReadout.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotMaskReadout.c (revision 26788) @@ -1,12 +1,40 @@ # include "psphotInternal.h" +bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // Generate the mask and weight images, including the user-defined analysis region of interest + if (!psphotSetMaskAndVarianceReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + // generate mask and variance if not defined, additional mask for restricted subregion -bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe) { +bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; - // ** 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 + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest + psAssert (file, "missing file?"); + + // find the currently selected readout + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + psAssert (readout, "missing readout?"); + + // 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); @@ -14,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); } @@ -73,33 +101,7 @@ } + // display the image, weight, mask (ch 1,2,3) + psphotVisualShowImage (readout); + return true; } - -bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe) { - - bool status = false; - - int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest - PS_ASSERT_PTR_NON_NULL (file, false); - - // find the currently selected readout - pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); - PS_ASSERT_PTR_NON_NULL (readout, false); - - // Generate the mask and weight images, including the user-defined analysis region of interest - if (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i); - return false; - } - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); - } - return true; -} Index: /branches/eam_branches/20091201/psphot/src/psphotMergeSources.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotMergeSources.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotMergeSources.c (revision 26788) @@ -5,10 +5,89 @@ PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply for PSF sources -// merge the externally supplied sources with the existing sources. mark them as having -// mode PM_SOURCE_MODE_EXTERNAL -bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources) { - - psArray *extSourcesCMF = NULL; +// for now, let's store the detections on the readout->analysis for each readout +bool psphotMergeSources (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotMergeSourcesReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +// add newly selected sources to the existing list of sources +bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *newSources = detections->newSources; + 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); + } + psArray *allSources = detections->allSources; + + for (int i = 0; i < newSources->n; i++) { + pmSource *source = newSources->data[i]; + psArrayAdd (allSources, 100, source); + } + + 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: @@ -17,9 +96,9 @@ 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]; - source->mode |= PM_SOURCE_MODE_EXTERNAL; + 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 @@ -27,10 +106,11 @@ source->peak->value = 1.0; - // drop the loaded source modelPSF - psFree (source->modelPSF); - source->modelPSF = NULL; - } - psphotMergeSources (sources, extSourcesCMF); - } + // drop the loaded source modelPSF + psFree (source->modelPSF); + source->modelPSF = NULL; + + psArrayAdd (detections->allSources, 100, source); + } + } } @@ -48,20 +128,21 @@ source->mode |= PM_SOURCE_MODE_EXTERNAL; - // drop the loaded source modelPSF - psFree (source->modelPSF); - source->modelPSF = NULL; - } - psphotMergeSources (sources, extSourcesTXT); - } + // drop the loaded source modelPSF + psFree (source->modelPSF); + source->modelPSF = NULL; + + psArrayAdd (detections->allSources, 100, source); + } + } } 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; @@ -71,16 +152,9 @@ } -// add newly selected sources to the existing list of sources -bool psphotMergeSources (psArray *oldSources, psArray *newSources) { - - for (int i = 0; i < newSources->n; i++) { - pmSource *source = newSources->data[i]; - psArrayAdd (oldSources, 100, source); - } - return true; -} - // extract the input sources corresponding to this readout +// 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 @@ -91,7 +165,14 @@ } - 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"); + return NULL; + } + + psArray *sources = detections->allSources; if (!sources) { psLogMsg ("psphot", 3, "no psf sources for this readout"); + return NULL; } @@ -99,6 +180,45 @@ } +// this function is used to fix sources which were loaded externally, but have passed from +// psphotDetectionsFromSources to psphotSourceStats and are now stored on +// detections->newSources. +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->newSources; + 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(); @@ -113,5 +233,5 @@ float snMin = psMetadataLookupF32(NULL, recipe, "MOMENTS_SN_MIN"); if (!isfinite(snMin)) { - return NULL; + return false; } @@ -152,5 +272,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; } @@ -194,2 +320,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/20091201/psphot/src/psphotModelBackground.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotModelBackground.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotModelBackground.c (revision 26788) @@ -32,5 +32,10 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotModelBackgroundReadout(psImage *model, // Model image +// +// 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 psMetadata *analysis, // Analysis metadata for outputs @@ -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); @@ -356,8 +354,9 @@ // find the currently selected readout pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - PS_ASSERT_PTR_NON_NULL (file, false); + psAssert (file, "missing file?"); pmFPA *inFPA = file->fpa; pmReadout *readout = pmFPAviewThisReadout(view, inFPA); + psAssert (readout, "missing readout?"); psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters Index: /branches/eam_branches/20091201/psphot/src/psphotOutput.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotOutput.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotOutput.c (revision 26788) @@ -126,22 +126,4 @@ } -bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { - - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest - PS_ASSERT (file, false); - - pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); - - // determine PHOTCODE from fpa & view, overwrite in readout->analysis - char *photcode = pmConceptsPhotcodeForView (file, view); - PS_ASSERT (photcode, false); - - psMetadataAddStr (readout->analysis, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "photcode from FPA concepts", photcode); - psLogMsg ("psphot", 3, "PHOTCODE is %s", photcode); - - psFree (photcode); - return true; -} - bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view) { @@ -161,5 +143,23 @@ } -bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources) { +bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) { + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest + PS_ASSERT (file, false); + + pmReadout *readout = pmFPAviewThisReadout (view, file->fpa); + + // determine PHOTCODE from fpa & view, overwrite in readout->analysis + char *photcode = pmConceptsPhotcodeForView (file, view); + PS_ASSERT (photcode, false); + + psMetadataAddStr (readout->analysis, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "photcode from FPA concepts", photcode); + psLogMsg ("psphot", 3, "PHOTCODE is %s", photcode); + + psFree (photcode); + return true; +} + +bool psphotSetHeaderNstars (psMetadata *header, psArray *sources) { int nSrc = 0; @@ -190,15 +190,14 @@ } - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NSTARS", PS_META_REPLACE, "Number of sources", nSrc); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT); - psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_CR", PS_META_REPLACE, "Number of cosmic rays", nCR); - return true; -} - -// 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 -psMetadata *psphotDefineHeader (psMetadata *recipe) { + psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS", PS_META_REPLACE, "Number of sources", nSrc); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT); + psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR", PS_META_REPLACE, "Number of cosmic rays", nCR); + return true; +} + +// these values are saved in an output header stub - they are added to either the PHU header +// (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/20091201/psphot/src/psphotRadiusChecks.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotRadiusChecks.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotRadiusChecks.c (revision 26788) @@ -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/20091201/psphot/src/psphotReadout.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadout.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReadout.c (revision 26788) @@ -3,5 +3,4 @@ // this should be called by every program that links against libpsphot bool psphotInit (void) { - psphotErrorRegister(); // register our error codes/messages psphotModelClassInit (); // load implementation-specific models @@ -12,8 +11,6 @@ bool psphotReadout(pmConfig *config, const pmFPAview *view) { - // measure the total elapsed time in psphotReadout. XXX the current threading plan - // for psphot envisions threading within psphotReadout, not multiple threads calling - // the same psphotReadout. In the current plan, this dtime is the elapsed time used - // jointly by the multiple threads, not the total time used by all threads. + // measure the total elapsed time in psphotReadout. dtime is the elapsed time used jointly + // by the multiple threads, not the total time used by all threads. psTimerStart ("psphotReadout"); @@ -26,4 +23,7 @@ return false; } + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + psAssert (breakPt, "configuration error: set BREAK_POINT"); // set the photcode for this image @@ -33,179 +33,165 @@ } - // 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"); - PS_ASSERT_PTR_NON_NULL (breakPt, false); - // 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); - } - - // display the image, weight, mask (ch 1,2,3) - psphotVisualShowImage (readout); + return psphotReadoutCleanup(config, view); + } // 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); - } - - // display the backsub and backgnd images - psphotVisualShowBackground (config, view, readout); - - // run a single-model test if desired (exits from here if test is run) - psphotModelTest (config, view, recipe); - - // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are saved in the recipe - pmPSF *psf = psphotLoadPSF (config, view, recipe); - - // several functions below behave differently if we have a PSF model already - bool havePSF = (psf != NULL); - + 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. - pmDetections *detections = psphotFindDetections (NULL, readout, recipe); - if (!detections) { + 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"); - 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); - } - - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, true); - if (!sources) return false; + return psphotReadoutCleanup (config, view); + } + + // construct sources and measure basic stats (saved on detections->newSources) + if (!psphotSourceStats (config, view, true)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); + } 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); + return psphotReadoutCleanup(config, view); + } + + // find blended neighbors of very saturated stars (detections->newSources) + if (!psphotDeblendSatstars (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + return psphotReadoutCleanup (config, view); + } + + // mark blended peaks PS_SOURCE_BLEND (detections->newSources) + if (!psphotBasicDeblend (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); + return psphotReadoutCleanup (config, view); + } + + // classify sources based on moments, brightness. if a PSF model has been loaded, the PSF + // clump defined for it is used not measured (detections->newSources) + if (!psphotRoughClass (config, view)) { // pass 1 + psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); + return psphotReadoutCleanup (config, view); } 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; + return psphotReadoutCleanup(config, view); + } + + // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) + if (!psphotImageQuality (config, view)) { // pass 1 + psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); + return psphotReadoutCleanup(config, view); + } + + // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, + // this step is skipped + if (!psphotChoosePSF (config, view)) { // pass 1 + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - psphotVisualShowPSFModel (readout, psf); + return psphotReadoutCleanup (config, view); + } // 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); + // XXX fix this in the new multi-input context + psphotLoadExtSources (config, view); // pass 1 + + // construct an initial model for each object, set the radius to fitRadius, set circular + // fit mask (detections->newSources) + psphotGuessModels (config, view); // pass 1 + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (config, view); // 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); + psphotFitSourcesLinear (config, view, false); // pass 1 (detections->allSources) + + // identify CRs and extended sources (only unmeasured sources are measured) + psphotSourceSize (config, view, true); // pass 1 (detections->allSources) 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); + psphotBlendFit (config, view); // pass 1 (detections->allSources) // replace all sources - psphotReplaceAllSources (sources, recipe); + psphotReplaceAllSources (config, view); // pass 1 (detections->allSources) // linear fit to include all sources (subtract again) - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); + // NOTE : apply to ALL sources (extended + psf) + psphotFitSourcesLinear (config, view, true); // pass 2 (detections->allSources) // 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 + if (!strcasecmp (breakPt, "PASS1")) goto pass1finish; + + // NOTE: possibly re-measure background model here with objects subtracted / or masked // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe); - - // find fainter sources (pass 2) - detections = psphotFindDetections (detections, readout, recipe); + psphotAddNoise (config, view); // pass 1 (detections->allSources) + + // find fainter sources + // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections + psphotFindDetections (config, view, false); // pass 2 (detections->peaks, detections->footprints) // remove noise for subtracted objects (ie, return to normal noise level) - psphotSubNoise (readout, sources, recipe); + // NOTE: this needs to operate only on the OLD sources + psphotSubNoise (config, view); // pass 1 (detections->allSources) // define new sources based on only the new peaks - psArray *newSources = psphotSourceStats (config, readout, detections, false); + // NOTE: new sources are saved on detections->newSources + psphotSourceStats (config, view, false); // pass 2 (detections->newSources) // set source type - if (!psphotRoughClass (readout, newSources, recipe, havePSF)) { + // NOTE: apply only to detections->newSources + if (!psphotRoughClass (config, view)) { // pass 2 (detections->newSources) psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // create full input models, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, newSources, psf); + // NOTE: apply only to detections->newSources + psphotGuessModels (config, view); // pass 2 (detections->newSources) // replace all sources so fit below applies to all at once - psphotReplaceAllSources (sources, recipe); + // NOTE: apply only to OLD sources (which have been subtracted) + psphotReplaceAllSources (config, view); // pass 2 // 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); + // NOTE: merge OLD and NEW + // XXX check on free of sources... + psphotMergeSources (config, view); // (detections->newSources + detections->allSources -> detections->allSources) + + // NOTE: apply to ALL sources + psphotFitSourcesLinear (config, view, true); // pass 3 (detections->allSources) pass1finish: // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, psf, 0); - - psphotExtendedSourceAnalysis (readout, sources, recipe); - - psphotExtendedSourceFits (readout, sources, recipe); + // NOTE: applies only to NEW (unmeasured) sources + psphotSourceSize (config, view, false); // pass 2 (detections->allSources) + + psphotExtendedSourceAnalysis (config, view); // pass 1 (detections->allSources) + psphotExtendedSourceFits (config, view); // pass 1 (detections->allSources) finish: @@ -215,13 +201,13 @@ // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { + if (!psphotApResid (config, view)) { // pass 1 (detections->allSources) 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); - - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + psphotMagnitudes(config, view); // pass 1 (detections->allSources) + + if (!psphotEfficiency(config, view)) { // pass 1 psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -232,11 +218,10 @@ // replace background in residual image - psphotSkyReplace (config, view); + psphotSkyReplace (config, view); // pass 1 // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (config, view); // pass 1 // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup(config, view); } - Index: /branches/eam_branches/20091201/psphot/src/psphotReadoutCleanup.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadoutCleanup.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReadoutCleanup.c (revision 26788) @@ -1,8 +1,8 @@ # include "psphotInternal.h" -// psphotReadoutCleanup is called on exit from psphotReadout. If the last raised error is -// not a DATA error, then there was a serious problem. Only in this case, or if the fail -// on the stats measurement, do we return false -bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view) +{ + bool status = true; // remove internal pmFPAfiles, if created @@ -12,17 +12,52 @@ } if (psErrorCodeLast() != PS_ERR_NONE) { - psFree (psf); - psFree (sources); - psFree (detections); return false; } + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i); + return false; + } + } + + // XXX move this to top of loop? + pmKapaClose (); + + return true; +} + +// psphotReadoutCleanup is called on exit from psphotReadout. If the last raised error is +// not a DATA error, then there was a serious problem. Only in this case, or if the fail +// on the stats measurement, do we return false +bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { + + bool status = true; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // when psphotReadoutCleanup is called, these are not necessarily defined + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psArray *sources = detections->allSources; + // XXX where do we free these, in here (psMetadataRemove?) + // 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"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } @@ -30,9 +65,6 @@ // 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"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } @@ -40,4 +72,5 @@ // Check to see if the image quality was measured + // XXX not sure we want / need this test if (!psf) { bool mdok; // Status of MD lookup @@ -45,24 +78,21 @@ if (!mdok || nIQ <= 0) { psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality"); - psFree (psf); - psFree (sources); - psFree (detections); return false; } } + // create an output header with stats results currently saved on readout->analysis + psMetadata *header = psphotDefineHeader (readout->analysis); // write NSTARS to the image header - psphotSetHeaderNstars (recipe, sources); - - // create an output header with stats results - psMetadata *header = psphotDefineHeader (recipe); + psphotSetHeaderNstars (header, sources); // save the results of the analysis + // this should happen way up stream (when needed?) psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER", PS_DATA_METADATA, "header stats", header); - if (sources) { - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources); - } + if (psf) { + // XXX this seems a little silly : we saved the psf on readout->analysis above, but now + // we are moving it to chip->analysis. // save the psf for possible output. if there was already an entry, it was loaded from external sources // the new one may have been updated or modified, so replace the existing entry. We @@ -79,12 +109,5 @@ } - // XXX move this to top of loop? - pmKapaClose (); - - psFree (detections); - psFree (psf); psFree (header); - psFree (sources); - return true; } Index: /branches/eam_branches/20091201/psphot/src/psphotReadoutFindPSF.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadoutFindPSF.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReadoutFindPSF.c (revision 26788) @@ -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; + // construct detections->newSources and measure basic stats (moments, local sky) + 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,23 @@ // 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); +# endif + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (config, view); + +# if 0 // 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/20091201/psphot/src/psphotReadoutKnownSources.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadoutKnownSources.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReadoutKnownSources.c (revision 26788) @@ -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,55 @@ // 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); + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (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/20091201/psphot/src/psphotReadoutMinimal.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadoutMinimal.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReadoutMinimal.c (revision 26788) @@ -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,96 +23,69 @@ } - // 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); + return psphotReadoutCleanup (config, view); } -#if 0 - if (!detections->peaks->n) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); + + // 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); } -#endif - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, false); - if (!sources) return false; + // find blended neighbors of very saturated stars + psphotDeblendSatstars (config, view); - // Do the efficiency before everything else, to ensure it's done - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + // 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 (use supplied psf shape parameters) + if (!psphotRoughClass (config, view)) { + psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); + return psphotReadoutCleanup (config, view); + } + + // construct an initial model for each object + psphotGuessModels (config, view); + + // linear PSF fit to source peaks + psphotFitSourcesLinear (config, view, false); + +// XXX eventually, add the extended source fits here +# if (0) + // measure source size for the remaining sources + psphotSourceSize (config, view); + + psphotExtendedSourceAnalysis (config, view); + + psphotExtendedSourceFits (config, view); +# endif + + // calculate source magnitudes + psphotMagnitudes(config, view); + + // 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(); } - if (detections->peaks->n == 0 || sources->n == 0) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); - } - - - - // 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 (use supplied psf shape parameters) - if (!psphotRoughClass (readout, sources, recipe, true)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // construct an initial model for each object - psphotGuessModels (config, readout, sources, psf); - - // 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); - -// XXX eventually, add the extended source fits here -# if (0) - // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, 0); - - psphotExtendedSourceAnalysis (readout, sources, recipe); - - psphotExtendedSourceFits (readout, sources, recipe); -# endif - - // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); - // 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/20091201/psphot/src/psphotReadoutStack.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReadoutStack.c (revision 26787) +++ (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/20091201/psphot/src/psphotReplaceUnfit.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotReplaceUnfit.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotReplaceUnfit.c (revision 26788) @@ -22,5 +22,27 @@ } -bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe) { +// for now, let's store the detections on the readout->analysis for each readout +bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; @@ -29,7 +51,20 @@ psTimerStart ("psphot.replace"); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels - assert (maskVal); + psAssert (maskVal, "missing mask value?"); for (int i = 0; i < sources->n; i++) { @@ -67,42 +102,2 @@ return true; } - -# if (0) -// add source, if the source has been subtracted; do not modify state -bool psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal) { - - // what is current state? (true : add; false : sub) - bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (state && useState) return true; - - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - return true; -} - -// sub source, if the source has been added; do not modify state -bool psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal) { - - // what is current state? (true : sub; false : add) - bool state = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (state && useState) return true; - - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - return true; -} - -// add or sub source to match recorded state: supply current state as true (add) or false (sub) -bool psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal) { - - // what is desired state? (true : add; false : sub) - bool newState = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED); - if (curState == newState) return true; - - if (curState && !newState) { - pmSourceSub (source, PM_MODEL_OP_FULL, maskVal); - } - if (newState && !curState) { - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - } - return true; -} -# endif Index: /branches/eam_branches/20091201/psphot/src/psphotRoughClass.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotRoughClass.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotRoughClass.c (revision 26788) @@ -7,12 +7,55 @@ } } -bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF); +// for now, let's store the detections on the readout->analysis for each readout +bool psphotRoughClass (pmConfig *config, const pmFPAview *view) +{ + bool status = true; -// 2006.02.02 : no leaks -bool psphotRoughClass (pmReadout *readout, psArray *sources, psMetadata *recipe, const bool havePSF) { + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status; psTimerStart ("psphot.rough"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + // if we have a PSF, use the existing PSF clump region below + bool havePSF = false; + if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) { + havePSF = true; + } + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->newSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification"); + return true; + } // we make this measurement on a NxM grid of regions across the readout @@ -75,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); @@ -103,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/20091201/psphot/src/psphotSetThreads.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSetThreads.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSetThreads.c (revision 26788) @@ -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/20091201/psphot/src/psphotSignificanceImage.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSignificanceImage.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSignificanceImage.c (revision 26788) @@ -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")); @@ -82,7 +81,6 @@ // 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/20091201/psphot/src/psphotSkyReplace.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSkyReplace.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSkyReplace.c (revision 26788) @@ -1,13 +1,33 @@ # include "psphotInternal.h" + +bool psphotSkyReplace (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} // XXX make this an option? // in order to successfully replace the sky, we must define a corresponding file... -bool psphotSkyReplace (pmConfig *config, const pmFPAview *view) { +bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) { psTimerStart ("psphot.skyreplace"); // find the currently selected readout - pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT"); - if (readout == NULL) psAbort("input not defined"); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); // select background pixels, from output background file, or create Index: /branches/eam_branches/20091201/psphot/src/psphotSourceFreePixels.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSourceFreePixels.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSourceFreePixels.c (revision 26788) @@ -1,7 +1,37 @@ # include "psphotInternal.h" -void psphotSourceFreePixels (psArray *sources) { +bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view) +{ + bool status = true; - if (!sources) return; + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSourceFreePixelsReadout (config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index) { + + bool status; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); for (int i = 0; i < sources->n; i++) { @@ -9,4 +39,4 @@ pmSourceFreePixels (source); } - return; + return true; } Index: /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSourceSize.c (revision 26788) @@ -20,22 +20,6 @@ bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options); bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options); -bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); -bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); -bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal); - -bool psphotMaskCosmicRayFootprintCheck (psArray *sources) { - - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - pmPeak *peak = source->peak; - pmFootprint *footprint = peak->footprint; - if (!footprint) continue; - for (int j = 0; j < footprint->spans->n; j++) { - pmSpan *sp = footprint->spans->data[j]; - psAssert (sp, "missing span"); - } - } - return true; -} +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal); +bool psphotMaskCosmicRayFootprintCheck (psArray *sources); // we need to call this function after sources have been fitted to the PSF model and @@ -45,6 +29,28 @@ // deviation from the psf model at the r = FWHM/2 position -// XXX use an internal flag to mark sources which have already been measured -bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first) +// for now, let's store the detections on the readout->analysis for each readout +bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) { + psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +// this function use an internal flag to mark sources which have already been measured +bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize) { bool status; @@ -52,4 +58,25 @@ psTimerStart ("psphot.size"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size"); + return true; + } + + pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); + psAssert (psf, "missing psf?"); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -91,127 +118,25 @@ // XXX move this to the code that generates the PSF? // XXX store the results on pmPSF? - psphotSourceSizePSF (&options, sources, psf); + + // XXX this should only be done on the first pass (ie, if we have newSources or allSources?) + if (getPSFsize) { + psphotSourceSizePSF (&options, sources, psf); + } // classify the sources based on ApResid and Moments (extended sources) + // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) psphotSourceClass(readout, sources, recipe, psf, &options); + // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) psphotSourceSizeCR (readout, sources, &options); - psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size")); + // XXX fix this (was source->n - first) + psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size")); psphotVisualPlotSourceSize (recipe, readout->analysis, sources); psphotVisualShowSourceSize (readout, sources); psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr); - - return true; -} - -// This attempt to mask the cosmic rays used the isophotal boundary -bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - // replace the source flux - pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); - - // flag this as a CR - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - // grab the matching footprint - pmFootprint *footprint = peak->footprint; - if (!footprint) { - psTrace("psphot.czw",2,"Using isophot CR mask code."); - - // if we have not footprint, use the old code to mask by isophot - psphotMaskCosmicRayIsophot (source, maskVal, crMask); - return true; - } - - if (!footprint->spans) { - psTrace("psphot.czw",2,"Using isophot CR mask code."); - - // if we have no footprint, use the old code to mask by isophot - psphotMaskCosmicRayIsophot (source, maskVal, crMask); - return true; - } - psphotMaskCosmicRayIsophot (source, maskVal, crMask); - // mask all of the pixels covered by the spans of the footprint - for (int j = 1; j < footprint->spans->n; j++) { - pmSpan *span1 = footprint->spans->data[j]; - - int iy = span1->y; - int xs = span1->x0; - int xe = span1->x1; - - for (int ix = xs; ix < xe; ix++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } - return true; -} - -bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { - - source->mode |= PM_SOURCE_MODE_CR_LIMIT; - pmPeak *peak = source->peak; - psAssert (peak, "NULL peak"); - - psImage *mask = source->maskView; - psImage *pixels = source->pixels; - psImage *variance = source->variance; - - // XXX This should be a recipe variable -# define SN_LIMIT 5.0 - - int xo = peak->x - pixels->col0; - int yo = peak->y - pixels->row0; - - // mark the pixels in this row to the left, then the right - for (int ix = xo; ix >= 0; ix--) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } - } - for (int ix = xo + 1; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); - if (SN > SN_LIMIT) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; - } - } - - // for each of the neighboring rows, mark the high pixels if they have a marked neighbor - // first go up: - for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } - // next go down: - for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { - // mark the pixels in this row to the left, then the right - for (int ix = 0; ix < pixels->numCols; ix++) { - float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); - if (SN < SN_LIMIT) continue; - - bool valid = false; - valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); - valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; - valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; - - if (!valid) continue; - mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; - } - } + psphotVisualShowSatStars (recipe, psf, sources); + return true; } @@ -254,6 +179,6 @@ float dMag = source->psfMag - apMag; - psVectorAppend (Ap, 100, dMag); - psVectorAppend (ApErr, 100, source->errMag); + psVectorAppend (Ap, dMag); + psVectorAppend (ApErr, source->errMag); } @@ -454,78 +379,4 @@ } -// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total -// deviation in sigmas. This is measured on the residual image - should we ignore negative -// deviations? NOTE: This function was an early attempt to classify extended objects, and is -// no longer used by psphot. -float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, - psImageMaskType maskVal, const pmModel *model, float Ro) -{ - psF32 *PAR = model->params->data.F32; // Model parameters - float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters - - // We treat the contour as an ellipse: - // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY - // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; - // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro - // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so: - // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2] - - // min/max value of x is where B^2 - 4AC = 0; solve this for x - float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); - if (Q < 0.0) { - // ellipse is imaginary - return NAN; - } - - int radius = sqrtf(Q) + 0.5; // Radius of ellipse - int nPts = 0; // Number of points in ellipse - float nSigma = 0.0; // - - for (int x = -radius; x <= radius; x++) { - // Polynomial coefficients - // XXX Should we be using the centre of the pixel as x or x+0.5? - float A = PS_SQR (1.0 / syy); - float B = x * sxy; - float C = PS_SQR (x / sxx) - Ro; - float T = PS_SQR(B) - 4*A*C; - if (T < 0.0) { - continue; - } - - // y position in source frame - float yP = (-B + sqrt (T)) / (2.0 * A); - float yM = (-B - sqrt (T)) / (2.0 * A); - - // Get the closest pixel positions (image frame) - int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; - int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; - - if (xPix < 0 || xPix >= image->numCols) { - continue; - } - - if (yPixM >= 0 && yPixM < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); - nSigma += dSigma; - nPts++; - } - - if (yPixM == yPixP) { - continue; - } - - if (yPixP >= 0 && yPixP < image->numRows && - !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { - float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); - nSigma += dSigma; - nPts++; - } - } - nSigma /= nPts; - return nSigma; -} - // given an object suspected to be a defect, generate a pixel mask using the Lapacian transform // if enough of the object is detected as 'sharp', consider the object a cosmic ray @@ -534,4 +385,10 @@ for (int i = 0; i < sources->n; i++) { pmSource *source = sources->data[i]; + + // skip source if it was already measured + if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) { + psTrace("psphot", 7, "Not calculating source size since it has already been measured\n"); + continue; + } // Integer position of peak @@ -582,5 +439,5 @@ // does no repair or recovery of the CR pixels, it only masks them out. My test code can be // found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c -bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal) { +bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) { // Get the actual images and information about the peak. @@ -760,4 +617,207 @@ } +bool psphotMaskCosmicRayFootprintCheck (psArray *sources) { + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + pmFootprint *footprint = peak->footprint; + if (!footprint) continue; + for (int j = 0; j < footprint->spans->n; j++) { + pmSpan *sp = footprint->spans->data[j]; + psAssert (sp, "missing span"); + } + } + return true; +} + +/**** ------ old versions of cosmic ray masking ----- ****/ + +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask); + +// This attempt to mask the cosmic rays used the isophotal boundary +bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + // replace the source flux + pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal); + + // flag this as a CR + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + // grab the matching footprint + pmFootprint *footprint = peak->footprint; + if (!footprint) { + psTrace("psphot.czw",2,"Using isophot CR mask code."); + + // if we have not footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + + if (!footprint->spans) { + psTrace("psphot.czw",2,"Using isophot CR mask code."); + + // if we have no footprint, use the old code to mask by isophot + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + return true; + } + psphotMaskCosmicRayIsophot (source, maskVal, crMask); + // mask all of the pixels covered by the spans of the footprint + for (int j = 1; j < footprint->spans->n; j++) { + pmSpan *span1 = footprint->spans->data[j]; + + int iy = span1->y; + int xs = span1->x0; + int xe = span1->x1; + + for (int ix = xs; ix < xe; ix++) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) { + + source->mode |= PM_SOURCE_MODE_CR_LIMIT; + pmPeak *peak = source->peak; + psAssert (peak, "NULL peak"); + + psImage *mask = source->maskView; + psImage *pixels = source->pixels; + psImage *variance = source->variance; + + // XXX This should be a recipe variable +# define SN_LIMIT 5.0 + + int xo = peak->x - pixels->col0; + int yo = peak->y - pixels->row0; + + // mark the pixels in this row to the left, then the right + for (int ix = xo; ix >= 0; ix--) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + for (int ix = xo + 1; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]); + if (SN > SN_LIMIT) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask; + } + } + + // for each of the neighboring rows, mark the high pixels if they have a marked neighbor + // first go up: + for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + // next go down: + for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) { + // mark the pixels in this row to the left, then the right + for (int ix = 0; ix < pixels->numCols; ix++) { + float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]); + if (SN < SN_LIMIT) continue; + + bool valid = false; + valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask); + valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0; + valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0; + + if (!valid) continue; + mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask; + } + } + return true; +} + +// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total +// deviation in sigmas. This is measured on the residual image - should we ignore negative +// deviations? NOTE: This function was an early attempt to classify extended objects, and is +// no longer used by psphot. +float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask, + psImageMaskType maskVal, const pmModel *model, float Ro) +{ + psF32 *PAR = model->params->data.F32; // Model parameters + float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters + + // We treat the contour as an ellipse: + // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY + // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0; + // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro + // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so: + // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2] + + // min/max value of x is where B^2 - 4AC = 0; solve this for x + float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0); + if (Q < 0.0) { + // ellipse is imaginary + return NAN; + } + + int radius = sqrtf(Q) + 0.5; // Radius of ellipse + int nPts = 0; // Number of points in ellipse + float nSigma = 0.0; // + + for (int x = -radius; x <= radius; x++) { + // Polynomial coefficients + // XXX Should we be using the centre of the pixel as x or x+0.5? + float A = PS_SQR (1.0 / syy); + float B = x * sxy; + float C = PS_SQR (x / sxx) - Ro; + float T = PS_SQR(B) - 4*A*C; + if (T < 0.0) { + continue; + } + + // y position in source frame + float yP = (-B + sqrt (T)) / (2.0 * A); + float yM = (-B - sqrt (T)) / (2.0 * A); + + // Get the closest pixel positions (image frame) + int xPix = x + PAR[PM_PAR_XPOS] - image->col0 + 0.5; + int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5; + + if (xPix < 0 || xPix >= image->numCols) { + continue; + } + + if (yPixM >= 0 && yPixM < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]); + nSigma += dSigma; + nPts++; + } + + if (yPixM == yPixP) { + continue; + } + + if (yPixP >= 0 && yPixP < image->numRows && + !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) { + float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]); + nSigma += dSigma; + nPts++; + } + } + nSigma /= nPts; + return nSigma; +} + // this was an old attempt to identify cosmic rays based on the peak curvature bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) { @@ -893,2 +953,3 @@ return true; } + Index: /branches/eam_branches/20091201/psphot/src/psphotSourceStats.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSourceStats.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSourceStats.c (revision 26788) @@ -1,7 +1,28 @@ # include "psphotInternal.h" -bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources); - -psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) { +// convert detections to sources and measure their basic properties (moments, local sky, sky +// variance) Note: this function only generates sources for the new peaks (peak->assigned) +bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow) +{ + bool status = true; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i); + return false; + } + } + return true; +} + +bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) { bool status = false; @@ -10,7 +31,22 @@ psTimerStart ("psphot.stats"); - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + 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 @@ -21,20 +57,39 @@ // determine properties (sky, moments) of initial sources - float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); - if (!status) return NULL; - + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + psAssert (status, "missing SKY_OUTER_RADIUS in recipe?"); + + // XXX this seems like an arbitrary number... 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 NULL; + 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; if (!peaks) { psError(PS_ERR_UNEXPECTED_NULL, false, "No peaks found!"); - return NULL; + return false; } // generate the array of sources, define the associated pixel sources = psArrayAllocEmpty (peaks->n); + + // if there are no peaks, we save the empty source array and return + if (!peaks->n) { + // save the new sources on the detection structure: + detections->newSources = sources; + return true; + } for (int i = 0; i < peaks->n; i++) { @@ -58,4 +113,159 @@ psArrayAdd (sources, 100, source); psFree (source); + } + + 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); + detections->newSources = sources; + return true; + } + + if (setWindow) { + if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); + psFree(sources); + return false; + } + } + + // 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; + 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 + + 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 + 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); + psFree(sources); + return false; + } + 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"); + 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; } @@ -67,9 +277,8 @@ } - if (setWindow) { - if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); - return NULL; - } + // 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; } @@ -94,4 +303,5 @@ 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); @@ -141,124 +351,4 @@ } -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 NULL; - - 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 NULL; - - 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); -} - bool psphotSourceStats_Threaded (psThreadJob *job) { @@ -268,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 @@ -349,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; @@ -367,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? @@ -398,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]); @@ -442,11 +523,15 @@ // 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); return true; } + +// 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/20091201/psphot/src/psphotSubtractBackground.c =================================================================== --- /branches/eam_branches/20091201/psphot/src/psphotSubtractBackground.c (revision 26787) +++ /branches/eam_branches/20091201/psphot/src/psphotSubtractBackground.c (revision 26788) @@ -4,5 +4,5 @@ // generate the median in NxN boxes, clipping heavily // linear interpolation to generate full-scale model -bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) +bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) { bool status = true; @@ -25,8 +25,4 @@ pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa); assert (model); - - // select the appropriate recipe information - psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); - assert (recipe); // user-defined masks to test for good/bad pixels (build from recipe list if not yet set) @@ -114,4 +110,8 @@ // the pmReadout selected in this function are all view on entries in config->files + // display the backsub and backgnd images + // move this inthe the subtract background loop + psphotVisualShowBackground (config, view, readout); + npass ++; return true; @@ -122,4 +122,8 @@ bool status = false; + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); @@ -127,5 +131,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i)) { + if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i); return false;