IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26788


Ignore:
Timestamp:
Feb 5, 2010, 1:38:43 PM (16 years ago)
Author:
eugene
Message:

updates to psphot APIs to enable stack photometry

Location:
branches/eam_branches/20091201/psphot
Files:
1 deleted
41 edited
1 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psphot

  • branches/eam_branches/20091201/psphot/doc/stack.txt

    r26542 r26788  
     1
     220100126:
     3
     4  * watch out for psphotSetMomentsWindow & MOMENTS_SX_MAX,etc
     5  * watch out for psphotSignificanceImage.c:,psphotEfficiency using the FWHM_MAJ from psphotChoosePSF
     6    * ppSimDetections.c : SIGMA_SMOOTH
     7ppSmooth/src/ppSmoothReadout.c:    psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
     8ppSmooth/src/ppSmoothReadout.c:    psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
     9
     10
     1120100120 : more stack processing mods:
     12
     13  there are a number of data collections generated during psphotReadout:
     14
     15  * detections (currently a structure containing multiple arrays)
     16  * sources (a psArray)
     17  * psf (the psf model structure)
     18 
     19  * new sources vs old sources:
     20    there is a sequence on the second pass in which we need to distinguish the sources
     21    from the first pass from those on the second pass:
     22
     23    - add noise (old sources only)
     24    - find detections
     25    - subtract noise (old sources only)
     26    - generate sources (new detections only)
     27    - source classifications (new sources only)
     28    - guess models (new sources only)
     29    - replace sources (old sources only)
     30    - merge sources (new + old -> sources)
     31
     32    currently we are distiguishing the old vs new based on different arrays.
     33    can we use the processing flags to distinguish the these cases and carry around
     34    only a single source list?
     35
     36    * detections->peaks holds the most recently detected set of peaks
     37      detections->oldPeaks holds the previous collection
     38
     39      detections->footprints holds the full set of merged footprints,
     40      including assigned peaks from the old and new set.
    141
    24220100107 : updates for stack processing

  • branches/eam_branches/20091201/psphot/src/Makefile.am

    r25984 r26788  
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotTest psphotMomentsStudy
    28 # bin_PROGRAMS = psphotPetrosianStudy
    29 # bin_PROGRAMS = psphot
     27bin_PROGRAMS = psphot psphotForced psphotMakePSF
     28# bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy
    3029
    3130psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    4140psphotMakePSF_LDADD = libpsphot.la
    4241
    43 psphotTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    44 psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    45 psphotTest_LDADD = libpsphot.la
    46 
    47 psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    48 psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    49 psphotMomentsStudy_LDADD = libpsphot.la
     42# psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     43# psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     44# psphotMomentsStudy_LDADD = libpsphot.la
    5045
    5146# psphotPetrosianStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    5247# psphotPetrosianStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    5348# psphotPetrosianStudy_LDADD = libpsphot.la
     49
     50# psphotTest_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     51# psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     52# psphotTest_LDADD = libpsphot.la
    5453
    5554# standard psphot for generic photometry
    … …  
    8281        psphotCleanup.c
    8382
    84 
    85 # psphot analysis of the detectability of specified positions
    86 psphotDetect_SOURCES =            \
    87         psphotDetect.c            \
    88         psphotDetectArguments.c   \
    89         psphotDetectParseCamera.c \
    90         psphotDetectImageLoop.c   \
    91         psphotDetectReadout.c     \
    92         psphotMosaicChip.c        \
    93         psphotCleanup.c
    94 
    95 psphotTest_SOURCES = \
    96         psphotTest.c
    97 
    98 psphotMomentsStudy_SOURCES = \
    99         psphotMomentsStudy.c
    100 
     83# # psphot analysis of the detectability of specified positions
     84# psphotDetect_SOURCES =            \
     85#         psphotDetect.c            \
     86#       psphotDetectArguments.c   \
     87#       psphotDetectParseCamera.c \
     88#       psphotDetectImageLoop.c   \
     89#       psphotDetectReadout.c     \
     90#       psphotMosaicChip.c        \
     91#       psphotCleanup.c
     92
     93# psphotTest_SOURCES = \
     94#         psphotTest.c
     95#
     96# psphotMomentsStudy_SOURCES = \
     97#         psphotMomentsStudy.c
     98#
    10199# psphotPetrosianStudy_SOURCES = \
    102100#         psphotPetrosianStudy.c
    … …  
    146144        psphotKernelFromPSF.c          \
    147145        psphotPSFConvModel.c           \
    148         psphotModelTest.c              \
    149146        psphotFitSet.c                 \
    150147        psphotSourceFreePixels.c       \
    … …  
    175172        psphotEfficiency.c
    176173
     174# XXX need to fix this for the new apis
     175#       psphotModelTest.c             
     176
    177177# re-instate these
    178178#       psphotIsophotal.c              \

  • branches/eam_branches/20091201/psphot/src/psphot.h

    r26635 r26788  
    1212
    1313#define PSPHOT_RECIPE_PSF_FAKE_ALLOW "PSF.FAKE.ALLOW" // Name for recipe component permitting fake PSFs
     14
     15// pmPCMData : PSF Convolved Model data storage structure
     16typedef struct {
     17    psImage *model;
     18    psArray *dmodels;
     19    psImage *modelConv;
     20    psArray *dmodelsConv;
     21} pmPCMData;
    1422
    1523// top-level psphot functions
    … …  
    2937bool            psphotReadoutMinimal(pmConfig *config, const pmFPAview *view);
    3038
    31 bool            psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources);
     39bool            psphotReadoutCleanup (pmConfig *config, const pmFPAview *view);
     40bool            psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     41
    3242bool            psphotDefineFiles (pmConfig *config, pmFPAfile *input);
    3343void            psphotFilesActivate(pmConfig *config, bool state);
    … …  
    3848// XXX test functions
    3949psArray        *psphotFakeSources (void);
    40 bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
     50bool            psphotMaskCosmicRayFootprintCheck (psArray *sources);
    4151
    4252// psphotReadout functions
     53bool            psphotAddPhotcode (pmConfig *config, const pmFPAview *view);
     54bool            psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
     55
     56bool            psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view);
     57bool            psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     58
    4359bool            psphotModelBackground (pmConfig *config, const pmFPAview *view);
    4460bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    4561
    46 // Create a background model for a readout, without saving the result as a pmFPAfile on
    47 // config->files.  Otherwise identical to psphotModelBackgroundFileIndex.
    48 psImage        *psphotModelBackgroundReadoutNoFile(pmReadout *readout, const pmConfig *config);
    49 
    50 // worker function for above background model functions
    51 bool psphotModelBackgroundReadout(psImage *model,  // Model image
    52                                   psImage *modelStdev, // Model stdev image
    53                                   psMetadata *analysis, // Analysis metadata for outputs
    54                                   pmReadout *readout, // Readout for which to generate a background model
    55                                   psImageBinning *binning, // Binning parameters
    56                                   const pmConfig *config // Configuration
    57     );
    58 
    59 psImageBinning *psphotBackgroundBinning(const psImage *, const pmConfig *);
    60 
    61 
    62 
    6362bool            psphotSubtractBackground (pmConfig *config, const pmFPAview *view);
    64 bool            psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
    65 
    66 
    67 pmDetections   *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe);
    68 
    69 bool            psphotRoughClass (pmReadout *readout, psArray *sources, psMetadata *recipe, const bool findPsfClump);
    70 bool            psphotBasicDeblend (psArray *sources, psMetadata *recipe);
    71 pmPSF          *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe);
    72 bool            psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf);
    73 bool            psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources);
    74 bool            psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final);
    75 bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
    76 
    77 bool            psphotReplaceAllSources (psArray *sources, psMetadata *recipe);
    78 bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
    79 
    80 bool            psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
     63bool            psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     64
     65bool            psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass);
     66bool            psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass);
     67
     68bool            psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow);
     69bool            psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow);
     70
     71bool            psphotDeblendSatstars (pmConfig *config, const pmFPAview *view);
     72bool            psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     73
     74bool            psphotBasicDeblend (pmConfig *config, const pmFPAview *view);
     75bool            psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     76
     77bool            psphotRoughClass (pmConfig *config, const pmFPAview *view);
     78bool            psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     79bool            psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *analysis, psMetadata *recipe, const bool havePSF);
     80
     81bool            psphotImageQuality (pmConfig *config, const pmFPAview *view);
     82bool            psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     83
     84bool            psphotChoosePSF (pmConfig *config, const pmFPAview *view);
     85bool            psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     86
     87bool            psphotGuessModels (pmConfig *config, const pmFPAview *view);
     88bool            psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     89
     90bool            psphotMergeSources (pmConfig *config, const pmFPAview *view);
     91bool            psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     92
     93bool            psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final);
     94bool            psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final);
     95
     96bool            psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize);
     97bool            psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize);
     98
     99bool            psphotBlendFit (pmConfig *config, const pmFPAview *view);
     100bool            psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    81101bool            psphotBlendFit_Threaded (psThreadJob *job);
    82102
    83 bool            psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout);
    84 psArray        *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow);
    85 bool            psphotSourceStats_Threaded (psThreadJob *job);
    86 
    87 bool            psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
    88 bool            psphotGuessModel_Threaded (psThreadJob *job);
    89 
    90 bool            psphotMagnitudes (pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf);
     103bool            psphotReplaceAllSources (pmConfig *config, const pmFPAview *view);
     104bool            psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     105
     106bool            psphotAddNoise (pmConfig *config, const pmFPAview *view);
     107bool            psphotSubNoise (pmConfig *config, const pmFPAview *view);
     108bool            psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add);
     109bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add);
     110
     111bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view);
     112bool            psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     113
     114bool            psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view);
     115bool            psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     116
     117bool            psphotApResid (pmConfig *config, const pmFPAview *view);
     118bool            psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
     119
     120bool            psphotMagnitudes (pmConfig *config, const pmFPAview *view);
     121bool            psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    91122bool            psphotMagnitudes_Threaded (psThreadJob *job);
     123
     124bool            psphotEfficiency (pmConfig *config, const pmFPAview *view);
     125bool            psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe);
    92126
    93127bool            psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources);
    94128bool            psphotPSFWeights_Threaded (psThreadJob *job);
    95129
    96 bool            psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf);
    97 
    98 bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
    99 bool            psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
    100 bool            psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe);
    101 bool            psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf, psMetadata *recipe, const psArray *realSources);
     130bool            psphotSkyReplace (pmConfig *config, const pmFPAview *view);
     131bool            psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     132
     133bool            psphotSourceFreePixels (pmConfig *config, const pmFPAview *view);
     134bool            psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index);
     135
     136// in psphotSourceStats.c:
     137bool            psphotSourceStats_Threaded (psThreadJob *job);
     138bool            psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout);
     139bool            psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources);
     140
     141// in psphotChoosePSF.c:
     142bool            psphotPSFstats (pmReadout *readout, pmPSF *psf);
     143bool            psphotMomentsStats (pmReadout *readout, psArray *sources);
     144
     145// in psphotGuessModel.c
     146bool            psphotGuessModel_Threaded (psThreadJob *job);
     147
     148// in psphotMergeSources.c:
     149bool            psphotLoadExtSources (pmConfig *config, const pmFPAview *view);
     150psArray        *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view);
     151bool            psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view);
     152bool            psphotCheckExtSources (pmConfig *config, const pmFPAview *view);
     153
     154// generate the detection structure for the supplied array of sources
     155bool            psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources);
     156
     157// generate the detection structure for the supplied array of sources
     158bool            psphotSetSourceParams (pmConfig *config, psArray *sources, pmPSF *psf);
     159
     160// in psphotModelBackground.c:
     161// Create a background model for a readout, without saving the result as a pmFPAfile on config->files.  Otherwise identical to psphotModelBackgroundFileIndex.
     162psImage        *psphotModelBackgroundReadoutNoFile(pmReadout *readout, const pmConfig *config);
     163psImageBinning *psphotBackgroundBinning(const psImage *image, const pmConfig *config);
     164
     165// in psphotReplaceUnfit.c:
     166bool            psphotRemoveAllSources (const psArray *sources, const psMetadata *recipe);
     167bool            psphotReplaceUnfitSources (psArray *sources, psImageMaskType maskVal);
    102168
    103169// thread-related:
    … …  
    113179psErrorCode     psphotCullPeaks(const psImage *img, const psImage *weight, const psMetadata *recipe, psArray *footprints);
    114180
    115 // generate the detection structure for the supplied array of sources
    116 pmDetections   *psphotDetectionsFromSources (pmConfig *config, psArray *sources);
    117 
    118 // used by ApResid
     181// in psphotApResid.c:
    119182pmTrend2D      *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag);
    120183bool            psphotApResidMags_Threaded (psThreadJob *job);
    … …  
    123186void            psphotModelClassInit (void);
    124187bool            psphotGrowthCurve (pmReadout *readout, pmPSF *psf, bool ignore, psImageMaskType maskVal);
    125 bool            psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
    126 bool            psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe);
    127 void            psphotSourceFreePixels (psArray *sources);
    128188
    129189// functions to set the correct source pixels
    130 bool            psphotInitRadiusPSF (const psMetadata *recipe, const pmModelType type);
     190bool            psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type);
     191
    131192bool            psphotCheckRadiusPSF (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
    132193bool            psphotCheckRadiusPSFBlend (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal, float dR);
    … …  
    134195bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
    135196float           psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal);
    136 
    137 // output functions
    138 bool            psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
    139 bool            psphotAddPhotcode (pmConfig *config, const pmFPAview *view);
    140197
    141198bool            psphotDumpMoments (psMetadata *recipe, psArray *sources);
    … …  
    169226bool            psphotFitSummary (void);
    170227
    171 bool            psphotMergeSources (psArray *oldSources, psArray *newSources);
    172 bool            psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources);
    173 psArray        *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view);
    174 
    175 pmPSF          *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe);
     228bool            psphotLoadPSF (pmConfig *config, const pmFPAview *view);
     229bool            psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index);
     230
    176231bool            psphotSetHeaderNstars (psMetadata *recipe, psArray *sources);
    177 bool            psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
    178 bool            psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe);
    179 bool            psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add);
    180232bool            psphotRadialPlot (int *kapa, const char *filename, pmSource *source);
    181233bool            psphotSourcePlots (pmReadout *readout, psArray *sources, psMetadata *recipe);
    182234bool            psphotMosaicSubimage (psImage *outImage, pmSource *source, int Xo, int Yo, int DX, int DY, bool normalize);
    183235
    184 bool            psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal);
    185 bool            psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal);
    186 bool            psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal);
    187 bool            psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe);
    188 bool            psphotSourceSize (pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first);
    189 
    190236bool            psphotMakeResiduals (psArray *sources, psMetadata *recipe, pmPSF *psf, psImageMaskType maskVal);
    191237
    … …  
    195241
    196242// functions related to extended source analysis
    197 bool  psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
    198 bool  psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax);
    199 float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
    200 bool  psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax);
    201 bool  psphotEllipticalProfile (pmSource *source);
    202 bool  psphotEllipticalContour (pmSource *source);
     243bool            psphotRadialProfile (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
     244bool            psphotRadialProfilesByAngles (pmSource *source, int Nsec, float Rmax);
     245float           psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
     246bool            psphotRadiiFromProfiles (pmSource *source, float fluxMin, float fluxMax);
     247bool            psphotEllipticalProfile (pmSource *source);
     248bool            psphotEllipticalContour (pmSource *source);
    203249
    204250// psphotVisual functions
    205 bool psphotVisualShowImage (pmReadout *readout);
    206 bool psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout);
    207 bool psphotVisualShowSignificance (psImage *image, float min, float max);
    208 bool psphotVisualShowPeaks (pmDetections *detections);
    209 bool psphotVisualShowFootprints (pmDetections *detections);
    210 bool psphotVisualShowMoments (psArray *sources);
    211 bool psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    212 bool psphotVisualShowRoughClass (psArray *sources);
    213 bool psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources);
    214 bool psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources);
    215 bool psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf);
    216 bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal);
    217 bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
    218 bool psphotVisualShowFlags (psArray *sources);
    219 bool psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
    220 bool psphotVisualShowResidualImage (pmReadout *readout);
    221 bool psphotVisualPlotApResid (psArray *sources, float mean, float error);
    222 bool psphotVisualPlotChisq (psArray *sources);
    223 bool psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    224 bool psphotVisualShowPetrosians (psArray *sources);
    225 
    226 // bool psphotPetrosianAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe);
    227 // bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise);
     251bool            psphotVisualShowImage (pmReadout *readout);
     252bool            psphotVisualShowBackground (pmConfig *config, const pmFPAview *view, pmReadout *readout);
     253bool            psphotVisualShowSignificance (psImage *image, float min, float max);
     254bool            psphotVisualShowPeaks (pmDetections *detections);
     255bool            psphotVisualShowFootprints (pmDetections *detections);
     256bool            psphotVisualShowMoments (psArray *sources);
     257bool            psphotVisualPlotMoments (psMetadata *recipe, psMetadata *analysis, psArray *sources);
     258bool            psphotVisualShowRoughClass (psArray *sources);
     259bool            psphotVisualShowPSFStars (psMetadata *recipe, pmPSF *psf, psArray *sources);
     260bool            psphotVisualShowSatStars (psMetadata *recipe, pmPSF *psf, psArray *sources);
     261bool            psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf);
     262bool            psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal);
     263bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
     264bool            psphotVisualShowFlags (psArray *sources);
     265bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
     266bool            psphotVisualShowResidualImage (pmReadout *readout);
     267bool            psphotVisualPlotApResid (psArray *sources, float mean, float error);
     268bool            psphotVisualPlotChisq (psArray *sources);
     269bool            psphotVisualPlotSourceSize (psMetadata *recipe, psMetadata *analysis, psArray *sources);
     270bool            psphotVisualShowPetrosians (psArray *sources);
     271bool            psphotVisualEraseOverlays (int channel, char *overlay);
    228272
    229273bool psphotPetrosian (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal);
    … …  
    231275bool psphotPetrosianStats (pmSource *source);
    232276
    233 // XXX old versions, currently disabled
     277// currently disabled:
    234278// bool            psphotIsophotal (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    235279// bool            psphotAnnuli (pmSource *source, psMetadata *recipe, psImageMaskType maskVal);
    … …  
    246290//                               float petFlux, float radiusForFlux);
    247291
    248 bool psphotImageQuality (psMetadata *recipe, psArray *sources);
    249 
    250292// structures & functions to support psf-convolved model fitting
    251 
    252 // pmPCMData : PSF Convolved Model data storage structure
    253 typedef struct {
    254     psImage *model;
    255     psArray *dmodels;
    256     psImage *modelConv;
    257     psArray *dmodelsConv;
    258 } pmPCMData;
    259 
    260293
    261294// psf-convolved model fitting
    … …  
    287320
    288321int psphotKapaOpen (void);
    289 bool psphotVisualEraseOverlays (int channel, char *overlay);
    290322bool psphotKapaClose (void);
    291323bool psphotImageBackgroundCellHistogram (psVector *values, float mean, float sigma, int ix, int iy);
    … …  
    301333int psphotSupplementStars (psArray *stars, psArray *sources, psImageBinning *binning, int ix, int iy);
    302334
    303 
    304335pmConfig *psphotForcedArguments(int argc, char **argv);
    305336bool psphotForcedImageLoop (pmConfig *config);

  • branches/eam_branches/20091201/psphot/src/psphotAddNoise.c

    r25755 r26788  
    11# include "psphotInternal.h"
    22
    3 bool psphotAddNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
    4   return psphotAddOrSubNoise (readout, sources, recipe, true);
     3bool psphotAddNoise (pmConfig *config, const pmFPAview *view) {
     4    return psphotAddOrSubNoise (config, view, true);
    55}
    66
    7 bool psphotSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe) {
    8   return psphotAddOrSubNoise (readout, sources, recipe, false);
     7bool psphotSubNoise (pmConfig *config, const pmFPAview *view) {
     8    return psphotAddOrSubNoise (config, view, false);
    99}
    1010
    11 bool psphotAddOrSubNoise (pmReadout *readout, const psArray *sources, const psMetadata *recipe, bool add) {
     11// for now, let's store the detections on the readout->analysis for each readout
     12bool psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, bool add)
     13{
     14    bool status = true;
     15
     16    // select the appropriate recipe information
     17    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     18    psAssert (recipe, "missing recipe?");
     19
     20    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     21    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     22
     23    // loop over the available readouts
     24    for (int i = 0; i < num; i++) {
     25        if (!psphotAddOrSubNoiseReadout (config, view, "PSPHOT.INPUT", i, recipe, add)) {
     26            psError (PSPHOT_ERR_CONFIG, false, "failed on to modify noise for PSPHOT.INPUT entry %d", i);
     27            return false;
     28        }
     29    }
     30    return true;
     31}
     32
     33bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool add) {
    1234
    1335    bool status = false;
    … …  
    1638    psEllipseAxes axes;
    1739
    18     PS_ASSERT (readout, false);
    19     PS_ASSERT (readout->parent, false);
    20     PS_ASSERT (readout->parent->concepts, false);
     40    // find the currently selected readout
     41    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     42    psAssert (file, "missing file?");
     43
     44    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     45    psAssert (readout, "missing readout?");
     46    psAssert (readout->parent, "missing cell?");
     47    psAssert (readout->parent->concepts, "missing concepts?");
     48
     49    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     50    psAssert (detections, "missing detections?");
     51
     52    psArray *sources = detections->allSources;
     53    psAssert (sources, "missing sources?");
    2154
    2255    psTimerStart ("psphot.noise");
    … …  
    2457    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    2558    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    26     assert (maskVal);
     59    psAssert (maskVal, "missing mask value?");
    2760
    2861    // increase variance by factor*(object noise):

  • branches/eam_branches/20091201/psphot/src/psphotApResid.c

    r26261 r26788  
    44// measure the aperture residual statistics and 2D variations
    55
    6 bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf)
     6// for now, let's store the detections on the readout->analysis for each readout
     7bool psphotApResid (pmConfig *config, const pmFPAview *view)
     8{
     9    bool status = true;
     10
     11    // select the appropriate recipe information
     12    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     13    psAssert (recipe, "missing recipe?");
     14
     15    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     16    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     17
     18    // loop over the available readouts
     19    for (int i = 0; i < num; i++) {
     20        if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     21            psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i);
     22            return false;
     23        }
     24    }
     25    return true;
     26}
     27
     28bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    729{
    830    int Nfail = 0;
    … …  
    1335    pmSource *source;
    1436
    15     PS_ASSERT_PTR_NON_NULL(config, false);
    16     PS_ASSERT_PTR_NON_NULL(readout, false);
    17     PS_ASSERT_PTR_NON_NULL(sources, false);
    18     PS_ASSERT_PTR_NON_NULL(psf, false);
    19 
    2037    psTimerStart ("psphot.apresid");
    2138
    22     // select the appropriate recipe information
    23     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    24     assert (recipe);
     39    // find the currently selected readout
     40    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     41    psAssert (file, "missing file?");
     42
     43    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     44    psAssert (readout, "missing readout?");
     45
     46    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     47    psAssert (detections, "missing detections?");
     48
     49    psArray *sources = detections->allSources;
     50    psAssert (sources, "missing sources?");
     51
     52    if (!sources->n) {
     53        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
     54        return true;
     55    }
     56
     57    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     58    psAssert (psf, "missing psf?");
    2559
    2660    // determine the number of allowed threads
    … …  
    3367    if (!measureAptrend) {
    3468        // save nan values since these were not calculated
    35         psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    36         psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    37         psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    38         psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     69        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     70        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     71        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     72        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    3973        return true;
    4074    }
    … …  
    288322
    289323    // save results for later output
    290     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
    291     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
    292     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
    293     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
     324    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
     325    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
     326    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
     327    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
    294328
    295329    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
    … …  
    308342escape:
    309343    // save nan values since these were not calculated
    310     psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
    311     psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
    312     psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
    313     psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
     344    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
     345    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
     346    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
     347    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
    314348
    315349    psFree (xPos);
    … …  
    318352    psFree (mag);
    319353    psFree (dMag);
    320     return false;
     354    return true;
     355    // this is a quality error, not a programming error
    321356}
    322357

  • branches/eam_branches/20091201/psphot/src/psphotBasicDeblend.c

    r20453 r26788  
    11# include "psphotInternal.h"
    22
    3 bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) {
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotBasicDeblend (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     9    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     10
     11    // loop over the available readouts
     12    for (int i = 0; i < num; i++) {
     13        if (!psphotBasicDeblendReadout (config, view, "PSPHOT.INPUT", i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on basic deblend analysis for PSPHOT.INPUT entry %d", i);
     15            return false;
     16        }
     17    }
     18    return true;
     19}
     20
     21bool psphotBasicDeblendReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) {
    422
    523    int N;
    … …  
    826    pmSource *source, *testSource;
    927
     28    psTimerStart ("psphot.deblend.basic");
     29
    1030    int Nblend = 0;
    1131
    12     psTimerStart ("psphot.deblend.basic");
     32    // find the currently selected readout
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest
     34    psAssert (file, "missing file?");
     35
     36    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     37    psAssert (readout, "missing readout?");
     38
     39    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     40    psAssert (detections, "missing detections?");
     41
     42    psArray *sources = detections->newSources;
     43    psAssert (sources, "missing sources?");
     44
     45    if (!sources->n) {
     46        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping basic deblend");
     47        return true;
     48    }
     49
     50    // select the appropriate recipe information
     51    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     52    psAssert (recipe, "missing recipe?");
    1353
    1454    float FRACTION = psMetadataLookupF32 (&status, recipe, "DEBLEND_PEAK_FRACTION");

  • branches/eam_branches/20091201/psphot/src/psphotBlendFit.c

    r25755 r26788  
    11# include "psphotInternal.h"
    22
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotBlendFit (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotBlendFitReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
    325// XXX I don't like this name
    4 bool psphotBlendFit (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
     26bool psphotBlendFitReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    527
    628    int Nfit = 0;
    … …  
    1234    psTimerStart ("psphot.fit.nonlinear");
    1335
    14     // select the appropriate recipe information
    15     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    16     assert (recipe);
     36    // find the currently selected readout
     37    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     38    psAssert (file, "missing file?");
     39
     40    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     41    psAssert (readout, "missing readout?");
     42
     43    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     44    psAssert (detections, "missing detections?");
     45
     46    psArray *sources = detections->allSources;
     47    psAssert (sources, "missing sources?");
     48
     49    if (!sources->n) {
     50        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
     51        return true;
     52    }
     53
     54    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     55    psAssert (psf, "missing psf?");
    1756
    1857    // determine the number of allowed threads
    … …  
    4079    psphotInitLimitsPSF (recipe, readout);
    4180    psphotInitLimitsEXT (recipe);
    42     psphotInitRadiusPSF (recipe, psf->type);
     81    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
    4382
    4483    // starts the timer, sets up the array of fitSets
    … …  
    4786    // source analysis is done in S/N order (brightest first)
    4887    sources = psArraySort (sources, pmSourceSortBySN);
     88    if (!sources->n) {
     89        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
     90        return true;
     91    }
    4992
    5093    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    … …  
    272315}
    273316
    274 # if (0)
    275 bool psphotBlendFit_Unthreaded (int *nfit, int *npsf, int *next, int *nfail, pmReadout *readout, psMetadata *recipe, psArray *sources, pmPSF *psf, psArray *newSources) {
    276 
    277     bool status = false;
    278     int Nfit = 0;
    279     int Npsf = 0;
    280     int Next = 0;
    281     int Nfail = 0;
    282 
    283     // bit-masks to test for good/bad pixels
    284     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    285     assert (maskVal);
    286 
    287     // bit-mask to mark pixels not used in analysis
    288     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    289     assert (markVal);
    290 
    291     // maskVal is used to test for rejected pixels, and must include markVal
    292     maskVal |= markVal;
    293 
    294     // S/N limit to perform full non-linear fits
    295     float FIT_SN_LIM = psMetadataLookupF32 (&status, recipe, "FULL_FIT_SN_LIM");
    296 
    297     // option to limit analysis to a specific region
    298     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
    299     psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
    300     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
    301 
    302     for (int i = 0; i < sources->n; i++) {
    303         pmSource *source = sources->data[i];
    304 
    305         // skip non-astronomical objects (very likely defects)
    306         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    307         if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
    308         if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    309         if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
    310 
    311         // skip DBL second sources (ie, added by psphotFitBlob)
    312         if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
    313 
    314         // limit selection to some SN limit
    315         if (source->peak->SN < FIT_SN_LIM) continue;
    316 
    317         // exclude sources outside optional analysis region
    318         if (source->peak->xf < AnalysisRegion.x0) continue;
    319         if (source->peak->yf < AnalysisRegion.y0) continue;
    320         if (source->peak->xf > AnalysisRegion.x1) continue;
    321         if (source->peak->yf > AnalysisRegion.y1) continue;
    322 
    323         // if model is NULL, we don't have a starting guess
    324         if (source->modelPSF == NULL) continue;
    325 
    326         // skip sources which are insignificant flux?
    327         // XXX this is somewhat ad-hoc
    328         if (source->modelPSF->params->data.F32[1] < 0.1) {
    329             psTrace ("psphot", 5, "skipping near-zero source: %f, %f : %f\n",
    330                      source->modelPSF->params->data.F32[1],
    331                      source->modelPSF->params->data.F32[2],
    332                      source->modelPSF->params->data.F32[3]);
    333             continue;
    334         }
    335 
    336         // replace object in image
    337         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    338             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    339         }
    340         Nfit ++;
    341 
    342         // try fitting PSFs or extended sources depending on source->mode
    343         // these functions subtract the resulting fitted source
    344         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
    345             if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
    346                 source->type = PM_SOURCE_TYPE_EXTENDED;
    347                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
    348                 Next ++;
    349                 continue;
    350             }
    351         } else {
    352             if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
    353                 source->type = PM_SOURCE_TYPE_STAR;
    354                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
    355                 Npsf ++;
    356                 continue;
    357             }
    358         }
    359 
    360         psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
    361         Nfail ++;
    362 
    363         // re-subtract the object, leave local sky
    364         pmSourceCacheModel (source, maskVal);
    365         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    366     }
    367 
    368     // change the value of a scalar on the array (wrap this and put it in psArray.h)
    369     *nfit  = Nfit;
    370     *npsf  = Npsf;
    371     *next  = Next;
    372     *nfail = Nfail;
    373 
    374     // moments are modified by the fit; re-display
    375     psphotVisualPlotMoments (recipe, sources);
    376     psphotVisualShowResidualImage (readout);
    377 
    378     return true;
    379 }
    380 # endif

  • branches/eam_branches/20091201/psphot/src/psphotChoosePSF.c

    r26262 r26788  
    11# include "psphotInternal.h"
    22
     3// generate a PSF model for inputs without PSF models already loaded
     4bool psphotChoosePSF (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
    325// try PSF models and select best option
    4 pmPSF *psphotChoosePSF (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     26bool psphotChoosePSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    527
    628    bool status;
    729
    830    psTimerStart ("psphot.choose.psf");
     31
     32    // find the currently selected readout
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     34    psAssert (file, "missing file?");
     35
     36    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     37    psAssert (readout, "missing readout?");
     38
     39    // do not generate a PSF if we already were supplied one
     40    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     41        psLogMsg ("psphot", PS_LOG_DETAIL, "psf model supplied for input file %d", index);
     42        return true;
     43    }
     44
     45    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     46    psAssert (detections, "missing detections?");
     47
     48    psArray *sources = detections->newSources;
     49    psAssert (sources, "missing sources?");
     50
     51    if (!sources->n) {
     52        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping PSF model");
     53        return true;
     54    }
    955
    1056    // bit-masks to test for good/bad pixels
    1157    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    12     assert (maskVal);
     58    psAssert (maskVal, "missing mask value?");
    1359
    1460    // bit-mask to mark pixels not used in analysis
    1561    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    16     assert (markVal);
     62    psAssert (markVal, "missing mark value?");
    1763
    1864    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    61107    // assert (status);
    62108
    63     // We have calculated a Gaussian window function, use that for both the PSF fit radius and
    64     // the aperture radius (scaling SIGMA)
    65     float gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
     109    // values on readout->analysis if we have calculated a Gaussian window function, use that
     110    // for both the PSF fit radius and the aperture radius (scaling SIGMA), otherwise base the value on the recipe value for MOMENTS_GAUSS_SIGMA
     111    float gaussSigma = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     112    if (!status) {
     113        gaussSigma = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     114    }
    66115    float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
    67116    float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
    … …  
    70119
    71120    // XXX use the same radii for standard analysis as for the PSF creation
    72     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
    73     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
     121    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_FIT_RADIUS", PS_META_REPLACE, "fit radius", options->fitRadius);
     122    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PSF_APERTURE", PS_META_REPLACE, "psf aperture", options->apRadius);
    74123
    75124    // dimensions of the field for which the PSF is defined
    … …  
    158207        bool status = true;
    159208        status &= psphotMakeFluxScale (readout->image, recipe, psf);
    160         status &= psphotPSFstats (readout, recipe, psf);
     209        status &= psphotPSFstats (readout, psf);
    161210        if (!status) {
    162211            psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF");
    … …  
    233282        bool status = true;
    234283        status &= psphotMakeFluxScale (readout->image, recipe, psf);
    235         status &= psphotPSFstats (readout, recipe, psf);
     284        status &= psphotPSFstats (readout, psf);
    236285        if (!status) {
    237286            psError(PSPHOT_ERR_PSF, false, "Failed to fit any models when choosing PSF");
    … …  
    290339    psFree (models);
    291340
    292     if (!psphotPSFstats (readout, recipe, psf)) {
     341    if (!psphotPSFstats (readout, psf)) {
    293342        psError(PSPHOT_ERR_PSF, false, "cannot measure PSF shape terms");
    294343        psFree(options);
    … …  
    302351
    303352    psFree (options);
    304     return (psf);
     353
     354    // save PSF on readout->analysis
     355    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) {
     356        psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
     357        return false;
     358    }
     359    psFree (psf); // XXX double-check
     360
     361    // move into psphotChoosePSF
     362    psphotVisualShowPSFModel (readout, psf);
     363
     364    return true;
    305365}
    306366
    307367// measure average parameters of the PSF model
    308 bool psphotPSFstats (pmReadout *readout, psMetadata *recipe, pmPSF *psf) {
     368bool psphotPSFstats (pmReadout *readout, pmPSF *psf) {
    309369
    310370    psEllipseShape shape;
    … …  
    312372
    313373    PS_ASSERT_PTR_NON_NULL(readout, false);
    314     PS_ASSERT_PTR_NON_NULL(recipe, false);
    315374    PS_ASSERT_PTR_NON_NULL(psf, false);
    316375
    … …  
    362421    }
    363422
    364     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
    365     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
    366     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
    367     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
     423    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", stats->sampleMean);
     424    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", stats->sampleStdev);
     425    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", stats->sampleLQ);
     426    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", stats->sampleUQ);
    368427
    369428    if (!psVectorStats (stats, fwhmMinor, NULL, NULL, 0)) {
    … …  
    371430        return false;
    372431    }
    373     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
    374     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
    375     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
    376     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
    377 
    378     psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
    379     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
    380     psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
     432    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", stats->sampleMean);
     433    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", stats->sampleStdev);
     434    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", stats->sampleLQ);
     435    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", stats->sampleUQ);
     436
     437    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           axes.theta);
     438    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", psf->nPSFstars);
     439    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", true);
    381440
    382441    psFree (fwhmMajor);
    … …  
    387446
    388447// determine approximate PSF shape parameters based on the moments
    389 bool psphotMomentsStats (pmReadout *readout, psMetadata *recipe, psArray *sources) {
     448bool psphotMomentsStats (pmReadout *readout, psArray *sources) {
    390449
    391450    // without the PSF model, we can only rely on the source->moments
    … …  
    401460
    402461    PS_ASSERT_PTR_NON_NULL(readout, false);
    403     PS_ASSERT_PTR_NON_NULL(recipe, false);
    404462    PS_ASSERT_PTR_NON_NULL(sources, false);
    405463
    … …  
    429487    FWHM_T /= FWHM_N;
    430488
    431     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X);
    432     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0);
    433     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0);
    434     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0);
    435 
    436     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y);
    437     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0);
    438     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0);
    439     psMetadataAddF32 (recipe, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0);
    440 
    441     psMetadataAddF32 (recipe, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           FWHM_T);
    442     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0);
    443     psMetadataAddBool(recipe, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false);
     489    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MAJ",   PS_META_REPLACE, "PSF FWHM Major axis (mean)", FWHM_X);
     490    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_SG",   PS_META_REPLACE, "PSF FWHM Major axis (sigma)", 0);
     491    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_LQ",   PS_META_REPLACE, "PSF FWHM Major axis (lower quartile)", 0);
     492    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MJ_UQ",   PS_META_REPLACE, "PSF FWHM Major axis (upper quartile)", 0);
     493    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FWHM_MIN",   PS_META_REPLACE, "PSF FWHM Minor axis (mean)", FWHM_Y);
     494    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_SG",   PS_META_REPLACE, "PSF FWHM Minor axis (sigma)", 0);
     495    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_LQ",   PS_META_REPLACE, "PSF FWHM Minor axis (lower quartile)", 0);
     496    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "FW_MN_UQ",   PS_META_REPLACE, "PSF FWHM Minor axis (upper quartile)", 0);
     497    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "ANGLE",    PS_META_REPLACE, "PSF angle",           FWHM_T);
     498    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "NPSFSTAR", PS_META_REPLACE, "Number of stars used to make PSF", 0);
     499    psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "PSFMODEL", PS_META_REPLACE, "Valid PSF Model?", false);
    444500
    445501    return true;

  • branches/eam_branches/20091201/psphot/src/psphotDeblendSatstars.c

    r25885 r26788  
    11# include "psphotInternal.h"
    22
    3 bool psphotDeblendSatstars (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotDeblendSatstars (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     9    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     10
     11    // loop over the available readouts
     12    for (int i = 0; i < num; i++) {
     13        if (!psphotDeblendSatstarsReadout (config, view, "PSPHOT.INPUT", i)) {
     14            psError (PSPHOT_ERR_CONFIG, false, "failed on saturated star deblend analysis for PSPHOT.INPUT entry %d", i);
     15            return false;
     16        }
     17    }
     18    return true;
     19}
     20
     21bool psphotDeblendSatstarsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int fileIndex) {
    422
    523    int N;
    624    pmSource *source;
     25    bool status;
    726
    827    psTimerStart ("psphot.deblend.sat");
    … …  
    1130    float SAT_MIN_RADIUS = 5.0;
    1231
    13     bool status;
     32    // find the currently selected readout
     33    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, fileIndex); // File of interest
     34    psAssert (file, "missing file?");
     35
     36    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     37    psAssert (readout, "missing readout?");
     38
     39    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     40    psAssert (detections, "missing detections?");
     41
     42    psArray *sources = detections->newSources;
     43    psAssert (sources, "missing sources?");
     44
     45    if (!sources->n) {
     46        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
     47        return true;
     48    }
     49
     50    // select the appropriate recipe information
     51    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     52    psAssert (recipe, "missing recipe?");
     53
    1454    pmCell *cell = readout->parent;
    1555    float SATURATION = 0.75*psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");

  • branches/eam_branches/20091201/psphot/src/psphotEfficiency.c

    r26705 r26788  
    66//#define TESTING
    77
     8
     9# if 0
    810
    911// Calculate the limiting magnitude for an image
    … …  
    146148}
    147149
     150# endif
     151
     152bool psphotEfficiency (pmConfig *config, const pmFPAview *view)
     153{
     154    bool status = true;
     155
     156    // select the appropriate recipe information
     157    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     158    psAssert (recipe, "missing recipe?");
     159
     160    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     161    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     162
     163    // loop over the available readouts
     164    for (int i = 0; i < num; i++) {
     165        if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     166            psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i);
     167            return false;
     168        }
     169    }
     170    return true;
     171}
    148172
    149173// Determine detection efficiency
    150 bool psphotEfficiency(pmConfig *config, pmReadout *readout, const pmFPAview *view, const pmPSF *psf,
    151                 psMetadata *recipe, const psArray *realSources)
     174bool psphotEfficiencyReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    152175{
     176# if 0
     177    bool status = true;
     178
     179    psTimerStart("psphot.fake");
     180
     181    // find the currently selected readout
     182    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     183    psAssert (file, "missing file?");
     184
     185    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     186    psAssert (readout, "missing readout?");
     187
     188    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     189    psAssert (detections, "missing detections?");
     190
     191    psArray *realSources = detections->allSources;
     192    psAssert (realSources, "missing sources?");
     193
     194    // XXX do we need to skip this step if we do not have a psf?
     195    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     196    psAssert (psf, "missing psf?");
     197
    153198    PM_ASSERT_READOUT_NON_NULL(readout, false);
    154199    PM_ASSERT_READOUT_IMAGE(readout, false);
    155200    PS_ASSERT_PTR_NON_NULL(psf, false);
    156201    PS_ASSERT_METADATA_NON_NULL(recipe, false);
    157     PS_ASSERT_ARRAY_NON_NULL(realSources, false);
    158 
    159     psTimerStart("psphot.fake");
    160202
    161203    // Collect recipe information
    … …  
    200242    // remove all sources, adding noise for subtracted sources
    201243    psphotRemoveAllSources(realSources, recipe);
    202     //    psphotAddNoise(readout, realSources, recipe);
     244    // psphotAddNoise(readout, realSources, recipe);
    203245
    204246
    … …  
    474516    psLogMsg("psphot", PS_LOG_INFO, "Detection efficiency: %lf sec\n", psTimerClear("psphot.fake"));
    475517
     518# endif
    476519
    477520    return true;

  • branches/eam_branches/20091201/psphot/src/psphotExtendedSourceAnalysis.c

    r25755 r26788  
    11# include "psphotInternal.h"
    22
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    // perform full non-linear fits / extended source analysis?
     13    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
     14        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
     15        return true;
     16    }
     17
     18    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     19    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     20
     21    // loop over the available readouts
     22    for (int i = 0; i < num; i++) {
     23        if (!psphotExtendedSourceAnalysisReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     24            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for PSPHOT.INPUT entry %d", i);
     25            return false;
     26        }
     27    }
     28    return true;
     29}
     30
    331// aperture-like measurements for extended sources
    4 bool psphotExtendedSourceAnalysis (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     32bool psphotExtendedSourceAnalysisReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    533
    634    bool status;
    … …  
    1139    int Nkron = 0;
    1240
     41    // find the currently selected readout
     42    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     43    psAssert (file, "missing file?");
     44
     45    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     46    psAssert (readout, "missing readout?");
     47
     48    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     49    psAssert (detections, "missing detections?");
     50
     51    psArray *sources = detections->allSources;
     52    psAssert (sources, "missing sources?");
     53
     54    if (!sources->n) {
     55        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     56        return true;
     57    }
     58
    1359    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    1460    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    1561    assert (maskVal);
    16 
    17     // perform full non-linear fits / extended source analysis?
    18     if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_ANALYSIS")) {
    19         psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
    20         return true;
    21     }
    2262
    2363    // XXX require petrosian analysis for non-linear fits?

  • branches/eam_branches/20091201/psphot/src/psphotExtendedSourceFits.c

    r25755 r26788  
    11# include "psphotInternal.h"
    22
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotExtendedSourceFits (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
     7
     8    // select the appropriate recipe information
     9    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     10    psAssert (recipe, "missing recipe?");
     11
     12    // perform full extended source non-linear fits?
     13    if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
     14        psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
     15        return true;
     16    }
     17
     18    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     19    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     20
     21    // loop over the available readouts
     22    for (int i = 0; i < num; i++) {
     23        if (!psphotExtendedSourceFitsReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     24            psError (PSPHOT_ERR_CONFIG, false, "failed on to fit extended sources for PSPHOT.INPUT entry %d", i);
     25            return false;
     26        }
     27    }
     28    return true;
     29}
     30
    331// non-linear model fitting for extended sources
    4 bool psphotExtendedSourceFits (pmReadout *readout, psArray *sources, psMetadata *recipe) {
     32bool psphotExtendedSourceFitsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    533
    634    bool status;
    … …  
    1240    bool savePics = false;
    1341
     42    // find the currently selected readout
     43    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     44    psAssert (file, "missing file?");
     45
     46    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     47    psAssert (readout, "missing readout?");
     48
     49    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     50    psAssert (detections, "missing detections?");
     51
     52    psArray *sources = detections->allSources;
     53    psAssert (sources, "missing sources?");
     54
     55    if (!sources->n) {
     56        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     57        return true;
     58    }
     59
    1460    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    1561    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    … …  
    2167    // maskVal is used to test for rejected pixels, and must include markVal
    2268    maskVal |= markVal;
    23 
    24     // perform full extended source non-linear fits?
    25     if (!psMetadataLookupBool (&status, recipe, "EXTENDED_SOURCE_FITS")) {
    26         psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source measurements\n");
    27         return true;
    28     }
    2969
    3070    // select the collection of desired models

  • branches/eam_branches/20091201/psphot/src/psphotFake.c

    r26705 r26788  
    211211
    212212    // Putting results on recipe because that appears to be the psphot standard, but it's not a good idea
    213     psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE,
    214                         "Efficiency fractions", frac);
    215     psMetadataAddVector(recipe, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE,
    216                         "Efficiency magnitudes", magOffsets);
    217     psMetadataAddF32(recipe, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE,
    218                      "Efficiency reference magnitude", magLim);
     213    psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "FAKE.EFF", PS_META_REPLACE, "Efficiency fractions", frac);
     214    psMetadataAddVector(readout->analysis, PS_LIST_TAIL, "FAKE.MAG", PS_META_REPLACE, "Efficiency magnitudes", magOffsets);
     215    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "FAKE.REF", PS_META_REPLACE, "Efficiency reference magnitude", magLim);
    219216
    220217    psFree(frac);

  • branches/eam_branches/20091201/psphot/src/psphotFindDetections.c

    r26634 r26788  
    11# include "psphotInternal.h"
    22
     3// we store the detections on the readout->analysis for each readout this function finds new
     4// peaks and new footprints.  any old peaks are saved on oldPeaks.  the resulting footprint set
     5// contains all footprints (old and new)
     6bool psphotFindDetections (pmConfig *config, const pmFPAview *view, bool firstPass)
     7{
     8    bool status = true;
     9
     10    // select the appropriate recipe information
     11    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     12    psAssert (recipe, "missing recipe?");
     13
     14    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     15    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     16
     17    // loop over the available readouts
     18    for (int i = 0; i < num; i++) {
     19        if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i, recipe, firstPass)) {
     20            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
     21            return false;
     22        }
     23    }
     24    return true;
     25}
     26
    327// smooth the image, search for peaks, optionally define footprints based on the peaks
    4 pmDetections *psphotFindDetections (pmDetections *detections, pmReadout *readout, psMetadata *recipe) {
     28bool psphotFindDetectionsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool firstPass) {
    529
    630    bool status;
    … …  
    933    int NMAX = 0;
    1034
     35    // find the currently selected readout
     36    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     37    psAssert (file, "missing file?");
     38
     39    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     40    psAssert (readout, "missing readout?");
     41
    1142    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    1243    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    13     assert (maskVal);
     44    psAssert (maskVal, "missing mask value?");
    1445
    1546    // Use the new pmFootprints approach?
    1647    const bool useFootprints = psMetadataLookupBool(NULL, recipe, "USE_FOOTPRINTS");
    1748
    18     // on first pass, detections have not yet been allocated
     49    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     50    // on the initial pass, detections have not yet been allocated or saved on readout->analysis
    1951    if (!detections) {
     52        // create the container
    2053        detections = pmDetectionsAlloc();
     54       
     55        // save detections on the readout->analysis
     56        if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) {
     57            psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
     58            return false;
     59        }
     60    } else {
     61        psMemIncrRefCounter(detections); // so we can free the detections below
     62    }
     63
     64    if (firstPass) {
    2165        pass = 1;
    2266        NSIGMA_PEAK = psMetadataLookupF32 (&status, recipe, "PEAKS_NSIGMA_LIMIT"); PS_ASSERT (status, NULL);
    … …  
    3074    float threshold = PS_SQR(NSIGMA_PEAK);
    3175
    32     // move the old peaks array (if it exists) to oldPeaks
    33     // XXX generically, we should be able to call this function an arbitrary number of times
     76    // move the old peaks array (if it exists) to oldPeaks XXX generically, we should be able
     77    // to call this function an arbitrary number of times the old peaks are saved so they can
     78    // be freed later -- the have to be freed after psphotFindFootprints is called below, since
     79    // they are also owned by the oldFootprints, which are in turn merged into the new
     80    // footprints.  (what about the source->peak entry?)
     81   
    3482    assert (detections->oldPeaks == NULL);
    3583    detections->oldPeaks = detections->peaks;
    … …  
    5199    // detect the peaks in the significance image
    52100    detections->peaks = psphotFindPeaks (significance, readout, recipe, threshold, NMAX);
    53     psMetadataAddF32  (recipe, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
     101    psMetadataAddF32  (readout->analysis, PS_LIST_TAIL, "PEAK_THRESHOLD", PS_META_REPLACE, "Peak Detection Threshold", threshold);
    54102    if (!detections->peaks) {
    55103        // we only get a NULL peaks array due to a programming or config error.
    … …  
    57105        psFree (detections);
    58106        psError (PSPHOT_ERR_CONFIG, false, "failed on peak search");
    59         return NULL;
     107        return false;
    60108    }
    61109
    … …  
    71119    psphotVisualShowFootprints (detections);
    72120
    73     return detections;
     121    psFree (detections);
     122
     123    return true;
    74124}
    75125
    76126// if we use the footprints, the output peaks list contains both old and new peaks,
    77127// otherwise it only contains the new peaks.
    78 
    79 # if 0
    80 // XXX where do we place the N sets of detections?
    81 bool psphotFindDetections (pmConfig *config, const pmFPAview *view)
    82 {
    83     bool status = true;
    84 
    85     int num = psMetadataAddS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    86     psAbort (!status, "programming error: must define PSPHOT.INPUT.NUM");
    87 
    88     // loop over the available readouts
    89     for (int i = 0; i < num; i++) {
    90         if (!psphotFindDetectionsReadout (config, view, "PSPHOT.INPUT", i)) {
    91             psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i);
    92             return false;
    93         }
    94     }
    95     return true;
    96 }
    97 # endif

  • branches/eam_branches/20091201/psphot/src/psphotFitSourcesLinear.c

    r25990 r26788  
    1212static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
    1313
    14 bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
     14// for now, let's store the detections on the readout->analysis for each readout
     15bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final)
     16{
     17    bool status = true;
     18
     19    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     20    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     21
     22    // loop over the available readouts
     23    for (int i = 0; i < num; i++) {
     24        if (!psphotFitSourcesLinearReadout (config, view, "PSPHOT.INPUT", i, final)) {
     25            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
     26            return false;
     27        }
     28    }
     29    return true;
     30}
     31
     32bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
    1533
    1634    bool status;
    … …  
    1937    float f;
    2038    // float r;
     39
     40    // select the appropriate recipe information
     41    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     42    assert (recipe);
     43
     44    // find the currently selected readout
     45    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     46    psAssert (file, "missing file?");
     47
     48    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     49    psAssert (readout, "missing readout?");
     50
     51    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     52    psAssert (detections, "missing detections?");
     53
     54    psArray *sources = detections->allSources;
     55    psAssert (sources, "missing sources?");
     56
     57    if (!sources->n) {
     58        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
     59        return true;
     60    }
     61
     62    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     63    psAssert (sources, "missing psf?");
    2164
    2265    psTimerStart ("psphot.linear");
    … …  
    248291    psphotVisualPlotChisq (sources);
    249292    // psphotVisualShowFlags (sources);
     293
     294    // We have to place this visualization here because the models are not realized until
     295    // psphotGuessModels or fitted until psphotFitSourcesLinear.
     296    psphotVisualShowPSFStars (recipe, psf, sources);
    250297
    251298    return true;

  • branches/eam_branches/20091201/psphot/src/psphotForcedReadout.c

    r26542 r26788  
    2525    }
    2626
    27     // find the currently selected readout
    28     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    29     PS_ASSERT_PTR_NON_NULL (readout, false);
    30 
    3127    // optional break-point for processing
    3228    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    … …  
    3430
    3531    // Generate the mask and weight images, including the user-defined analysis region of interest
    36     psphotSetMaskAndVariance (config, view, recipe);
     32    psphotSetMaskAndVariance (config, view);
    3733    if (!strcasecmp (breakPt, "NOTHING")) {
    38         return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
     34        return psphotReadoutCleanup(config, view);
    3935    }
    40 
    41     // display the image, weight, mask (ch 1,2,3)
    42     psphotVisualShowImage (readout);
    4336
    4437    // generate a background model (median, smoothed image)
    4538    if (!psphotModelBackground (config, view)) {
    46         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     39        return psphotReadoutCleanup (config, view);
    4740    }
    4841    if (!psphotSubtractBackground (config, view)) {
    49         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     42        return psphotReadoutCleanup (config, view);
    5043    }
    5144    if (!strcasecmp (breakPt, "BACKMDL")) {
    52         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     45        return psphotReadoutCleanup (config, view);
    5346    }
    5447
    55     // display the backsub and backgnd images
    56     psphotVisualShowBackground (config, view, readout);
    57 
    58     // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    59     pmPSF *psf = psphotLoadPSF (config, view, recipe);
    60     if (!psf) {
    61         psError(PSPHOT_ERR_CONFIG, false, "unable to load psf model");
    62         return false;
     48    if (!psphotLoadPSF (config, view)) {
     49        // this only happens if we had a programming error in psphotLoadPSF
     50        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
     51        return psphotReadoutCleanup (config, view);
    6352    }
    6453
    6554    // include externally-supplied sources
    66     psArray *sources = psArrayAllocEmpty(100);
    67     psphotLoadExtSources (config, view, sources);
     55    psphotLoadExtSources (config, view);
    6856
    6957    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    70     psphotGuessModels (config, readout, sources, psf);
     58    psphotGuessModels (config, view);
    7159
    7260    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    73     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
     61    psphotFitSourcesLinear (config, view, false);
    7462
    7563    // identify CRs and extended sources
    … …  
    7967
    8068    // calculate source magnitudes
    81     psphotMagnitudes(config, readout, view, sources, psf);
     69    psphotMagnitudes(config, view);
    8270
    8371    // XXX do I want to do this?
    … …  
    9179
    9280    // drop the references to the image pixels held by each source
    93     psphotSourceFreePixels (sources);
     81    psphotSourceFreePixels (config, view);
    9482
    9583    // create the exported-metadata and free local data
    96     return psphotReadoutCleanup(config, readout, recipe, NULL, psf, sources);
     84    return psphotReadoutCleanup(config, view);
    9785}

  • branches/eam_branches/20091201/psphot/src/psphotGuessModels.c

    r25755 r26788  
    77// 3) define the threaded function to work with sources for a given cell
    88
    9 // construct an initial PSF model for each object
    10 bool psphotGuessModels (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf) {
     9// for now, let's store the detections on the readout->analysis for each readout
     10bool psphotGuessModels (pmConfig *config, const pmFPAview *view)
     11{
     12    bool status = true;
     13
     14    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     15    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     16
     17    // loop over the available readouts
     18    for (int i = 0; i < num; i++) {
     19        if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) {
     20            psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i);
     21            return false;
     22        }
     23    }
     24    return true;
     25}
     26
     27// construct an initial PSF model for each object (new sources only)
     28bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    1129
    1230    bool status;
    1331
    1432    psTimerStart ("psphot.models");
     33
     34    // find the currently selected readout
     35    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     36    psAssert (file, "missing file?");
     37
     38    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     39    psAssert (readout, "missing readout?");
     40
     41    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     42    psAssert (detections, "missing detections?");
     43
     44    psArray *sources = detections->newSources;
     45    psAssert (sources, "missing sources?");
     46
     47    if (!sources->n) {
     48        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping model guess");
     49        return true;
     50    }
     51
     52    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     53    psAssert (sources, "missing psf?");
    1554
    1655    // select the appropriate recipe information
    … …  
    3675
    3776    // setup the PSF fit radius details
    38     psphotInitRadiusPSF (recipe, psf->type);
     77    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
    3978
    4079    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)

  • branches/eam_branches/20091201/psphot/src/psphotImageQuality.c

    r25755 r26788  
    44// XXX Lots of code duplication here
    55
     6// for now, let's store the detections on the readout->analysis for each readout
     7bool psphotImageQuality (pmConfig *config, const pmFPAview *view)
     8{
     9    bool status = true;
     10
     11    // select the appropriate recipe information
     12    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     13    psAssert (recipe, "missing recipe?");
     14
     15    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     16    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     17
     18    // loop over the available readouts
     19    for (int i = 0; i < num; i++) {
     20        if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     21            psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i);
     22            return false;
     23        }
     24    }
     25    return true;
     26}
     27
    628// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
    7 bool psphotImageQuality(psMetadata *recipe, psArray *sources)
     29bool psphotImageQualityReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    830{
     31    bool status = true;
     32
     33    // find the currently selected readout
     34    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     35    psAssert (file, "missing file?");
     36
     37    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     38    psAssert (readout, "missing readout?");
     39
     40    // if we have a PSF, skip this analysis
     41    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     42        return true;
     43    }
     44
     45    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     46    psAssert (detections, "missing detections?");
     47
     48    psArray *sources = detections->newSources;
     49    psAssert (sources, "missing sources?");
     50
     51    if (!sources->n) {
     52        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping image quality");
     53        return true;
     54    }
     55
    956    psVector *FWHM_MAJOR = psVectorAllocEmpty(100, PS_TYPE_F32);
    1057    psVector *FWHM_MINOR = psVectorAllocEmpty(100, PS_TYPE_F32);
    … …  
    81128
    82129    if (num == 0) {
    83       psWarning("Unable to find sources from which to measure image quality");
    84       return false;
    85     }
    86 
    87     psMetadataAddS32(recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE,
    88                      "Number of stars used for IQ measurements", M2->n);
     130        psLogMsg ("psphot", PS_LOG_INFO, "no valid sources for image quality, skipping");
     131        return true;
     132    }
     133
     134    psMetadataAddS32(readout->analysis, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, "Number of stars used for IQ measurements", M2->n);
    89135
    90136// XXX make this a recipe option
    … …  
    98144
    99145    float fwhm_major = stats->sampleMean;
    100     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE,
    101                      "FWHM of Major Axis from moments", stats->sampleMean);
    102     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE,
    103                      "FWHM scatter (Major) from moments", stats->sampleStdev);
     146    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE, "FWHM of Major Axis from moments", stats->sampleMean);
     147    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->sampleStdev);
    104148
    105149    if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) {
    … …  
    108152    }
    109153    float fwhm_minor = stats->sampleMean;
    110     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE,
    111                      "FWHM of Minor Axis from moments", stats->sampleMean);
    112     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE,
    113                      "FWHM scatter (Minor) from moments", stats->sampleStdev);
     154    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->sampleMean);
     155    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->sampleStdev);
    114156
    115157    if (!psVectorStats(stats, M2, NULL, NULL, 0)) {
    … …  
    119161    float vM2 = stats->sampleMean;
    120162    float dM2 = stats->sampleStdev;
    121     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE,
    122                      "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
    123     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE,
    124                      "Stdev of M_2", dM2);
    125     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE,
    126                      "Lower Quartile of M_2", stats->sampleLQ);
    127     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE,
    128                      "Upper Quartile of M_2", stats->sampleUQ);
     163    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
     164    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2);
     165    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->sampleLQ);
     166    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->sampleUQ);
    129167
    130168    if (!psVectorStats(stats, M2c, NULL, NULL, 0)) {
    … …  
    132170        goto FAIL;
    133171    }
    134     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE,
    135                      "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
    136     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE,
    137                      "Stdev of M_2c", stats->sampleStdev);
    138     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE,
    139                      "Lower Quartile of M_2c", stats->sampleLQ);
    140     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE,
    141                      "Upper Quartile of M_2c", stats->sampleUQ);
     172    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
     173    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->sampleStdev);
     174    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->sampleLQ);
     175    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->sampleUQ);
    142176
    143177    if (!psVectorStats(stats, M2s, NULL, NULL, 0)) {
    … …  
    145179        goto FAIL;
    146180    }
    147     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE,
    148                      "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
    149     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE,
    150                      "Stdev of M_2s", stats->sampleStdev);
    151     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE,
    152                      "Lower Quartile of M_2s", stats->sampleLQ);
    153     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE,
    154                      "Upper Quartile of M_2s", stats->sampleUQ);
     181    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
     182    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->sampleStdev);
     183    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->sampleLQ);
     184    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->sampleUQ);
    155185
    156186    if (!psVectorStats(stats, M3, NULL, NULL, 0)) {
    … …  
    160190    float vM3 = stats->sampleMean;
    161191    float dM3 = stats->sampleStdev;
    162     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE,
    163                      "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
    164     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE,
    165                      "Stdev of M_3", dM3);
    166     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE,
    167                      "Lower Quartile of M_3", stats->sampleLQ);
    168     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE,
    169                      "Upper Quartile of M_3", stats->sampleUQ);
     192    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
     193    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3);
     194    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->sampleLQ);
     195    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->sampleUQ);
    170196
    171197    if (!psVectorStats(stats, M4, NULL, NULL, 0)) {
    … …  
    175201    float vM4 = stats->sampleMean;
    176202    float dM4 = stats->sampleStdev;
    177     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE,
    178                      "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
    179     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE,
    180                      "Stdev of M_4", dM4);
    181     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE,
    182                      "Lower Quartile of M_4", stats->sampleLQ);
    183     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE,
    184                      "Upper Quartile of M_4", stats->sampleUQ);
     203    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
     204    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4);
     205    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->sampleLQ);
     206    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->sampleUQ);
    185207
    186208#else
    … …  
    192214    }
    193215    float fwhm_major = stats->robustMedian;
    194     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE,
    195                      "FWHM of Major Axis from moments", stats->robustMedian);
    196     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE,
    197                      "FWHM scatter (Major) from moments", stats->robustStdev);
     216    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE, "FWHM of Major Axis from moments", stats->robustMedian);
     217    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->robustStdev);
    198218
    199219    if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) {
    … …  
    202222    }
    203223    float fwhm_minor = stats->robustMedian;
    204     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE,
    205                      "FWHM of Minor Axis from moments", stats->robustMedian);
    206     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE,
    207                      "FWHM scatter (Minor) from moments", stats->robustStdev);
     224    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->robustMedian);
     225    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->robustStdev);
    208226
    209227    if (!psVectorStats(stats, M2, NULL, NULL, 0)) {
    … …  
    213231    float vM2 = stats->robustMedian;
    214232    float dM2 = stats->robustStdev;
    215     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE,
    216                      "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
    217     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE,
    218                      "Stdev of M_2", dM2);
    219     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE,
    220                      "Lower Quartile of M_2", stats->robustLQ);
    221     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE,
    222                      "Upper Quartile of M_2", stats->robustUQ);
     233    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
     234    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2);
     235    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->robustLQ);
     236    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->robustUQ);
    223237
    224238    if (!psVectorStats(stats, M2c, NULL, NULL, 0)) {
    … …  
    226240        goto FAIL;
    227241    }
    228     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE,
    229                      "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
    230     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE,
    231                      "Stdev of M_2c", stats->robustStdev);
    232     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE,
    233                      "Lower Quartile of M_2c", stats->robustLQ);
    234     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE,
    235                      "Upper Quartile of M_2c", stats->robustUQ);
     242    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
     243    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->robustStdev);
     244    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->robustLQ);
     245    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->robustUQ);
    236246
    237247    if (!psVectorStats(stats, M2s, NULL, NULL, 0)) {
    … …  
    239249        goto FAIL;
    240250    }
    241     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE,
    242                      "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
    243     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE,
    244                      "Stdev of M_2s", stats->robustStdev);
    245     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE,
    246                      "Lower Quartile of M_2s", stats->robustLQ);
    247     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE,
    248                      "Upper Quartile of M_2s", stats->robustUQ);
     251    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
     252    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->robustStdev);
     253    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->robustLQ);
     254    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->robustUQ);
    249255
    250256    if (!psVectorStats(stats, M3, NULL, NULL, 0)) {
    … …  
    254260    float vM3 = stats->robustMedian;
    255261    float dM3 = stats->robustStdev;
    256     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE,
    257                      "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
    258     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE,
    259                      "Stdev of M_3", dM3);
    260     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE,
    261                      "Lower Quartile of M_3", stats->robustLQ);
    262     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE,
    263                      "Upper Quartile of M_3", stats->robustUQ);
     262    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
     263    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3);
     264    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->robustLQ);
     265    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->robustUQ);
    264266
    265267    if (!psVectorStats(stats, M4, NULL, NULL, 0)) {
    … …  
    269271    float vM4 = stats->robustMedian;
    270272    float dM4 = stats->robustStdev;
    271     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE,
    272                      "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
    273     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE,
    274                      "Stdev of M_4", dM4);
    275     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE,
    276                      "Lower Quartile of M_4", stats->robustLQ);
    277     psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE,
    278                      "Upper Quartile of M_4", stats->robustUQ);
     273    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
     274    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4);
     275    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->robustLQ);
     276    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->robustUQ);
    279277#endif
    280278

  • branches/eam_branches/20091201/psphot/src/psphotLoadPSF.c

    r18002 r26788  
    11# include "psphotInternal.h"
    22
     3// NOTE : pmPSF_IO.c functions must load the psf model onto the chip->analysis metadata because
     4// the I/O operation likely occurs before the readout exists.  This implementation assumes that
     5// a single psf model is valid for the entire set of readouts (not valid for a time series of readouts)
     6
     7// XXX for now (2010.01.27), the supporting programs do not define multiple PSPHOT.PSF.LOAD
     8// files to go with multiple PSPHOT.INPUT files.  as a result, the implementation below is
     9// currently going to work for the case of a single input file, but will fail if we try with a
     10// stack of images.
     11
    312// load an externally supplied psf model
    4 pmPSF *psphotLoadPSF (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {
     13bool psphotLoadPSFReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     14
     15    bool status;
     16
     17    // find the currently selected readout
     18    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     19    psAssert (file, "missing file?");
    520
    621    // find the currently selected chip
    7     pmChip *chip = pmFPAfileThisChip (config->files, view, "PSPHOT.PSF.LOAD");
    8     if (!chip) return NULL;
     22    pmChip *chip = pmFPAviewThisChip (view, file->fpa);
     23    if (!chip) return false;
    924
    1025    // find the currently selected readout
    11     pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.PSF.LOAD");
    12     if (!readout) return NULL;
     26    pmReadout *readout = pmFPAviewThisReadout (view, file->fpa);
     27    if (!readout) return false;
    1328
    1429    // check if a PSF model is supplied by the user
    15     pmPSF *psf = psMetadataLookupPtr (NULL, chip->analysis, "PSPHOT.PSF");
     30    pmPSF *psf = psMetadataLookupPtr (&status, chip->analysis, "PSPHOT.PSF");
    1631    if (psf == NULL) {
    1732        psLogMsg ("psphot", 3, "no psf supplied for this chip");
    18         return NULL;
     33        return true;
    1934    }
    2035
    21     if (!psphotPSFstats (readout, recipe, psf)) psAbort("cannot measure PSF shape terms");
     36    if (!psphotPSFstats (readout, psf)) {
     37        psAbort("cannot measure PSF shape terms");
     38    }
    2239
    2340    psLogMsg ("psphot", 3, "using externally supplied PSF model for this readout");
    2441
    25     // we return a psf which can be free (and which may be removed from the metadata)
    26     psMemIncrRefCounter (psf);
    27     return psf;
     42    // save PSF on readout->analysis
     43    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot psf model", psf)) {
     44        psError (PSPHOT_ERR_UNKNOWN, false, "problem saving sources on readout");
     45        return false;
     46    }
     47
     48    return true;
    2849}
     50
     51bool psphotLoadPSF (pmConfig *config, const pmFPAview *view) {
     52
     53    bool status = false;
     54
     55    // XXX PSPHOT.PSF.LOAD vs PSPHOT.INPUT -- see note at top
     56    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     57    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     58
     59    // loop over the available readouts
     60    for (int i = 0; i < num; i++) {
     61
     62        // Generate the mask and weight images, including the user-defined analysis region of interest
     63        if (!psphotLoadPSFReadout (config, view, "PSPHOT.PSF.LOAD", i)) {
     64            psError (PSPHOT_ERR_CONFIG, false, "failed to load PSF model for PSPHOT.PSF.LOAD entry %d", i);
     65            return false;
     66        }
     67    }
     68    return true;
     69}

  • branches/eam_branches/20091201/psphot/src/psphotMagnitudes.c

    r25755 r26788  
    11# include "psphotInternal.h"
    22
    3 bool psphotMagnitudes(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources, const pmPSF *psf) {
     3bool psphotMagnitudes (pmConfig *config, const pmFPAview *view)
     4{
     5    bool status = true;
     6
     7    // select the appropriate recipe information
     8    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     9    psAssert (recipe, "missing recipe?");
     10
     11    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     12    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     13
     14    // loop over the available readouts
     15    for (int i = 0; i < num; i++) {
     16        if (!psphotMagnitudesReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     17            psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for PSPHOT.INPUT entry %d", i);
     18            return false;
     19        }
     20    }
     21    return true;
     22}
     23
     24bool psphotMagnitudesReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    425
    526    bool status = false;
    … …  
    829    psTimerStart ("psphot.mags");
    930
    10     // select the appropriate recipe information
    11     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    12     assert (recipe);
     31    // find the currently selected readout
     32    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     33    psAssert (file, "missing file?");
     34
     35    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     36    psAssert (readout, "missing readout?");
     37
     38    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     39    psAssert (detections, "missing detections?");
     40
     41    psArray *sources = detections->allSources;
     42    psAssert (sources, "missing sources?");
     43
     44    if (!sources->n) {
     45        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     46        return true;
     47    }
     48
     49    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     50    psAssert (psf, "missing psf?");
    1351
    1452    // determine the number of allowed threads
    … …  
    64102
    65103            psArrayAdd(job->args, 1, cells->data[j]); // sources
    66             psArrayAdd(job->args, 1, (pmPSF*)psf);    // Casting away const
     104            psArrayAdd(job->args, 1, psf);
    67105            psArrayAdd(job->args, 1, binning);
    68106            psArrayAdd(job->args, 1, backModel);
    … …  
    179217}
    180218
    181 # if (0)
    182 bool psphotMagnitudes_Unthreaded (int *nap, psArray *sources, pmPSF *psf, psImageBinning *binning, pmReadout *backModel, pmReadout *backStdev, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
    183 
    184     bool status;
    185     int Nap = 0;
    186 
    187     for (int i = 0; i < sources->n; i++) {
    188         pmSource *source = (pmSource *) sources->data[i];
    189         status = pmSourceMagnitudes (source, psf, photMode, maskVal);
    190         if (status && isfinite(source->apMag)) Nap ++;
    191 
    192         if (backModel) {
    193             psAssert (binning, "if backModel is defined, so should binning be");
    194             source->sky = psImageUnbinPixel(source->peak->x, source->peak->y, backModel->image, binning);
    195             if (isnan(source->sky) && false) {
    196                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.sky");
    197             }
    198         } else {
    199             source->sky = NAN;
    200         }
    201 
    202         if (backStdev) {
    203             psAssert (binning, "if backStdev is defined, so should binning be");
    204             source->skyErr = psImageUnbinPixel(source->peak->x, source->peak->y, backStdev->image, binning);
    205             if (isnan(source->skyErr) && false) {
    206                 psLogMsg ("psphot.magnitudes", PS_LOG_DETAIL, "error setting pmSource.skyErr");
    207             }
    208         } else {
    209             source->skyErr = NAN;
    210         }
    211     }
    212 
    213     // change the value of a scalar on the array (wrap this and put it in psArray.h)
    214     *nap = Nap;
    215 
    216     return true;
    217 }
    218 # endif
    219 
    220219bool psphotPSFWeights(pmConfig *config, pmReadout *readout, const pmFPAview *view, psArray *sources) {
    221220

  • branches/eam_branches/20091201/psphot/src/psphotMakePSFReadout.c

    r26542 r26788  
    2424    }
    2525
    26     // find the currently selected readout
    27     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    28     PS_ASSERT_PTR_NON_NULL (readout, false);
    29 
    3026    // optional break-point for processing
    3127    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    … …  
    3329
    3430    // Generate the mask and weight images, including the user-defined analysis region of interest
    35     psphotSetMaskAndVariance (config, view, recipe);
     31    psphotSetMaskAndVariance (config, view);
    3632    if (!strcasecmp (breakPt, "NOTHING")) {
    37         return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
     33        return psphotReadoutCleanup(config, view);
    3834    }
    39 
    40     // display the image, weight, mask (ch 1,2,3)
    41     psphotVisualShowImage (readout);
    4235
    4336    // generate a background model (median, smoothed image)
    4437    if (!psphotModelBackground (config, view)) {
    45         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     38        return psphotReadoutCleanup (config, view);
    4639    }
    4740    if (!psphotSubtractBackground (config, view)) {
    48         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     41        return psphotReadoutCleanup (config, view);
    4942    }
    5043    if (!strcasecmp (breakPt, "BACKMDL")) {
    51         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     44        return psphotReadoutCleanup (config, view);
    5245    }
    5346
    54     // display the backsub and backgnd images
    55     psphotVisualShowBackground (config, view, readout);
     47    psphotLoadExtSources (config, view);
    5648
    57     pmDetections *detections = NULL;
    58 
    59     // If sources have been supplied, then these should be used to measure the PSF
    60     // include externally-supplied sources
    61 
    62     // XXX sources loaded from a text file have no valid SN values, but psphotChoosePSF
    63     // selected the top PSF_MAX_NSTARS to generate the PSF, excluding an arbitrary subset.
    64     psArray *sources = psArrayAllocEmpty(100);
    65     psphotLoadExtSources (config, view, sources);
    66     if (sources->n) {
    67         // the user wants to make the psf from these stars; define them as psf stars:
    68         for (int i = 0; i < sources->n; i++) {
    69             pmSource *source = sources->data[i];
    70             source->mode |= PM_SOURCE_MODE_PSFSTAR;
    71         }
    72         // force psphotChoosePSF to use all loaded sources
    73         psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF_MAX_NSTARS", PS_META_REPLACE, "fit radius", sources->n);
    74 
    75         // measure stats of externally specified sources
    76         if (!psphotSourceStatsUpdate (sources, config, readout)) {
    77             psError(PSPHOT_ERR_CONFIG, false, "failure to measure stats of existing sources");
    78             return false;
    79         }
    80 
    81     } else {
    82         // find the detections (by peak and/or footprint) in the image.
    83         detections = psphotFindDetections (NULL, readout, recipe);
    84         if (!detections) {
    85             psLogMsg ("psphot", 3, "unable to find detections in this image");
    86             return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources);
    87         }
    88 
    89         // construct sources and measure basic stats
    90         psFree (sources);
    91         sources = psphotSourceStats (config, readout, detections, true);
    92         if (!sources) return false;
    93 
    94         // find blended neighbors of very saturated stars
    95         // XXX merge this with Basic Deblend?
    96         psphotDeblendSatstars (readout, sources, recipe);
    97 
    98         // mark blended peaks PS_SOURCE_BLEND
    99         if (!psphotBasicDeblend (sources, recipe)) {
    100             psLogMsg ("psphot", 3, "failed on deblend analysis");
    101             return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources);
    102         }
    103 
    104         // classify sources based on moments, brightness
    105         if (!psphotRoughClass (readout, sources, recipe, false)) {
    106             psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    107             return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources);
    108         }
     49    // If sources have been supplied, then these should be used to measure the PSF include
     50    // externally-supplied sources; if not, we need to generate a set of possible PSF sources.
     51    // This function updates the SN entries for the loaded sources or generates a set of
     52    // detections from the image, if no external ones have been supplied.  Sources loaded from
     53    // a text file have no valid SN values, but psphotChoosePSF needs to select the top
     54    // PSF_MAX_NSTARS to generate the PSF.
     55    if (!psphotCheckExtSources (config, view)) {
     56        psLogMsg ("psphot", 3, "failure to select possible PSF sources (external or internal)");
     57        return psphotReadoutCleanup (config, view);
    10958    }
    11059
    111     // use bright stellar objects to measure PSF
    112     // XXX if we do not have enough stars to generate the PSF, build one
    113     // from the SEEING guess and model class
    114     pmPSF *psf = psphotChoosePSF (readout, sources, recipe);
    115     if (psf == NULL) {
     60    // Use bright stellar objects to measure PSF. If we do not have enough stars to generate
     61    // the PSF, build one from the SEEING guess and model class
     62    if (!psphotChoosePSF (config, view)) {
    11663        psLogMsg ("psphot", 3, "failure to construct a psf model");
    117         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     64        return psphotReadoutCleanup (config, view);
    11865    }
    119     psphotVisualShowPSFModel (readout, psf);
    120     return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    121 }
    12266
    12367    // measure aperture photometry corrections
    124     // XXX isn't this part of the PSF??
    125     // if (!psphotApResid (config, readout, sources, psf)) {
    126     //     psLogMsg ("psphot", 3, "failed on psphotApResid");
    127     //     return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    128     // }
     68# if 0
     69    if (!psphotApResid (config, view)) {
     70        psLogMsg ("psphot", 3, "failed on psphotApResid");
     71        return psphotReadoutCleanup (config, view);
     72    }
     73# endif
     74
     75    return psphotReadoutCleanup (config, view);
     76}

  • branches/eam_branches/20091201/psphot/src/psphotMaskReadout.c

    r26542 r26788  
    11# include "psphotInternal.h"
    22
     3bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view) {
     4
     5    bool status = false;
     6
     7    // select the appropriate recipe information
     8    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     9    psAssert (recipe, "missing recipe?");
     10
     11    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     12    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     13
     14    // loop over the available readouts
     15    for (int i = 0; i < num; i++) {
     16
     17        // Generate the mask and weight images, including the user-defined analysis region of interest
     18        if (!psphotSetMaskAndVarianceReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     19            psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i);
     20            return false;
     21        }
     22    }
     23    return true;
     24}
     25
    326// generate mask and variance if not defined, additional mask for restricted subregion
    4 bool psphotSetMaskAndVarianceReadout (pmConfig *config, pmReadout *readout, psMetadata *recipe) {
     27bool psphotSetMaskAndVarianceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    528
    629    bool status;
    730
    8     // ** Interpret the mask values:
    9     // XXX drop the write to recipe and move config into psphotRoughClass?
    10     // XXX alternatively, define a function to set the psphot recipe masks
     31    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest
     32    psAssert (file, "missing file?");
     33
     34    // find the currently selected readout
     35    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     36    psAssert (readout, "missing readout?");
     37
     38    // save maskSat and maskBad on the psphot recipe (mostly for psphotRoughClass)
    1139    psImageMaskType maskSat  = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
    1240    psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.SAT", PS_META_REPLACE, "user-defined mask", maskSat);
    … …  
    1442    psImageMaskType maskBad  = pmConfigMaskGet("LOW", config); // Mask value for low pixels
    1543    if (!maskBad) {
    16         // XXX: for backward compatability look up old name
     44        // for backward compatability look up old name
    1745        maskBad  = pmConfigMaskGet("BAD", config);
    1846    }
    … …  
    73101    }
    74102
     103    // display the image, weight, mask (ch 1,2,3)
     104    psphotVisualShowImage (readout);
     105
    75106    return true;
    76107}
    77 
    78 bool psphotSetMaskAndVariance (pmConfig *config, const pmFPAview *view, psMetadata *recipe) {
    79 
    80     bool status = false;
    81 
    82     int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    83     psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    84 
    85     // loop over the available readouts
    86     for (int i = 0; i < num; i++) {
    87 
    88         pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest
    89         PS_ASSERT_PTR_NON_NULL (file, false);
    90 
    91         // find the currently selected readout
    92         pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
    93         PS_ASSERT_PTR_NON_NULL (readout, false);
    94 
    95         // Generate the mask and weight images, including the user-defined analysis region of interest
    96         if (!psphotSetMaskAndVarianceReadout (config, readout, recipe)) {
    97             psError (PSPHOT_ERR_CONFIG, false, "failed to generate mask for PSPHOT.INPUT entry %d", i);
    98             return false;
    99         }
    100 
    101         // display the image, weight, mask (ch 1,2,3)
    102         psphotVisualShowImage (readout);
    103     }
    104     return true;
    105 }

  • branches/eam_branches/20091201/psphot/src/psphotMergeSources.c

    r26645 r26788  
    55                         PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply for PSF sources
    66
    7 // merge the externally supplied sources with the existing sources.  mark them as having
    8 // mode PM_SOURCE_MODE_EXTERNAL
    9 bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view, psArray *sources) {
    10 
    11     psArray *extSourcesCMF = NULL;
     7// for now, let's store the detections on the readout->analysis for each readout
     8bool psphotMergeSources (pmConfig *config, const pmFPAview *view)
     9{
     10    bool status = true;
     11
     12    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     13    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     14
     15    // loop over the available readouts
     16    for (int i = 0; i < num; i++) {
     17        if (!psphotMergeSourcesReadout (config, view, "PSPHOT.INPUT", i)) {
     18            psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i);
     19            return false;
     20        }
     21    }
     22    return true;
     23}
     24
     25// add newly selected sources to the existing list of sources
     26bool psphotMergeSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     27
     28    bool status;
     29
     30    // find the currently selected readout
     31    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     32    psAssert (file, "missing file?");
     33
     34    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     35    psAssert (readout, "missing readout?");
     36
     37    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     38    psAssert (detections, "missing detections?");
     39
     40    psArray *newSources = detections->newSources;
     41    psAssert (newSources, "missing sources?");
     42
     43    // XXX TEST:
     44    if (detections->allSources) {
     45        psphotMaskCosmicRayFootprintCheck(detections->allSources);
     46    }
     47    if (detections->newSources) {
     48        psphotMaskCosmicRayFootprintCheck(detections->newSources);
     49    }
     50
     51    if (!detections->allSources) {
     52        detections->allSources = psArrayAllocEmpty(newSources->n);
     53    }
     54    psArray *allSources = detections->allSources;
     55
     56    for (int i = 0; i < newSources->n; i++) {
     57        pmSource *source = newSources->data[i];
     58        psArrayAdd (allSources, 100, source);
     59    }
     60
     61    psFree (detections->newSources);
     62    detections->newSources = NULL;
     63
     64    return true;
     65}
     66
     67// Merge the externally supplied sources with the existing sources.  Mark them as having mode
     68// PM_SOURCE_MODE_EXTERNAL.
     69
     70// XXX This function needs to be updated to loop over set of input files.  At the moment, we
     71// only expect a single entry for PSPHOT.INPUT.CMF and PSPHOT.SOURCES.TEXT, so we can only
     72// associate input sources with a single entry for PSPHOT.INPUT
     73bool psphotLoadExtSources (pmConfig *config, const pmFPAview *view) {
     74
     75    bool status;
     76    pmDetections *extCMF = NULL;
    1277    psArray *extSourcesTXT = NULL;
     78
     79    // find the currently selected readout
     80    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     81    psAssert (file, "missing file?");
     82
     83    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     84    psAssert (readout, "missing readout?");
     85
     86    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     87    psAssert (detections, "missing detections?");
     88
     89    // XXX allSources of newSources?
     90    psArray *sources = detections->allSources;
     91    psAssert (sources, "missing sources?");
    1392
    1493    // load data from input CMF file:
    … …  
    1796        if (!readoutCMF) goto loadTXT;
    1897
    19         extSourcesCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.SOURCES");
    20         if (extSourcesCMF) {
    21             for (int i = 0; i < extSourcesCMF->n; i++) {
    22                 pmSource *source = extSourcesCMF->data[i];
    23                 source->mode |= PM_SOURCE_MODE_EXTERNAL;
     98        extCMF = psMetadataLookupPtr (NULL, readoutCMF->analysis, "PSPHOT.DETECTIONS");
     99        if (extCMF) {
     100            for (int i = 0; i < extCMF->allSources->n; i++) {
     101                pmSource *source = extCMF->allSources->data[i];
     102                source->mode |= PM_SOURCE_MODE_EXTERNAL;
    24103
    25104                // the supplied peak flux needs to be re-normalized
    … …  
    27106                source->peak->value = 1.0;
    28107
    29                 // drop the loaded source modelPSF
    30                 psFree (source->modelPSF);
    31                 source->modelPSF = NULL;
    32             }
    33             psphotMergeSources (sources, extSourcesCMF);
    34         }
     108                // drop the loaded source modelPSF
     109                psFree (source->modelPSF);
     110                source->modelPSF = NULL;
     111
     112                psArrayAdd (detections->allSources, 100, source);
     113            }
     114        }
    35115    }
    36116
    … …  
    48128                source->mode |= PM_SOURCE_MODE_EXTERNAL;
    49129
    50                 // drop the loaded source modelPSF
    51                 psFree (source->modelPSF);
    52                 source->modelPSF = NULL;
    53             }
    54             psphotMergeSources (sources, extSourcesTXT);
    55         }
     130                // drop the loaded source modelPSF
     131                psFree (source->modelPSF);
     132                source->modelPSF = NULL;
     133
     134                psArrayAdd (detections->allSources, 100, source);
     135            }
     136        }
    56137    }
    57138
    58139finish:
    59140
    60     if (!extSourcesTXT && !extSourcesTXT) {
     141    if (!extCMF && !extSourcesTXT) {
    61142        psLogMsg ("psphot", 3, "no external sources for this readout");
    62143        return true;
    63144    }
    64145
    65     int nCMF = extSourcesCMF ? extSourcesCMF->n : 0;
     146    int nCMF = extCMF        ? extCMF->allSources->n        : 0;
    66147    int nTXT = extSourcesTXT ? extSourcesTXT->n : 0;
    67148
    … …  
    71152}
    72153
    73 // add newly selected sources to the existing list of sources
    74 bool psphotMergeSources (psArray *oldSources, psArray *newSources) {
    75 
    76     for (int i = 0; i < newSources->n; i++) {
    77         pmSource *source = newSources->data[i];
    78         psArrayAdd (oldSources, 100, source);
    79     }
    80     return true;
    81 }
    82 
    83154// extract the input sources corresponding to this readout
     155// XXX this function needs to be updated to work with the new context of pshot inputs
    84156psArray *psphotLoadPSFSources (pmConfig *config, const pmFPAview *view) {
     157
     158    bool status;
    85159
    86160    // find the currently selected readout
    … …  
    91165    }
    92166
    93     psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES");
     167    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     168    if (!detections) {
     169        psLogMsg ("psphot", 3, "no psf sources for this readout");
     170        return NULL;
     171    }
     172
     173    psArray *sources = detections->allSources;
    94174    if (!sources) {
    95175        psLogMsg ("psphot", 3, "no psf sources for this readout");
     176        return NULL;
    96177    }
    97178
    … …  
    99180}
    100181
     182// this function is used to fix sources which were loaded externally, but have passed from
     183// psphotDetectionsFromSources to psphotSourceStats and are now stored on
     184// detections->newSources.
     185bool psphotRepairLoadedSources (pmConfig *config, const pmFPAview *view) {
     186
     187    // find the currently selected readout
     188    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     189    psAssert (file, "missing file?");
     190
     191    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     192    psAssert (readout, "missing readout?");
     193
     194    pmDetections *detections = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.DETECTIONS");
     195    if (!detections) {
     196        psError(PSPHOT_ERR_CONFIG, false, "missing detections");
     197        return false;
     198    }
     199
     200    psArray *sources = detections->newSources;
     201    psAssert (sources, "missing sources?");
     202
     203    // peak flux is wrong : set based on previous image
     204    // use the peak measured in the moments analysis:
     205    for (int i = 0; i < sources->n; i++) {
     206      pmSource *source = sources->data[i];
     207      source->peak->flux = source->moments->Peak;
     208    }
     209
     210    return true;
     211}
     212
    101213// generate the detection structure for the supplied array of sources
    102 pmDetections *psphotDetectionsFromSources (pmConfig *config, psArray *sources) {
     214// XXX this currently assumes there is a single input file
     215bool psphotDetectionsFromSources (pmConfig *config, const pmFPAview *view, psArray *sources) {
     216
     217    // find the currently selected readout
     218    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     219    psAssert (file, "missing file?");
     220
     221    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     222    psAssert (readout, "missing readout?");
    103223
    104224    pmDetections *detections = pmDetectionsAlloc();
    … …  
    113233    float snMin = psMetadataLookupF32(NULL, recipe, "MOMENTS_SN_MIN");
    114234    if (!isfinite(snMin)) {
    115         return NULL;
     235        return false;
    116236    }
    117237
    … …  
    152272    psLogMsg ("psphot", 3, "%ld PSF sources loaded", detections->peaks->n);
    153273
    154     return detections;
     274    // save detections on the readout->analysis
     275    if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detectinos", detections)) {
     276        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
     277        return false;
     278    }
     279    psFree (detections);
     280    return true;
    155281}
    156282
    … …  
    194320    return true;
    195321}
     322
     323bool psphotCheckExtSources (pmConfig *config, const pmFPAview *view) {
     324
     325    bool status;
     326
     327    psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
     328    psAssert (recipe, "missing recipe");
     329
     330    // find the currently selected readout
     331    pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest
     332    psAssert (file, "missing file?");
     333
     334    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     335    psAssert (readout, "missing readout?");
     336
     337    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     338    psAssert (detections, "missing detections?");
     339
     340    // XXX allSources of newSources?
     341    psArray *sources = detections->allSources;
     342    psAssert (sources, "missing sources?");
     343
     344    if (sources->n) {
     345        // the user wants to make the psf from these stars; define them as psf stars:
     346        for (int i = 0; i < sources->n; i++) {
     347            pmSource *source = sources->data[i];
     348            source->mode |= PM_SOURCE_MODE_PSFSTAR;
     349        }
     350        // force psphotChoosePSF to use all loaded sources
     351        psMetadataAddS32 (recipe, PS_LIST_TAIL, "PSF_MAX_NSTARS", PS_META_REPLACE, "max number of sources for PSF model", sources->n);
     352
     353        // measure stats of externally specified sources
     354        if (!psphotSourceStatsUpdate (sources, config, readout)) {
     355            psError(PSPHOT_ERR_CONFIG, false, "failure to measure stats of existing sources");
     356            return false;
     357        }
     358    } else {
     359
     360        // find the detections (by peak and/or footprint) in the image.
     361        if (!psphotFindDetections (config, view, true)) {
     362            psError(PSPHOT_ERR_CONFIG, false, "unable to find detections in this image");
     363            return psphotReadoutCleanup (config, view);
     364        }
     365
     366        // construct sources and measure basic stats
     367        psphotSourceStats (config, view, true);
     368
     369        // find blended neighbors of very saturated stars
     370        psphotDeblendSatstars (config, view);
     371
     372        // mark blended peaks PS_SOURCE_BLEND
     373        if (!psphotBasicDeblend (config, view)) {
     374            psLogMsg ("psphot", 3, "failed on deblend analysis");
     375            return psphotReadoutCleanup (config, view);
     376        }
     377
     378        // classify sources based on moments, brightness
     379        if (!psphotRoughClass (config, view)) {
     380            psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     381            return psphotReadoutCleanup (config, view);
     382        }
     383    }
     384
     385    return true;
     386}

  • branches/eam_branches/20091201/psphot/src/psphotModelBackground.c

    r26542 r26788  
    3232// generate the median in NxN boxes, clipping heavily
    3333// linear interpolation to generate full-scale model
    34 bool psphotModelBackgroundReadout(psImage *model,  // Model image
     34//
     35// NOTE that the 'analysis' metedata pass in here is used to store the binning information.
     36// This may be the analysis for this readout, but it may be the analysis for the pmFPAfile
     37// corresponding to the model.  Other information about the background model is saved on the
     38// readout->analysis
     39static bool psphotModelBackgroundReadout(psImage *model,  // Model image
    3540                                  psImage *modelStdev, // Model stdev image
    3641                                  psMetadata *analysis, // Analysis metadata for outputs
    … …  
    140145
    141146    // we save the binning structure for use in psphotMagnitudes
    142     psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING",
    143                      PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
     147    psMetadataAddPtr(analysis, PS_LIST_TAIL, "PSPHOT.BACKGROUND.BINNING", PS_DATA_UNKNOWN | PS_META_REPLACE, "Background binning", binning);
    144148
    145149    psF32 **modelData = model->data.F32;
    … …  
    296300    psLogMsg ("psphot", PS_LOG_INFO, "built background image: %f sec\n", psTimerMark ("psphot.background"));
    297301
    298     psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
    299     psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     302    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_MEAN", PS_META_REPLACE, "sky mean", Value);
     303    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
    300304    psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
    301305
    … …  
    306310                                      PS_STAT_MAX);
    307311    psImageStats (statsBck, model, NULL, 0);
    308     psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MN",
    309                       PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
    310     psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_SIG",
    311                       PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
    312     psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MAX",
    313                       PS_META_REPLACE, "sky model maximum value", statsBck->max);
    314     psMetadataAddF32 (recipe, PS_LIST_TAIL, "MSKY_MIN",
    315                       PS_META_REPLACE, "sky model minimum value", statsBck->min);
    316     psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NX",
    317                       PS_META_REPLACE, "sky model size (x)",      model->numCols);
    318     psMetadataAddS32 (recipe, PS_LIST_TAIL, "MSKY_NY",
    319                       PS_META_REPLACE, "sky model size (y)",      model->numRows);
     312    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MN", PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     313    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_SIG", PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     314    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MAX", PS_META_REPLACE, "sky model maximum value", statsBck->max);
     315    psMetadataAddF32 (readout->analysis, PS_LIST_TAIL, "MSKY_MIN", PS_META_REPLACE, "sky model minimum value", statsBck->min);
     316    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NX", PS_META_REPLACE, "sky model size (x)",      model->numCols);
     317    psMetadataAddS32 (readout->analysis, PS_LIST_TAIL, "MSKY_NY", PS_META_REPLACE, "sky model size (y)",      model->numRows);
    320318    psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
    321319              statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
    … …  
    356354    // find the currently selected readout
    357355    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
    358     PS_ASSERT_PTR_NON_NULL (file, false);
     356    psAssert (file, "missing file?");
    359357
    360358    pmFPA *inFPA = file->fpa;
    361359    pmReadout *readout = pmFPAviewThisReadout(view, inFPA);
     360    psAssert (readout, "missing readout?");
    362361
    363362    psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters

  • branches/eam_branches/20091201/psphot/src/psphotOutput.c

    r26542 r26788  
    126126}
    127127
    128 bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
    129 
    130     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
    131     PS_ASSERT (file, false);
    132 
    133     pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
    134 
    135     // determine PHOTCODE from fpa & view, overwrite in readout->analysis
    136     char *photcode = pmConceptsPhotcodeForView (file, view);
    137     PS_ASSERT (photcode, false);
    138 
    139     psMetadataAddStr (readout->analysis, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "photcode from FPA concepts", photcode);
    140     psLogMsg ("psphot", 3, "PHOTCODE is %s", photcode);
    141 
    142     psFree (photcode);
    143     return true;
    144 }
    145 
    146128bool psphotAddPhotcode (pmConfig *config, const pmFPAview *view) {
    147129
    … …  
    161143}
    162144
    163 bool psphotSetHeaderNstars (psMetadata *recipe, psArray *sources) {
     145bool psphotAddPhotcodeReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
     146
     147    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
     148    PS_ASSERT (file, false);
     149
     150    pmReadout  *readout = pmFPAviewThisReadout (view, file->fpa);
     151
     152    // determine PHOTCODE from fpa & view, overwrite in readout->analysis
     153    char *photcode = pmConceptsPhotcodeForView (file, view);
     154    PS_ASSERT (photcode, false);
     155
     156    psMetadataAddStr (readout->analysis, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "photcode from FPA concepts", photcode);
     157    psLogMsg ("psphot", 3, "PHOTCODE is %s", photcode);
     158
     159    psFree (photcode);
     160    return true;
     161}
     162
     163bool psphotSetHeaderNstars (psMetadata *header, psArray *sources) {
    164164
    165165    int nSrc = 0;
    … …  
    190190    }
    191191
    192     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NSTARS",   PS_META_REPLACE, "Number of sources", nSrc);
    193     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT);
    194     psMetadataAddS32 (recipe, PS_LIST_TAIL, "NDET_CR",  PS_META_REPLACE, "Number of cosmic rays", nCR);
    195     return true;
    196 }
    197 
    198 // these values are saved in an output header stub - they are added to either the
    199 // PHU header (CMP) or the MEF table header (CMF)
    200 // XXX these are currently carried by the recipe -- this will not work in a Stack context
    201 // XXX they should be place in the readout->analysis of the given image
    202 psMetadata *psphotDefineHeader (psMetadata *recipe) {
     192    psMetadataAddS32 (header, PS_LIST_TAIL, "NSTARS",   PS_META_REPLACE, "Number of sources", nSrc);
     193    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_EXT", PS_META_REPLACE, "Number of extended sources", nEXT);
     194    psMetadataAddS32 (header, PS_LIST_TAIL, "NDET_CR",  PS_META_REPLACE, "Number of cosmic rays", nCR);
     195    return true;
     196}
     197
     198// these values are saved in an output header stub - they are added to either the PHU header
     199// (CMP) or the MEF table header (CMF).  before the header is created, each readout has these
     200// values stored on readout->analysis
     201psMetadata *psphotDefineHeader (psMetadata *analysis) {
    203202
    204203    bool status = true;
    … …  
    207206
    208207    // write necessary information to output header
    209     psMetadataItemSupplement (&status, header, recipe, "ZERO_PT");
    210     psMetadataItemSupplement (&status, header, recipe, "PHOTCODE");
    211 
    212     psMetadataItemSupplement (&status, header, recipe, "APMIFIT");
    213     psMetadataItemSupplement (&status, header, recipe, "DAPMIFIT");
    214     psMetadataItemSupplement (&status, header, recipe, "NAPMIFIT");
     208    psMetadataItemSupplement (&status, header, analysis, "ZERO_PT");
     209    psMetadataItemSupplement (&status, header, analysis, "PHOTCODE");
     210
     211    psMetadataItemSupplement (&status, header, analysis, "APMIFIT");
     212    psMetadataItemSupplement (&status, header, analysis, "DAPMIFIT");
     213    psMetadataItemSupplement (&status, header, analysis, "NAPMIFIT");
    215214
    216215    // PSF model parameters (shape values for image center)
    217     psMetadataItemSupplement (&status, header, recipe, "NPSFSTAR");
    218     psMetadataItemSupplement (&status, header, recipe, "APLOSS");
    219 
    220     psMetadataItemSupplement (&status, header, recipe, "FWHM_MAJ");
    221     psMetadataItemSupplement (&status, header, recipe, "FW_MJ_SG");
    222     psMetadataItemSupplement (&status, header, recipe, "FW_MJ_LQ");
    223     psMetadataItemSupplement (&status, header, recipe, "FW_MJ_UQ");
    224 
    225     psMetadataItemSupplement (&status, header, recipe, "FWHM_MIN");
    226     psMetadataItemSupplement (&status, header, recipe, "FW_MN_SG");
    227     psMetadataItemSupplement (&status, header, recipe, "FW_MN_LQ");
    228     psMetadataItemSupplement (&status, header, recipe, "FW_MN_UQ");
    229 
    230     psMetadataItemSupplement (&status, header, recipe, "ANGLE");
     216    psMetadataItemSupplement (&status, header, analysis, "NPSFSTAR");
     217    psMetadataItemSupplement (&status, header, analysis, "APLOSS");
     218
     219    psMetadataItemSupplement (&status, header, analysis, "FWHM_MAJ");
     220    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_SG");
     221    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_LQ");
     222    psMetadataItemSupplement (&status, header, analysis, "FW_MJ_UQ");
     223
     224    psMetadataItemSupplement (&status, header, analysis, "FWHM_MIN");
     225    psMetadataItemSupplement (&status, header, analysis, "FW_MN_SG");
     226    psMetadataItemSupplement (&status, header, analysis, "FW_MN_LQ");
     227    psMetadataItemSupplement (&status, header, analysis, "FW_MN_UQ");
     228
     229    psMetadataItemSupplement (&status, header, analysis, "ANGLE");
    231230
    232231    // Image Quality measurements
    233     psMetadataItemSupplement (&status, header, recipe, "IQ_NSTAR");
    234 
    235     psMetadataItemSupplement (&status, header, recipe, "IQ_FW1");
    236     psMetadataItemSupplement (&status, header, recipe, "IQ_FW1_E");
    237     psMetadataItemSupplement (&status, header, recipe, "IQ_FW2");
    238     psMetadataItemSupplement (&status, header, recipe, "IQ_FW2_E");
    239 
    240     psMetadataItemSupplement (&status, header, recipe, "IQ_M2");
    241     psMetadataItemSupplement (&status, header, recipe, "IQ_M2_ER");
    242     psMetadataItemSupplement (&status, header, recipe, "IQ_M2_LQ");
    243     psMetadataItemSupplement (&status, header, recipe, "IQ_M2_UQ");
    244 
    245     psMetadataItemSupplement (&status, header, recipe, "IQ_M2C");
    246     psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_E");
    247     psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_L");
    248     psMetadataItemSupplement (&status, header, recipe, "IQ_M2C_U");
    249 
    250     psMetadataItemSupplement (&status, header, recipe, "IQ_M2S");
    251     psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_E");
    252     psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_L");
    253     psMetadataItemSupplement (&status, header, recipe, "IQ_M2S_U");
    254 
    255     psMetadataItemSupplement (&status, header, recipe, "IQ_M3");
    256     psMetadataItemSupplement (&status, header, recipe, "IQ_M3_ER");
    257     psMetadataItemSupplement (&status, header, recipe, "IQ_M3_LQ");
    258     psMetadataItemSupplement (&status, header, recipe, "IQ_M3_UQ");
    259 
    260     psMetadataItemSupplement (&status, header, recipe, "IQ_M4");
    261     psMetadataItemSupplement (&status, header, recipe, "IQ_M4_ER");
    262     psMetadataItemSupplement (&status, header, recipe, "IQ_M4_LQ");
    263     psMetadataItemSupplement (&status, header, recipe, "IQ_M4_UQ");
     232    psMetadataItemSupplement (&status, header, analysis, "IQ_NSTAR");
     233
     234    psMetadataItemSupplement (&status, header, analysis, "IQ_FW1");
     235    psMetadataItemSupplement (&status, header, analysis, "IQ_FW1_E");
     236    psMetadataItemSupplement (&status, header, analysis, "IQ_FW2");
     237    psMetadataItemSupplement (&status, header, analysis, "IQ_FW2_E");
     238
     239    psMetadataItemSupplement (&status, header, analysis, "IQ_M2");
     240    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_ER");
     241    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_LQ");
     242    psMetadataItemSupplement (&status, header, analysis, "IQ_M2_UQ");
     243
     244    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C");
     245    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_E");
     246    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_L");
     247    psMetadataItemSupplement (&status, header, analysis, "IQ_M2C_U");
     248
     249    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S");
     250    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_E");
     251    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_L");
     252    psMetadataItemSupplement (&status, header, analysis, "IQ_M2S_U");
     253
     254    psMetadataItemSupplement (&status, header, analysis, "IQ_M3");
     255    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_ER");
     256    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_LQ");
     257    psMetadataItemSupplement (&status, header, analysis, "IQ_M3_UQ");
     258
     259    psMetadataItemSupplement (&status, header, analysis, "IQ_M4");
     260    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_ER");
     261    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_LQ");
     262    psMetadataItemSupplement (&status, header, analysis, "IQ_M4_UQ");
    264263
    265264    // XXX these need to be defined from elsewhere
    266265    psMetadataAdd (header, PS_LIST_TAIL, "FSATUR",   PS_DATA_F32 | PS_META_REPLACE, "SATURATION MAG",      0.0);
    267266    psMetadataAdd (header, PS_LIST_TAIL, "FLIMIT",   PS_DATA_F32 | PS_META_REPLACE, "COMPLETENESS MAG",    0.0);
    268     psMetadataItemSupplement (&status, header, recipe, "NSTARS");
    269 
    270     psMetadataItemSupplement (&status, header, recipe, "NDET_EXT");
    271     psMetadataItemSupplement (&status, header, recipe, "NDET_CR");
     267    psMetadataItemSupplement (&status, header, analysis, "NSTARS");
     268
     269    psMetadataItemSupplement (&status, header, analysis, "NDET_EXT");
     270    psMetadataItemSupplement (&status, header, analysis, "NDET_CR");
    272271
    273272    // sky background model statistics
    274     psMetadataItemSupplement (&status, header, recipe, "MSKY_MN");
    275     psMetadataItemSupplement (&status, header, recipe, "MSKY_SIG");
    276     psMetadataItemSupplement (&status, header, recipe, "MSKY_MIN");
    277     psMetadataItemSupplement (&status, header, recipe, "MSKY_MAX");
    278     psMetadataItemSupplement (&status, header, recipe, "MSKY_NX");
    279     psMetadataItemSupplement (&status, header, recipe, "MSKY_NY");
     273    psMetadataItemSupplement (&status, header, analysis, "MSKY_MN");
     274    psMetadataItemSupplement (&status, header, analysis, "MSKY_SIG");
     275    psMetadataItemSupplement (&status, header, analysis, "MSKY_MIN");
     276    psMetadataItemSupplement (&status, header, analysis, "MSKY_MAX");
     277    psMetadataItemSupplement (&status, header, analysis, "MSKY_NX");
     278    psMetadataItemSupplement (&status, header, analysis, "MSKY_NY");
    280279
    281280    psMetadataAddF32 (header, PS_LIST_TAIL, "DT_PHOT", PS_META_REPLACE, "elapsed psphot time", psTimerMark ("psphotReadout"));

  • branches/eam_branches/20091201/psphot/src/psphotRadiusChecks.c

    r25755 r26788  
    44static float PSF_FIT_NSIGMA;
    55static float PSF_FIT_PADDING;
     6static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
    67static float PSF_FIT_RADIUS = 0;        // radius to use in fitting (ignored if <= 0,
    78                                        // and a per-object radius is calculated)
    89
    9 static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
    10 
    11 
    12 bool psphotInitRadiusPSF(const psMetadata *recipe, const pmModelType type) {
     10bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
    1311
    1412    bool status = true;
    1513
    16     PSF_FIT_NSIGMA  = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA");
     14    PSF_FIT_NSIGMA = psMetadataLookupF32(&status, recipe, "PSF_FIT_NSIGMA");
    1715    PSF_FIT_PADDING = psMetadataLookupF32(&status, recipe, "PSF_FIT_PADDING");
    18     PSF_FIT_RADIUS  =  psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
    19     PSF_APERTURE    =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
     16
     17    PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
     18    if (!status) {
     19        PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
     20    }
     21
     22    PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
     23    if (!status) {
     24        PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
     25    }
    2026
    2127    return true;

  • branches/eam_branches/20091201/psphot/src/psphotReadout.c

    r26596 r26788  
    33// this should be called by every program that links against libpsphot
    44bool psphotInit (void) {
    5 
    65    psphotErrorRegister();              // register our error codes/messages
    76    psphotModelClassInit ();            // load implementation-specific models
    … …  
    1211bool psphotReadout(pmConfig *config, const pmFPAview *view) {
    1312
    14     // measure the total elapsed time in psphotReadout.  XXX the current threading plan
    15     // for psphot envisions threading within psphotReadout, not multiple threads calling
    16     // the same psphotReadout.  In the current plan, this dtime is the elapsed time used
    17     // jointly by the multiple threads, not the total time used by all threads.
     13    // measure the total elapsed time in psphotReadout.  dtime is the elapsed time used jointly
     14    // by the multiple threads, not the total time used by all threads.
    1815    psTimerStart ("psphotReadout");
    1916
    … …  
    2623        return false;
    2724    }
     25    // optional break-point for processing
     26    char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
     27    psAssert (breakPt, "configuration error: set BREAK_POINT");
    2828
    2929    // set the photcode for this image
    … …  
    3333    }
    3434
    35     // find the currently selected readout
    36     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    37     PS_ASSERT_PTR_NON_NULL (readout, false);
    38 
    39     // optional break-point for processing
    40     char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT");
    41     PS_ASSERT_PTR_NON_NULL (breakPt, false);
    42 
    4335    // Generate the mask and weight images, including the user-defined analysis region of interest
    44     psphotSetMaskAndVariance (config, view, recipe);
     36    psphotSetMaskAndVariance (config, view);
    4537    if (!strcasecmp (breakPt, "NOTHING")) {
    46         return psphotReadoutCleanup(config, readout, recipe, NULL, NULL, NULL);
    47     }
    48 
    49     // display the image, weight, mask (ch 1,2,3)
    50     psphotVisualShowImage (readout);
     38        return psphotReadoutCleanup(config, view);
     39    }
    5140
    5241    // generate a background model (median, smoothed image)
    5342    if (!psphotModelBackground (config, view)) {
    54         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     43        return psphotReadoutCleanup (config, view);
    5544    }
    5645    if (!psphotSubtractBackground (config, view)) {
    57         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     46        return psphotReadoutCleanup (config, view);
    5847    }
    5948    if (!strcasecmp (breakPt, "BACKMDL")) {
    60         return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
    61     }
    62 
    63     // display the backsub and backgnd images
    64     psphotVisualShowBackground (config, view, readout);
    65 
    66     // run a single-model test if desired (exits from here if test is run)
    67     psphotModelTest (config, view, recipe);
    68 
    69     // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    70     pmPSF *psf = psphotLoadPSF (config, view, recipe);
    71 
    72     // several functions below behave differently if we have a PSF model already
    73     bool havePSF = (psf != NULL);
    74 
     49        return psphotReadoutCleanup (config, view);
     50    }
     51
     52    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are determined and saved on
     53    // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT
     54    if (!psphotLoadPSF (config, view)) {
     55        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
     56        return psphotReadoutCleanup (config, view);
     57    }
     58       
    7559    // find the detections (by peak and/or footprint) in the image.
    76     pmDetections *detections = psphotFindDetections (NULL, readout, recipe);
    77     if (!detections) {
     60    if (!psphotFindDetections (config, view, true)) { // pass 1
    7861        // this only happens if we had an error in psphotFindDetections
    7962        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    80         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    81     }
    82     if (!detections->peaks->n) {
    83         psLogMsg ("psphot", 3, "unable to find detections in this image");
    84         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    85     }
    86 
    87     // construct sources and measure basic stats
    88     psArray *sources = psphotSourceStats (config, readout, detections, true);
    89     if (!sources) return false;
     63        return psphotReadoutCleanup (config, view);
     64    }
     65
     66    // construct sources and measure basic stats (saved on detections->newSources)
     67    if (!psphotSourceStats (config, view, true)) { // pass 1
     68        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     69        return psphotReadoutCleanup (config, view);
     70    }
    9071    if (!strcasecmp (breakPt, "PEAKS")) {
    91         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    92     }
    93 
    94     // find blended neighbors of very saturated stars
    95     // XXX merge this with Basic Deblend?
    96     psphotDeblendSatstars (readout, sources, recipe);
    97 
    98     // mark blended peaks PS_SOURCE_BLEND
    99     if (!psphotBasicDeblend (sources, recipe)) {
    100         psLogMsg ("psphot", 3, "failed on deblend analysis");
    101         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    102     }
    103 
    104     // classify sources based on moments, brightness
    105     if (!psphotRoughClass (readout, sources, recipe, havePSF)) {
    106         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    107         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     72        return psphotReadoutCleanup(config, view);
     73    }
     74
     75    // find blended neighbors of very saturated stars (detections->newSources)
     76    if (!psphotDeblendSatstars (config, view)) {
     77        psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
     78        return psphotReadoutCleanup (config, view);
     79    }
     80
     81    // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
     82    if (!psphotBasicDeblend (config, view)) {
     83        psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis");
     84        return psphotReadoutCleanup (config, view);
     85    }
     86
     87    // classify sources based on moments, brightness.  if a PSF model has been loaded, the PSF
     88    // clump defined for it is used not measured (detections->newSources)
     89    if (!psphotRoughClass (config, view)) { // pass 1
     90        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications");
     91        return psphotReadoutCleanup (config, view);
    10892    }
    10993    if (!strcasecmp (breakPt, "MOMENTS")) {
    110         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    111     }
    112 
    113     if (!havePSF && !psphotImageQuality (recipe, sources)) {
    114         psLogMsg("psphot", 3, "failed to measure image quality");
    115         return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
    116     }
    117 
    118     // if we were not supplied a PSF, choose one here
    119     if (psf == NULL) {
    120         // use bright stellar objects to measure PSF
    121         // XXX if we do not have enough stars to generate the PSF, build one
    122         // from the SEEING guess and model class
    123         psf = psphotChoosePSF (readout, sources, recipe);
    124         if (psf == NULL) {
    125             psLogMsg ("psphot", 3, "failure to construct a psf model");
    126             return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    127         }
    128         havePSF = true;
     94        return psphotReadoutCleanup(config, view);
     95    }
     96
     97    // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     98    if (!psphotImageQuality (config, view)) { // pass 1
     99        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     100        return psphotReadoutCleanup(config, view);
     101    }
     102
     103    // use bright stellar objects to measure PSF if we were supplied a PSF for any input file,
     104    // this step is skipped
     105    if (!psphotChoosePSF (config, view)) { // pass 1
     106        psLogMsg ("psphot", 3, "failure to construct a psf model");
     107        return psphotReadoutCleanup (config, view);
    129108    }
    130109    if (!strcasecmp (breakPt, "PSFMODEL")) {
    131         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    132     }
    133     psphotVisualShowPSFModel (readout, psf);
     110        return psphotReadoutCleanup (config, view);
     111    }
    134112
    135113    // include externally-supplied sources
    136     psphotLoadExtSources (config, view, sources);
    137 
    138     // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    139     psphotGuessModels (config, readout, sources, psf);
     114    // XXX fix this in the new multi-input context
     115    psphotLoadExtSources (config, view); // pass 1
     116
     117    // construct an initial model for each object, set the radius to fitRadius, set circular
     118    // fit mask (detections->newSources)
     119    psphotGuessModels (config, view); // pass 1
     120
     121    // merge the newly selected sources into the existing list
     122    // NOTE: merge OLD and NEW
     123    psphotMergeSources (config, view);
    140124
    141125    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    142     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    143 
    144     // We have to place these visualizations here because the models are not realized until
    145     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    146     psphotVisualShowPSFStars (recipe, psf, sources);
    147 
    148     // identify CRs and extended sources
    149     psphotSourceSize (config, readout, sources, recipe, psf, 0);
     126    psphotFitSourcesLinear (config, view, false); // pass 1 (detections->allSources)
     127
     128    // identify CRs and extended sources (only unmeasured sources are measured)
     129    psphotSourceSize (config, view, true); // pass 1 (detections->allSources)
    150130    if (!strcasecmp (breakPt, "ENSEMBLE")) {
    151131        goto finish;
    152132    }
    153     psphotVisualShowSatStars (recipe, psf, sources);
    154133
    155134    // non-linear PSF and EXT fit to brighter sources
    156135    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    157     psphotBlendFit (config, readout, sources, psf);
     136    psphotBlendFit (config, view); // pass 1 (detections->allSources)
    158137
    159138    // replace all sources
    160     psphotReplaceAllSources (sources, recipe);
     139    psphotReplaceAllSources (config, view); // pass 1 (detections->allSources)
    161140
    162141    // linear fit to include all sources (subtract again)
    163     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     142    // NOTE : apply to ALL sources (extended + psf)
     143    psphotFitSourcesLinear (config, view, true); // pass 2 (detections->allSources)
    164144
    165145    // if we only do one pass, skip to extended source analysis
    166     if (!strcasecmp (breakPt, "PASS1")) {
    167         goto pass1finish;
    168     }
    169     // NOTE: possibly re-measure background model here with objects subtracted
     146    if (!strcasecmp (breakPt, "PASS1")) goto pass1finish;
     147
     148    // NOTE: possibly re-measure background model here with objects subtracted / or masked
    170149
    171150    // add noise for subtracted objects
    172     psphotAddNoise (readout, sources, recipe);
    173 
    174     // find fainter sources (pass 2)
    175     detections = psphotFindDetections (detections, readout, recipe);
     151    psphotAddNoise (config, view); // pass 1 (detections->allSources)
     152
     153    // find fainter sources
     154    // NOTE: finds new peaks and new footprints, OLD and FULL set are saved on detections
     155    psphotFindDetections (config, view, false); // pass 2 (detections->peaks, detections->footprints)
    176156
    177157    // remove noise for subtracted objects (ie, return to normal noise level)
    178     psphotSubNoise (readout, sources, recipe);
     158    // NOTE: this needs to operate only on the OLD sources
     159    psphotSubNoise (config, view); // pass 1 (detections->allSources)
    179160
    180161    // define new sources based on only the new peaks
    181     psArray *newSources = psphotSourceStats (config, readout, detections, false);
     162    // NOTE: new sources are saved on detections->newSources
     163    psphotSourceStats (config, view, false); // pass 2 (detections->newSources)
    182164
    183165    // set source type
    184     if (!psphotRoughClass (readout, newSources, recipe, havePSF)) {
     166    // NOTE: apply only to detections->newSources
     167    if (!psphotRoughClass (config, view)) { // pass 2 (detections->newSources)
    185168        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    186         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     169        return psphotReadoutCleanup (config, view);
    187170    }
    188171
    189172    // create full input models, set the radius to fitRadius, set circular fit mask
    190     psphotGuessModels (config, readout, newSources, psf);
     173    // NOTE: apply only to detections->newSources
     174    psphotGuessModels (config, view); // pass 2 (detections->newSources)
    191175
    192176    // replace all sources so fit below applies to all at once
    193     psphotReplaceAllSources (sources, recipe);
     177    // NOTE: apply only to OLD sources (which have been subtracted)
     178    psphotReplaceAllSources (config, view); // pass 2
    194179
    195180    // merge the newly selected sources into the existing list
    196     psphotMergeSources (sources, newSources);
    197     psFree (newSources);
    198 
    199     // linear fit to all sources
    200     psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE);
     181    // NOTE: merge OLD and NEW
     182    // XXX check on free of sources...
     183    psphotMergeSources (config, view); // (detections->newSources + detections->allSources -> detections->allSources)
     184
     185    // NOTE: apply to ALL sources
     186    psphotFitSourcesLinear (config, view, true); // pass 3 (detections->allSources)
    201187
    202188pass1finish:
    203189
    204190    // measure source size for the remaining sources
    205     psphotSourceSize (config, readout, sources, recipe, psf, 0);
    206 
    207     psphotExtendedSourceAnalysis (readout, sources, recipe);
    208 
    209     psphotExtendedSourceFits (readout, sources, recipe);
     191    // NOTE: applies only to NEW (unmeasured) sources
     192    psphotSourceSize (config, view, false); // pass 2 (detections->allSources)
     193
     194    psphotExtendedSourceAnalysis (config, view); // pass 1 (detections->allSources)
     195    psphotExtendedSourceFits (config, view); // pass 1 (detections->allSources)
    210196
    211197finish:
    … …  
    215201
    216202    // measure aperture photometry corrections
    217     if (!psphotApResid (config, readout, sources, psf)) {
     203    if (!psphotApResid (config, view)) { // pass 1 (detections->allSources)
    218204        psLogMsg ("psphot", 3, "failed on psphotApResid");
    219         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     205        return psphotReadoutCleanup (config, view);
    220206    }
    221207
    222208    // calculate source magnitudes
    223     psphotMagnitudes(config, readout, view, sources, psf);
    224 
    225     if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
     209    psphotMagnitudes(config, view); // pass 1 (detections->allSources)
     210
     211    if (!psphotEfficiency(config, view)) { // pass 1
    226212        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    227213        psErrorClear();
    … …  
    232218
    233219    // replace background in residual image
    234     psphotSkyReplace (config, view);
     220    psphotSkyReplace (config, view); // pass 1
    235221
    236222    // drop the references to the image pixels held by each source
    237     psphotSourceFreePixels (sources);
     223    psphotSourceFreePixels (config, view); // pass 1
    238224
    239225    // create the exported-metadata and free local data
    240     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     226    return psphotReadoutCleanup(config, view);
    241227}
    242 

  • branches/eam_branches/20091201/psphot/src/psphotReadoutCleanup.c

    r24203 r26788  
    11# include "psphotInternal.h"
    22
    3 // psphotReadoutCleanup is called on exit from psphotReadout.  If the last raised error is
    4 // not a DATA error, then there was a serious problem.  Only in this case, or if the fail
    5 // on the stats measurement, do we return false
    6 bool psphotReadoutCleanup (pmConfig *config, pmReadout *readout, psMetadata *recipe, pmDetections *detections, pmPSF *psf, psArray *sources) {
     3// for now, let's store the detections on the readout->analysis for each readout
     4bool psphotReadoutCleanup (pmConfig *config, const pmFPAview *view)
     5{
     6    bool status = true;
    77
    88    // remove internal pmFPAfiles, if created
    … …  
    1212    }
    1313    if (psErrorCodeLast() != PS_ERR_NONE) {
    14         psFree (psf);
    15         psFree (sources);
    16         psFree (detections);
    1714        return false;
    1815    }
    1916
     17    // select the appropriate recipe information
     18    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     19    psAssert (recipe, "missing recipe?");
     20
     21    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     22    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     23
     24    // loop over the available readouts
     25    for (int i = 0; i < num; i++) {
     26        if (!psphotReadoutCleanupReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     27            psError (PSPHOT_ERR_CONFIG, false, "failed on psphotReadoutCleanup for PSPHOT.INPUT entry %d", i);
     28            return false;
     29        }
     30    }
     31
     32    // XXX move this to top of loop?
     33    pmKapaClose ();
     34
     35    return true;
     36}
     37
     38// psphotReadoutCleanup is called on exit from psphotReadout.  If the last raised error is
     39// not a DATA error, then there was a serious problem.  Only in this case, or if the fail
     40// on the stats measurement, do we return false
     41bool psphotReadoutCleanupReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
     42
     43    bool status = true;
     44
     45    // find the currently selected readout
     46    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     47    psAssert (file, "missing file?");
     48
     49    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     50    psAssert (readout, "missing readout?");
     51
     52    // when psphotReadoutCleanup is called, these are not necessarily defined
     53    pmPSF        *psf        = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     54    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     55    psArray      *sources    = detections->allSources;
     56    // XXX where do we free these, in here (psMetadataRemove?)
     57
    2058    // use the psf-model to measure FWHM stats
    2159    if (psf) {
    22         if (!psphotPSFstats (readout, recipe, psf)) {
     60        if (!psphotPSFstats (readout, psf)) {
    2361            psError(PSPHOT_ERR_PROG, false, "Failed to measure PSF shape parameters");
    24             psFree (psf);
    25             psFree (sources);
    26             psFree (detections);
    2762            return false;
    2863        }
    … …  
    3065    // otherwise, use the source moments to measure FWHM stats
    3166    if (!psf && sources) {
    32         if (!psphotMomentsStats (readout, recipe, sources)) {
     67        if (!psphotMomentsStats (readout, sources)) {
    3368            psError(PSPHOT_ERR_PROG, false, "Failed to measure Moment shape parameters");
    34             psFree (psf);
    35             psFree (sources);
    36             psFree (detections);
    3769            return false;
    3870        }
    … …  
    4072
    4173    // Check to see if the image quality was measured
     74    // XXX not sure we want / need this test
    4275    if (!psf) {
    4376        bool mdok;                      // Status of MD lookup
    … …  
    4578        if (!mdok || nIQ <= 0) {
    4679            psError(PSPHOT_ERR_DATA, false, "Unable to measure image quality");
    47             psFree (psf);
    48             psFree (sources);
    49             psFree (detections);
    5080            return false;
    5181        }
    5282    }
    5383
     84    // create an output header with stats results currently saved on readout->analysis
     85    psMetadata *header = psphotDefineHeader (readout->analysis);
    5486
    5587    // write NSTARS to the image header
    56     psphotSetHeaderNstars (recipe, sources);
    57 
    58     // create an output header with stats results
    59     psMetadata *header = psphotDefineHeader (recipe);
     88    psphotSetHeaderNstars (header, sources);
    6089
    6190    // save the results of the analysis
     91    // this should happen way up stream (when needed?)
    6292    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.HEADER",  PS_DATA_METADATA, "header stats", header);
    63     if (sources) {
    64         psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "psphot sources", sources);
    65     }
     93
    6694    if (psf) {
     95        // XXX this seems a little silly : we saved the psf on readout->analysis above, but now
     96        // we are moving it to chip->analysis.
    6797        // save the psf for possible output.  if there was already an entry, it was loaded from external sources
    6898        // the new one may have been updated or modified, so replace the existing entry.  We
    … …  
    79109    }
    80110
    81     // XXX move this to top of loop?
    82     pmKapaClose ();
    83 
    84     psFree (detections);
    85     psFree (psf);
    86111    psFree (header);
    87     psFree (sources);
    88 
    89112    return true;
    90113}

  • branches/eam_branches/20091201/psphot/src/psphotReadoutFindPSF.c

    r26611 r26788  
    77    psTimerStart ("psphotReadout");
    88
    9     // select the current recipe
    10     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    11     if (!recipe) {
    12         psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
    13         return false;
    14     }
    15 
    169    // set the photcode for the PSPHOT.INPUT
    1710    if (!psphotAddPhotcode(config, view)) {
    … …  
    2013    }
    2114
    22     // find the currently selected readout
    23     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    24     PS_ASSERT_PTR_NON_NULL (readout, false);
    25 
    2615    // Generate the mask and variance images, including the user-defined analysis region of interest
    27     psphotSetMaskAndVariance (config, view, recipe);
    28 
    29     // display the image, weight, mask (ch 1,2,3)
    30     psphotVisualShowImage (readout);
     16    psphotSetMaskAndVariance (config, view);
    3117
    3218    // Note that in this implementation, we do NOT model the background and we do not
    … …  
    3420
    3521    // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF)
    36     pmDetections *detections = psphotDetectionsFromSources (config, inSources);
    37     if (!detections || !detections->peaks) {
     22    // XXX we assume a single set of input sources is supplied
     23    if (!psphotDetectionsFromSources (config, view, inSources)) {
    3824        psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars");
    39         return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL);
     25        return psphotReadoutCleanup(config, view);
    4026    }
    4127
    42     // construct sources and measure basic stats (moments, local sky)
    43     psArray *sources = psphotSourceStats(config, readout, detections, true);
    44     if (!sources) return false;
     28    // construct detections->newSources and measure basic stats (moments, local sky)
     29    if (!psphotSourceStats(config, view, true)) {
     30        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     31        return false;
     32    }
    4533
    46     // peak flux is wrong : set based on previous image
    47     // use the peak measured in the moments analysis:
    48     for (int i = 0; i < sources->n; i++) {
    49       pmSource *source = sources->data[i];
    50       source->peak->flux = source->moments->Peak;
     34    // peak flux is wrong : use the peak measured in the moments analysis:
     35    if (!psphotRepairLoadedSources(config, view)) {
     36        psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources");
     37        return false;
    5138    }
    5239
    5340    // classify sources based on moments, brightness (psf is not known)
    54     if (!psphotRoughClass (readout, sources, recipe, false)) {
    55         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    56         return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources);
     41    if (!psphotRoughClass (config, view)) {
     42        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class");
     43        return psphotReadoutCleanup (config, view);
    5744    }
    5845
    59     if (!psphotImageQuality (recipe, sources)) {
    60         psLogMsg("psphot", 3, "failed to measure image quality");
    61         return psphotReadoutCleanup(config, readout, recipe, detections, NULL, sources);
     46    if (!psphotImageQuality (config, view)) {
     47        psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality");
     48        return psphotReadoutCleanup(config, view);
    6249    }
    6350
    64     pmPSF *psf = psphotChoosePSF(readout, sources, recipe);
    65     if (!psf) {
     51    if (!psphotChoosePSF(config, view)) {
    6652        psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
    67         psFree(sources);
    68         return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL);
     53        return psphotReadoutCleanup(config, view);
    6954    }
    70     psphotVisualShowPSFModel(readout, psf);
    7155
    7256# if 0
    … …  
    7559    // fits from that analysis, or run the linear PSF fit for all objects currently in hand
    7660    // construct an initial model for each object, set the radius to fitRadius, set circular fit mask
    77     psphotGuessModels (config, readout, sources, psf);
     61    psphotGuessModels (config, view);
     62# endif
    7863
     64    // merge the newly selected sources into the existing list
     65    // NOTE: merge OLD and NEW
     66    psphotMergeSources (config, view);
     67
     68# if 0
    7969    // measure aperture photometry corrections
    80     if (!psphotApResid (config, readout, sources, psf)) {
     70    if (!psphotApResid (config, view)) {
    8171        psLogMsg ("psphot", 3, "failed on psphotApResid");
    82         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     72        return psphotReadoutCleanup (config, view);
    8373    }
    8474# endif
    8575
    8676    // drop the references to the image pixels held by each source
    87     psphotSourceFreePixels(sources);
    88     psFree(sources);
     77    psphotSourceFreePixels(config, view);
    8978
    9079    // create the exported-metadata and free local data
    91     return psphotReadoutCleanup(config, readout, recipe, detections, psf, NULL);
     80    return psphotReadoutCleanup(config, view);
    9281}

  • branches/eam_branches/20091201/psphot/src/psphotReadoutKnownSources.c

    r26542 r26788  
    11# include "psphotInternal.h"
    22
    3 // in this psphotReadout-variant, we are only measuring the photometry for known sources,
    4 // using a PSF generated from this observation those sources
     3// in this psphotReadout-variant, we are only measuring the photometry for known sources, using
     4// a PSF generated for this observation from those sources
    55bool psphotReadoutKnownSources(pmConfig *config, const pmFPAview *view, psArray *inSources) {
    66
    77    psTimerStart ("psphotReadout");
    8 
    9     // select the current recipe
    10     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    11     if (!recipe) {
    12         psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
    13         return false;
    14     }
    158
    169    // set the photcode for this image
    … …  
    2013    }
    2114
    22     // find the currently selected readout
    23     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    24     PS_ASSERT_PTR_NON_NULL (readout, false);
    25 
    2615    // Generate the mask and weight images, including the user-defined analysis region of interest
    27     psphotSetMaskAndVariance (config, view, recipe);
    28 
    29     // display the image, weight, mask (ch 1,2,3)
    30     psphotVisualShowImage (readout);
     16    psphotSetMaskAndVariance (config, view);
    3117
    3218    // Note that in this implementation, we do NOT model the background and we do not
    … …  
    3420
    3521    // include externally-supplied sources (supplied as PSPHOT.INPUT.CMF)
    36     pmDetections *detections = psphotDetectionsFromSources (config, inSources);
    37     if (!detections || !detections->peaks) {
     22    if (!psphotDetectionsFromSources (config, view, inSources)) {
    3823        psError(PSPHOT_ERR_ARGUMENTS, true, "Can't find PSF stars");
    39         return psphotReadoutCleanup(config, readout, recipe, detections, NULL, NULL);
     24        return psphotReadoutCleanup(config, view);
    4025    }
    4126
    4227    // construct sources and measure basic stats
    43     psArray *sources = psphotSourceStats (config, readout, detections, true);
    44     if (!sources) return false;
     28    if (!psphotSourceStats (config, view, true)) {
     29        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     30        return false;
     31    }
    4532
    46     // peak flux is wrong : set based on previous image
    47     // use the peak measured in the moments analysis:
    48     for (int i = 0; i < sources->n; i++) {
    49         pmSource *source = sources->data[i];
    50         source->peak->flux = source->moments->Peak;
     33    // peak flux is wrong : use the peak measured in the moments analysis:
     34    if (!psphotRepairLoadedSources(config, view)) {
     35        psError(PSPHOT_ERR_UNKNOWN, false, "failure to repair sources");
     36        return false;
    5137    }
    5238
    5339    // classify sources based on moments, brightness (psf is not known)
    54     if (!psphotRoughClass (readout, sources, recipe, false)) {
    55         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    56         return psphotReadoutCleanup (config, readout, recipe, detections, NULL, sources);
     40    if (!psphotRoughClass (config, view)) {
     41        psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough source class");
     42        return psphotReadoutCleanup(config, view);
    5743    }
    5844
    59     pmPSF *psf = psphotChoosePSF (readout, sources, recipe);
    60     if (!psf) {
    61         psLogMsg ("psphot", 3, "failure to construct a psf model");
    62         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     45    if (!psphotChoosePSF (config, view)) {
     46        psError(PSPHOT_ERR_PSF, false, "Failed to construct a psf model");
     47        return psphotReadoutCleanup(config, view);
    6348    }
    64     psphotVisualShowPSFModel (readout, psf);
    6549
    6650    // construct an initial model for each object
    67     psphotGuessModels (config, readout, sources, psf);
     51    psphotGuessModels (config, view);
     52
     53    // merge the newly selected sources into the existing list
     54    // NOTE: merge OLD and NEW
     55    psphotMergeSources (config, view);
    6856
    6957    // linear PSF fit to source peaks
    70     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    71 
    72     // We have to place these visualizations here because the models are not realized until
    73     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    74     psphotVisualShowPSFStars (recipe, psf, sources);
    75     psphotVisualShowSatStars (recipe, psf, sources);
     58    psphotFitSourcesLinear (config, view, false);
    7659
    7760    // measure aperture photometry corrections
    78     if (!psphotApResid (config, readout, sources, psf)) {
     61    if (!psphotApResid (config, view)) {
    7962        psLogMsg ("psphot", 3, "failed on psphotApResid");
    80         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
     63        return psphotReadoutCleanup(config, view);
    8164    }
    8265
    8366    // calculate source magnitudes
    84     psphotMagnitudes(config, readout, view, sources, psf);
     67    psphotMagnitudes(config, view);
    8568
    8669    // drop the references to the image pixels held by each source
    87     psphotSourceFreePixels (sources);
     70    psphotSourceFreePixels (config, view);
    8871
    8972    // create the exported-metadata and free local data
    90     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     73    return psphotReadoutCleanup(config, view);
    9174}

  • branches/eam_branches/20091201/psphot/src/psphotReadoutMinimal.c

    r26648 r26788  
    1717    pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
    1818
    19     // select the current recipe
    20     psMetadata *recipe = psMetadataLookupPtr (NULL, config->recipes, PSPHOT_RECIPE);
    21     if (!recipe) {
    22         psError(PSPHOT_ERR_CONFIG, false, "missing recipe %s", PSPHOT_RECIPE);
    23         return false;
    24     }
    25 
    2619    // set the photcode for this image
    2720    if (!psphotAddPhotcode(config, view)) {
    … …  
    3023    }
    3124
    32     // find the currently selected readout
    33     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    34     PS_ASSERT_PTR_NON_NULL (readout, false);
     25    // Generate the mask and weight images, including the user-defined analysis region of interest
     26    psphotSetMaskAndVariance (config, view);
    3527
    36     // Generate the mask and weight images, including the user-defined analysis region of interest
    37     psphotSetMaskAndVariance (config, view, recipe);
    38 
    39     // display the image, weight, mask (ch 1,2,3)
    40     psphotVisualShowImage (readout);
    41 
    42     // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved in the recipe
    43     pmPSF *psf = psphotLoadPSF (config, view, recipe);
    44     if (!psf) {
     28    // load the psf model, if suppled.  FWHM_X,FWHM_Y,etc are saved on readout->analysis
     29    if (!psphotLoadPSF (config, view)) {
    4530      psError (PSPHOT_ERR_CONFIG, false, "missing psf model");
    46       return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL);
     31      return psphotReadoutCleanup (config, view);
    4732    }
    4833
    49     // find the detections (by peak and/or footprint) in the image.
    50     pmDetections *detections = pmDetectionsAlloc(); // New detections; allocated to ensure pass=2
    51     detections = psphotFindDetections(detections, readout, recipe);
    52     if (!detections) {
     34    // find the detections (by peak and/or footprint) in the image. (final pass)
     35    if (!psphotFindDetections(config, view, false)) {
    5336        psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis");
    54         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
     37        return psphotReadoutCleanup (config, view);
    5538    }
    56 #if 0
    57     if (!detections->peaks->n) {
    58         psLogMsg ("psphot", 3, "unable to find detections in this image");
    59         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
     39
     40    // construct sources and measure basic stats (saved on detections->newSources)
     41    if (!psphotSourceStats (config, view, false)) { // pass 1
     42        psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources");
     43        return psphotReadoutCleanup (config, view);
    6044    }
    61 #endif
    6245
    63     // construct sources and measure basic stats
    64     psArray *sources = psphotSourceStats (config, readout, detections, false);
    65     if (!sources) return false;
     46    // find blended neighbors of very saturated stars
     47    psphotDeblendSatstars (config, view);
    6648
    67     // Do the efficiency before everything else, to ensure it's done
    68     if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) {
     49    // mark blended peaks PS_SOURCE_BLEND
     50    if (!psphotBasicDeblend (config, view)) {
     51        psLogMsg ("psphot", 3, "failed on deblend analysis");
     52        return psphotReadoutCleanup (config, view);
     53    }
     54
     55    // classify sources based on moments, brightness (use supplied psf shape parameters)
     56    if (!psphotRoughClass (config, view)) {
     57        psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
     58        return psphotReadoutCleanup (config, view);
     59    }
     60
     61    // construct an initial model for each object
     62    psphotGuessModels (config, view);
     63
     64    // linear PSF fit to source peaks
     65    psphotFitSourcesLinear (config, view, false);
     66
     67// XXX eventually, add the extended source fits here
     68# if (0)
     69    // measure source size for the remaining sources
     70    psphotSourceSize (config, view);
     71
     72    psphotExtendedSourceAnalysis (config, view);
     73
     74    psphotExtendedSourceFits (config, view);
     75# endif
     76
     77    // calculate source magnitudes
     78    psphotMagnitudes(config, view);
     79
     80    // XXX ensure this is measured if the analysis succeeds (even if quality is low)
     81    if (!psphotEfficiency(config, view)) {
    6982        psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
    7083        psErrorClear();
    7184    }
    7285
    73     if (detections->peaks->n == 0 || sources->n == 0) {
    74         psLogMsg ("psphot", 3, "unable to find detections in this image");
    75         return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    76     }
    77 
    78 
    79 
    80     // find blended neighbors of very saturated stars
    81     // XXX merge this with Basic Deblend?
    82     psphotDeblendSatstars (readout, sources, recipe);
    83 
    84     // mark blended peaks PS_SOURCE_BLEND
    85     if (!psphotBasicDeblend (sources, recipe)) {
    86         psLogMsg ("psphot", 3, "failed on deblend analysis");
    87         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    88     }
    89 
    90     // classify sources based on moments, brightness (use supplied psf shape parameters)
    91     if (!psphotRoughClass (readout, sources, recipe, true)) {
    92         psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image");
    93         return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources);
    94     }
    95 
    96     // construct an initial model for each object
    97     psphotGuessModels (config, readout, sources, psf);
    98 
    99     // linear PSF fit to source peaks
    100     psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE);
    101 
    102     // We have to place these visualizations here because the models are not realized until
    103     // psphotGuessModels or fitted until psphotFitSourcesLinear.
    104     psphotVisualShowPSFStars (recipe, psf, sources);
    105     psphotVisualShowSatStars (recipe, psf, sources);
    106 
    107 // XXX eventually, add the extended source fits here
    108 # if (0)
    109     // measure source size for the remaining sources
    110     psphotSourceSize (config, readout, sources, recipe, 0);
    111 
    112     psphotExtendedSourceAnalysis (readout, sources, recipe);
    113 
    114     psphotExtendedSourceFits (readout, sources, recipe);
    115 # endif
    116 
    117     // calculate source magnitudes
    118     psphotMagnitudes(config, readout, view, sources, psf);
    119 
    12086    // drop the references to the image pixels held by each source
    121     psphotSourceFreePixels (sources);
     87    psphotSourceFreePixels (config, view);
    12288
    12389    // create the exported-metadata and free local data
    124     return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     90    return psphotReadoutCleanup(config, view);
    12591}

  • branches/eam_branches/20091201/psphot/src/psphotReplaceUnfit.c

    r25755 r26788  
    2222}
    2323
    24 bool psphotReplaceAllSources (psArray *sources, psMetadata *recipe) {
     24// for now, let's store the detections on the readout->analysis for each readout
     25bool psphotReplaceAllSources (pmConfig *config, const pmFPAview *view)
     26{
     27    bool status = true;
     28
     29    // select the appropriate recipe information
     30    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     31    psAssert (recipe, "missing recipe?");
     32
     33    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     34    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     35
     36    // loop over the available readouts
     37    for (int i = 0; i < num; i++) {
     38        if (!psphotReplaceAllSourcesReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     39            psError (PSPHOT_ERR_CONFIG, false, "failed to replace all sources for PSPHOT.INPUT entry %d", i);
     40            return false;
     41        }
     42    }
     43    return true;
     44}
     45
     46bool psphotReplaceAllSourcesReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    2547
    2648    bool status;
    … …  
    2951    psTimerStart ("psphot.replace");
    3052
     53    // find the currently selected readout
     54    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     55    psAssert (file, "missing file?");
     56
     57    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     58    psAssert (readout, "missing readout?");
     59
     60    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     61    psAssert (detections, "missing detections?");
     62
     63    psArray *sources = detections->allSources;
     64    psAssert (sources, "missing sources?");
     65
    3166    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    3267    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    33     assert (maskVal);
     68    psAssert (maskVal, "missing mask value?");
    3469
    3570    for (int i = 0; i < sources->n; i++) {
    … …  
    67102    return true;
    68103}
    69 
    70 # if (0)
    71 // add source, if the source has been subtracted; do not modify state
    72 bool psphotAddWithTest (pmSource *source, bool useState, psImageMaskType maskVal) {
    73 
    74     // what is current state? (true : add; false : sub)
    75     bool state = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    76     if (state && useState) return true;
    77 
    78     pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    79     return true;
    80 }
    81 
    82 // sub source, if the source has been added; do not modify state
    83 bool psphotSubWithTest (pmSource *source, bool useState, psImageMaskType maskVal) {
    84 
    85     // what is current state? (true : sub; false : add)
    86     bool state = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    87     if (state && useState) return true;
    88 
    89     pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    90     return true;
    91 }
    92 
    93 // add or sub source to match recorded state: supply current state as true (add) or false (sub)
    94 bool psphotSetState (pmSource *source, bool curState, psImageMaskType maskVal) {
    95 
    96     // what is desired state? (true : add; false : sub)
    97     bool newState = !(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
    98     if (curState == newState) return true;
    99 
    100     if (curState && !newState) {
    101         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    102     }
    103     if (newState && !curState) {
    104         pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    105     }
    106     return true;
    107 }
    108 # endif

  • branches/eam_branches/20091201/psphot/src/psphotRoughClass.c

    r25989 r26788  
    77        } }
    88
    9 bool psphotRoughClassRegion (int nRegion, psRegion *region, psArray *sources, psMetadata *target, psMetadata *recipe, const bool havePSF);
     9// for now, let's store the detections on the readout->analysis for each readout
     10bool psphotRoughClass (pmConfig *config, const pmFPAview *view)
     11{
     12    bool status = true;
    1013
    11 // 2006.02.02 : no leaks
    12 bool psphotRoughClass (pmReadout *readout, psArray *sources, psMetadata *recipe, const bool havePSF) {
     14    // select the appropriate recipe information
     15    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     16    psAssert (recipe, "missing recipe?");
     17
     18    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     19    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     20
     21    // loop over the available readouts
     22    for (int i = 0; i < num; i++) {
     23        if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
     24            psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i);
     25            return false;
     26        }
     27    }
     28    return true;
     29}
     30
     31bool psphotRoughClassReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe) {
    1332
    1433    bool status;
    1534
    1635    psTimerStart ("psphot.rough");
     36
     37    // find the currently selected readout
     38    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     39    psAssert (file, "missing file?");
     40
     41    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     42    psAssert (readout, "missing readout?");
     43
     44    // if we have a PSF, use the existing PSF clump region below
     45    bool havePSF = false;
     46    if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
     47        havePSF = true;
     48    }
     49
     50    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     51    psAssert (detections, "missing detections?");
     52
     53    psArray *sources = detections->newSources;
     54    psAssert (sources, "missing sources?");
     55
     56    if (!sources->n) {
     57        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
     58        return true;
     59    }
    1760
    1861    // we make this measurement on a NxM grid of regions across the readout
    … …  
    75118        // determine the PSF parameters from the source moment values
    76119        // XXX why not save the psfClump as a PTR?
    77         psfClump = pmSourcePSFClump (region, sources, recipe);
     120
     121        float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     122        float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
     123
     124        float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");
     125        if (!status) {
     126            PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
     127            psAssert (status, "missing PSF_CLUMP_GRID_SCALE");
     128        }
     129        float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");
     130        if (!status) {
     131            MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
     132            psAssert (status, "missing MOMENTS_SX_MAX");
     133        }
     134        float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");
     135        if (!status) {
     136            MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
     137            psAssert (status, "missing MOMENTS_SY_MAX");
     138        }
     139
     140        psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
     141
    78142        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X",  PS_META_REPLACE, "psf clump center", psfClump.X);
    79143        psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.Y",  PS_META_REPLACE, "psf clump center", psfClump.Y);
    … …  
    103167    psLogMsg ("psphot", 3, "psf clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY);
    104168
     169    // get basic parameters, or set defaults
     170    float PSF_SN_LIM = psMetadataLookupF32 (&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     171    float PSF_CLUMP_NSIGMA = psMetadataLookupF32 (&status, recipe, "PSF_CLUMP_NSIGMA"); psAssert (status, "missing PSF_CLUMP_NSIGMA");
     172
    105173    // group into STAR, COSMIC, EXTENDED, SATURATED, etc.
    106     if (!pmSourceRoughClass (region, sources, recipe, psfClump, maskSat)) {
     174    if (!pmSourceRoughClass (region, sources, PSF_SN_LIM, PSF_CLUMP_NSIGMA, psfClump, maskSat)) {
    107175        psError(PSPHOT_ERR_PROG, false, "programming error calling pmSourceRoughClass");
    108176        return false;

  • branches/eam_branches/20091201/psphot/src/psphotSetThreads.c

    r25755 r26788  
    2525    psFree(task);
    2626
    27     task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 5);
     27    task = psThreadTaskAlloc("PSPHOT_SOURCE_STATS", 10);
    2828    task->function = &psphotSourceStats_Threaded;
    2929    psThreadTaskAdd(task);

  • branches/eam_branches/20091201/psphot/src/psphotSignificanceImage.c

    r26705 r26788  
    2222    }
    2323
     24    // if we have already determined the PSF model, then we have a better idea how to smooth this image
    2425    bool statusMajor, statusMinor;
    25     float fwhmMajor = psMetadataLookupF32(&statusMajor, recipe, "FWHM_MAJ");
    26     float fwhmMinor = psMetadataLookupF32(&statusMinor, recipe, "FWHM_MIN");
     26    float fwhmMajor = psMetadataLookupF32(&statusMajor, readout->analysis, "FWHM_MAJ");
     27    float fwhmMinor = psMetadataLookupF32(&statusMinor, readout->analysis, "FWHM_MIN");
    2728    if (statusMajor && statusMinor) {
    2829        // if we know the FHWM, use that to set the smoothing kernel (XXX allow an optional override?)
    … …  
    4344    }
    4445    // record the actual smoothing sigma
    45     psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections",
    46                      SIGMA_SMTH);
     46    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGMA_SMOOTH", PS_META_REPLACE, "Smoothing sigma for detections", SIGMA_SMTH);
    4747
    4848    // smooth the image, applying the mask as we go
    … …  
    5959    // effective per-pixel variance is maintained.
    6060    psImage *smooth_wt = psImageCopy(NULL, readout->variance, PS_TYPE_F32);
    61     psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2,
    62                       NSIGMA_SMTH, minGauss);
     61    psImageSmoothMask_Threaded(smooth_wt, smooth_wt, readout->mask, maskVal, SIGMA_SMTH * M_SQRT1_2, NSIGMA_SMTH, minGauss);
    6362    psLogMsg("psphot", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("psphot.smooth"));
    6463
    … …  
    8281    // record the effective area and significance scaling factor
    8382    float effArea = 8.0 * M_PI * PS_SQR(SIGMA_SMTH);
    84     psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
    85     psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE,
    86                      "Signicance scale factor", factor);
     83    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
     84    psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
    8785
    8886    // XXX multithread this if needed

  • branches/eam_branches/20091201/psphot/src/psphotSkyReplace.c

    r21183 r26788  
    11# include "psphotInternal.h"
     2
     3bool psphotSkyReplace (pmConfig *config, const pmFPAview *view)
     4{
     5    bool status = true;
     6
     7    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     8    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     9
     10    // loop over the available readouts
     11    for (int i = 0; i < num; i++) {
     12        if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) {
     13            psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i);
     14            return false;
     15        }
     16    }
     17    return true;
     18}
    219
    320// XXX make this an option?
    421// in order to  successfully replace the sky, we must define a corresponding file...
    5 bool psphotSkyReplace (pmConfig *config, const pmFPAview *view) {
     22bool psphotSkyReplaceReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index) {
    623
    724    psTimerStart ("psphot.skyreplace");
    825
    926    // find the currently selected readout
    10     pmReadout *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
    11     if (readout == NULL) psAbort("input not defined");
     27    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     28    psAssert (file, "missing file?");
     29
     30    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     31    psAssert (readout, "missing readout?");
    1232
    1333    // select background pixels, from output background file, or create

  • branches/eam_branches/20091201/psphot/src/psphotSourceFreePixels.c

    r12792 r26788  
    11# include "psphotInternal.h"
    22
    3 void psphotSourceFreePixels (psArray *sources) {
     3bool psphotSourceFreePixels (pmConfig *config, const pmFPAview *view)
     4{
     5    bool status = true;
    46
    5     if (!sources) return;
     7    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     8    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     9
     10    // loop over the available readouts
     11    for (int i = 0; i < num; i++) {
     12        if (!psphotSourceFreePixelsReadout (config, view, "PSPHOT.INPUT", i)) {
     13            psError (PSPHOT_ERR_CONFIG, false, "failed to free source pixels for PSPHOT.INPUT entry %d", i);
     14            return false;
     15        }
     16    }
     17    return true;
     18}
     19
     20bool psphotSourceFreePixelsReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index) {
     21
     22    bool status;
     23
     24    // find the currently selected readout
     25    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     26    psAssert (file, "missing file?");
     27
     28    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     29    psAssert (readout, "missing readout?");
     30
     31    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     32    psAssert (detections, "missing detections?");
     33
     34    psArray *sources = detections->allSources;
     35    psAssert (sources, "missing sources?");
    636
    737    for (int i = 0; i < sources->n; i++) {
    … …  
    939        pmSourceFreePixels (source);
    1040    }
    11     return;
     41    return true;
    1242}

  • branches/eam_branches/20091201/psphot/src/psphotSourceSize.c

    r26646 r26788  
    2020bool psphotSourceClassRegion (psRegion *region, pmPSFClump *psfClump, psArray *sources, psMetadata *recipe, pmPSF *psf, psphotSourceSizeOptions *options);
    2121bool psphotSourceSizeCR (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options);
    22 bool psphotMaskCosmicRay (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    23 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
    24 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
    25 
    26 bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
    27 
    28     for (int i = 0; i < sources->n; i++) {
    29         pmSource *source = sources->data[i];
    30         pmPeak *peak = source->peak;
    31         pmFootprint *footprint = peak->footprint;
    32         if (!footprint) continue;
    33         for (int j = 0; j < footprint->spans->n; j++) {
    34             pmSpan *sp = footprint->spans->data[j];
    35             psAssert (sp, "missing span");
    36         }
    37     }
    38     return true;
    39 }
     22bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal);
     23bool psphotMaskCosmicRayFootprintCheck (psArray *sources);
    4024
    4125// we need to call this function after sources have been fitted to the PSF model and
    … …  
    4529// deviation from the psf model at the r = FWHM/2 position
    4630
    47 // XXX use an internal flag to mark sources which have already been measured
    48 bool psphotSourceSize(pmConfig *config, pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, long first)
     31// for now, let's store the detections on the readout->analysis for each readout
     32bool psphotSourceSize (pmConfig *config, const pmFPAview *view, bool getPSFsize)
     33{
     34    bool status = true;
     35
     36    // select the appropriate recipe information
     37    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     38    psAssert (recipe, "missing recipe?");
     39
     40    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     41    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     42
     43    // loop over the available readouts
     44    for (int i = 0; i < num; i++) {
     45        if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) {
     46            psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i);
     47            return false;
     48        }
     49    }
     50    return true;
     51}
     52
     53// this function use an internal flag to mark sources which have already been measured
     54bool psphotSourceSizeReadout(pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool getPSFsize)
    4955{
    5056    bool status;
    … …  
    5258
    5359    psTimerStart ("psphot.size");
     60
     61    // find the currently selected readout
     62    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     63    psAssert (file, "missing file?");
     64
     65    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     66    psAssert (readout, "missing readout?");
     67
     68    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     69    psAssert (detections, "missing detections?");
     70
     71    psArray *sources = detections->allSources;
     72    psAssert (sources, "missing sources?");
     73
     74    if (!sources->n) {
     75        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
     76        return true;
     77    }
     78
     79    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
     80    psAssert (psf, "missing psf?");
    5481
    5582    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    91118    // XXX move this to the code that generates the PSF?
    92119    // XXX store the results on pmPSF?
    93     psphotSourceSizePSF (&options, sources, psf);
     120
     121    // XXX this should only be done on the first pass (ie, if we have newSources or allSources?)
     122    if (getPSFsize) {
     123        psphotSourceSizePSF (&options, sources, psf);
     124    }
    94125
    95126    // classify the sources based on ApResid and Moments (extended sources)
     127    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    96128    psphotSourceClass(readout, sources, recipe, psf, &options);
    97129
     130    // NOTE: only sources not already measured !(source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED)
    98131    psphotSourceSizeCR (readout, sources, &options);
    99132
    100     psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n - first, psTimerMark ("psphot.size"));
     133    // XXX fix this (was source->n  - first)
     134    psLogMsg ("psphot.size", PS_LOG_INFO, "measure source sizes for %ld sources: %f sec\n", sources->n, psTimerMark ("psphot.size"));
    101135
    102136    psphotVisualPlotSourceSize (recipe, readout->analysis, sources);
    103137    psphotVisualShowSourceSize (readout, sources);
    104138    psphotVisualPlotApResid (sources, options.ApResid, options.ApSysErr);
    105 
    106     return true;
    107 }
    108 
    109 // This attempt to mask the cosmic rays used the isophotal boundary
    110 bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    111 
    112     // replace the source flux
    113     pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    114 
    115     // flag this as a CR
    116     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    117     pmPeak *peak = source->peak;
    118     psAssert (peak, "NULL peak");
    119 
    120     // grab the matching footprint
    121     pmFootprint *footprint = peak->footprint;
    122     if (!footprint) {
    123       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    124      
    125         // if we have not footprint, use the old code to mask by isophot
    126         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    127         return true;
    128     }
    129 
    130     if (!footprint->spans) {
    131       psTrace("psphot.czw",2,"Using isophot CR mask code.");
    132      
    133         // if we have no footprint, use the old code to mask by isophot
    134         psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    135         return true;
    136     }
    137     psphotMaskCosmicRayIsophot (source, maskVal, crMask);
    138     // mask all of the pixels covered by the spans of the footprint
    139     for (int j = 1; j < footprint->spans->n; j++) {
    140         pmSpan *span1 = footprint->spans->data[j];
    141 
    142         int iy = span1->y;
    143         int xs = span1->x0;
    144         int xe = span1->x1;
    145 
    146         for (int ix = xs; ix < xe; ix++) {
    147             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    148         }
    149     }
    150     return true;
    151 }
    152 
    153 bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
    154 
    155     source->mode |= PM_SOURCE_MODE_CR_LIMIT;
    156     pmPeak *peak = source->peak;
    157     psAssert (peak, "NULL peak");
    158 
    159     psImage *mask   = source->maskView;
    160     psImage *pixels = source->pixels;
    161     psImage *variance = source->variance;
    162 
    163     // XXX This should be a recipe variable
    164 # define SN_LIMIT 5.0
    165 
    166     int xo = peak->x - pixels->col0;
    167     int yo = peak->y - pixels->row0;
    168 
    169     // mark the pixels in this row to the left, then the right
    170     for (int ix = xo; ix >= 0; ix--) {
    171         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    172         if (SN > SN_LIMIT) {
    173             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    174         }
    175     }
    176     for (int ix = xo + 1; ix < pixels->numCols; ix++) {
    177         float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
    178         if (SN > SN_LIMIT) {
    179             mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
    180         }
    181     }
    182 
    183     // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
    184     // first go up:
    185     for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
    186         // mark the pixels in this row to the left, then the right
    187         for (int ix = 0; ix < pixels->numCols; ix++) {
    188             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    189             if (SN < SN_LIMIT) continue;
    190 
    191             bool valid = false;
    192             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
    193             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
    194             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
    195 
    196             if (!valid) continue;
    197             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    198         }
    199     }
    200     // next go down:
    201     for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
    202         // mark the pixels in this row to the left, then the right
    203         for (int ix = 0; ix < pixels->numCols; ix++) {
    204             float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
    205             if (SN < SN_LIMIT) continue;
    206 
    207             bool valid = false;
    208             valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
    209             valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
    210             valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
    211 
    212             if (!valid) continue;
    213             mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
    214         }
    215     }
     139    psphotVisualShowSatStars (recipe, psf, sources);
     140
    216141    return true;
    217142}
    … …  
    254179        float dMag = source->psfMag - apMag;
    255180
    256         psVectorAppend (Ap, 100, dMag);
    257         psVectorAppend (ApErr, 100, source->errMag);
     181        psVectorAppend (Ap, dMag);
     182        psVectorAppend (ApErr, source->errMag);
    258183    }
    259184
    … …  
    454379}
    455380
    456 // given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
    457 // deviation in sigmas.  This is measured on the residual image - should we ignore negative
    458 // deviations?  NOTE: This function was an early attempt to classify extended objects, and is
    459 // no longer used by psphot.
    460 float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
    461                          psImageMaskType maskVal, const pmModel *model, float Ro)
    462 {
    463     psF32 *PAR = model->params->data.F32; // Model parameters
    464     float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
    465 
    466     // We treat the contour as an ellipse:
    467     // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
    468     // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
    469     // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
    470     // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
    471     // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
    472 
    473     // min/max value of x is where B^2 - 4AC = 0; solve this for x
    474     float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
    475     if (Q < 0.0) {
    476         // ellipse is imaginary
    477         return NAN;
    478     }
    479 
    480     int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
    481     int nPts = 0;                       // Number of points in ellipse
    482     float nSigma = 0.0;                 //
    483 
    484     for (int x = -radius; x <= radius; x++) {
    485         // Polynomial coefficients
    486         // XXX Should we be using the centre of the pixel as x or x+0.5?
    487         float A = PS_SQR (1.0 / syy);
    488         float B = x * sxy;
    489         float C = PS_SQR (x / sxx) - Ro;
    490         float T = PS_SQR(B) - 4*A*C;
    491         if (T < 0.0) {
    492             continue;
    493         }
    494 
    495         // y position in source frame
    496         float yP = (-B + sqrt (T)) / (2.0 * A);
    497         float yM = (-B - sqrt (T)) / (2.0 * A);
    498 
    499         // Get the closest pixel positions (image frame)
    500         int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
    501         int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    502         int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
    503 
    504         if (xPix < 0 || xPix >= image->numCols) {
    505             continue;
    506         }
    507 
    508         if (yPixM >= 0 && yPixM < image->numRows &&
    509             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
    510             float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
    511             nSigma += dSigma;
    512             nPts++;
    513         }
    514 
    515         if (yPixM == yPixP) {
    516             continue;
    517         }
    518 
    519         if (yPixP >= 0 && yPixP < image->numRows &&
    520             !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
    521             float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
    522             nSigma += dSigma;
    523             nPts++;
    524         }
    525     }
    526     nSigma /= nPts;
    527     return nSigma;
    528 }
    529 
    530381// given an object suspected to be a defect, generate a pixel mask using the Lapacian transform
    531382// if enough of the object is detected as 'sharp', consider the object a cosmic ray
    … …  
    534385    for (int i = 0; i < sources->n; i++) {
    535386        pmSource *source = sources->data[i];
     387
     388        // skip source if it was already measured
     389        if (source->tmpFlags & PM_SOURCE_TMPF_SIZE_MEASURED) {
     390            psTrace("psphot", 7, "Not calculating source size since it has already been measured\n");
     391            continue;
     392        }
    536393
    537394        // Integer position of peak
    … …  
    582439// does no repair or recovery of the CR pixels, it only masks them out.  My test code can be
    583440// found at /data/ipp031.0/watersc1/psphot.20091209/algo_check.c
    584 bool psphotMaskCosmicRayCZW (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
     441bool psphotMaskCosmicRay (pmReadout *readout, pmSource *source, psImageMaskType maskVal) {
    585442
    586443    // Get the actual images and information about the peak.
    … …  
    760617}
    761618
     619bool psphotMaskCosmicRayFootprintCheck (psArray *sources) {
     620
     621    for (int i = 0; i < sources->n; i++) {
     622        pmSource *source = sources->data[i];
     623        pmPeak *peak = source->peak;
     624        pmFootprint *footprint = peak->footprint;
     625        if (!footprint) continue;
     626        for (int j = 0; j < footprint->spans->n; j++) {
     627            pmSpan *sp = footprint->spans->data[j];
     628            psAssert (sp, "missing span");
     629        }
     630    }
     631    return true;
     632}
     633
     634/**** ------ old versions of cosmic ray masking ----- ****/
     635
     636bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask);
     637
     638// This attempt to mask the cosmic rays used the isophotal boundary
     639bool psphotMaskCosmicRay_V1 (psImage *mask, pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     640
     641    // replace the source flux
     642    pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     643
     644    // flag this as a CR
     645    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     646    pmPeak *peak = source->peak;
     647    psAssert (peak, "NULL peak");
     648
     649    // grab the matching footprint
     650    pmFootprint *footprint = peak->footprint;
     651    if (!footprint) {
     652      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     653     
     654        // if we have not footprint, use the old code to mask by isophot
     655        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     656        return true;
     657    }
     658
     659    if (!footprint->spans) {
     660      psTrace("psphot.czw",2,"Using isophot CR mask code.");
     661     
     662        // if we have no footprint, use the old code to mask by isophot
     663        psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     664        return true;
     665    }
     666    psphotMaskCosmicRayIsophot (source, maskVal, crMask);
     667    // mask all of the pixels covered by the spans of the footprint
     668    for (int j = 1; j < footprint->spans->n; j++) {
     669        pmSpan *span1 = footprint->spans->data[j];
     670
     671        int iy = span1->y;
     672        int xs = span1->x0;
     673        int xe = span1->x1;
     674
     675        for (int ix = xs; ix < xe; ix++) {
     676            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     677        }
     678    }
     679    return true;
     680}
     681
     682bool psphotMaskCosmicRayIsophot (pmSource *source, psImageMaskType maskVal, psImageMaskType crMask) {
     683
     684    source->mode |= PM_SOURCE_MODE_CR_LIMIT;
     685    pmPeak *peak = source->peak;
     686    psAssert (peak, "NULL peak");
     687
     688    psImage *mask   = source->maskView;
     689    psImage *pixels = source->pixels;
     690    psImage *variance = source->variance;
     691
     692    // XXX This should be a recipe variable
     693# define SN_LIMIT 5.0
     694
     695    int xo = peak->x - pixels->col0;
     696    int yo = peak->y - pixels->row0;
     697
     698    // mark the pixels in this row to the left, then the right
     699    for (int ix = xo; ix >= 0; ix--) {
     700        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     701        if (SN > SN_LIMIT) {
     702            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     703        }
     704    }
     705    for (int ix = xo + 1; ix < pixels->numCols; ix++) {
     706        float SN = pixels->data.F32[yo][ix] / sqrt(variance->data.F32[yo][ix]);
     707        if (SN > SN_LIMIT) {
     708            mask->data.PS_TYPE_IMAGE_MASK_DATA[yo][ix] |= crMask;
     709        }
     710    }
     711
     712    // for each of the neighboring rows, mark the high pixels if they have a marked neighbor
     713    // first go up:
     714    for (int iy = PS_MIN(yo, mask->numRows-2); iy >= 0; iy--) {
     715        // mark the pixels in this row to the left, then the right
     716        for (int ix = 0; ix < pixels->numCols; ix++) {
     717            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     718            if (SN < SN_LIMIT) continue;
     719
     720            bool valid = false;
     721            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix] & crMask);
     722            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix-1] & crMask) : 0;
     723            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy+1][ix+1] & crMask) : 0;
     724
     725            if (!valid) continue;
     726            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     727        }
     728    }
     729    // next go down:
     730    for (int iy = PS_MIN(yo+1, mask->numRows-1); iy < pixels->numRows; iy++) {
     731        // mark the pixels in this row to the left, then the right
     732        for (int ix = 0; ix < pixels->numCols; ix++) {
     733            float SN = pixels->data.F32[iy][ix] / sqrt(variance->data.F32[iy][ix]);
     734            if (SN < SN_LIMIT) continue;
     735
     736            bool valid = false;
     737            valid |= (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix] & crMask);
     738            valid |= (ix > 0) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix-1] & crMask) : 0;
     739            valid |= (ix <= mask->numCols) ? (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy-1][ix+1] & crMask) : 0;
     740
     741            if (!valid) continue;
     742            mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] |= crMask;
     743        }
     744    }
     745    return true;
     746}
     747
     748// given the PSF ellipse parameters, navigate around the 1sigma contour, return the total
     749// deviation in sigmas.  This is measured on the residual image - should we ignore negative
     750// deviations?  NOTE: This function was an early attempt to classify extended objects, and is
     751// no longer used by psphot.
     752float psphotModelContour(const psImage *image, const psImage *variance, const psImage *mask,
     753                         psImageMaskType maskVal, const pmModel *model, float Ro)
     754{
     755    psF32 *PAR = model->params->data.F32; // Model parameters
     756    float sxx = PAR[PM_PAR_SXX], sxy = PAR[PM_PAR_SXY], syy = PAR[PM_PAR_SYY]; // Ellipse parameters
     757
     758    // We treat the contour as an ellipse:
     759    // Ro = (x / SXX)^2 + (y / SYY)^2 + x y SXY
     760    // y^2 (1/SYY^2) + y (x SXY) + (x / SXX)^2 - Ro = 0;
     761    // This is a quadratic, Ay^2 + By + C with A = 1/SYY^2, B = x*SXY, C = (x / SXX)^2 - Ro
     762    // The solution is y = [-B +/- sqrt (B^2 - 4 A C)] / [2 A], so:
     763    // y = [-(x SXY) +/- sqrt ((x SXY)^2 - 4 (1/SYY^2) ((x/SXX)^2 - Ro))] * [SYY^2 / 2]
     764
     765    // min/max value of x is where B^2 - 4AC = 0; solve this for x
     766    float Q = Ro * PS_SQR(sxx) / (1.0 - PS_SQR(sxx * syy * sxy) / 4.0);
     767    if (Q < 0.0) {
     768        // ellipse is imaginary
     769        return NAN;
     770    }
     771
     772    int radius = sqrtf(Q) + 0.5;        // Radius of ellipse
     773    int nPts = 0;                       // Number of points in ellipse
     774    float nSigma = 0.0;                 //
     775
     776    for (int x = -radius; x <= radius; x++) {
     777        // Polynomial coefficients
     778        // XXX Should we be using the centre of the pixel as x or x+0.5?
     779        float A = PS_SQR (1.0 / syy);
     780        float B = x * sxy;
     781        float C = PS_SQR (x / sxx) - Ro;
     782        float T = PS_SQR(B) - 4*A*C;
     783        if (T < 0.0) {
     784            continue;
     785        }
     786
     787        // y position in source frame
     788        float yP = (-B + sqrt (T)) / (2.0 * A);
     789        float yM = (-B - sqrt (T)) / (2.0 * A);
     790
     791        // Get the closest pixel positions (image frame)
     792        int xPix  = x  + PAR[PM_PAR_XPOS] - image->col0 + 0.5;
     793        int yPixM = yM + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     794        int yPixP = yP + PAR[PM_PAR_YPOS] - image->row0 + 0.5;
     795
     796        if (xPix < 0 || xPix >= image->numCols) {
     797            continue;
     798        }
     799
     800        if (yPixM >= 0 && yPixM < image->numRows &&
     801            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixM][xPix] & maskVal))) {
     802            float dSigma = image->data.F32[yPixM][xPix] / sqrtf(variance->data.F32[yPixM][xPix]);
     803            nSigma += dSigma;
     804            nPts++;
     805        }
     806
     807        if (yPixM == yPixP) {
     808            continue;
     809        }
     810
     811        if (yPixP >= 0 && yPixP < image->numRows &&
     812            !(mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[yPixP][xPix] & maskVal))) {
     813            float dSigma = image->data.F32[yPixP][xPix] / sqrtf(variance->data.F32[yPixP][xPix]);
     814            nSigma += dSigma;
     815            nPts++;
     816        }
     817    }
     818    nSigma /= nPts;
     819    return nSigma;
     820}
     821
    762822// this was an old attempt to identify cosmic rays based on the peak curvature
    763823bool psphotSourcePeakCurvature (pmReadout *readout, psArray *sources, psphotSourceSizeOptions *options) {
    … …  
    893953    return true;
    894954}
     955

  • branches/eam_branches/20091201/psphot/src/psphotSourceStats.c

    r26596 r26788  
    11# include "psphotInternal.h"
    22
    3 bool psphotSetMomentsWindow (psMetadata *recipe, psMetadata *analysis, psArray *sources);
    4 
    5 psArray *psphotSourceStats (pmConfig *config, pmReadout *readout, pmDetections *detections, bool setWindow) {
     3// convert detections to sources and measure their basic properties (moments, local sky, sky
     4// variance) Note: this function only generates sources for the new peaks (peak->assigned)
     5bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow)
     6{
     7    bool status = true;
     8
     9    // select the appropriate recipe information
     10    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     11    psAssert (recipe, "missing recipe?");
     12
     13    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
     14    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
     15
     16    // loop over the available readouts
     17    for (int i = 0; i < num; i++) {
     18        if (!psphotSourceStatsReadout (config, view, "PSPHOT.INPUT", i, recipe, setWindow)) {
     19            psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for PSPHOT.INPUT entry %d", i);
     20            return false;
     21        }
     22    }
     23    return true;
     24}
     25
     26bool psphotSourceStatsReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe, bool setWindow) {
    627
    728    bool status = false;
    … …  
    1031    psTimerStart ("psphot.stats");
    1132
    12     // select the appropriate recipe information
    13     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    14     assert (recipe);
     33    // find the currently selected readout
     34    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest
     35    psAssert (file, "missing file?");
     36
     37    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     38    psAssert (readout, "missing readout?");
     39
     40    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     41    psAssert (detections, "missing detections?");
     42    psAssert (!detections->newSources, "new sources already defined?");
     43
     44    // XXX TEST:
     45    if (detections->allSources) {
     46        psphotMaskCosmicRayFootprintCheck(detections->allSources);
     47    }
     48    if (detections->newSources) {
     49        psphotMaskCosmicRayFootprintCheck(detections->newSources);
     50    }
    1551
    1652    // determine the number of allowed threads
    … …  
    2157
    2258    // determine properties (sky, moments) of initial sources
    23     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    24     if (!status) return NULL;
    25 
     59    float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     60    psAssert (status, "missing SKY_OUTER_RADIUS in recipe?");
     61
     62    // XXX this seems like an arbitrary number...
    2663    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
    2764
    2865    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    29     if (!status) return NULL;
     66    psAssert (status, "missing BREAK_POINT?");
     67
     68    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS"); psAssert (status, "missing SKY_INNER_RADIUS");
     69    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
     70
     71    // bit-masks to test for good/bad pixels
     72    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
     73    psAssert (maskVal, "missing MASK.PSPHOT");
     74
     75    // bit-mask to mark pixels not used in analysis
     76    psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
     77    psAssert (markVal, "missing MARK.PSPHOT");
    3078
    3179    psArray *peaks = detections->peaks;
    3280    if (!peaks) {
    3381        psError(PS_ERR_UNEXPECTED_NULL, false, "No peaks found!");
    34         return NULL;
     82        return false;
    3583    }
    3684
    3785    // generate the array of sources, define the associated pixel
    3886    sources = psArrayAllocEmpty (peaks->n);
     87
     88    // if there are no peaks, we save the empty source array and return
     89    if (!peaks->n) {
     90        // save the new sources on the detection structure:
     91        detections->newSources = sources;
     92        return true;
     93    }
    3994
    4095    for (int i = 0; i < peaks->n; i++) {
    … …  
    58113        psArrayAdd (sources, 100, source);
    59114        psFree (source);
     115    }
     116
     117    if (!strcasecmp (breakPt, "PEAKS")) {
     118        psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
     119        psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
     120        psphotVisualShowMoments (sources);
     121        detections->newSources = sources;
     122        return true;
     123    }
     124
     125    if (setWindow) {
     126        if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
     127            psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
     128            psFree(sources);
     129            return false;
     130        }
     131    }
     132
     133    // if we have measured the window, we will be saving the modified version of these recipe values on readout->analysis
     134    float SIGMA = psMetadataLookupF32 (&status, readout->analysis, "MOMENTS_GAUSS_SIGMA");
     135    if (!status) {
     136        SIGMA = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
     137    }
     138    float RADIUS = psMetadataLookupF32 (&status, readout->analysis, "PSF_MOMENTS_RADIUS");
     139    if (!status) {
     140        RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     141    }
     142
     143    // threaded measurement of the source magnitudes
     144    int Nfail = 0;
     145    int Nmoments = 0;
     146    int Nfaint = 0;
     147
     148    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
     149    int Cx = 1, Cy = 1;
     150    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
     151
     152    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
     153
     154    for (int i = 0; i < cellGroups->n; i++) {
     155
     156        psArray *cells = cellGroups->data[i];
     157
     158        for (int j = 0; j < cells->n; j++) {
     159
     160            // allocate a job -- if threads are not defined, this just runs the job
     161            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
     162
     163            psArrayAdd(job->args, 1, cells->data[j]); // sources
     164
     165            PS_ARRAY_ADD_SCALAR(job->args, INNER,   PS_TYPE_F32);
     166            PS_ARRAY_ADD_SCALAR(job->args, MIN_SN,  PS_TYPE_F32);
     167            PS_ARRAY_ADD_SCALAR(job->args, RADIUS,  PS_TYPE_F32);
     168            PS_ARRAY_ADD_SCALAR(job->args, SIGMA,   PS_TYPE_F32);
     169
     170            PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
     171            PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
     172
     173            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
     174            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
     175            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
     176
     177            if (!psThreadJobAddPending(job)) {
     178                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
     179                psFree (job);
     180                psFree(sources);
     181                return false;
     182            }
     183            psFree(job);
     184        }
     185
     186        // wait for the threads to finish and manage results
     187        if (!psThreadPoolWait (false)) {
     188            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
     189            psFree(sources);
     190            return false;
     191        }
     192
     193        // we have only supplied one type of job, so we can assume the types here
     194        psThreadJob *job = NULL;
     195        while ((job = psThreadJobGetDone()) != NULL) {
     196            if (job->args->n < 1) {
     197                fprintf (stderr, "error with job\n");
     198            } else {
     199                psScalar *scalar = NULL;
     200                scalar = job->args->data[7];
     201                Nmoments += scalar->data.S32;
     202                scalar = job->args->data[8];
     203                Nfail += scalar->data.S32;
     204                scalar = job->args->data[9];
     205                Nfaint += scalar->data.S32;
     206            }
     207            psFree(job);
     208        }
     209    }
     210
     211    psFree (cellGroups);
     212
     213    psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
     214
     215    psphotVisualShowMoments (sources);
     216
     217    // save the new sources on the detection structure:
     218    detections->newSources = sources;
     219
     220
     221    if (detections->allSources) {
     222        psphotMaskCosmicRayFootprintCheck(detections->allSources);
     223    }
     224    if (detections->newSources) {
     225        psphotMaskCosmicRayFootprintCheck(detections->newSources);
     226    }
     227
     228    return true;
     229}
     230
     231// this function is currently only called by psphotCheckExtSources
     232bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
     233
     234    bool status = false;
     235
     236    psTimerStart ("psphot.stats");
     237
     238    // select the appropriate recipe information
     239    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     240    assert (recipe);
     241
     242    // determine the number of allowed threads
     243    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     244    if (!status) {
     245        nThreads = 0;
     246    }
     247
     248    // determine properties (sky, moments) of initial sources
     249    float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
     250    if (!status) return false;
     251
     252    OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
     253
     254    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
     255    if (!status) return false;
     256
     257    for (int i = 0; i < sources->n; i++) {
     258
     259        pmSource *source = sources->data[i];
     260        if (!source->peak) continue; // XXX how can we have a peak-less source?
     261
     262        // allocate space for moments
     263        if (!source->moments) {
     264            source->moments = pmMomentsAlloc();
     265        }
     266
     267        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
     268        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
     269        source->peak->assigned = true;
    60270    }
    61271
    … …  
    67277    }
    68278
    69     if (setWindow) {
    70         if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
    71             psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
    72             return NULL;
    73         }
     279    // XXX how else could we get the window size in?
     280    if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
     281        psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
     282        return NULL;
    74283    }
    75284
    … …  
    94303            psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
    95304
     305            // XXX: this must match the above
    96306            psArrayAdd(job->args, 1, cells->data[j]); // sources
    97307            psArrayAdd(job->args, 1, recipe);
    … …  
    141351}
    142352
    143 bool psphotSourceStatsUpdate (psArray *sources, pmConfig *config, pmReadout *readout) {
    144 
    145     bool status = false;
    146 
    147     psTimerStart ("psphot.stats");
    148 
    149     // select the appropriate recipe information
    150     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    151     assert (recipe);
    152 
    153     // determine the number of allowed threads
    154     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    155     if (!status) {
    156         nThreads = 0;
    157     }
    158 
    159     // determine properties (sky, moments) of initial sources
    160     float OUTER    = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS");
    161     if (!status) return NULL;
    162 
    163     OUTER = PS_MAX(OUTER, 20.0); // XXX Guarantee that we can encompass the max moments radius
    164 
    165     char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    166     if (!status) return NULL;
    167 
    168     for (int i = 0; i < sources->n; i++) {
    169 
    170         pmSource *source = sources->data[i];
    171         if (!source->peak) continue; // XXX how can we have a peak-less source?
    172 
    173         // allocate space for moments
    174         if (!source->moments) {
    175             source->moments = pmMomentsAlloc();
    176         }
    177 
    178         // allocate image, weight, mask arrays for each peak (square of radius OUTER)
    179         pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    180         source->peak->assigned = true;
    181     }
    182 
    183     if (!strcasecmp (breakPt, "PEAKS")) {
    184         psLogMsg ("psphot", PS_LOG_INFO, "%ld sources : %f sec\n", sources->n, psTimerMark ("psphot.stats"));
    185         psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n");
    186         psphotVisualShowMoments (sources);
    187         return sources;
    188     }
    189 
    190     // XXX how else could we get the window size in?
    191     if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) {
    192         psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!");
    193         return NULL;
    194     }
    195 
    196     // threaded measurement of the source magnitudes
    197     int Nfail = 0;
    198     int Nmoments = 0;
    199     int Nfaint = 0;
    200 
    201     // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
    202     int Cx = 1, Cy = 1;
    203     psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
    204 
    205     psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
    206 
    207     for (int i = 0; i < cellGroups->n; i++) {
    208 
    209         psArray *cells = cellGroups->data[i];
    210 
    211         for (int j = 0; j < cells->n; j++) {
    212 
    213             // allocate a job -- if threads are not defined, this just runs the job
    214             psThreadJob *job = psThreadJobAlloc ("PSPHOT_SOURCE_STATS");
    215 
    216             psArrayAdd(job->args, 1, cells->data[j]); // sources
    217             psArrayAdd(job->args, 1, recipe);
    218             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nmoments
    219             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
    220             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfaint
    221 
    222             if (!psThreadJobAddPending(job)) {
    223                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    224                 psFree (job);
    225                 return NULL;
    226             }
    227             psFree(job);
    228         }
    229 
    230         // wait for the threads to finish and manage results
    231         if (!psThreadPoolWait (false)) {
    232             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    233             return NULL;
    234         }
    235 
    236         // we have only supplied one type of job, so we can assume the types here
    237         psThreadJob *job = NULL;
    238         while ((job = psThreadJobGetDone()) != NULL) {
    239             if (job->args->n < 1) {
    240                 fprintf (stderr, "error with job\n");
    241             } else {
    242                 psScalar *scalar = NULL;
    243                 scalar = job->args->data[2];
    244                 Nmoments += scalar->data.S32;
    245                 scalar = job->args->data[3];
    246                 Nfail += scalar->data.S32;
    247                 scalar = job->args->data[4];
    248                 Nfaint += scalar->data.S32;
    249             }
    250             psFree(job);
    251         }
    252     }
    253 
    254     psFree (cellGroups);
    255 
    256     psLogMsg ("psphot", PS_LOG_INFO, "%ld sources, %d moments, %d faint, %d failed: %f sec\n", sources->n, Nmoments, Nfaint, Nfail, psTimerMark ("psphot.stats"));
    257 
    258     psphotVisualShowMoments (sources);
    259 
    260     return (sources);
    261 }
    262 
    263353bool psphotSourceStats_Threaded (psThreadJob *job) {
    264354
    … …  
    268358
    269359    psArray *sources                = job->args->data[0];
    270     psMetadata *recipe              = job->args->data[1];
    271 
    272     float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    273     if (!status) return false;
    274     float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    275     if (!status) return false;
    276     float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
    277     if (!status) return false;
    278     float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
    279     if (!status) return false;
    280 
    281     // bit-masks to test for good/bad pixels
    282     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT");
    283     assert (maskVal);
    284 
    285     // bit-mask to mark pixels not used in analysis
    286     psImageMaskType markVal = psMetadataLookupImageMask(&status, recipe, "MARK.PSPHOT");
    287     assert (markVal);
     360
     361    float INNER                     = PS_SCALAR_VALUE(job->args->data[1],F32);
     362    float MIN_SN                    = PS_SCALAR_VALUE(job->args->data[2],F32);
     363    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[3],F32);
     364    float SIGMA                     = PS_SCALAR_VALUE(job->args->data[4],F32);
     365
     366    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
     367    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
    288368
    289369    // maskVal is used to test for rejected pixels, and must include markVal
    … …  
    349429
    350430    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    351     scalar = job->args->data[2];
     431    scalar = job->args->data[7];
    352432    scalar->data.S32 = Nmoments;
    353433
    354     scalar = job->args->data[3];
     434    scalar = job->args->data[8];
    355435    scalar->data.S32 = Nfail;
    356436
    357     scalar = job->args->data[4];
     437    scalar = job->args->data[9];
    358438    scalar->data.S32 = Nfaint;
    359439
    … …  
    367447    bool status;
    368448
    369     float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    370     if (!status) return false;
     449    float MIN_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN"); psAssert (status, "missing MOMENTS_SN_MIN");
     450    float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
     451    psF32 MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
    371452
    372453    // XXX move this to a config file?
    … …  
    398479
    399480        // choose a grid scale that is a fixed fraction of the psf sigma^2
    400         psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(sigma[i]));
    401         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
    402         psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(sigma[i]));
     481        float PSF_CLUMP_GRID_SCALE = 0.1*PS_SQR(sigma[i]);
     482        float MOMENTS_SX_MAX = 2.0*PS_SQR(sigma[i]);
     483        float MOMENTS_SY_MAX = 2.0*PS_SQR(sigma[i]);
    403484
    404485        // determine the PSF parameters from the source moment values
    405         pmPSFClump psfClump = pmSourcePSFClump (NULL, sources, recipe);
     486        pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);
    406487        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]);
    407488
    … …  
    442523
    443524    // choose a grid scale that is a fixed fraction of the psf sigma^2
    444     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
    445     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    446     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
    447     psMetadataAddF32(recipe, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
    448     psMetadataAddF32(recipe, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
     525    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_CLUMP_GRID_SCALE", PS_META_REPLACE, "clump grid", 0.1*PS_SQR(Sigma));
     526    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SX_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     527    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_SY_MAX", PS_META_REPLACE, "moments limit", 2.0*PS_SQR(Sigma));
     528    psMetadataAddF32(analysis, PS_LIST_TAIL, "MOMENTS_GAUSS_SIGMA", PS_META_REPLACE, "moments limit", Sigma);
     529    psMetadataAddF32(analysis, PS_LIST_TAIL, "PSF_MOMENTS_RADIUS", PS_META_REPLACE, "moments limit", 4.0*Sigma);
    449530
    450531    psLogMsg ("psphot", 3, "using window function with sigma = %f\n", Sigma);
    451532    return true;
    452533}
     534
     535// if we use the footprints, the output peaks list contains both old and new peaks,
     536// otherwise it only contains the new peaks.
     537

  • branches/eam_branches/20091201/psphot/src/psphotSubtractBackground.c

    r26542 r26788  
    44// generate the median in NxN boxes, clipping heavily
    55// linear interpolation to generate full-scale model
    6 bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index)
     6bool psphotSubtractBackgroundReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
    77{
    88    bool status = true;
    … …  
    2525    pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa);
    2626    assert (model);
    27 
    28     // select the appropriate recipe information
    29     psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
    30     assert (recipe);
    3127
    3228    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    … …  
    114110    // the pmReadout selected in this function are all view on entries in config->files
    115111
     112    // display the backsub and backgnd images
     113    // move this inthe the subtract background loop
     114    psphotVisualShowBackground (config, view, readout);
     115
    116116    npass ++;
    117117    return true;
    … …  
    122122    bool status = false;
    123123
     124    // select the appropriate recipe information
     125    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
     126    psAssert (recipe, "missing recipe?");
     127
    124128    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
    125129    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
    … …  
    127131    // loop over the available readouts
    128132    for (int i = 0; i < num; i++) {
    129         if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i)) {
     133        if (!psphotSubtractBackgroundReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
    130134            psError (PSPHOT_ERR_CONFIG, false, "failed to subtract background for PSPHOT.INPUT entry %d", i);
    131135            return false;
Note: See TracChangeset for help on using the changeset viewer.