Index: /trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-simple.mdc (revision 27656) +++ /trunk/ippconfig/recipes/filerules-simple.mdc (revision 27657) @@ -167,5 +167,5 @@ PSPHOT.PSF.RAW.SAVE OUTPUT {OUTPUT}.psf PSF NONE FPA TRUE NONE PSPHOT.PSF.SKY.SAVE OUTPUT {OUTPUT}.psf PSF NONE FPA TRUE NONE - + # outputs for psphotStack: PSPHOT.CHISQ.IMAGE OUTPUT {OUTPUT}.chisq.im.fits IMAGE NONE FPA TRUE SIMPLE Index: /trunk/ppSim/src/ppSimInsertGalaxies.c =================================================================== --- /trunk/ppSim/src/ppSimInsertGalaxies.c (revision 27656) +++ /trunk/ppSim/src/ppSimInsertGalaxies.c (revision 27657) @@ -74,8 +74,13 @@ int dY = PM_CELL_TO_CHIP (0.0, y0Cell, yParityCell, binning); - // psMetadataLookupPtr (readout->analysis, "PSPHOT.SOURCES", 0); - - psArray *sources = psArrayAllocEmpty (galaxies->n); - + pmDetections *detections = psMetadataLookupPtr (&mdok, readout->analysis, "PSPHOT.DETECTIONS"); + if (!detections) { + detections = pmDetectionsAlloc(); + detections->allSources = psArrayAllocEmpty (galaxies->n); + psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_DATA_ARRAY | PS_META_REPLACE, "psphot detections", detections); + } else { + psMemIncrRefCounter (detections); + } + psArray *sources = sources = detections->allSources; // add sources to the readout image & weight @@ -181,8 +186,7 @@ } - // NOTE: readout must be part of the pmFPAfile named "PPSIM.OUTPUT" - // psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, "psphot sources", sources); + // XXX many leaks in here, i think + psFree (detections); - // XXX many leaks in here, i think return true; } Index: /trunk/ppSim/src/ppSimInsertStars.c =================================================================== --- /trunk/ppSim/src/ppSimInsertStars.c (revision 27656) +++ /trunk/ppSim/src/ppSimInsertStars.c (revision 27657) @@ -145,7 +145,13 @@ fclose (outfile); - // NOTE: the pmFPAfile "PPSIM.OUTPUT" points at these sources - psMetadataAdd (readout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY, "fake sources", sources); - psFree(sources); + pmDetections *detections = pmDetectionsAlloc(); + detections->allSources = sources; + + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "fake sources", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree(detections); // XXX many leaks in here, i think Index: /trunk/ppSim/src/ppSimPhotomReadout.c =================================================================== --- /trunk/ppSim/src/ppSimPhotomReadout.c (revision 27656) +++ /trunk/ppSim/src/ppSimPhotomReadout.c (revision 27657) @@ -6,5 +6,8 @@ PS_ASSERT_PTR_NON_NULL (readout, NULL); - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; return sources; } @@ -151,6 +154,4 @@ // create the exported-metadata and free local data - // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? - // (or don't supply the sources, and do this with a different function) psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); @@ -167,6 +168,14 @@ psAssert (forceReadout, "no forceReadout?"); pmChipSetDataStatus (forceChip, true); - psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "forced photometry ", forceSources); - psFree (forceSources); + + pmDetections *detectionsForce = pmDetectionsAlloc(); + detectionForce->allSources = forceSources; + + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "force sources", detectionsForce)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree(detectionsForce); pmCell *fakeCell = pmFPAfileThisCell (config->files, view, "PPSIM.FAKE.SOURCES"); psAssert (fakeCell, "no cell?"); @@ -179,5 +188,14 @@ psAssert (fakeReadout, "no fakeReadout?"); pmChipSetDataStatus (fakeChip, true); - psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + pmDetections *detectionsFake = pmDetectionsAlloc(); + detectionFake->allSources = fakeSources; + + // save detections on the readout->analysis + if (!psMetadataAddPtr (readout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "fake sources", detectionsFake)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + psFree(detectionsFake); return true; Index: /trunk/ppSim/src/ppSimPhotomReadoutFake.c =================================================================== --- /trunk/ppSim/src/ppSimPhotomReadoutFake.c (revision 27656) +++ /trunk/ppSim/src/ppSimPhotomReadoutFake.c (revision 27657) @@ -111,5 +111,10 @@ psAssert (fakeReadout, "no fakeReadout?"); pmChipSetDataStatus (fakeChip, true); - psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "fake photometry ", fakeSources); + + pmDetections *detections = pmDetectionsAlloc(); + detection->allSources = fakeSources; + + psMetadataAddArray (fakeReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE, "fake photometry ", detections); + psFree(detections); return true; Index: /trunk/ppSim/src/ppSimPhotomReadoutForce.c =================================================================== --- /trunk/ppSim/src/ppSimPhotomReadoutForce.c (revision 27656) +++ /trunk/ppSim/src/ppSimPhotomReadoutForce.c (revision 27657) @@ -98,6 +98,4 @@ // create the exported-metadata and free local data - // XXX this places the sources on readout->analysis as PSPHOT.SOURCES. modify? - // (or don't supply the sources, and do this with a different function) psphotReadoutCleanup(config, readout, recipe, NULL, psf, NULL); @@ -111,5 +109,9 @@ psAssert (forceReadout, "no forceReadout?"); pmChipSetDataStatus (forceChip, true); - psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_META_REPLACE, "force photometry ", forceSources); + + pmDetections *detections = pmDetectionsAlloc(); + detection->allSources = forceSources; + psMetadataAddArray (forceReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE, "force photometry ", detections); + psFree(detections); return true; Index: /trunk/ppSim/src/ppSimUtils.c =================================================================== --- /trunk/ppSim/src/ppSimUtils.c (revision 27656) +++ /trunk/ppSim/src/ppSimUtils.c (revision 27657) @@ -270,8 +270,13 @@ psArray *ppSimSelectSources (pmConfig *config, const pmFPAview *view, const char *filename) { + bool status; + pmReadout *readout = pmFPAfileThisReadout (config->files, view, filename); PS_ASSERT_PTR_NON_NULL (readout, NULL); - psArray *sources = psMetadataLookupPtr (NULL, readout->analysis, "PSPHOT.SOURCES"); + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; return sources; } Index: /trunk/psModules/src/camera/pmFPAfile.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.c (revision 27656) +++ /trunk/psModules/src/camera/pmFPAfile.c (revision 27657) @@ -111,5 +111,6 @@ file->save = false; - file->index = fileNum++; + file->fileIndex = fileNum++; + file->fileID = 0; file->imageId = 0; @@ -372,6 +373,14 @@ // Number of the file in list psString num = NULL; // Number to use - psStringAppend(&num, "%d", file->index); + psStringAppend(&num, "%d", file->fileIndex); psStringSubstitute(&newRule, num, "{FILE.INDEX}"); + psFree(num); + } + + if (strstr(newRule, "{FILE.ID}")) { + // Number of the file in list + psString num = NULL; // Number to use + psStringAppend(&num, "%03d", file->fileID); + psStringSubstitute(&newRule, num, "{FILE.ID}"); psFree(num); } @@ -638,5 +647,5 @@ psFree(iter); - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to find instance %d of file %s", num, name); + psLogMsg("psModules.camera", PS_LOG_MINUTIA, "Unable to find instance %d of file %s", num, name); return NULL; } Index: /trunk/psModules/src/camera/pmFPAfile.h =================================================================== --- /trunk/psModules/src/camera/pmFPAfile.h (revision 27656) +++ /trunk/psModules/src/camera/pmFPAfile.h (revision 27657) @@ -111,5 +111,6 @@ psString formatName; // name of the camera format - int index; // Index of file + int fileIndex; // Index of file + int fileID; // internal sequence number psS64 imageId, sourceId; // Image and source identifiers Index: /trunk/psModules/src/camera/pmFPAfileDefine.c =================================================================== --- /trunk/psModules/src/camera/pmFPAfileDefine.c (revision 27656) +++ /trunk/psModules/src/camera/pmFPAfileDefine.c (revision 27657) @@ -1271,6 +1271,11 @@ file->name = psStringCopy (name); + // free a previously existing readout + psFree(file->readout); file->readout = readout; - psMetadataAddPtr(files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN, "", file); + + // allow for multiple entries + // XXX handle replace vs multiple? + psMetadataAddPtr(files, PS_LIST_TAIL, name, PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", file); psFree(file); // we free this copy of file, but 'files' still has a copy @@ -1313,14 +1318,16 @@ const char *name, // name of internal/external file const pmFPA *fpa, // use this fpa to generate - const psImageBinning *binning) { + const psImageBinning *binning, + int index) { pmReadout *readout = NULL; - bool status = true; - pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, name, index); // if the file does not exist, it is not being used as an I/O file: define an internal version if (file == NULL) { - readout = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32); - return readout; + // XXX currently, we do not guarantee that the defined file lands on entry 'index' + psAssert (binning, "internal files must be supplied a psImageBinning for the output images size"); + readout = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32); + return readout; } Index: /trunk/psModules/src/camera/pmFPAfileDefine.h =================================================================== --- /trunk/psModules/src/camera/pmFPAfileDefine.h (revision 27656) +++ /trunk/psModules/src/camera/pmFPAfileDefine.h (revision 27657) @@ -172,5 +172,7 @@ const char *name, // name of internal/external file const pmFPA *fpa, // use this fpa to generate - const psImageBinning *binning); + const psImageBinning *binning, + int index + ); /// @} Index: /trunk/psModules/src/objects/Makefile.am =================================================================== --- /trunk/psModules/src/objects/Makefile.am (revision 27656) +++ /trunk/psModules/src/objects/Makefile.am (revision 27657) @@ -20,4 +20,5 @@ pmModelUtils.c \ pmSource.c \ + pmPhotObj.c \ pmSourceMasks.c \ pmSourceMoments.c \ @@ -80,4 +81,5 @@ pmModelUtils.h \ pmSource.h \ + pmPhotObj.h \ pmSourceMasks.h \ pmSourceDiffStats.h \ Index: /trunk/psModules/src/objects/pmPeaks.h =================================================================== --- /trunk/psModules/src/objects/pmPeaks.h (revision 27656) +++ /trunk/psModules/src/objects/pmPeaks.h (revision 27657) @@ -63,5 +63,5 @@ bool assigned; ///< is peak assigned to a source? pmPeakType type; ///< Description of peak. - pmFootprint *footprint; ///< reference to containing footprint + pmFootprint *footprint; ///< reference to containing footprint } pmPeak; Index: /trunk/psModules/src/objects/pmPhotObj.c =================================================================== --- /trunk/psModules/src/objects/pmPhotObj.c (revision 27656) +++ /trunk/psModules/src/objects/pmPhotObj.c (revision 27657) @@ -17,4 +17,5 @@ #include #include "pmPhotObj.h" +#include "pmSource.h" static void pmPhotObjFree (pmPhotObj *tmp) @@ -38,2 +39,27 @@ } +bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source) { + + psAssert (source, "programming error: NULL source"); + if (!source->peak) { + psError(PS_ERR_UNKNOWN, true, "source missing peak"); + return false; + } + if (!finite(source->peak->xf)) { + psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate"); + return false; + } + if (!finite(source->peak->yf)) { + psError(PS_ERR_UNKNOWN, true, "NAN peak coordinate"); + return false; + } + + // XXX we should probably use the fitted position if it exists + if (!object->sources) { + object->sources = psArrayAllocEmpty(1); + object->x = source->peak->xf; + object->y = source->peak->yf; + } + psArrayAdd (object->sources, 1, source); + return true; +} Index: /trunk/psModules/src/objects/pmPhotObj.h =================================================================== --- /trunk/psModules/src/objects/pmPhotObj.h (revision 27656) +++ /trunk/psModules/src/objects/pmPhotObj.h (revision 27657) @@ -34,8 +34,15 @@ */ typedef struct { - int seq; ///< ID for output (generated on write OR set on read) - psArray *sources; + int id; ///< ID for output (generated on write OR set on read) + psArray *sources; + int flags; + float x; + float y; } pmPhotObj; + +bool pmPhotObjAddSource(pmPhotObj *object, pmSource *source); +pmPhotObj *pmPhotObjAlloc(void); /// @} # endif /* PM_PHOT_OBJ_H */ + Index: /trunk/psModules/src/objects/pmSource.c =================================================================== --- /trunk/psModules/src/objects/pmSource.c (revision 27656) +++ /trunk/psModules/src/objects/pmSource.c (revision 27657) @@ -180,4 +180,5 @@ source->type = in->type; source->mode = in->mode; + source->imageID = in->imageID; return(source); @@ -1058,4 +1059,19 @@ psF32 fA = (A->peak == NULL) ? 0 : A->peak->y; psF32 fB = (B->peak == NULL) ? 0 : B->peak->y; + + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); +} + +// sort by X (ascending) +int pmSourceSortByX (const void **a, const void **b) +{ + pmSource *A = *(pmSource **)a; + pmSource *B = *(pmSource **)b; + + psF32 fA = (A->peak == NULL) ? 0 : A->peak->x; + psF32 fB = (B->peak == NULL) ? 0 : B->peak->x; psF32 diff = fA - fB; Index: /trunk/psModules/src/objects/pmSource.h =================================================================== --- /trunk/psModules/src/objects/pmSource.h (revision 27656) +++ /trunk/psModules/src/objects/pmSource.h (revision 27657) @@ -43,4 +43,5 @@ PM_SOURCE_TMPF_SIZE_MEASURED = 0x0004, PM_SOURCE_TMPF_SIZE_CR_CANDIDATE = 0x0008, + PM_SOURCE_TMPF_MOMENTS_MEASURED = 0x0010, } pmSourceTmpF; @@ -90,4 +91,5 @@ pmSourceExtendedPars *extpars; ///< extended source parameters pmSourceDiffStats *diffStats; ///< extra parameters for difference detections + int imageID; }; @@ -245,4 +247,5 @@ int pmSourceSortBySN (const void **a, const void **b); int pmSourceSortByY (const void **a, const void **b); +int pmSourceSortByX (const void **a, const void **b); int pmSourceSortBySeq (const void **a, const void **b); Index: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 27656) +++ /trunk/psModules/src/psmodules.h (revision 27657) @@ -124,4 +124,5 @@ #include #include +#include #include #include Index: /trunk/psphot/src/Makefile.am =================================================================== --- /trunk/psphot/src/Makefile.am (revision 27656) +++ /trunk/psphot/src/Makefile.am (revision 27657) @@ -17,13 +17,13 @@ # Force recompilation of psphotVersion.c, since it gets the version information psphotVersion.c: psphotVersionDefinitions.h -psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE - -$(RM) psphotVersionDefinitions.h - $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h -FORCE: ; +# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE +# -$(RM) psphotVersionDefinitions.h +# $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h +# FORCE: ; libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotMakePSF +bin_PROGRAMS = psphot psphotForced psphotMakePSF psphotStack # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -39,4 +39,8 @@ psphotMakePSF_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotMakePSF_LDADD = libpsphot.la + +psphotStack_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotStack_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotStack_LDADD = libpsphot.la # psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -80,4 +84,18 @@ psphotMosaicChip.c \ psphotCleanup.c + +# a psphot-variant for stack photometry +psphotStack_SOURCES = \ + psphotStack.c \ + psphotStackArguments.c \ + psphotStackParseCamera.c \ + psphotStackImageLoop.c \ + psphotStackReadout.c \ + psphotStackChisqImage.c \ + psphotFitSourcesLinearStack.c \ + psphotSourceMatch.c \ + psphotCleanup.c + + # # psphot analysis of the detectability of specified positions Index: /trunk/psphot/src/psphot.h =================================================================== --- /trunk/psphot/src/psphot.h (revision 27656) +++ /trunk/psphot/src/psphot.h (revision 27657) @@ -341,3 +341,27 @@ bool psphotMakePSFReadout(pmConfig *config, const pmFPAview *view); +/**** psphotStack prototypes ****/ + +pmConfig *psphotStackArguments(int argc, char **argv); +bool psphotStackParseCamera (pmConfig *config); +bool psphotStackImageLoop (pmConfig *config); +bool psphotStackReadout (pmConfig *config, const pmFPAview *view); +bool psphotStackChisqImage (pmConfig *config, const pmFPAview *view); +bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration + const pmFPAview *view, + pmReadout **chiReadout, + char *filename, + int index); + +bool psphotStackRemoveChisqFromInputs (pmConfig *config); +bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num); + +psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view); +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index); +bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS); + +bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final); +int pmPhotObjSortBySN (const void **a, const void **b); +int pmPhotObjSortByX (const void **a, const void **b); + #endif Index: /trunk/psphot/src/psphotApResid.c =================================================================== --- /trunk/psphot/src/psphotApResid.c (revision 27656) +++ /trunk/psphot/src/psphotApResid.c (revision 27657) @@ -16,6 +16,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotChoosePSF.c =================================================================== --- /trunk/psphot/src/psphotChoosePSF.c (revision 27656) +++ /trunk/psphot/src/psphotChoosePSF.c (revision 27657) @@ -13,6 +13,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotChoosePSFReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to choose a psf model for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotEfficiency.c =================================================================== --- /trunk/psphot/src/psphotEfficiency.c (revision 27656) +++ /trunk/psphot/src/psphotEfficiency.c (revision 27657) @@ -160,6 +160,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotEfficiencyReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed to measure detection efficiency for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27656) +++ /trunk/psphot/src/psphotFitSourcesLinearStack.c (revision 27657) @@ -1,54 +1,16 @@ # include "psphotInternal.h" +// XXX need array of covar factors for each image +// XXX define the 'good' / 'bad' flags? -// for now, let's store the detections on the readout->analysis for each readout -bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final) -{ - bool status = true; +# define COVAR_FACTOR 1.0 - int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); - psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); - - // loop over the available readouts - for (int i = 0; i < num; i++) { - if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) { - psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i); - return false; - } - } - return true; -} - -bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) { +bool psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) { bool status; - float x; - float y; float f; - // float r; // select the appropriate recipe information psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); assert (recipe); - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest - psAssert (file, "missing file?"); - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); - psAssert (readout, "missing readout?"); - - pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psAssert (detections, "missing detections?"); - - psArray *sources = detections->allSources; - psAssert (sources, "missing sources?"); - - if (!sources->n) { - psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit"); - return true; - } - - pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); - psAssert (sources, "missing psf?"); psTimerStart ("psphot.linear"); @@ -65,36 +27,17 @@ maskVal |= markVal; - // source analysis is done in spatial order - sources = psArraySort (sources, pmSourceSortByY); + // analysis is done in spatial order (to speed up overlap search) + // sort by first element in each source list + objects = psArraySort (objects, pmPhotObjSortByX); // storage array for fitSources - psArray *fitSources = psArrayAllocEmpty (sources->n); + psArray *fitSources = psArrayAllocEmpty (objects->n); - // option to limit analysis to a specific region - char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION"); - psRegion AnalysisRegion = psRegionFromString (region); - AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion); - if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined"); + bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); + psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - bool CONSTANT_PHOTOMETRIC_WEIGHTS = - psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS"); - if (!status) { - psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS"); - } - int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER"); - if (!status) { - SKY_FIT_ORDER = 0; - } - bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR"); - if (!status) { - SKY_FIT_LINEAR = false; - } - - // XXX test: choose a larger-than expected radius: - float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix - psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor); - - // XXX do not apply covarFactor for the moment... - // covarFactor = 1.0; + // XXX store a local static array of covar factors for each of the images (by image ID) + // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix + // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor); // select the sources which will be used for the fitting analysis @@ -104,6 +47,6 @@ if (!object->sources) continue; - // check an element of the group to see if we should use it - if (!object->flags & PM_PHOT_OBJ_BAD) continue; + // XXX check an element of the group to see if we should use it + // if (!object->flags & PM_PHOT_OBJ_BAD) continue; for (int j = 0; j < object->sources->n; j++) { @@ -123,5 +66,5 @@ } } - psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n); + psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n); if (fitSources->n == 0) { @@ -143,10 +86,10 @@ // diagonal elements of the sparse matrix (auto-cross-product) - f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseMatrixElement (sparse, i, i, f); // the formal error depends on the weighting scheme if (CONSTANT_PHOTOMETRIC_WEIGHTS) { - float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor); + float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR); errors->data.F32[i] = 1.0 / sqrt(var); } else { @@ -155,5 +98,5 @@ // find the image x model value - f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseVectorElement (sparse, i, f); @@ -162,15 +105,15 @@ pmSource *SRCj = fitSources->data[j]; - // XXX I need to know if this source is on the same image as SRCi -- - if (!sameImge) { continue; } + // we only need to generate dot terms for source on the same image + if (SRCj->imageID != SRCi->imageID) { continue; } // skip over disjoint source images, break after last possible overlap - if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break; - if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; - if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue; - if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; + if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue; // source(i) is above source(j) + if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) continue; // source(i) is below source(j) + if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue; // source(i) is right of source(j) + if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) break; // source(i) is left of source(j) [no other source(j) can overlap source(i)] // got an overlap; calculate cross-product and add to output array - f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor); + f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR); psSparseMatrixElement (sparse, j, i, f); } @@ -190,10 +133,4 @@ psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); - // XXXX **** philosophical question: - // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit: - // should retain the chisq and errors from the intermediate non-linear fit? - // the non-linear fit provides better values for the position errors, and for - // extended sources, the shape errors - // adjust I0 for fitSources and subtract for (int i = 0; i < fitSources->n; i++) { @@ -208,5 +145,4 @@ model->params->data.F32[PM_PAR_I0] = norm->data.F32[i]; model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i]; - // XXX is the value of 'errors' modified by the sky fit? // subtract object @@ -221,5 +157,5 @@ if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor); + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); @@ -229,100 +165,23 @@ psFree (fitSources); psFree (norm); - psFree (skyfit); psFree (errors); - psFree (border); psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear")); - - psphotVisualShowResidualImage (readout); - psphotVisualPlotChisq (sources); - // psphotVisualShowFlags (sources); - - // We have to place this visualization here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); return true; } -// XXX do we need this? -// XXX disallow the simultaneous sky fit and remove this code... +// sort by X (ascending) +int pmPhotObjSortByX (const void **a, const void **b) +{ + pmPhotObj *objA = *(pmPhotObj **)a; + pmPhotObj *objB = *(pmPhotObj **)b; -// Calculate the weight terms for the sky fit component of the matrix. This function operates -// on the pixels which correspond to all of the sources of interest. These elements fill in -// the border matrix components in the sparse matrix equation. -static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) { + psF32 fA = objA->x; + psF32 fB = objB->x; - // generate the image-wide weight terms - // turn on MARK for all image pixels - psRegion fullArray = psRegionSet (0, 0, 0, 0); - fullArray = psRegionForImage (readout->mask, fullArray); - psImageMaskRegion (readout->mask, fullArray, "OR", markVal); - - // turn off MARK for all object pixels - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; - pmModel *model = pmSourceGetModel (NULL, source); - if (model == NULL) continue; - float x = model->params->data.F32[PM_PAR_XPOS]; - float y = model->params->data.F32[PM_PAR_YPOS]; - psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal)); - } - - // accumulate the image statistics from the masked regions - psF32 **image = readout->image->data.F32; - psF32 **variance = readout->variance->data.F32; - psImageMaskType **mask = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA; - - double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy; - w = x = y = x2 = xy = y2 = fo = fx = fy = 0; - - int col0 = readout->image->col0; - int row0 = readout->image->row0; - - for (int j = 0; j < readout->image->numRows; j++) { - for (int i = 0; i < readout->image->numCols; i++) { - if (mask[j][i]) continue; - if (constant_weights) { - wt = 1.0; - } else { - wt = variance[j][i]; - } - f = image[j][i]; - w += 1/wt; - fo += f/wt; - if (SKY_FIT_ORDER == 0) continue; - - xc = i + col0; - yc = j + row0; - x += xc/wt; - y += yc/wt; - x2 += xc*xc/wt; - xy += xc*yc/wt; - y2 += yc*yc/wt; - fx += f*xc/wt; - fy += f*yc/wt; - } - } - - // turn off MARK for all image pixels - psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // set the Border T elements - psSparseBorderElementG (border, 0, fo); - psSparseBorderElementT (border, 0, 0, w); - if (SKY_FIT_ORDER > 0) { - psSparseBorderElementG (border, 0, fx); - psSparseBorderElementG (border, 0, fy); - psSparseBorderElementT (border, 1, 0, x); - psSparseBorderElementT (border, 2, 0, y); - psSparseBorderElementT (border, 0, 1, x); - psSparseBorderElementT (border, 1, 1, x2); - psSparseBorderElementT (border, 2, 1, xy); - psSparseBorderElementT (border, 0, 2, y); - psSparseBorderElementT (border, 1, 2, xy); - psSparseBorderElementT (border, 2, 2, y2); - } - - return true; + psF32 diff = fA - fB; + if (diff > FLT_EPSILON) return (+1); + if (diff < FLT_EPSILON) return (-1); + return (0); } Index: /trunk/psphot/src/psphotGuessModels.c =================================================================== --- /trunk/psphot/src/psphotGuessModels.c (revision 27656) +++ /trunk/psphot/src/psphotGuessModels.c (revision 27657) @@ -15,6 +15,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotGuessModelsReadout (config, view, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotImageLoop.c =================================================================== --- /trunk/psphot/src/psphotImageLoop.c (revision 27656) +++ /trunk/psphot/src/psphotImageLoop.c (revision 27657) @@ -91,15 +91,15 @@ } - status = true; - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); - status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); - if (!status) { - psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); - psFree (view); - return false; - } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } } - // save output which is saved at the chip level if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); Index: /trunk/psphot/src/psphotImageQuality.c =================================================================== --- /trunk/psphot/src/psphotImageQuality.c (revision 27656) +++ /trunk/psphot/src/psphotImageQuality.c (revision 27657) @@ -16,6 +16,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotImageQualityReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on to measure image quality for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotMagnitudes.c =================================================================== --- /trunk/psphot/src/psphotMagnitudes.c (revision 27656) +++ /trunk/psphot/src/psphotMagnitudes.c (revision 27657) @@ -12,6 +12,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image // find the currently selected readout Index: /trunk/psphot/src/psphotMergeSources.c =================================================================== --- /trunk/psphot/src/psphotMergeSources.c (revision 27656) +++ /trunk/psphot/src/psphotMergeSources.c (revision 27657) @@ -76,7 +76,8 @@ pmDetections *extCMF = NULL; psArray *extSourcesTXT = NULL; - - // find the currently selected readout - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 0); // File of interest + int index = 0; + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", index); // File of interest psAssert (file, "missing file?"); @@ -118,4 +119,5 @@ psFree (source->modelPSF); source->modelPSF = NULL; + source->imageID = index; psArrayAdd (detections->newSources, 100, source); @@ -140,4 +142,5 @@ psFree (source->modelPSF); source->modelPSF = NULL; + source->imageID = index; psArrayAdd (detections->newSources, 100, source); Index: /trunk/psphot/src/psphotModelBackground.c =================================================================== --- /trunk/psphot/src/psphotModelBackground.c (revision 27656) +++ /trunk/psphot/src/psphotModelBackground.c (revision 27657) @@ -371,6 +371,6 @@ psImageBinning *binning = psphotBackgroundBinning(readout->image, config); // Image binning parameters - pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning); - pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning); + pmReadout *model = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL", inFPA, binning, index); + pmReadout *modelStdev = pmFPAGenerateReadout(config, view, "PSPHOT.BACKMDL.STDEV", inFPA, binning, index); if (!psphotModelBackgroundReadout(model->image, modelStdev->image, model->analysis, readout, binning, config)) { Index: /trunk/psphot/src/psphotReadout.c =================================================================== --- /trunk/psphot/src/psphotReadout.c (revision 27656) +++ /trunk/psphot/src/psphotReadout.c (revision 27657) @@ -34,5 +34,7 @@ // Generate the mask and weight images, including the user-defined analysis region of interest - psphotSetMaskAndVariance (config, view); + if (!psphotSetMaskAndVariance (config, view)) { + return psphotReadoutCleanup(config, view); + } if (!strcasecmp (breakPt, "NOTHING")) { return psphotReadoutCleanup(config, view); Index: /trunk/psphot/src/psphotReadoutCleanup.c =================================================================== --- /trunk/psphot/src/psphotReadoutCleanup.c (revision 27656) +++ /trunk/psphot/src/psphotReadoutCleanup.c (revision 27657) @@ -53,5 +53,5 @@ pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF"); pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); - psArray *sources = detections->allSources; + psArray *sources = detections ? detections->allSources : NULL; // XXX where do we free these, in here (psMetadataRemove?) @@ -73,5 +73,5 @@ // Check to see if the image quality was measured // XXX not sure we want / need this test - if (!psf) { + if (0 && !psf) { bool mdok; // Status of MD lookup int nIQ = psMetadataLookupS32(&mdok, recipe, "IQ_NSTAR"); // Number of stars for IQ measurement Index: /trunk/psphot/src/psphotRoughClass.c =================================================================== --- /trunk/psphot/src/psphotRoughClass.c (revision 27656) +++ /trunk/psphot/src/psphotRoughClass.c (revision 27657) @@ -19,6 +19,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotRoughClassReadout (config, view, "PSPHOT.INPUT", i, recipe)) { psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotSkyReplace.c =================================================================== --- /trunk/psphot/src/psphotSkyReplace.c (revision 27656) +++ /trunk/psphot/src/psphotSkyReplace.c (revision 27657) @@ -8,6 +8,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotSkyReplaceReadout (config, view, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to replace sky for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotSourceMatch.c =================================================================== --- /trunk/psphot/src/psphotSourceMatch.c (revision 27656) +++ /trunk/psphot/src/psphotSourceMatch.c (revision 27657) @@ -1,50 +1,256 @@ -# include "psphotInternal.h" - -// we need a couple of functions to distinguish coincident sources: -// XXX these are similar (identical?) to the goals of pmSourceMatch.c - -# define NEXT1 { i++; continue; } -# define NEXT2 { j++; continue; } - -psphotSourceMerge () { - - float dx, dy; - - // sort the source list by X - pmSourceSortByX (sources1); - pmSourceSortByX (sources2); - - int i, j; - for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { - - pmSource *src1 = sources1->data[i]; - pmSource *src2 = sources2->data[j]; - - if (!src1) NEXT1; - if (!src1->peak) NEXT1; - if (!finite(src1->peak->xf)) NEXT1; - if (!finite(src1->peak->yf)) NEXT1; - - if (!src2) NEXT2; - if (!src2->peak) NEXT2; - if (!finite(src2->peak->xf)) NEXT2; - if (!finite(src2->peak->yf)) NEXT2; - - dx = src1->peak->xf - src2->peak->xf; - if (dx < -1.02*RADIUS) NEXT1; - if (dx > +1.02*RADIUS) NEXT2; - - // we are within match range, look for matches: - for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { - - dx = src1->peak->xf - src2->peak->xf; - dy = src1->peak->yf - src2->peak->yf; - - dr = dx*dx + dy*dy; - if (dr > RADIUS2) continue; - - // add to group? - } - i++; - } -} +# include "psphotInternal.h" + +bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects); + +psArray *psphotMatchSources (pmConfig *config, const pmFPAview *view) +{ + bool status = true; + + psArray *objects = psArrayAllocEmpty(100); + + int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (!psphotMatchSourcesReadout (objects, config, view, "PSPHOT.INPUT", i)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to merge sources for PSPHOT.INPUT entry %d", i); + psFree (objects); + return NULL; + } + } + + psphotMatchSourcesGenerate (config, view, objects); + + return objects; +} + +bool psphotMatchSourcesReadout (psArray *objects, pmConfig *config, const pmFPAview *view, char *filename, int index) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int RADIUS = psMetadataLookupF32 (&status, recipe, "PSPHOT.STACK.MATCH.RADIUS"); + psAssert (status, "programming error: must define PSPHOT.STACK.MATCH.RADIUS"); + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + psAssert (detections->newSources, "new sources not defined?"); + psAssert (!detections->allSources, "all sources already defined?"); + + // XXX TEST: + if (detections->newSources) { + psphotMatchSourcesToObjects(objects, detections->newSources, RADIUS); + } + + return true; +} + +# define NEXT1 { i++; continue; } +# define NEXT2 { j++; continue; } +bool psphotMatchSourcesToObjects (psArray *objects, psArray *sources, float RADIUS) { + + float dx, dy; + + float RADIUS2 = RADIUS*RADIUS; + + // sort the source list by X + sources = psArraySort (sources, pmSourceSortByX); + objects = psArraySort (objects, pmPhotObjSortByX); + + psVector *found = psVectorAlloc(sources->n, PS_TYPE_U8); + psVectorInit (found, 0); + + // match sources to existing objects + + psLogMsg ("psphot", PS_LOG_DETAIL, "attempt to match sources (%ld vs %ld)", sources->n, objects->n); + + int i, j; + for (i = j = 0; (i < sources->n) && (j < objects->n); ) { + + pmSource *src = sources->data[i]; + pmPhotObj *obj = objects->data[j]; + + if (!src) NEXT1; + if (!src->peak) NEXT1; + if (!finite(src->peak->xf)) NEXT1; + if (!finite(src->peak->yf)) NEXT1; + + if (!obj) NEXT2; + if (!finite(obj->x)) NEXT2; + if (!finite(obj->y)) NEXT2; + + dx = src->peak->xf - obj->x; + if (dx < -1.02*RADIUS) NEXT1; + if (dx > +1.02*RADIUS) NEXT2; + + // we are within match range, look for matches: + int Jmin = -1; + float Rmin = RADIUS2; + for (int J = j; (dx > -1.02*RADIUS) && (J < objects->n); J++) { + + obj = objects->data[J]; + + dx = src->peak->xf - obj->x; + dy = src->peak->yf - obj->y; + + float dr = dx*dx + dy*dy; + if (dr > RADIUS2) continue; + if (dr < Rmin) { + Rmin = dr; + Jmin = J; + } + } + + // no match, try next source + if (Jmin == -1) { + i++; + continue; + } + obj = objects->data[Jmin]; + + // add to object + pmPhotObjAddSource (obj, src); + found->data.U8[i] = 1; + i++; + } + + // add missed sources to new objects + + for (i = 0; i < sources->n; i++) { + + if (found->data.U8[i]) continue; + + pmSource *src = sources->data[i]; + + pmPhotObj *obj = pmPhotObjAlloc(); + pmPhotObjAddSource(obj, src); + psArrayAdd (objects, 100, obj); + psFree (obj); + } + psLogMsg ("psphot", PS_LOG_DETAIL, "matched sources (%ld vs %ld)", sources->n, objects->n); + + psFree (found); + return true; +} + +bool psphotMatchSourcesGenerate (pmConfig *config, const pmFPAview *view, psArray *objects) { + + bool status = false; + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + // determine properties (sky, moments) of initial sources + float OUTER = psMetadataLookupF32 (&status, recipe, "SKY_OUTER_RADIUS"); + psAssert (status, "missing SKY_OUTER_RADIUS in recipe?"); + + int nImages = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + // generate look-up arrays for detections and readouts + psArray *detArrays = psArrayAlloc(nImages); + psArray *readouts = psArrayAlloc(nImages); + + for (int i = 0; i < nImages; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + detArrays->data[i] = psMemIncrRefCounter(detections); + readouts->data[i] = psMemIncrRefCounter(readout); + } + + // vector to track if source for an image is found + psVector *found = psVectorAlloc(nImages, PS_TYPE_U8); + + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + + // mark the images for which sources have been found + psVectorInit (found, 0); + for (int j = 0; j < obj->sources->n; j++) { + + pmSource *src = obj->sources->data[j]; + int index = src->imageID; + psAssert (index >= 0, "invalid index"); + psAssert (index < found->n, "invalid index"); + + found->data.U8[index] = 1; + } + + // generate new sources for the image that are missing + for (int index = 0; index < found->n; index++) { + if (found->data.U8[index]) continue; + + pmDetections *detections = detArrays->data[index]; + pmReadout *readout = readouts->data[index]; + int row0 = readout->image->row0; + int col0 = readout->image->col0; + + // XXX the peak type is not really used in psphot + // PM_PEAK_LONE is certainly not true, but irrelevant + float peakFlux = readout->image->data.F32[(int)(obj->y-row0-0.5)][(int)(obj->x-col0-0.5)]; + pmPeak *peak = pmPeakAlloc(obj->x, obj->y, peakFlux, PM_PEAK_LONE); + peak->flux = peakFlux; + peak->SN = 1.0; + peak->xf = obj->x; + peak->yf = obj->y; + peak->dx = NAN; + peak->dy = NAN; + + // XXX assign to a footprint? + + // create a new source + pmSource *source = pmSourceAlloc(); + source->imageID = index; + + // add the peak + source->peak = psMemIncrRefCounter(peak); + + // allocate space for moments + source->moments = pmMomentsAlloc(); + + // allocate image, weight, mask arrays for each peak (square of radius OUTER) + pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER); + + peak->assigned = true; + pmPhotObjAddSource(obj, source); + psArrayAdd (detections->newSources, 100, source); + psFree (source); + psFree (peak); + } + } + + // how many sources do we have now? + int nSources = 0; + for (int i = 0; i < objects->n; i++) { + pmPhotObj *obj = objects->data[i]; + nSources += obj->sources->n; + } + psLogMsg ("psphot", PS_LOG_DETAIL, "total of %d sources for %d images", nSources, nImages); + + + psFree (found); + psFree (detArrays); + psFree (readouts); + return true; +} Index: unk/psphot/src/psphotSourceMerge.c =================================================================== --- /trunk/psphot/src/psphotSourceMerge.c (revision 27656) +++ (revision ) @@ -1,49 +1,0 @@ - -// we need a couple of functions to distinguish coincident sources: -// XXX these are similar (identical?) to the goals of pmSourceMatch.c - -# define NEXT1 { i++; continue; } -# define NEXT2 { j++; continue; } - -psphotSourceMerge () { - - float dx, dy; - - // sort the source list by X - pmSourceSortByX (sources1); - pmSourceSortByX (sources2); - - int i, j; - for (i = j = 0; (i < sources1->n) && (j < sources2->n); ) { - - pmSource *src1 = sources1->data[i]; - pmSource *src2 = sources2->data[j]; - - if (!src1) NEXT1; - if (!src1->peak) NEXT1; - if (!finite(src1->peak->xf)) NEXT1; - if (!finite(src1->peak->yf)) NEXT1; - - if (!src2) NEXT2; - if (!src2->peak) NEXT2; - if (!finite(src2->peak->xf)) NEXT2; - if (!finite(src2->peak->yf)) NEXT2; - - dx = src1->peak->xf - src2->peak->xf; - if (dx < -1.02*RADIUS) NEXT1; - if (dx > +1.02*RADIUS) NEXT2; - - // we are within match range, look for matches: - for (int J = j; (dx > -1.02*radius) && (J < sources2->n); J++) { - - dx = src1->peak->xf - src2->peak->xf; - dy = src1->peak->yf - src2->peak->yf; - - dr = dx*dx + dy*dy; - if (dr > RADIUS2) continue; - - // add to group? - } - i++; - } -} Index: /trunk/psphot/src/psphotSourceSize.c =================================================================== --- /trunk/psphot/src/psphotSourceSize.c (revision 27656) +++ /trunk/psphot/src/psphotSourceSize.c (revision 27657) @@ -44,6 +44,11 @@ psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + // loop over the available readouts for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image if (!psphotSourceSizeReadout (config, view, "PSPHOT.INPUT", i, recipe, getPSFsize)) { psError (PSPHOT_ERR_CONFIG, false, "failed on source size analysis for PSPHOT.INPUT entry %d", i); Index: /trunk/psphot/src/psphotSourceStats.c =================================================================== --- /trunk/psphot/src/psphotSourceStats.c (revision 27656) +++ /trunk/psphot/src/psphotSourceStats.c (revision 27657) @@ -2,5 +2,7 @@ // convert detections to sources and measure their basic properties (moments, local sky, sky -// variance) Note: this function only generates sources for the new peaks (peak->assigned) +// variance) Note: this function only generates sources for the new peaks (peak->assigned). +// The new sources are added to any existing sources on detections->newSources. The sources +// on detections->allSources are ignored. bool psphotSourceStats (pmConfig *config, const pmFPAview *view, bool setWindow) { @@ -40,5 +42,4 @@ pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); psAssert (detections, "missing detections?"); - psAssert (!detections->newSources, "new sources already defined?"); // XXX TEST: @@ -84,10 +85,11 @@ // generate the array of sources, define the associated pixel - sources = psArrayAllocEmpty (peaks->n); + if (!detections->newSources) { + detections->newSources = psArrayAllocEmpty (peaks->n); + } + sources = detections->newSources; // if there are no peaks, we save the empty source array and return if (!peaks->n) { - // save the new sources on the detection structure: - detections->newSources = sources; return true; } @@ -100,4 +102,5 @@ // create a new source pmSource *source = pmSourceAlloc(); + source->imageID = index; // add the peak @@ -119,5 +122,4 @@ psLogMsg ("psphot", PS_LOG_INFO, "break point PEAKS, skipping MOMENTS\n"); psphotVisualShowMoments (sources); - detections->newSources = sources; return true; } @@ -126,5 +128,5 @@ if (!psphotSetMomentsWindow(recipe, readout->analysis, sources)) { psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); - psFree(sources); + psFree(detections->newSources); return false; } @@ -178,5 +180,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS"); psFree (job); - psFree(sources); + psFree(detections->newSources); return false; } @@ -187,5 +189,5 @@ if (!psThreadPoolWait (false)) { psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS"); - psFree(sources); + psFree(detections->newSources); return false; } @@ -215,8 +217,4 @@ psphotVisualShowMoments (sources); - // save the new sources on the detection structure: - detections->newSources = sources; - - if (detections->allSources) { psphotMaskCosmicRayFootprintCheck(detections->allSources); @@ -377,4 +375,7 @@ pmSource *source = sources->data[i]; + if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue; + source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED; + // skip faint sources for moments measurement if (source->peak->SN < MIN_SN) { Index: /trunk/psphot/src/psphotStackArguments.c =================================================================== --- /trunk/psphot/src/psphotStackArguments.c (revision 27656) +++ /trunk/psphot/src/psphotStackArguments.c (revision 27657) @@ -2,14 +2,13 @@ static void usage(pmConfig *config, int exitCode); -static void writeHelpInfo(pmConfig* config, FILE* ofile); +static void writeHelpInfo(FILE* ofile); pmConfig *psphotStackArguments(int argc, char **argv) { int N; - bool status; // print help info - if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(argv[0], config, stdout); - if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(argv[0], config, stdout); + if (psArgumentGet(argc, argv, "-help")) writeHelpInfo(stdout); + if (psArgumentGet(argc, argv, "-h")) writeHelpInfo(stdout); // load config data from default locations @@ -29,15 +28,4 @@ // Number of threads is handled PSARGUMENTS_INSTANTIATE_THREADSARG( psphot, config, argc, argv ) - - // photcode : used in output to supplement header data (argument or recipe?) - if ((N = psArgumentGet (argc, argv, "-photcode"))) { - if (argc <= N+1) { - psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); - exit(PS_EXIT_CONFIG_ERROR); - } - psArgumentRemove (N, &argc, argv); - psMetadataAddStr (options, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", argv[N]); - psArgumentRemove (N, &argc, argv); - } // visual : interactive display mode @@ -90,14 +78,14 @@ } -static void writeHelpInfo(pmConfig* config, FILE* ofile) +static void writeHelpInfo(FILE* ofile) { fprintf(ofile, - "Usage: psphotStack -input (input.mdc) outroot\n" + "Usage: psphotStack -input (INPUTS.mdc) (OUTROOT)\n" "\n" - "where:\n" - " FileNameList is a text file containing filenames, one per line\n" - " MaskFileNameList is a text file of mask filenames, one per line\n" - " VarFileNameList is a text file of variance filenames, one per line\n" - " OutFileBaseName is the 'root name' for output files\n" + "where INPUTS.mdc contains various METADATAs, each with:\n" + "\tIMAGE(STR): Image filename\n" + "\tMASK(STR): Mask filename\n" + "\tVARIANCE(STR): Variance map filename\n" + "OUTROOT is the 'root name' for output files\n" "\n" "additional options:\n" @@ -136,6 +124,4 @@ " -trace-levels\n" " print current trace levels\n"); - psFree(config); - pmConfigDone(); psLibFinalize(); exit(PS_EXIT_SUCCESS); Index: /trunk/psphot/src/psphotStackChisqImage.c =================================================================== --- /trunk/psphot/src/psphotStackChisqImage.c (revision 27656) +++ /trunk/psphot/src/psphotStackChisqImage.c (revision 27657) @@ -12,6 +12,8 @@ psTimerStart ("psphot.chisq.image"); - // create a holder for the image - pmReadout *chiImage = pmReadoutAlloc(); + pmFPAfile *chisqFile = pmFPAfileSelectSingle(config->files, "PSPHOT.CHISQ.IMAGE", 0); + psAssert (chisqFile, "missing chisq image FPA?"); + + pmReadout *chiReadout = NULL; int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); @@ -20,5 +22,5 @@ // loop over the available readouts for (int i = 0; i < num; i++) { - if (!psphotStackChisqImageAddReadout(config, view, chiImage, "PSPHOT.INPUT", i)) { + if (!psphotStackChisqImageAddReadout(config, view, &chiReadout, "PSPHOT.INPUT", i)) { psError (PSPHOT_ERR_CONFIG, false, "failed to model background for PSPHOT.INPUT entry %d", i); return false; @@ -26,5 +28,16 @@ } + num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.CHISQ.NUM", PS_META_REPLACE, "", num); + num++; + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", num); + // need to save the resulting image somewhere (PSPHOT.INPUT?) + if (!psMetadataAddPtr(config->files, PS_LIST_TAIL, "PSPHOT.INPUT", PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "", chisqFile)) { + psError(PM_ERR_CONFIG, false, "could not add chisqFPA to config files"); + return false; + } psLogMsg ("psphot", PS_LOG_INFO, "built chisq image: %f sec\n", psTimerMark ("psphot.chisq.image")); @@ -34,6 +47,6 @@ bool psphotStackChisqImageAddReadout(const pmConfig *config, // Configuration - pmFPAview *view, - pmReadout *chiReadout, + const pmFPAview *view, + pmReadout **chiReadout, char *filename, int index) @@ -49,10 +62,24 @@ psImage *inImage = inReadout->image; + psAssert (inImage, "missing image?"); + psImage *inVariance = inReadout->variance; + psAssert (inVariance, "missing variance?"); + psImage *inMask = inReadout->mask; + psAssert (inMask, "missing mask?"); - psImage *chiImage = chiReadout->image; - psImage *chiVariance = chiReadout->variance; - psImage *chiMask = chiReadout->mask; + if (*chiReadout == NULL) { + *chiReadout = pmFPAGenerateReadout(config, view, "PSPHOT.CHISQ.IMAGE", input->fpa, NULL, 0); + } + + psImage *chiImage = (*chiReadout)->image; + psAssert (chiImage, "missing chi image"); + + psImage *chiVariance = (*chiReadout)->variance; + psAssert (chiVariance, "missing chi variance"); + + psImage *chiMask = (*chiReadout)->mask; + psAssert (chiMask, "missing chi mask"); // select the appropriate recipe information @@ -83,2 +110,72 @@ return true; } + +bool psphotStackRemoveChisqFromInputs (pmConfig *config) { + + bool status = false; + + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + psAssert (status, "programming error: must define PSPHOT.CHISQ.NUM"); + + int inputNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM"); + psAssert (status, "programming error: must define PSPHOT.INPUT.NUM"); + + pmFPAfileRemoveSingle (config->files, "PSPHOT.INPUT", chisqNum); + + inputNum --; + psMetadataAddS32(config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "", inputNum); + + return true; +} + +bool pmFPAfileRemoveSingle(psMetadata *files, const char *name, int num) +{ + PS_ASSERT_PTR_NON_NULL(files, NULL); + PS_ASSERT_INT_NONNEGATIVE(num, NULL); + + psList* mdList = files->list; + psHash* mdHash = files->hash; + + // Generate a REGEX to select only items that match 'name' + psString regex = NULL; // Regular expression + if (name) { + if (!psMetadataLookup(files, name)) { + // No files match the requested name + return false; + } + psStringAppend(®ex, "^%s$", name); + } + + psMetadataIterator *iter = psMetadataIteratorAlloc(files, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + + bool found = false; + int i = 0; // Counter + for (i = 0; !found && (item = psMetadataGetAndIncrement(iter)); i++) { + if (i == num) found = true; + } + psFree(iter); + if (!found) { + return false; + } + + char *key = item->name; + + // look up the name via hash to see if we have a multi or not + psMetadataItem* hashItem = psHashLookup(mdHash, name); + if (hashItem == NULL) { + psError(PS_ERR_UNKNOWN, false, _("Failed to remove metadata item, %s, from metadata table."), key); + return false; + } + + if (hashItem->type == PS_DATA_METADATA_MULTI) { + // multiple entries with same key, remove just the specified one + psListRemoveData(hashItem->data.list, item); + } else { + psHashRemove(mdHash, key); + } + psListRemoveData(mdList, item); + + return true; +} Index: /trunk/psphot/src/psphotStackImageLoop.c =================================================================== --- /trunk/psphot/src/psphotStackImageLoop.c (revision 27656) +++ /trunk/psphot/src/psphotStackImageLoop.c (revision 27657) @@ -7,5 +7,5 @@ } -bool psphotImageLoop (pmConfig *config) { +bool psphotStackImageLoop (pmConfig *config) { bool status; @@ -14,4 +14,10 @@ pmReadout *readout; + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + pmFPAview *view = pmFPAviewAlloc (0); @@ -19,6 +25,44 @@ if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); - // XXX for now, we assume there is only a single chip in the PHU: - psphotStackReadout (config, view); + // for psphot, we force data to be read at the chip level + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process); + if (! chip->process || ! chip->file_exists) { continue; } + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Chip in psphotStack."); + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for Cell in psphotStack."); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphot", 6, "Readout %d: %x %x\n", view->readout, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + + // XXX for now, we assume there is only a single chip in the PHU: + if (!psphotStackReadout (config, view)) { + psError(psErrorCodeLast(), false, "failure in psphotReadout for chip %d, cell %d, readout %d\n", view->chip, view->cell, view->readout); + psFree (view); + return false; + } + + } + // drop all versions of the internal files + status = true; + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV"); + status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files"); + psFree (view); + return false; + } + } + // save output which is saved at the chip level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed output for Chip in psphot."); + } + // save output which is saved at the fpa level + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed ouput for FPA in psphot."); // fail if we failed to handle an error @@ -32,4 +76,4 @@ the easiest way to implement this is to assume we can pre-load the full set of images up front. with 5 filters and 6000^2 (image, mask, var = 10 byte per pixel), we need 1.8GB, which is not too bad. - */ +*/ Index: /trunk/psphot/src/psphotStackParseCamera.c =================================================================== --- /trunk/psphot/src/psphotStackParseCamera.c (revision 27656) +++ /trunk/psphot/src/psphotStackParseCamera.c (revision 27657) @@ -15,7 +15,7 @@ } - psMetadata *item == NULL; - int nInputs = 0; - while ((item = psMetadataGetIndex(inputs, nInputs)) != NULL) { + int nInputs = inputs->list->n; + for (int i = 0; i < nInputs; i++) { + psMetadataItem *item = psMetadataGet(inputs, i); if (item->type != PS_DATA_METADATA) { psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type METADATA", item->name); @@ -40,5 +40,5 @@ psString mask = psMetadataLookupStr(&status, input, "MASK"); // Name of mask if (mask && strlen(mask) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.INPUT.MASK", mask, PM_FPA_FILE_MASK)) { + if (!defineFile(config, imageFile, "PSPHOT.MASK", mask, PM_FPA_FILE_MASK)) { psError(PS_ERR_UNKNOWN, false, "Unable to define file from mask %d (%s)", i, mask); return false; @@ -46,16 +46,59 @@ } - psString variance = psMetadataLookupStr(&mdok, input, "VARIANCE"); // Name of variance map + psString variance = psMetadataLookupStr(&status, input, "VARIANCE"); // Name of variance map if (variance && strlen(variance) > 0) { - if (!defineFile(config, imageFile, "PSPHOT.INPUT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { + if (!defineFile(config, imageFile, "PSPHOT.VARIANCE", variance, PM_FPA_FILE_VARIANCE)) { psError(PS_ERR_UNKNOWN, false, "Unable to define file from variance %d (%s)", i, variance); return false; } } - nInputs ++; + // the output sources are carried on the input->fpa structures + pmFPAfile *outsources = pmFPAfileDefineOutputFromFile (config, imageFile, "PSPHOT.STACK.OUTPUT"); + if (!outsources) { + psError(PSPHOT_ERR_CONFIG, false, "Cannot find a rule for PSPHOT.STACK.OUTPUT"); + return false; + } + outsources->save = true; + outsources->fileID = i; // this is used to generate output names } psMetadataRemoveKey(config->arguments, "FILENAMES"); psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", nInputs); + if (!psphotSetMaskBits (config)) { + psError (PS_ERR_UNKNOWN, false, "failed to set mask bit values"); + return NULL; + } + + // generate an pmFPAimage for the chisqImage + pmFPAfile *chisqImage = pmFPAfileDefineOutput(config, NULL, "PSPHOT.CHISQ.IMAGE"); + if (!chisqImage) { + psError(PSPHOT_ERR_CONFIG, false, "Trouble defining PSPHOT.CHISQ.IMAGE"); + return false; + } + chisqImage->save = true; + + pmFPAfile *chisqMask = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.MASK"); + if (!chisqMask) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.MASK")); + return NULL; + } + if (chisqMask->type != PM_FPA_FILE_MASK) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.MASK is not of type MASK"); + return NULL; + } + chisqMask->save = true; + + pmFPAfile *chisqVariance = pmFPAfileDefineOutput(config, chisqImage->fpa, "PSPHOT.CHISQ.VARIANCE"); + if (!chisqVariance) { + psError(PS_ERR_IO, false, _("Unable to generate output file from PSPHOT.CHISQ.VARIANCE")); + return NULL; + } + if (chisqVariance->type != PM_FPA_FILE_VARIANCE) { + psError(PS_ERR_IO, true, "PSPHOT.CHISQ.VARIANCE is not of type VARIANCE"); + return NULL; + } + chisqVariance->save = true; + +# if (0) // define the additional input/output files associated with psphot // XXX figure out which files are needed by psphotStack @@ -64,4 +107,5 @@ return false; } +# endif psTrace("psphot", 1, "Done with psphotStackParseCamera...\n"); Index: /trunk/psphot/src/psphotStackReadout.c =================================================================== --- /trunk/psphot/src/psphotStackReadout.c (revision 27656) +++ /trunk/psphot/src/psphotStackReadout.c (revision 27657) @@ -13,6 +13,9 @@ return false; } + // optional break-point for processing + char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); + PS_ASSERT_PTR_NON_NULL (breakPt, false); - // set the photcode for this image (XXX currently goes into recipe, should it go into analysis?) + // set the photcode for each image if (!psphotAddPhotcode (config, view)) { psError (PSPHOT_ERR_CONFIG, false, "trouble defining the photcode"); @@ -20,173 +23,129 @@ } - // generate a background model (median, smoothed image) - if (!psphotSetMaskAndVariance (config, view, recipe)) { - return psphotReadoutCleanup (config, NULL, recipe, NULL, NULL, NULL); + // Generate the mask and weight images + if (!psphotSetMaskAndVariance (config, view)) { + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "NOTHING")) { - return psphotReadoutCleanup(config, NULL, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } - // optional break-point for processing - char *breakPt = psMetadataLookupStr (NULL, recipe, "BREAK_POINT"); - PS_ASSERT_PTR_NON_NULL (breakPt, false); - // generate a background model (median, smoothed image) + // XXX I think this is not defined correctly for an array of images. if (!psphotModelBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!psphotSubtractBackground (config, view)) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); } if (!strcasecmp (breakPt, "BACKMDL")) { - return psphotReadoutCleanup (config, readout, recipe, NULL, NULL, NULL); + return psphotReadoutCleanup (config, view); + } + + // load the psf model, if suppled. FWHM_X,FWHM_Y,etc are determined and saved on + // readout->analysis XXX this function currently only works with a single PSPHOT.INPUT + if (!psphotLoadPSF (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model"); + return psphotReadoutCleanup (config, view); + } + + if (!psphotStackChisqImage(config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failure to generate chisq image"); + return psphotReadoutCleanup (config, view); + } + if (!strcasecmp (breakPt, "CHISQ")) { + return psphotReadoutCleanup (config, view); } // find the detections (by peak and/or footprint) in the image. - pmDetections *detections = psphotFindDetections (NULL, readout, recipe); - if (!detections) { - psLogMsg ("psphot", 3, "unable to find detections in this image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL); + // This finds the detections on Chisq image as well as the individuals + if (!psphotFindDetections (config, view, true)) { // pass 1 + // this only happens if we had an error in psphotFindDetections + psError (PSPHOT_ERR_UNKNOWN, false, "failure in peak analysis"); + return psphotReadoutCleanup (config, view); } - // construct sources and measure basic stats - psArray *sources = psphotSourceStats (config, readout, detections, true); - if (!sources) return false; - if (!strcasecmp (breakPt, "PEAKS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + // construct sources and measure basic stats (saved on detections->newSources) + // only run this on detections from the input images, not chisq image + if (!psphotSourceStats (config, view, true)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); } - // find blended neighbors of very saturated stars - // XXX merge this with Basic Deblend? - psphotDeblendSatstars (readout, sources, recipe); + // *** generate the objects (which unify the sources from the different images) + psArray *objects = psphotMatchSources (config, view); - // mark blended peaks PS_SOURCE_BLEND - if (!psphotBasicDeblend (sources, recipe)) { - psLogMsg ("psphot", 3, "failed on deblend analysis"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + // construct sources for the newly-generated sources (from other images) + if (!psphotSourceStats (config, view, false)) { // pass 1 + psError(PSPHOT_ERR_UNKNOWN, false, "failure to generate sources"); + return psphotReadoutCleanup (config, view); } + // find blended neighbors of very saturated stars (detections->newSources) + // if (!psphotDeblendSatstars (config, view)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis"); + // return psphotReadoutCleanup (config, view); + // } + + // mark blended peaks PS_SOURCE_BLEND (detections->newSources) + // if (!psphotBasicDeblend (config, view)) { + // psError (PSPHOT_ERR_UNKNOWN, false, "failed on deblend analysis"); + // return psphotReadoutCleanup (config, view); + // } + // classify sources based on moments, brightness - if (!psphotRoughClass (readout, sources, recipe, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + // only run this on detections from the input images, not chisq image + if (!psphotRoughClass (config, view)) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to determine rough classifications"); + return psphotReadoutCleanup (config, view); + } + // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources) + // only run this on detections from the input images, not chisq image + if (!psphotImageQuality (config, view)) { // pass 1 + psError (PSPHOT_ERR_UNKNOWN, false, "failed to measure image quality"); + return psphotReadoutCleanup(config, view); } if (!strcasecmp (breakPt, "MOMENTS")) { - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } - if (!havePSF && !psphotImageQuality (recipe, sources)) { - psLogMsg("psphot", 3, "failed to measure image quality"); - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + // use bright stellar objects to measure PSF if we were supplied a PSF for any input file, + // this step is skipped + if (!psphotChoosePSF (config, view)) { // pass 1 + psLogMsg ("psphot", 3, "failure to construct a psf model"); + return psphotReadoutCleanup (config, view); + } + if (!strcasecmp (breakPt, "PSFMODEL")) { + return psphotReadoutCleanup (config, view); } - // if we were not supplied a PSF, choose one here - if (psf == NULL) { - // use bright stellar objects to measure PSF - // XXX if we do not have enough stars to generate the PSF, build one - // from the SEEING guess and model class - psf = psphotChoosePSF (readout, sources, recipe); - if (psf == NULL) { - psLogMsg ("psphot", 3, "failure to construct a psf model"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - havePSF = true; - } - if (!strcasecmp (breakPt, "PSFMODEL")) { - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - psphotVisualShowPSFModel (readout, psf); - // include externally-supplied sources - psphotLoadExtSources (config, view, sources); + // XXX fix this in the new multi-input context + // psphotLoadExtSources (config, view); // pass 1 // construct an initial model for each object, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, sources, psf); + psphotGuessModels (config, view); + + // merge the newly selected sources into the existing list + // NOTE: merge OLD and NEW + psphotMergeSources (config, view); // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) - psphotFitSourcesLinear (readout, sources, recipe, psf, FALSE); - - // We have to place these visualizations here because the models are not realized until - // psphotGuessModels or fitted until psphotFitSourcesLinear. - psphotVisualShowPSFStars (recipe, psf, sources); + psphotFitSourcesLinearStack (config, objects, FALSE); + psFree (objects); // identify CRs and extended sources - psphotSourceSize (config, readout, sources, recipe, psf, 0); - if (!strcasecmp (breakPt, "ENSEMBLE")) { - goto finish; - } - psphotVisualShowSatStars (recipe, psf, sources); - - // non-linear PSF and EXT fit to brighter sources - // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) - psphotBlendFit (config, readout, sources, psf); - - // replace all sources - psphotReplaceAllSources (sources, recipe); - - // linear fit to include all sources (subtract again) - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); - - // if we only do one pass, skip to extended source analysis - if (!strcasecmp (breakPt, "PASS1")) { - goto pass1finish; - } - // NOTE: possibly re-measure background model here with objects subtracted - - // add noise for subtracted objects - psphotAddNoise (readout, sources, recipe); - - // find fainter sources (pass 2) - detections = psphotFindDetections (detections, readout, recipe); - - // remove noise for subtracted objects (ie, return to normal noise level) - psphotSubNoise (readout, sources, recipe); - - // define new sources based on only the new peaks - psArray *newSources = psphotSourceStats (config, readout, detections, false); - - // set source type - if (!psphotRoughClass (readout, newSources, recipe, havePSF)) { - psLogMsg ("psphot", 3, "failed to find a valid PSF clump for image"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); - } - - // create full input models, set the radius to fitRadius, set circular fit mask - psphotGuessModels (config, readout, newSources, psf); - - // replace all sources so fit below applies to all at once - psphotReplaceAllSources (sources, recipe); - - // merge the newly selected sources into the existing list - psphotMergeSources (sources, newSources); - psFree (newSources); - - // linear fit to all sources - psphotFitSourcesLinear (readout, sources, recipe, psf, TRUE); - -pass1finish: - - // measure source size for the remaining sources - psphotSourceSize (config, readout, sources, recipe, psf, 0); - - psphotExtendedSourceAnalysis (readout, sources, recipe); - - psphotExtendedSourceFits (readout, sources, recipe); - -finish: - - // plot positive sources - // psphotSourcePlots (readout, sources, recipe); + psphotSourceSize (config, view, TRUE); // measure aperture photometry corrections - if (!psphotApResid (config, readout, sources, psf)) { + if (!psphotApResid (config, view)) { psLogMsg ("psphot", 3, "failed on psphotApResid"); - return psphotReadoutCleanup (config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } // calculate source magnitudes - psphotMagnitudes(config, readout, view, sources, psf); + psphotMagnitudes(config, view); - if (!psphotEfficiency(config, readout, view, psf, recipe, sources)) { + if (!psphotEfficiency(config, view)) { psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources"); psErrorClear(); @@ -200,8 +159,11 @@ // drop the references to the image pixels held by each source - psphotSourceFreePixels (sources); + psphotSourceFreePixels (config, view); + + // remove chisq image from config->file:PSPHOT.INPUT (why?) + psphotStackRemoveChisqFromInputs(config); // create the exported-metadata and free local data - return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources); + return psphotReadoutCleanup (config, view); } Index: /trunk/psphot/src/psphotSubtractBackground.c =================================================================== --- /trunk/psphot/src/psphotSubtractBackground.c (revision 27656) +++ /trunk/psphot/src/psphotSubtractBackground.c (revision 27657) @@ -23,5 +23,11 @@ pmFPAfile *modelFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest assert (modelFile); - pmReadout *model = pmFPAviewThisReadout (view, modelFile->fpa); + + pmReadout *model = NULL; + if (modelFile->mode == PM_FPA_MODE_INTERNAL) { + model = modelFile->readout; + } else { + model = pmFPAviewThisReadout (view, modelFile->fpa); + } assert (model);