Index: /branches/bills_branches/bills_201312/psphot/src/Makefile.am =================================================================== --- /branches/bills_branches/bills_201312/psphot/src/Makefile.am (revision 36406) +++ /branches/bills_branches/bills_201312/psphot/src/Makefile.am (revision 36407) @@ -25,5 +25,5 @@ libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) -bin_PROGRAMS = psphot psphotForced psphotFullForce psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff +bin_PROGRAMS = psphot psphotForced psphotFullForce psphotFullForceSummary psphotMinimal psphotMakePSF psphotStack psphotModelTest psmakecff # bin_PROGRAMS = psphotPetrosianStudy psphotTest psphotMomentsStudy @@ -39,4 +39,8 @@ psphotFullForce_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) psphotFullForce_LDADD = libpsphot.la + +psphotFullForceSummary_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) +psphotFullForceSummary_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS) +psphotFullForceSummary_LDADD = libpsphot.la psphotMinimal_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS) @@ -99,4 +103,9 @@ psphotMosaicChip.c \ psphotCleanup.c + +# combine full force results from several inputs +psphotFullForceSummary_SOURCES = \ + psphotFullForceSummary.c \ + psphotFullForceSummaryReadout.c # forced photometry of specified positions given a specified psf Index: /branches/bills_branches/bills_201312/psphot/src/psphot.h =================================================================== --- /branches/bills_branches/bills_201312/psphot/src/psphot.h (revision 36406) +++ /branches/bills_branches/bills_201312/psphot/src/psphot.h (revision 36407) @@ -361,4 +361,6 @@ bool psphotModelTestReadout(pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotFullForceSummaryReadout (pmConfig * config, const pmFPAview *view); + int psphotFileruleCount(const pmConfig *config, const char *filerule); bool psphotFileruleCountSet(const pmConfig *config, const char *filerule, int num); Index: /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummary.c =================================================================== --- /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummary.c (revision 36407) +++ /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummary.c (revision 36407) @@ -0,0 +1,273 @@ +# ifdef HAVE_CONFIG_H +# include +# endif + +#include +#include +#include +#include "psphot.h" + +// For simplicilty, this program's (simple) functions are all contained in this file. +static pmConfig* psphotFullForceSummaryArguments(int, char**); +static bool psphotFullForceSummaryParseCamera(pmConfig *); +static bool psphotFullForceSummaryImageLoop(pmConfig*); + +int main (int argc, char **argv) { + + psMemInit(); + psTimerStart ("complete"); + pmErrorRegister(); // register psModule's error codes/messages + psphotInit(); + + // load command-line arguments, options, and system config data + pmConfig *config = psphotFullForceSummaryArguments (argc, argv); + assert(config); + +// psphotVersionPrint(); + + if (!psphotFullForceSummaryParseCamera (config)) { + psErrorStackPrint(stderr, "Error setting up the camera\n"); + exit (1); + } + + if (!psphotFullForceSummaryImageLoop (config)) { + psErrorStackPrint(stderr, "Error in the psphot image loop\n"); + exit(1); + } + + // XXX:check for memory leaks + exit (0); +} + +// all functions which return to this level must raise one of the top-level error codes if they +// exit with an error. these error codes are used to specify the program exit status + +void usage() { + fprintf(stderr, "usage: psphotFullForceSummary -inputs \n"); + exit (1); +} + +static pmConfig* psphotFullForceSummaryArguments(int argc, char **argv) { + + pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE); + int N; + + if (config == NULL) { + psErrorStackPrint(stderr, "Can't read site configuration"); + exit(PS_EXIT_CONFIG_ERROR); + } + if ((N = psArgumentGet(argc, argv, "-input"))) { + if (argc <= N+1) { + psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1); + usage(); + } + psArgumentRemove(N, &argc, argv); + + unsigned int numBad = 0; // Number of bad lines + psMetadata *inputs = psMetadataConfigRead(NULL, &numBad, argv[N], false); // Input file info + if (!inputs || numBad > 0) { + psErrorStackPrint(stderr, "Unable to cleanly read MDC file with inputs."); + exit(PS_EXIT_CONFIG_ERROR); + } + psMetadataAddMetadata(config->arguments, PS_LIST_TAIL, "INPUTS", 0, "Metadata with input details", inputs); + psFree(inputs); + + psArgumentRemove(N, &argc, argv); + } else { + psErrorStackPrint(stderr, "-input must be supplied."); + usage(); + } + + if (argc < 2) { + psErrorStackPrint(stderr, "Output is required."); + usage(); + } + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", PS_DATA_STRING, "", argv[1]); + + return config; +} + +static bool psphotFullForceSummaryParseCamera(pmConfig *config) { + bool status = false; + + psMetadata *inputs = psMetadataLookupMetadata(&status, config->arguments, "INPUTS"); // The inputs info + if (!inputs) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find inputs."); + return false; + } + + int nInputs = inputs->list->n; + + for (int i = 0; i < nInputs; i++) { + psMetadataItem *item = psMetadataGet(inputs, i); + + if (item->type != PS_DATA_STRING) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Component %s of the input metadata is not of type STRING", item->name); + return false; + + } + psString sourcesFilename = item->data.str; + + psArray *dummy = psArrayAlloc(1); // dummy array of filenames + dummy->data[0] = psStringCopy(sourcesFilename); + + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "FILENAMES", PS_META_REPLACE, + "Filenames for file rule definition", dummy); + psFree(dummy); + + bool found = false; + pmFPAfile *file = pmFPAfileDefineFromArgs(&found, config, "PSPHOT.INPUT.CMF", "FILENAMES"); + if (!file || !found) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file %s from %s", "PSPHOT.INPUT.CMF", sourcesFilename); + return false; + } + if (file->type != PM_FPA_FILE_CMF) { + psError(PS_ERR_IO, true, "%s is not of type %s", sourcesFilename, pmFPAfileStringFromType(PM_FPA_FILE_CMF)); + return false; + } + } + psMetadataRemoveKey(config->arguments, "FILENAMES"); + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.CMF.NUM", PS_META_REPLACE, "number of inputs", + nInputs); + + + pmFPA *outputFPA = pmFPAConstruct(config->camera, config->cameraName); + if (!outputFPA) { + psError(psErrorCodeLast(), false, "Unable to construct an FPA from camera configuration."); + return false; + } + pmFPAfile *output = pmFPAfileDefineOutput(config, outputFPA, "PSPHOT.OUTPUT"); + psFree(outputFPA); // Drop reference + if (!output) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT")); + return false; + } + if (output->type != PM_FPA_FILE_CMF) { + psError(PSPHOT_ERR_CONFIG, true, "PSPHOT.OUTPUT is not of type CMF"); + return false; + } + output->save = true; + + return true; +} + +# define ESCAPE(MESSAGE) { \ + psError(PSPHOT_ERR_DATA, false, MESSAGE); \ + psFree (view); \ + return false; \ + } + +bool psphotFullForceSummaryImageLoop (pmConfig *config) { + + bool status; + pmChip *chip; + pmCell *cell; + pmReadout *readout; + + pmFPAfile *input = psMetadataLookupPtr (&status, config->files, "PSPHOT.INPUT.CMF"); + if (!status) { + psError(PSPHOT_ERR_PROG, false, "Can't find input data!"); + return false; + } + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT"); + if (!output) { + psError(PSPHOT_ERR_PROG, false, "Can't find output data!"); + return false; + } + + pmFPAview *view = pmFPAviewAlloc (0); + if (!pmFPAAddSourceFromView(output->fpa, view, output->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to insert HDU into FPA for writing.\n"); + psFree(view); + return NULL; + } +#ifdef notdef + pmHDU *lastHDU = NULL; // Last HDU updated +#endif + + // files associated with the science image + if (!pmFPAfileIOChecks (config, view, PM_FPA_BEFORE)) ESCAPE ("failed input for fpa in psphot."); + + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphotFullForceSummary", 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 psphotFullForceSummary."); + + // We read the WCS from the first input + { + pmHDU *hduLow = pmHDUGetLowest(input->fpa, chip, NULL); + if (hduLow && !pmAstromReadWCS(input->fpa, chip, hduLow->header, 1.0)) { + psWarning("Unable to read WCS astrometry from header."); + psErrorClear(); + pmHDU *hduHigh = pmHDUGetHighest(input->fpa, chip, NULL); + if (hduHigh && hduHigh != hduLow && + !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) { + psWarning("Unable to read WCS astrometry from primary header."); + psErrorClear(); + } + } + } + + // Copy the transformations from the input + output->fpa->fromTPA = psMemIncrRefCounter(input->fpa->fromTPA); + output->fpa->toTPA = psMemIncrRefCounter(input->fpa->toTPA); + output->fpa->toSky = psMemIncrRefCounter(input->fpa->toSky); + pmChip *outputChip = pmFPAviewThisChip(view, output->fpa); + outputChip->toFPA = psMemIncrRefCounter(chip->toFPA); + outputChip->fromFPA = psMemIncrRefCounter(chip->fromFPA); + if (output->fpa->hdu->header == NULL) { + output->fpa->hdu->header = psMetadataAlloc(); + } + // XXX: how come psphot and psphotStack don't have to do this? + if (!pmAstromWriteWCS(output->fpa->hdu->header, output->fpa, outputChip, 0.001)) { + ESCAPE("failure to copy WCS to header"); + } + + // there is now only a single chip (multiple readouts?). loop over it and process + while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphotFullForceSummary", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + + // process each of the readouts + while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) { + psLogMsg ("psphotFullForceSummary", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (! readout->data_exists) { continue; } + +#ifdef notdef + pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell); + if (hdu && hdu != lastHDU) { + // XXX: need to do this + // psphotVersionHeaderFull(hdu->header); + lastHDU = hdu; + } +#endif + if (!psphotFullForceSummaryReadout(config, view)) { + ESCAPE ("failure in psphotFullForceSummaryReadout"); + } + } + } + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for Chip in psphotFullForceSummary."); + } + + + // If these keywords are not set we get a warning message on output. + int numCols = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS1"); + int numRows = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS2"); + psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "", numCols); + psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "", numRows); + // XXX: What other keywords and concepts should be set (or copied)? + + if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot."); + + // fail if we encountered an unhandled error + if (psErrorCodeLast() != PS_ERR_NONE) psAbort ("failed to handle an error!"); + + psFree (view); + return true; +} Index: /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummaryReadout.c (revision 36407) +++ /branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummaryReadout.c (revision 36407) @@ -0,0 +1,55 @@ +#include "psphotInternal.h" + +bool psphotFullForceSummaryReadout (pmConfig *config, const pmFPAview *view) { + + bool status; + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, "PSPHOT.INPUT.CMF"); + + pmFPAfile *output = psMetadataLookupPtr (&status, config->files, "PSPHOT.OUTPUT"); + pmCell *outputCell = pmFPAviewThisCell(view, output->fpa); + pmReadout *outputReadout = pmFPAviewThisReadout(view, output->fpa); + if (!outputReadout) { + outputReadout = pmReadoutAlloc(outputCell); + } + + + // Create objects from the various input's sources + // loop over the available readouts + pmReadout *savedReadout = NULL; + for (int index = 0; index < num; index++) { + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT.CMF", index); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + if (savedReadout == NULL) { + // save this readout for help in setting some concepts below + savedReadout = readout; + } + + // find detections + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + + if (index == 0) { + // As a test, just copy readout zero's detections to the output + if (!psMetadataAddPtr (outputReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) { + psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout"); + return false; + } + + outputReadout->data_exists = true; + outputReadout->parent->data_exists = true; + outputReadout->parent->parent->data_exists = true; + } + } + + // Loop over objects and compute the summaries + // save the summaries on the output readout->analysis + + + return true; +}