Index: trunk/ppStack/src/Makefile.am =================================================================== --- trunk/ppStack/src/Makefile.am (revision 34372) +++ trunk/ppStack/src/Makefile.am (revision 34800) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppStack +bin_PROGRAMS = ppStack ppStackMedian if HAVE_SVNVERSION @@ -25,4 +25,7 @@ ppStack_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) $(PPSTATS_LIBS) $(PPSTACK_LIBS) +ppStackMedian_CFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PSPHOT_CFLAGS) $(PPSTATS_CFLAGS) $(PPSTACK_CFLAGS) +ppStackMedian_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PSPHOT_LIBS) $(PPSTATS_LIBS) $(PPSTACK_LIBS) + ppStack_SOURCES = \ ppStack.c \ @@ -43,4 +46,35 @@ ppStackCombinePrepare.c \ ppStackCombineInitial.c \ + ppStackCombineAlternate.c \ + ppStackReject.c \ + ppStackCombineFinal.c \ + ppStackCleanup.c \ + ppStackPhotometry.c \ + ppStackFinish.c \ + ppStackTarget.c \ + ppStackUpdateHeader.c \ + ppStackJPEGs.c \ + ppStackStats.c \ + ppStackErrorCodes.c + +ppStackMedian_SOURCES = \ + ppStackMedian.c \ + ppStackArguments.c \ + ppStackCamera.c \ + ppStackFiles.c \ + ppStackMedianLoop.c \ + ppStackPSF.c \ + ppStackReadout.c \ + ppStackVersion.c \ + ppStackMatch.c \ + ppStackSources.c \ + ppStackThread.c \ + ppStackOptions.c \ + ppStackSetup.c \ + ppStackPrepare.c \ + ppStackConvolve.c \ + ppStackCombinePrepare.c \ + ppStackCombineInitial.c \ + ppStackCombineAlternate.c \ ppStackReject.c \ ppStackCombineFinal.c \ Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 34372) +++ trunk/ppStack/src/ppStack.h (revision 34800) @@ -38,5 +38,8 @@ PPSTACK_FILES_STACK, // Stack files PPSTACK_FILES_UNCONV, // Unconvolved stack files - PPSTACK_FILES_PHOT // Files for photometry + PPSTACK_FILES_PHOT, // Files for photometry + PPSTACK_FILES_BKG, // Files for bkg + PPSTACK_FILES_MEDIAN_IN, // Files for median only stacks. + PPSTACK_FILES_MEDIAN_OUT // Files for median only stacks. } ppStackFileList; @@ -111,4 +114,7 @@ bool ppStackReadoutFinalThread(psThreadJob *job // Job to process ); +// Perform median stacking for background +bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); +bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); // Return software version Index: trunk/ppStack/src/ppStackArguments.c =================================================================== --- trunk/ppStack/src/ppStackArguments.c (revision 34372) +++ trunk/ppStack/src/ppStackArguments.c (revision 34800) @@ -324,4 +324,5 @@ } + psTrace("ppStack", 1, "Done parsing arguments\n"); Index: trunk/ppStack/src/ppStackCamera.c =================================================================== --- trunk/ppStack/src/ppStackCamera.c (revision 34372) +++ trunk/ppStack/src/ppStackCamera.c (revision 34800) @@ -156,5 +156,6 @@ psString psf = psMetadataLookupStr(&mdok, input, "PSF"); // Name of PSF psString sources = psMetadataLookupStr(&mdok, input, "SOURCES"); // Name of sources - + psString bkgmodel = psMetadataLookupStr(&mdok, input, "BKGMODEL"); // Name of warped background model + pmFPAfile *imageFile = defineFile(config, NULL, "PPSTACK.INPUT", image, PM_FPA_FILE_IMAGE); // File for image @@ -215,4 +216,21 @@ } + // Grab bkgmodel information here + if (!bkgmodel || strlen(bkgmodel) == 0) { + // We have no background models. + } + else { + pmFPAfile *inputBKG = defineFile(config,NULL,"PPSTACK.INPUT.BKGMODEL",bkgmodel, + PM_FPA_FILE_IMAGE); + if (!inputBKG) { + psError(psErrorCodeLast(), false, + "Unable to define file from bkgmodel %d (%s)",i,bkgmodel); + return(false); + } + }// End bkgmodel + + + + i++; } @@ -460,4 +478,25 @@ jpeg2->save = true; + // Output background + pmFPAfile *outBkg = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGMODEL"); + if (!outBkg) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGMODEL")); + return(false); + } + if (outBkg->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGMODEL is not of type IMAGE"); + return(false); + } + outBkg->save = true; + pmFPAfile *outBkgRest = pmFPAfileDefineOutput(config,NULL,"PPSTACK.OUTPUT.BKGREST"); + if (!outBkgRest) { + psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT.BKGREST")); + return(false); + } + if (outBkgRest->type != PM_FPA_FILE_IMAGE) { + psError(PPSTACK_ERR_CONFIG, true, "PPSTACK.OUTPUT.BKGREST is not of type IMAGE"); + return(false); + } + // For photometry, we operate on the chip-mosaicked image // we create a copy of the mosaicked image for psphot so we can write out a clean image @@ -496,4 +535,6 @@ } + // Define output file here. + return true; } Index: trunk/ppStack/src/ppStackCleanup.c =================================================================== --- trunk/ppStack/src/ppStackCleanup.c (revision 34372) +++ trunk/ppStack/src/ppStackCleanup.c (revision 34800) @@ -22,10 +22,24 @@ options->outRO = NULL; - options->expRO->data_exists = false; - options->expRO->parent->data_exists = false; - options->expRO->parent->parent->data_exists = false; - psFree(options->expRO); - options->expRO = NULL; + if (options->expRO) { + options->expRO->data_exists = false; + if (options->expRO->parent) { + options->expRO->parent->data_exists = false; + options->expRO->parent->parent->data_exists = false; + } + psFree(options->expRO); + options->expRO = NULL; + } + if (options->bkgRO) { + options->bkgRO->data_exists = false; + if (options->bkgRO->parent) { + options->bkgRO->parent->data_exists = false; + options->bkgRO->parent->parent->data_exists = false; + } + psFree(options->bkgRO); + options->bkgRO = NULL; + } + for (int i = 0; i < options->num; i++) { pmCellFreeData(options->cells->data[i]); @@ -62,7 +76,6 @@ bool ppStackCleanup (pmConfig *config, ppStackOptions *options) { - psExit exitValue = ppStackExitCode(PS_EXIT_SUCCESS); // Exit code - + // Ensure everything closes if (config) { @@ -72,4 +85,6 @@ ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true); ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); + ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_IN, true); + ppStackFileActivation(config, PPSTACK_FILES_MEDIAN_OUT, true); if (!ppStackFilesIterateUp(config)) { psError(psErrorCodeLast(), false, "Unable to close files."); Index: trunk/ppStack/src/ppStackCombineAlternate.c =================================================================== --- trunk/ppStack/src/ppStackCombineAlternate.c (revision 34800) +++ trunk/ppStack/src/ppStackCombineAlternate.c (revision 34800) @@ -0,0 +1,115 @@ +#include "ppStack.h" + +// This is the doomsday switch. +// #define TESTING // Enable test output +bool ppStackCombineMedian(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *outRO = options->outRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(outRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.median.fits", + outRO->parent->parent->parent->hdu->header, + outRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, false, i); + } + psFree(inputs); + outRO->data_exists = true; + outRO->parent->data_exists = true; + outRO->parent->parent->data_exists = true; + + return(true); +} + + +bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_BKGMED"); + + + pmReadout *bkgRO = options->bkgRO; + + psArray *inputs = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + inputs->data[i] = ro; + pmFPAfileClose(file,view); + } + if (!pmStackSimpleMedianCombine(bkgRO,inputs)) { + psFree(inputs); + return(false); + } +#if 0 + if (!ppStackWriteImage("/tmp/test_forced.bkgmdl.fits", + bkgRO->parent->parent->parent->hdu->header, + bkgRO->image, + config)) { + fprintf(stderr,"Failed to write image because fail.\n"); + } +#endif + for (int i = 0; i < options->num; i++) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); + pmFPAview *view = ppStackFilesIterateDown(config); + bool success = pmFPAfileClose(file,view); + if (!success) { + psTrace("ppStack",5,"I failed at closing a file.\n"); + } + + psFitsClose(file->fits); + file->fits = NULL; + file->header = NULL; + file->state = PM_FPA_STATE_CLOSED; + file->wrote_phu = false; + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, false, i); + } + psFree(inputs); + bkgRO->data_exists = true; + bkgRO->parent->data_exists = true; + bkgRO->parent->parent->data_exists = true; + + return(true); +} + Index: trunk/ppStack/src/ppStackCombinePrepare.c =================================================================== --- trunk/ppStack/src/ppStackCombinePrepare.c (revision 34372) +++ trunk/ppStack/src/ppStackCombinePrepare.c (revision 34800) @@ -1,5 +1,5 @@ #include "ppStack.h" -bool ppStackCombinePrepare(const char *outName, const char *expName, +bool ppStackCombinePrepare(const char *outName, const char *expName, const char *bkgName, ppStackFileList files, ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) @@ -27,6 +27,46 @@ options->outRO = pmReadoutAlloc(cell); // Output readout - pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell - options->expRO = pmReadoutAlloc(expCell); // Output readout + if (expName) { + pmCell *expCell = pmFPAfileThisCell(config->files, view, expName); // Exposure cell + options->expRO = pmReadoutAlloc(expCell); //Output readout + } +/* else { */ +/* options->expRO = NULL; */ +/* } */ + pmCell *bkgCell; + int bkg_r0,bkg_c0; + int bkg_nC,bkg_nR; + if (bkgName) { + ppStackFileActivationSingle(config, PPSTACK_FILES_BKG, true, 0); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", 0); + pmFPAview *view = ppStackFilesIterateDown(config); + pmReadout *ro = pmFPAviewThisReadout(view,file->fpa); + + bkg_r0 = ro->image->row0; + bkg_c0 = ro->image->col0; + bkg_nC = ro->image->numCols; + bkg_nR = ro->image->numRows; + bkgCell = pmFPAfileThisCell(config->files, view, bkgName); // Bkg cell + + options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout + // if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) { + options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL); + if (!options->bkgRO->parent->parent->parent->hdu) { + fprintf(stderr,"failed to generate a HDU for this thing.\n"); + } + options->bkgRO->parent->parent->parent->hdu->header = psMetadataCopy(options->bkgRO->parent->parent->parent->hdu->header, + ro->parent->parent->parent->hdu->header); + + options->bkgRO->parent->concepts = psMetadataCopy(options->bkgRO->parent->concepts, + ro->parent->concepts); + options->bkgRO->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->concepts, + ro->parent->parent->concepts); + options->bkgRO->parent->parent->parent->concepts = psMetadataCopy(options->bkgRO->parent->parent->parent->concepts, + ro->parent->parent->parent->concepts); + + } + else { + options->bkgRO = NULL; + } psFree(view); @@ -42,7 +82,19 @@ } - if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) { + if (expName) { + if (!pmReadoutStackDefineOutput(options->expRO, col0, row0, numCols, numRows, true, true, 0)) { psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); return false; + } + } + + if (bkgName) { + if (!pmReadoutStackDefineOutput(options->bkgRO, bkg_c0, bkg_r0, bkg_nC, bkg_nR, false, false, 0)) { + psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); + return false; + } + + + } Index: trunk/ppStack/src/ppStackFiles.c =================================================================== --- trunk/ppStack/src/ppStackFiles.c (revision 34372) +++ trunk/ppStack/src/ppStackFiles.c (revision 34800) @@ -14,4 +14,14 @@ /// Files required for the convolution static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; + +/// Files required for the background +static char *filesBkg[] = { "PPSTACK.INPUT.BKGMODEL", NULL }; + +/// Files required for median only stacking +static char *filesMedianIn[] = { "PPSTACK.INPUT", + NULL }; +/// Files required for median only stacking +static char *filesMedianOut[] = { "PPSTACK.OUTPUT", + NULL }; /// Regular (convolved) stack files @@ -19,4 +29,5 @@ "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT", "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", + "PPSTACK.OUTPUT.BKGMODEL", NULL }; /// Unconvolved stack files @@ -41,4 +52,7 @@ case PPSTACK_FILES_UNCONV: return filesUnconv; case PPSTACK_FILES_PHOT: return filesPhot; + case PPSTACK_FILES_BKG: return filesBkg; + case PPSTACK_FILES_MEDIAN_IN: return filesMedianIn; + case PPSTACK_FILES_MEDIAN_OUT: return filesMedianOut; default: psAbort("Unrecognised file list: %x", list); Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 34372) +++ trunk/ppStack/src/ppStackLoop.c (revision 34800) @@ -63,5 +63,5 @@ // Prepare for combination - if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", PPSTACK_FILES_STACK, stack, options, config)) { + if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.BKGMODEL", PPSTACK_FILES_STACK, stack, options, config)) { psError(psErrorCodeLast(), false, "Unable to prepare for combination."); psFree(stack); @@ -84,5 +84,5 @@ pmCellFreeData(options->cells->data[i]); } - psFree(stack); + // psFree(stack); // Pixel rejection @@ -142,4 +142,11 @@ } + // Generate median background stack here. + if (!ppStackCombineBackground(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to generate median of background images."); + psFree(stack); + return false; + } + ppStackFileActivation(config,PPSTACK_FILES_BKG ,false); // Photometry psTrace("ppStack", 1, "Photometering stacked image....\n"); @@ -171,5 +178,5 @@ return false; } - psFree(stack); + // psFree(stack); psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("cleanup"); @@ -186,5 +193,5 @@ // Prepare for combination - if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", PPSTACK_FILES_UNCONV, + if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, stack, options, config)) { psError(psErrorCodeLast(), false, "Unable to prepare for combination."); Index: trunk/ppStack/src/ppStackLoop.h =================================================================== --- trunk/ppStack/src/ppStackLoop.h (revision 34372) +++ trunk/ppStack/src/ppStackLoop.h (revision 34800) @@ -4,4 +4,9 @@ // Loop over the inputs, doing the combination bool ppStackLoop( + pmConfig *config, // Configuration + ppStackOptions *options // Options for stacking + ); +// Median only loop. +bool ppStackMedianLoop( pmConfig *config, // Configuration ppStackOptions *options // Options for stacking @@ -32,4 +37,5 @@ const char *outName, // Name of output file const char *expName, // Name of exposure file + const char *bkgName, // Name of background file ppStackFileList files, // Files of interest ppStackThreadData *stack, // Stack Index: trunk/ppStack/src/ppStackMedian.c =================================================================== --- trunk/ppStack/src/ppStackMedian.c (revision 34800) +++ trunk/ppStack/src/ppStackMedian.c (revision 34800) @@ -0,0 +1,51 @@ +#include "ppStack.h" + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + psTimerStart("PPSTACK"); + psTimerStart("PPSTACK_STEPS"); + + pmErrorRegister(); + ppStackErrorRegister(); + psphotErrorRegister(); + + ppStackOptions *options = NULL; // Options for stacking + + pmConfig *config = pmConfigRead(&argc, argv, PPSTACK_RECIPE); // Configuration + if (!config) { + ppStackCleanup(config, options); + } + + ppStackVersionPrint(); + + if (!pmModelClassInit()) { + psError(PPSTACK_ERR_PROG, false, "Unable to initialise model classes."); + ppStackCleanup(config, options); + } + + if (!psphotInit()) { + psError(PPSTACK_ERR_PROG, false, "Error initialising psphot."); + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsSetup(argc, argv, config)) { + ppStackCleanup(config, options); + } + + if (!ppStackCamera(config)) { + ppStackCleanup(config, options); + } + + if (!ppStackArgumentsParse(config)) { + ppStackCleanup(config, options); + } + + options = ppStackOptionsAlloc(); + if (!ppStackMedianLoop(config, options)) { + ppStackCleanup(config, options); + } + + ppStackCleanup(config, options); +} + Index: trunk/ppStack/src/ppStackMedianLoop.c =================================================================== --- trunk/ppStack/src/ppStackMedianLoop.c (revision 34800) +++ trunk/ppStack/src/ppStackMedianLoop.c (revision 34800) @@ -0,0 +1,204 @@ +#include "ppStack.h" + +// static functions are defined below +/* static int stackSummary(const ppStackOptions *options, const char *place); */ + +bool ppStackMedianLoop(pmConfig *config, ppStackOptions *options) +{ + assert(config); + + psTimerStart("PPSTACK_TOTAL"); + psTimerStart("PPSTACK_STEPS"); + + // Setup + psTrace("ppStack", 1, "Setup....\n"); + if (!ppStackSetup(options, config)) { + psError(psErrorCodeLast(), false, "Unable to setup."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("setup"); + + // Duplicate code that is in ppStackConvolve + options->cells = psArrayAlloc(options->num); // Cells for convolved images --- a handle for reading again + pmFPAview *view = ppStackFilesIterateDown(config); + for (int i = 0; i < options->num; i++) { + ppStackFileActivationSingle(config, PPSTACK_FILES_MEDIAN_IN, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + + // pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!view) { + return false; + } + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout + pmCell *inCell = readout->parent; // Input cell */ + options->cells->data[i] = psMemIncrRefCounter(inCell); + + } + psFree(view); + + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, options->convolve); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Preparation for stacking + if (!ppStackCombinePrepare("PPSTACK.OUTPUT",NULL,NULL,PPSTACK_FILES_MEDIAN_OUT,stack,options,config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } + + if (!ppStackCombineMedian(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unabel to combine data."); + psFree(stack); + return false; + } +/* for (int x = 0; x < options->outRO->image->numCols; x++) { */ +/* for (int y = 0; y < options->outRO->image->numRows; y++) { */ +/* options->outRO->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; */ +/* } */ +/* } */ + psFree(options->outRO->mask); + // Do this before performing photometry so that the cmf header gets all of the information. + options->zpInput = NULL; + options->expTimeInput = NULL; + options->airmassInput = NULL; + if (!ppStackUpdateHeader(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to update header."); + psFree(stack); + return false; + } +/* // Generate median background stack here. */ +/* if (!ppStackCombineBackground(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to generate median of background images."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +// ppStackFileActivation(config,PPSTACK_FILES_BKG ,false); +/* // Photometry */ +/* psTrace("ppStack", 1, "Photometering stacked image....\n"); */ +/* if (!ppStackPhotometry(options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform photometry."); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("photometry"); */ + + // Create JPEGS +/* if (!ppStackJPEGs(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to make jpegs."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Assemble Stats */ +/* if (!ppStackStats(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to assemble statistics."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + + // Clean up + + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_MEDIAN_OUT, PPSTACK_FILES_MEDIAN_OUT, false)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + // psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); + +/* // Unconvolved stack --- it's cheap to calculate, compared to everything else! */ +/* // XXX unconvolved stack is currently using the convolved mask! oops! */ +/* if (options->convolve) { */ +/* // Start threading */ +/* ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false); */ +/* if (!stack) { */ +/* psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); */ +/* return false; */ +/* } */ + +/* // Prepare for combination */ +/* if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, */ +/* stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to prepare for combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ + +/* // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been */ +/* psTrace("ppStack", 2, "Stack of unconvolved images....\n"); */ +/* if (!ppStackCombineFinal(stack, options->origCovars, options, config, false, true, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); */ +/* ppStackMemDump("unconv"); */ + +/* // Update Header */ +/* if (!ppStackUpdateHeader(stack, options, config)) { */ +/* psError(psErrorCodeLast(), false, "Unable to update header."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* // Clean up unconvolved stack */ +/* psTrace("ppStack", 2, "Cleaning up after unconvolved stack....\n"); */ +/* if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_UNCONV, PPSTACK_FILES_NONE, false)) { */ +/* psError(psErrorCodeLast(), false, "Unable to clean up."); */ +/* psFree(stack); */ +/* return false; */ +/* } */ +/* psFree(stack); */ +/* } */ + psFree(options->cells); options->cells = NULL; + + // Finish up + psTrace("ppStack", 1, "Finishing up....\n"); + if (!ppStackFinish(options, config)) { + psError(psErrorCodeLast(), false, "Unable to finish up."); + return false; + } + ppStackMemDump("finish"); + + return true; +} + +/* /// Print a summary of the inputs, and return the number of good inputs */ +/* static int stackSummary(const ppStackOptions *options, const char *place) */ +/* { */ +/* int numGood = 0; // Number of good inputs */ +/* psString summary = NULL; // Summary of images */ +/* for (int i = 0; i < options->num; i++) { */ +/* psString reason = NULL; // Reason for rejecting */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) { */ +/* psStringAppend(&reason, " Good."); */ +/* numGood++; */ +/* } else { */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CAL) { */ +/* psStringAppend(&reason, " Calibration failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_PSF) { */ +/* psStringAppend(&reason, " PSF measurement failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_MATCH) { */ +/* psStringAppend(&reason, " PSF matching failed."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_CHI2) { */ +/* psStringAppend(&reason, " PSF matching chi^2 deviant."); */ +/* } */ +/* if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_REJECT) { */ +/* psStringAppend(&reason, " Rejection exceeded threshold."); */ +/* } */ +/* } */ +/* psStringAppend(&summary, "Image %d: %s\n", i, reason); */ +/* psFree(reason); */ +/* } */ +/* psLogMsg("ppStack", PS_LOG_INFO, "Summary of images for %s:\n%s", place, summary); */ +/* psFree(summary); */ + +/* return numGood; */ +/* } */ Index: trunk/ppStack/src/ppStackOptions.h =================================================================== --- trunk/ppStack/src/ppStackOptions.h (revision 34372) +++ trunk/ppStack/src/ppStackOptions.h (revision 34800) @@ -48,4 +48,7 @@ // Rejection psArray *rejected; // Rejected pixels + // Background + pmReadout *bkgRO; // Output background readout + psArray *bkgImages; // Input background images } ppStackOptions; Index: trunk/ppStack/src/ppStackPrepare.c =================================================================== --- trunk/ppStack/src/ppStackPrepare.c (revision 34372) +++ trunk/ppStack/src/ppStackPrepare.c (revision 34800) @@ -132,4 +132,5 @@ psVectorInit(options->inputMask, 0); options->exposures = psVectorAlloc(options->num, PS_TYPE_F32); + psVectorInit(options->exposures, NAN); @@ -150,4 +151,8 @@ options->exposures->data.F32[i] = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); + if ((options->exposures->data.F32[i] == 0)|| + (!(isfinite(options->exposures->data.F32[i])))){ + options->exposures->data.F32[i] = psMetadataLookupF32(NULL,recipe,"DEFAULT.EXPTIME"); + } options->sumExposure += options->exposures->data.F32[i]; Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 34372) +++ trunk/ppStack/src/ppStackReadout.c (revision 34800) @@ -194,5 +194,4 @@ - bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, Index: trunk/ppStack/src/ppStackSetup.c =================================================================== --- trunk/ppStack/src/ppStackSetup.c (revision 34372) +++ trunk/ppStack/src/ppStackSetup.c (revision 34800) @@ -89,5 +89,8 @@ options->origMasks = psArrayAlloc(num); options->origVariances = psArrayAlloc(num); + options->bkgImages = psArrayAlloc(num); pmFPAview *view = pmFPAviewAlloc(0); + int nullMasks = 0; + int nullVariances = 0; for (int i = 0; i < num; i++) { { @@ -97,11 +100,34 @@ { // We want the convolved mask, since that defines the area that has been tested for outliers + if (options->convolve) { options->origMasks->data[i] = psMemIncrRefCounter(options->convMasks->data[i]); + } + else { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.MASK", i); + options->origMasks->data[i] = pmFPAfileName(file, view, config); + } + if (!(options->origMasks->data[i])) { + nullMasks++; + } } { pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i); options->origVariances->data[i] = pmFPAfileName(file, view, config); + if (!(options->origVariances->data[i])) { + nullVariances++; + } } +/* { */ +/* pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.BKGMODEL", i); */ +/* options->bkgImages->data[i] = pmFPAfileName(file, view, config); */ +/* } */ } + if (nullMasks == num) { + psFree(options->origMasks); + } + if (nullVariances == num) { + psFree(options->origVariances); + } + psFree(view); Index: trunk/ppStack/src/ppStackThread.c =================================================================== --- trunk/ppStack/src/ppStackThread.c (revision 34372) +++ trunk/ppStack/src/ppStackThread.c (revision 34800) @@ -40,4 +40,5 @@ psFree(stack->maskFits); psFree(stack->varianceFits); + psFree(stack->bkgFits); return; } @@ -67,5 +68,5 @@ PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, covariances, NULL); } - + ppStackThreadData *stack = psAlloc(sizeof(ppStackThreadData)); // Thread data, to return psMemSetDeallocator(stack, (psFreeFunc)stackThreadDataFree); @@ -77,4 +78,5 @@ stack->maskFits = psArrayAlloc(numInputs); stack->varianceFits = psArrayAlloc(numInputs); + stack->bkgFits = psArrayAlloc(numInputs); for (int i = 0; i < numInputs; i++) { if (!cells->data[i]) { @@ -95,8 +97,12 @@ psFree(resolved); \ } - + IMAGE_OPEN(imageNames, stack->imageFits, i); - IMAGE_OPEN(maskNames, stack->maskFits, i); - IMAGE_OPEN(varianceNames, stack->varianceFits, i); + if (maskNames) { + IMAGE_OPEN(maskNames, stack->maskFits, i); + } + if (varianceNames) { + IMAGE_OPEN(varianceNames, stack->varianceFits, i); + } } Index: trunk/ppStack/src/ppStackThread.h =================================================================== --- trunk/ppStack/src/ppStackThread.h (revision 34372) +++ trunk/ppStack/src/ppStackThread.h (revision 34800) @@ -25,4 +25,5 @@ psArray *maskFits; // FITS file pointers for masks psArray *varianceFits; // FITS file pointers for variances + psArray *bkgFits; // FITS file pointers for background models } ppStackThreadData; Index: trunk/ppStack/src/ppStackUpdateHeader.c =================================================================== --- trunk/ppStack/src/ppStackUpdateHeader.c (revision 34372) +++ trunk/ppStack/src/ppStackUpdateHeader.c (revision 34800) @@ -6,9 +6,10 @@ pmReadout *outRO = options->outRO; // Output readout + pmReadout *expRO = options->expRO; // Propagate WCS bool wcsDone = false; // Have we done the WCS? for (int i = 0; i < options->num && !wcsDone; i++) { - if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + if ((options->inputMask)&&(options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i])) { continue; } @@ -123,5 +124,23 @@ snprintf (field, 64, "AIR_%04d", i); psMetadataAddF32(hdu->header, PS_LIST_TAIL, field, PS_META_REPLACE, "input image airmass", value); - } + } + + // Copy information into expRO, because it should be there too. + if ((expRO)&&(expRO->parent)) { + pmHDU *expROhdu = pmHDUFromCell(expRO->parent); + if (!expROhdu->header) { + expROhdu->header = psMetadataAlloc(); + } + expRO->parent->parent->parent->hdu->header = psMetadataCopy(expRO->parent->parent->parent->hdu->header, + outRO->parent->parent->parent->hdu->header); + + expRO->parent->concepts = psMetadataCopy(expRO->parent->concepts, + outRO->parent->concepts); + expRO->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->concepts, + outRO->parent->concepts); + expRO->parent->parent->parent->concepts = psMetadataCopy(expRO->parent->parent->parent->concepts, + outRO->parent->parent->parent->concepts); + } + return true; }