Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 23340) +++ /trunk/ppStack/src/Makefile.am (revision 23341) @@ -17,15 +17,28 @@ ppStackArguments.c \ ppStackCamera.c \ + ppStackFiles.c \ ppStackLoop.c \ ppStackPSF.c \ ppStackReadout.c \ - ppStackPhotometry.c \ ppStackVersion.c \ ppStackMatch.c \ ppStackSources.c \ - ppStackThread.c + ppStackThread.c \ + ppStackOptions.c \ + ppStackSetup.c \ + ppStackPrepare.c \ + ppStackConvolve.c \ + ppStackCombineInitial.c \ + ppStackReject.c \ + ppStackCombineFinal.c \ + ppStackCleanup.c \ + ppStackPhotometry.c \ + ppStackFinish.c noinst_HEADERS = \ - ppStack.h + ppStack.h \ + ppStackLoop.h \ + ppStackOptions.h \ + ppStackThread.h CLEANFILES = *~ Index: /trunk/ppStack/src/ppStack.c =================================================================== --- /trunk/ppStack/src/ppStack.c (revision 23340) +++ /trunk/ppStack/src/ppStack.c (revision 23341) @@ -9,4 +9,5 @@ #include "ppStack.h" +#include "ppStackLoop.h" #define TIMER_NAME "PPSTACK" // Name of timer Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 23340) +++ /trunk/ppStack/src/ppStack.h (revision 23341) @@ -1,4 +1,4 @@ -#ifndef PP_STACK_H -#define PP_STACK_H +#ifndef PPSTACK_H +#define PPSTACK_H #define PPSTACK_RECIPE "PPSTACK" // Name of the recipe @@ -18,47 +18,12 @@ } ppStackMask; -// Thread for stacking chunks -// -// Each input file contributes a readout, into which is read a chunk from that file -typedef struct { - psArray *readouts; // Input readouts to read and stack - bool read; // Has the scan been read? - bool busy; // Is the scan being processed? - int firstScan; // First row of the chunk to be read for this group - int lastScan; // Last row of the chunk to be read for this group -} ppStackThread; +// List of files +typedef enum { + PPSTACK_FILES_PREPARE, // Files for preparation + PPSTACK_FILES_CONVOLVE, // Files for convolution + PPSTACK_FILES_COMBINE, // Files for combination + PPSTACK_FILES_PHOT // Files for photometry +} ppStackFileList; -// Allocator -ppStackThread *ppStackThreadAlloc(psArray *readouts // Inputs readouts to read and stack - ); - -// Data for threads -typedef struct { - psArray *threads; // Threads doing stacking - int lastScan; // Last row that's been read - psArray *imageFits; // FITS file pointers for images - psArray *maskFits; // FITS file pointers for masks - psArray *varianceFits; // FITS file pointers for variances -} ppStackThreadData; - -// Set up thread data -ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, // Array of input cells - const psArray *imageNames, // Names of images to read - const psArray *maskNames, // Names of masks to read - const psArray *varianceNames, // Names of variance maps to read - const psArray *covariances, // Covariance matrices (already read) - const pmConfig *config // Configuration - ); - -// Read chunk into the first available file thread -ppStackThread *ppStackThreadRead(bool *status, // Status of read - ppStackThreadData *stack, // Stacks available for reading - pmConfig *config, // Configuration - int numChunk, // Chunk number (only for interest) - int overlap // Overlap between subsequent scans - ); - -// Initialise the threads -void ppStackThreadInit(void); // Setup command-line arguments @@ -73,8 +38,4 @@ // Parse cameras bool ppStackCamera(pmConfig *config // Configuration - ); - -// Loop over the inputs, doing the combination -bool ppStackLoop(pmConfig *config // Configuration ); @@ -127,10 +88,4 @@ ); -// Perform photometry on stack -bool ppStackPhotometry(pmConfig *config, // Configuration - const pmReadout *readout, // Readout to be photometered - const pmFPAview *view // View to readout - ); - // Return software version psString ppStackVersion(void); @@ -171,3 +126,46 @@ ); +/// Dump memory debugging information +void ppStackMemDump( + const char *name ///< Stage name, for inclusion in the output file name + ); + +/// Activate/deactivate a list of files +void ppStackFileActivation( + pmConfig *config, // Configuration + ppStackFileList list, // Files to turn on/off + bool state // Activation state + ); + +// Activate/deactivate a single element for a list +void ppStackFileActivationSingle( + pmConfig *config, // Configuration + ppStackFileList list, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ); + +/// Iterate down the hierarchy, loading files +/// +/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells. +pmFPAview *ppStackFilesIterateDown( + pmConfig *config // Configuration + ); + +/// Iterate up the hierarchy, writing files +/// +/// We can get away with this simplistic treatment of the FPA hierarchy because we're working on skycells. +bool ppStackFilesIterateUp( + pmConfig *config // Configuration + ); + +/// Write an image to a FITS file +bool ppStackWriteImage( + const char *name, // Name of image + psMetadata *header, // Header + const psImage *image, // Image + pmConfig *config // Configuration + ); + + #endif Index: /trunk/ppStack/src/ppStackCleanup.c =================================================================== --- /trunk/ppStack/src/ppStackCleanup.c (revision 23341) +++ /trunk/ppStack/src/ppStackCleanup.c (revision 23341) @@ -0,0 +1,130 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + + +bool ppStackCleanup(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + bool mdok; // Status of MD lookup + bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files? + +#if 0 + // Ensure masked regions really look masked + { + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + if (!pmReadoutMaskApply(options->outRO, maskBad)) { + psWarning("Unable to apply mask"); + } + } +#endif + + // Close up + bool wcsDone = false; // Have we done the WCS? + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + + ppStackThread *thread = stack->threads->data[0]; // Representative stack + pmReadout *inRO = thread->readouts->data[i]; // Template readout + if (inRO && !wcsDone) { + // Copy astrometry over + wcsDone = true; + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU + pmHDU *outHDU = pmHDUFromCell(options->outRO->parent); // Output HDU + pmChip *outChip = options->outRO->parent->parent; // Output chip + pmFPA *outFPA = outChip->parent; // Output FPA + if (!outHDU || !inHDU) { + psWarning("Unable to find HDU at FPA level to copy astrometry."); + } else { + if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA."); + wcsDone = false; + } else { + if (!outHDU->header) { + outHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { + psErrorClear(); + psWarning("Unable to write WCS astrometry to output FPA."); + wcsDone = false; + } + } + } + } + + if (tempDelete) { + psString imageResolved = pmConfigConvertFilename(options->imageNames->data[i], + config, false, false); + psString maskResolved = pmConfigConvertFilename(options->maskNames->data[i], + config, false, false); + psString varianceResolved = pmConfigConvertFilename(options->varianceNames->data[i], + config, false, false); + if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || + unlink(varianceResolved) == -1) { + psWarning("Unable to delete temporary files for image %d", i); + } + psFree(imageResolved); + psFree(maskResolved); + psFree(varianceResolved); + } + } + + ppStackMemDump("cleanup"); + + // Generate binned JPEGs + { + int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level + int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level + + // Target cells + pmFPAview *view = pmFPAviewAlloc(0); // View to cells of interest + view->chip = view->cell = 0; + pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1"); + pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2"); + psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask + + pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts + if (!pmReadoutRebin(ro1, options->outRO, maskValue, bin1, bin1) || + !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { + psError(PS_ERR_UNKNOWN, false, "Unable to bin output."); + psFree(ro1); + psFree(ro2); + return false; + } + psFree(ro1); + psFree(ro2); + } + + // Put version information into the header + pmHDU *hdu = pmHDUFromCell(options->outRO->parent); + if (!hdu) { + psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); + return false; + } + if (!hdu->header) { + hdu->header = psMetadataAlloc(); + } + ppStackVersionHeader(hdu->header); + + return true; +} Index: /trunk/ppStack/src/ppStackCombineFinal.c =================================================================== --- /trunk/ppStack/src/ppStackCombineFinal.c (revision 23341) +++ /trunk/ppStack/src/ppStackCombineFinal.c (revision 23341) @@ -0,0 +1,55 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +bool ppStackCombineFinal(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + stack->lastScan = 0; // Reset read + bool status; // Status of read + for (int numChunk = 0; true; numChunk++) { + ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0); + if (!status) { + // Something went wrong + psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); + return false; + } + if (!thread) { + // Nothing more to read + break; + } + + // call: ppStackReadoutFinal(config, outRO, readouts, rejected) + psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start + psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, options); + psArrayAdd(job->args, 1, config); + if (!psThreadJobAddPending(job)) { + psFree(job); + return false; + } + psFree(job); + } + + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do final combination."); + return false; + } + +#ifdef TESTING + pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); + ppStackWriteImage("combined_final.fits", NULL, outRO->image, config); +#endif + + return true; +} Index: /trunk/ppStack/src/ppStackCombineInitial.c =================================================================== --- /trunk/ppStack/src/ppStackCombineInitial.c (revision 23341) +++ /trunk/ppStack/src/ppStackCombineInitial.c (revision 23341) @@ -0,0 +1,135 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +bool ppStackCombineInitial(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_FINAL"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + int overlap = 2 * psMetadataLookupS32(NULL, ppsub, + "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans + + pmFPAview *view = NULL; // View to readout + + int row0, col0; // Offset for readout + int numCols, numRows; // Size of readout + ppStackThread *thread = stack->threads->data[0]; // Representative thread + if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) { + psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); + return false; + } + + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + view = ppStackFilesIterateDown(config); + if (!view) { + return false; + } + + pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell + options->outRO = pmReadoutAlloc(outCell); // Output readout + psFree(view); + + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) { + psError(PS_ERR_UNKNOWN, false, "Unable to prepare output."); + return false; + } + + psFree(options->cells); options->cells = NULL; + + bool status; // Status of read + int numChunk; // Number of chunks + for (numChunk = 0; true; numChunk++) { + ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap); + if (!status) { + // Something went wrong + psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); + return false; + } + if (!thread) { + // Nothing more to read + break; + } + + // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels) + psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start + psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, options); + psArrayAdd(job->args, 1, config); + if (!psThreadJobAddPending(job)) { + psFree(job); + return false; + } + psFree(job); + } + + if (!psThreadPoolWait(false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination."); + return false; + } + + // Harvest the jobs, gathering the inspection lists + options->inspect = psArrayAlloc(options->num); + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + options->inspect->data[i] = psArrayAllocEmpty(numChunk); + } + psThreadJob *job; // Completed job + while ((job = psThreadJobGetDone())) { + psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0, + "Job has incorrect type: %s", job->type); + psArray *results = job->results; // Results of job + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, results->data[i]); + } + psFree(job); + } + + ppStackMemDump("initial"); + +#ifdef TESTING + ppStackWriteImage("combined_initial.fits", NULL, outRO->image, config); + pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); +#endif + + // Sum covariance matrices + for (int i = 0; i < options->num; i++) { + if (options->inputMask->data.U8[i]) { + psFree(options->covariances->data[i]); + options->covariances->data[i] = NULL; + } + } + options->outRO->covariance = psImageCovarianceSum(options->covariances); + psFree(options->covariances); options->covariances = NULL; + psFree(options->matchChi2); options->matchChi2 = NULL; + + + if (options->stats) { + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_FINAL", 0, + "Time to make final stack", psTimerMark("PPSTACK_FINAL")); + } + + return true; +} Index: /trunk/ppStack/src/ppStackConvolve.c =================================================================== --- /trunk/ppStack/src/ppStackConvolve.c (revision 23341) +++ /trunk/ppStack/src/ppStackConvolve.c (revision 23341) @@ -0,0 +1,238 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +// Update the value of a concept +#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \ + psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \ + psAssert(item, "Concept should be present"); \ + psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \ + item->data.F32 = VALUE; \ +} + + + +bool ppStackConvolve(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + int num = options->num; // Number of inputs + options->cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again + options->subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking + options->subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking + int numGood = 0; // Number of good frames + options->numCols = 0; + options->numRows = 0; + options->matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching + psVectorInit(options->matchChi2, NAN); + options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2) + psVectorInit(options->weightings, NAN); + options->covariances = psArrayAlloc(num); // Covariance matrices + + psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging + psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + + for (int i = 0; i < num; i++) { + if (options->inputMask->data.U8[i]) { + continue; + } + psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num); + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest + pmFPAview *view = ppStackFilesIterateDown(config); + if (!view) { + psFree(rng); + psFree(fpaList); + psFree(cellList); + return false; + } + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout + psFree(view); + + if (options->numCols == 0 && options->numRows == 0) { + options->numCols = readout->image->numCols; + options->numRows = readout->image->numRows; + } else if (options->numCols != readout->image->numCols || + options->numRows != readout->image->numRows) { + psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d", + readout->image->numCols, readout->image->numRows, options->numCols, options->numRows); + psFree(rng); + psFree(fpaList); + psFree(cellList); + return false; + } + + // Background subtraction, scaling and normalisation is performed automatically by the image matching + psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction + psTimerStart("PPSTACK_MATCH"); + + if (!ppStackMatch(readout, ®ions, &kernels, &options->matchChi2->data.F32[i], + &options->weightings->data.F32[i], options->sourceLists->data[i], + options->psf, rng, config)) { + psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH; + psErrorClear(); + continue; + } + options->covariances->data[i] = psMemIncrRefCounter(readout->covariance); + + if (options->stats) { + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK, + "Time to match PSF", psTimerMark("PPSTACK_MATCH")); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK, + "Mean deviation for stamps", kernels->mean); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK, + "RMS deviation for stamps", kernels->rms); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK, + "Number of stamps", kernels->numStamps); + float deconv = psMetadataLookupF32(NULL, readout->analysis, + PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction + psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK, + "Deconvolution fraction for kernel", deconv); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK, + "Weighting for image", options->weightings->data.F32[i]); + } + psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH")); + + options->subRegions->data[i] = regions; + options->subKernels->data[i] = kernels; + + // Write the temporary convolved files + pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image + assert(hdu); + ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config); + psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask + pmConfigMaskWriteHeader(config, maskHeader); + ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config); + psFree(maskHeader); + psImageCovarianceTransfer(readout->variance, readout->covariance); + ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config); +#ifdef TESTING + { + psString name = NULL; + psStringAppend(&name, "covariance_%d.fits", i); + ppStackWriteImage(name, hdu->header, readout->covariance->image, config); + pmStackVisualPlotTestImage(readout->covariance->image, name); + psFree(name); + } +#endif + + pmCell *inCell = readout->parent; // Input cell + + psListAdd(cellList, PS_LIST_TAIL, inCell); + psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent); + + options->cells->data[i] = psMemIncrRefCounter(inCell); + if (!ppStackFilesIterateUp(config)) { + psFree(fpaList); + psFree(cellList); + psFree(rng); + return false; + } + numGood++; + + ppStackMemDump("match"); + } + psFree(rng); + + psFree(options->sourceLists); options->sourceLists = NULL; + psFree(options->psf); options->psf = NULL; + + + if (numGood == 0) { + psError(PS_ERR_UNKNOWN, false, "No good images to combine."); + psFree(fpaList); + psFree(cellList); + return false; + } + + // Update concepts for output + // XXX This should probably be placed later, so that it's not influenced + // by images that get rejected later. + { + pmFPAview view; // View for output + view.chip = view.cell = view.readout = 0; + pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell + pmFPA *outFPA = outCell->parent->parent; // Output FPA + pmConceptsAverageFPAs(outFPA, fpaList); + pmConceptsAverageCells(outCell, cellList, NULL, NULL, true); + psFree(fpaList); + psFree(cellList); + + UPDATE_CONCEPT(outFPA, "FPA.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure); + UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN); + } + + + // Reject images out-of-hand on the basis of their match chi^2 + { + psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort + for (int i = 0; i < num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { + continue; + } + values->data.F32[values->n++] = options->matchChi2->data.F32[i]; + } + assert(values->n == numGood); + if (!psVectorSortInPlace(values)) { + psError(PS_ERR_UNKNOWN, false, "Unable to sort vector."); + psFree(values); + return false; + } + float median = numGood % 2 ? values->data.F32[numGood / 2] : + 0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]); + + float rms = 0.74 * (values->data.F32[numGood * 3 / 4] - + values->data.F32[numGood / 4]); // Estimated RMS from interquartile range + psFree(values); + + float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs) + if (isfinite(rej)) { + float thresh = median + rej * rms; // Threshold for rejection + psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f", + median, rej, rms, thresh); + + int numRej = 0; // Number rejected + numGood = 0; // Number of good images + for (int i = 0; i < num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { + continue; + } + if (options->matchChi2->data.F32[i] > thresh) { + numRej++; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2; + psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f", + i, options->matchChi2->data.F32[i]); + } else { + psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f", + i, options->matchChi2->data.F32[i]); + numGood++; + } + } + } + } + + if (numGood == 0) { + psError(PS_ERR_UNKNOWN, false, "No good images to combine."); + return false; + } + + return true; +} Index: /trunk/ppStack/src/ppStackFiles.c =================================================================== --- /trunk/ppStack/src/ppStackFiles.c (revision 23341) +++ /trunk/ppStack/src/ppStackFiles.c (revision 23341) @@ -0,0 +1,184 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" + + +// Here follows lists of files for activation/deactivation at various stages. Each must be NULL-terminated. + +/// Files required in preparation for convolution +static char *filesPrepare[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL }; + +/// Files required for the convolution +static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; + +/// Output files for the combination +static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", + "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; + +/// Files for photometry +static char *filesPhot[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", + "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", + "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", + "PSPHOT.INPUT.CMF", NULL }; + +static char **stackFiles(ppStackFileList list) +{ + switch (list) { + case PPSTACK_FILES_PREPARE: return filesPrepare; + case PPSTACK_FILES_CONVOLVE: return filesConvolve; + case PPSTACK_FILES_COMBINE: return filesCombine; + case PPSTACK_FILES_PHOT: return filesPhot; + default: + psAbort("Unrecognised file list: %x", list); + } + return NULL; +} + + +void ppStackMemDump(const char *name) +{ + return; + + static int num = 0; // Counter, to make files unique and give an idea of sequence + + psString filename = NULL; // Name of file + psStringAppend(&filename, "memdump_%s_%03d.txt", name, num); + FILE *memFile = fopen(filename, "w"); + psFree(filename); + + psMemBlock **leaks = NULL; + int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); + fprintf(memFile, "# MemBlock Size Source\n"); + unsigned long total = 0; // Total memory used + for (int i = 0; i < numLeaks; i++) { + psMemBlock *mb = leaks[i]; + fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, + mb->file, mb->lineno); + total += mb->userMemorySize; + } + fclose(memFile); + psFree(leaks); + + fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0)); + num++; +} + + + +// Activate/deactivate a list of files +void ppStackFileActivation(pmConfig *config, // Configuration + ppStackFileList list, // Files to turn on/off + bool state // Activation state + ) +{ + assert(config); + + char **files = stackFiles(list); // Files to turn on/off + for (int i = 0; files[i] != NULL; i++) { + pmFPAfileActivate(config->files, state, files[i]); + } + return; +} + +// Activate/deactivate a single element for a list +void ppStackFileActivationSingle(pmConfig *config, // Configuration + ppStackFileList list, // Files to turn on/off + bool state, // Activation state + int num // Number of file in sequence + ) +{ + assert(config); + + char **files = stackFiles(list); // Files to turn on/off + for (int i = 0; files[i] != NULL; i++) { + pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file + } + return; +} + +// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells +pmFPAview *ppStackFilesIterateDown(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + return view; +} + +// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells +bool ppStackFilesIterateUp(pmConfig *config // Configuration + ) +{ + assert(config); + + pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + return false; + } + psFree(view); + return true; +} + +// Write an image to a FITS file +bool ppStackWriteImage(const char *name, // Name of image + psMetadata *header, // Header + const psImage *image, // Image + pmConfig *config // Configuration + ) +{ + assert(name); + assert(image); + + psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name + psFits *fits = psFitsOpen(resolved, "w"); + if (!fits) { + psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved); + psFree(resolved); + return false; + } + if (!psFitsWriteImage(fits, header, image, 0, NULL)) { + psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved); + psFitsClose(fits); + psFree(resolved); + return false; + } + psFitsClose(fits); + psFree(resolved); + return true; +} Index: /trunk/ppStack/src/ppStackFinish.c =================================================================== --- /trunk/ppStack/src/ppStackFinish.c (revision 23341) +++ /trunk/ppStack/src/ppStackFinish.c (revision 23341) @@ -0,0 +1,63 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +bool ppStackFinish(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psThreadPoolFinalize(); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + // Statistics on output + if (options->stats) { + psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout + view->chip = view->cell = view->readout = 0; + + ppStatsFPA(options->stats, options->outRO->parent->parent->parent, view, maskBad, config); + + psFree(view); + } + + ppStackMemDump("stats"); + + psFree(options->outRO); options->outRO = NULL; + + // Write out the output files + if (!ppStackFilesIterateUp(config)) { + return false; + } + + + // Write out summary statistics + if (options->stats) { + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_STACK", 0, + "Total time", psTimerClear("PPSTACK_TOTAL")); + + const char *statsMDC = psMetadataConfigFormat(options->stats); + if (!statsMDC || strlen(statsMDC) == 0) { + psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); + } else { + fprintf(options->statsFile, "%s", statsMDC); + } + psFree((void *)statsMDC); + fclose(options->statsFile); options->statsFile = NULL; + } + + return true; +} Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 23340) +++ /trunk/ppStack/src/ppStackLoop.c (revision 23341) @@ -4,177 +4,9 @@ #include -#include -#include -#include #include #include -#include #include "ppStack.h" - -//#define TESTING - -#define WCS_TOLERANCE 0.001 // Tolerance for WCS -#define PIXELS_BUFFER 1024 // Initial size of pixel lists - -// Here follows lists of files for activation/deactivation at various stages. Each must be NULL-terminated. - -// Files required in preparation for convolution -static char *prepareFiles[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL }; - -// Files required for the convolution -static char *convolveFiles[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; - -// Output files for the combination -static char *combineFiles[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", - "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; - -// Files for photometry -static char *photFiles[] = { "PSPHOT.INPUT", "PSPHOT.OUTPUT", "PSPHOT.RESID", "PSPHOT.BACKMDL", - "PSPHOT.BACKMDL.STDEV", "PSPHOT.BACKGND", "PSPHOT.BACKSUB", - "SOURCE.PLOT.MOMENTS", "SOURCE.PLOT.PSFMODEL", "SOURCE.PLOT.APRESID", - "PSPHOT.INPUT.CMF", NULL }; - -static void memDump(const char *name) -{ - return; - - static int num = 0; // Counter, to make files unique and give an idea of sequence - - psString filename = NULL; // Name of file - psStringAppend(&filename, "memdump_%s_%03d.txt", name, num); - FILE *memFile = fopen(filename, "w"); - psFree(filename); - - psMemBlock **leaks = NULL; - int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); - fprintf(memFile, "# MemBlock Size Source\n"); - unsigned long total = 0; // Total memory used - for (int i = 0; i < numLeaks; i++) { - psMemBlock *mb = leaks[i]; - fprintf(memFile, "%12lu\t%12zd\t%s:%d\n", mb->id, mb->userMemorySize, - mb->file, mb->lineno); - total += mb->userMemorySize; - } - fclose(memFile); - psFree(leaks); - - fprintf(stderr, "Memdump %s %d: Memory use: %ld, sbrk: %p\n", name, num, total, sbrk(0)); - num++; -} - - - -// Activate/deactivate a list of files -static void fileActivation(pmConfig *config, // Configuration - char **files, // Files to turn on/off - bool state // Activation state - ) -{ - assert(config); - assert(files); - - for (int i = 0; files[i] != NULL; i++) { - pmFPAfileActivate(config->files, state, files[i]); - } - return; -} - -// Activate/deactivate a single element for a list -static void fileActivationSingle(pmConfig *config, // Configuration - char **files, // Files to turn on/off - bool state, // Activation state - int num // Number of file in sequence - ) -{ - assert(config); - assert(files); - - for (int i = 0; files[i] != NULL; i++) { - pmFPAfileActivateSingle(config->files, state, files[i], num); // Activated file - } - return; -} - -// Iterate down the hierarchy, loading files; we can get away with this because we're working on skycells -static pmFPAview *filesIterateDown(pmConfig *config // Configuration - ) -{ - assert(config); - - pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return NULL; - } - view->chip = 0; - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return NULL; - } - view->cell = 0; - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return NULL; - } - view->readout = 0; - if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { - return NULL; - } - return view; -} - -// Iterate up the hierarchy, writing files; we can get away with this because we're working on skycells -static bool filesIterateUp(pmConfig *config // Configuration - ) -{ - assert(config); - - pmFPAview *view = pmFPAviewAlloc(0);// Pointer into FPA hierarchy - view->chip = view->cell = view->readout = 0; - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - view->readout = -1; - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - view->cell = -1; - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - view->chip = -1; - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - return false; - } - psFree(view); - return true; -} - -// Write an image to a FITS file -static bool writeImage(const char *name, // Name of image - psMetadata *header, // Header - const psImage *image, // Image - pmConfig *config // Configuration - ) -{ - assert(name); - assert(image); - - psString resolved = pmConfigConvertFilename(name, config, true, true); // Resolved file name - psFits *fits = psFitsOpen(resolved, "w"); - if (!fits) { - psError(PS_ERR_IO, false, "Unable to open FITS file %s to write image.", resolved); - psFree(resolved); - return false; - } - if (!psFitsWriteImage(fits, header, image, 0, NULL)) { - psError(PS_ERR_IO, false, "Unable to write FITS image %s.", resolved); - psFitsClose(fits); - psFree(resolved); - return false; - } - psFitsClose(fits); - psFree(resolved); - return true; -} - +#include "ppStackLoop.h" bool ppStackLoop(pmConfig *config) @@ -183,1080 +15,126 @@ psTimerStart("PPSTACK_TOTAL"); + psTimerStart("PPSTACK_STEPS"); - psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe - psAssert(recipe, "We've thrown an error on this before."); + ppStackOptions *options = ppStackOptionsAlloc(); // Options for stacking - bool mdok; // Status of MD lookup - bool tempDelete = psMetadataLookupBool(&mdok, recipe, "TEMP.DELETE"); // Delete temporary files? - const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images - if (!tempDir) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe"); + // Setup + psTrace("ppStack", 1, "Setup....\n"); + if (!ppStackSetup(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to setup."); + psFree(options); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 0: Setup: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("setup"); + + + // Preparation for stacking + psTrace("ppStack", 1, "Preparation for stacking: merging sources, determining target PSF....\n"); + if (!ppStackPrepare(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to prepare for stacking."); + psFree(options); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 1: Load Sources and Generate Target PSF: %f sec", + psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("prepare"); + + + // Convolve inputs + psTrace("ppStack", 1, "Convolving inputs to target PSF....\n"); + if (!ppStackConvolve(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve images."); + psFree(options); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec", + psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("convolve"); + + + // Start threading + ppStackThreadInit(); + ppStackThreadData *stack = ppStackThreadDataSetup(options, config); + if (!stack) { + psError(PS_ERR_IO, false, "Unable to initialise stack threads."); + psFree(options); return false; } - psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, "OUTPUT")); // Name for temporary files - const char *tempName = basename(outputName); - if (!tempName) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files."); + // Initial combination + psTrace("ppStack", 1, "Initial stack of convolved images....\n"); + if (!ppStackCombineInitial(stack, options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform initial combination."); + psFree(stack); + psFree(options); return false; } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec", + psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("convolve"); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 4: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS")); - const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images - const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks - const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp variance maps - if (!tempImage || !tempMask || !tempVariance) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe"); + + // Pixel rejection + psTrace("ppStack", 1, "Reject pixels....\n"); + if (!ppStackReject(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject pixels."); + psFree(stack); + psFree(options); return false; } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Pixel Rejection: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("reject"); - float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution - float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject - // before rejecting entire image - float poorFrac = psMetadataLookupF32(&mdok, recipe, "POOR.FRACTION"); // Fraction for "poor" - - const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics - psMetadata *stats = NULL; // Container for statistics - FILE *statsFile = NULL; // File stream for statistics - if (statsName && strlen(statsName) > 0) { - psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename - statsFile = fopen(resolved, "w"); - if (!statsFile) { - psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); - psFree(resolved); - return false; - } else { - stats = psMetadataAlloc(); - } - psFree(resolved); - } - - pmFPAfile *output = psMetadataLookupPtr(NULL, config->files, "PPSTACK.OUTPUT"); // Output file - if (!output) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Can't find output data!"); - return false; - } - int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs - - psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe - int overlap = 2 * psMetadataLookupS32(NULL, ppsub, - "KERNEL.SIZE"); // Overlap by kernel size between consecutive scans - - if (!pmConfigMaskSetBits(NULL, NULL, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); - return false; - } - - memDump("start"); - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 0 : Initialization and Configuration : %f sec", psTimerClear("PPSTACK_STEPS")); - - // Preparation iteration: Load the sources, and get a target PSF model - psTrace("ppStack", 1, "Determining target PSF....\n"); - psArray *sourceLists = psArrayAlloc(num); // Individual lists of sources for matching - pmPSF *targetPSF = NULL; // Target PSF - float sumExposure = NAN; // Sum of exposure times - psVector *inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs - psVectorInit(inputMask, 0); - if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { - pmFPAfileActivate(config->files, false, NULL); - fileActivation(config, prepareFiles, true); - pmFPAview *view = filesIterateDown(config); - if (!view) { - return false; - } - - // Generate list of PSFs - psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, - "^PPSTACK.INPUT$"); // Iterator - psMetadataItem *fileItem; // Item from iteration - psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope - int numCols = 0, numRows = 0; // Size of image - int index = 0; // Index for file - while ((fileItem = psMetadataGetAndIncrement(fileIter))) { - assert(fileItem->type == PS_DATA_UNKNOWN); - pmFPAfile *inputFile = fileItem->data.V; // An input file - pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF - pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF - if (!psf) { - psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); - psFree(view); - psFree(sourceLists); - psFree(fileIter); - psFree(psfs); - psFree(inputMask); - return false; - } - psfs->data[index] = psMemIncrRefCounter(psf); - psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF"); - - pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest - pmHDU *hdu = pmHDUFromCell(cell); - assert(hdu && hdu->header); - int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns - int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows - if (naxis1 <= 0 || naxis2 <= 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF."); - psFree(view); - psFree(sourceLists); - psFree(fileIter); - psFree(psfs); - psFree(inputMask); - return false; - } - if (numCols == 0 && numRows == 0) { - numCols = naxis1; - numRows = naxis2; - } - - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources - psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources - if (!sources) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout."); - return NULL; - } - sourceLists->data[index] = psMemIncrRefCounter(sources); - - // Re-do photometry if we don't trust the source lists - if (psMetadataLookupBool(NULL, recipe, "PHOT")) { - psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num); - pmFPAfileActivate(config->files, false, NULL); - fileActivationSingle(config, convolveFiles, true, index); - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File - pmFPAview *view = filesIterateDown(config); - if (!view) { - psFree(sourceLists); - psFree(targetPSF); - psFree(inputMask); - return false; - } - - pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest - - if (!ppStackInputPhotometry(ro, sources, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); - psFree(sourceLists); - psFree(targetPSF); - psFree(inputMask); - return false; - } - - psFree(view); - if (!filesIterateUp(config)) { - psFree(sourceLists); - psFree(targetPSF); - psFree(inputMask); - return false; - } - pmFPAfileActivate(config->files, false, NULL); - fileActivation(config, prepareFiles, true); - } - - index++; - } - psFree(fileIter); - - // Zero point calibration - sumExposure = ppStackSourcesTransparency(sourceLists, inputMask, view, config); - if (!isfinite(sumExposure) || sumExposure <= 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences"); - psFree(sourceLists); - psFree(targetPSF); - psFree(inputMask); - return false; - } - - // Generate target PSF - targetPSF = ppStackPSF(config, numCols, numRows, psfs, inputMask); - psFree(psfs); - if (!targetPSF) { - psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); - psFree(sourceLists); - psFree(view); - psFree(inputMask); - return false; - } - psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, - "Target PSF for stack", targetPSF); - - pmChip *outChip = pmFPAviewThisChip(view, output->fpa); // Output chip - psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, - "Target PSF", targetPSF); - outChip->data_exists = true; - - psFree(view); - if (!filesIterateUp(config)) { - return false; - } - } - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 1 : Load Sources and Generate Target PSF : %f sec", psTimerClear("PPSTACK_STEPS")); - memDump("psf"); - - psArray *imageNames = psArrayAlloc(num); - psArray *maskNames = psArrayAlloc(num); - psArray *varianceNames = psArrayAlloc(num); - for (int i = 0; i < num; i++) { - psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images - psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage); - psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask); - psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance); - psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName); - imageNames->data[i] = imageName; - maskNames->data[i] = maskName; - varianceNames->data[i] = varianceName; - } - // Free the outputName that we grabbed above to set up tempName. - psFree(outputName); - - // Generate convolutions and write them to disk - psTrace("ppStack", 1, "Convolving inputs to target PSF....\n"); - psArray *cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again - psArray *subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking - psArray *subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - int numGood = 0; // Number of good frames - int numCols = 0, numRows = 0; // Size of image - psVector *matchChi2 = psVectorAlloc(num, PS_TYPE_F32); // chi^2 for stamps when matching - psVectorInit(matchChi2, NAN); - psVector *weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2) - psVectorInit(weightings, NAN); - psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging - psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging - psArray *covariances = psArrayAlloc(num); // Covariance matrices - for (int i = 0; i < num; i++) { - if (inputMask->data.U8[i]) { - continue; - } - psTrace("ppStack", 2, "Convolving input %d of %d to target PSF....\n", i, num); - pmFPAfileActivate(config->files, false, NULL); - fileActivationSingle(config, convolveFiles, true, i); - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest - pmFPAview *view = filesIterateDown(config); - if (!view) { - psFree(sourceLists); - psFree(targetPSF); - psFree(rng); - psFree(inputMask); - psFree(matchChi2); - psFree(fpaList); - psFree(cellList); - psFree(covariances); - return false; - } - - pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout - psFree(view); - - if (numCols == 0 && numRows == 0) { - numCols = readout->image->numCols; - numRows = readout->image->numRows; - } else if (numCols != readout->image->numCols || numRows != readout->image->numRows) { - psError(PS_ERR_UNKNOWN, true, "Sizes of input images don't match: %dx%d vs %dx%d", - readout->image->numCols, readout->image->numRows, numCols, numRows); - psFree(sourceLists); - psFree(targetPSF); - psFree(rng); - psFree(inputMask); - psFree(matchChi2); - psFree(fpaList); - psFree(cellList); - psFree(covariances); - return false; - } - - // Background subtraction, scaling and normalisation is performed automatically by the image matching - psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction - psTimerStart("PPSTACK_MATCH"); - - if (!ppStackMatch(readout, ®ions, &kernels, &matchChi2->data.F32[i], &weightings->data.F32[i], - sourceLists->data[i], targetPSF, rng, config)) { - psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); - inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH; - psErrorClear(); - continue; - } - covariances->data[i] = psMemIncrRefCounter(readout->covariance); - - if (stats) { - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK, - "Time to match PSF", psTimerMark("PPSTACK_MATCH")); - psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK, - "Mean deviation for stamps", kernels->mean); - psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK, - "RMS deviation for stamps", kernels->rms); - psMetadataAddF32(stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK, - "Number of stamps", kernels->numStamps); - float deconv = psMetadataLookupF32(NULL, readout->analysis, - PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction - psMetadataAddF32(stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK, - "Deconvolution fraction for kernel", deconv); - psMetadataAddF32(stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK, - "Weighting for image", weightings->data.F32[i]); - } - psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH")); - - subRegions->data[i] = regions; - subKernels->data[i] = kernels; - - // Write the temporary convolved files - pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image - assert(hdu); - writeImage(imageNames->data[i], hdu->header, readout->image, config); - psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask - pmConfigMaskWriteHeader(config, maskHeader); - writeImage(maskNames->data[i], maskHeader, readout->mask, config); - psFree(maskHeader); - psImageCovarianceTransfer(readout->variance, readout->covariance); - writeImage(varianceNames->data[i], hdu->header, readout->variance, config); -#ifdef TESTING - { - psString name = NULL; - psStringAppend(&name, "covariance_%d.fits", i); - writeImage(name, hdu->header, readout->covariance->image, config); - pmStackVisualPlotTestImage(readout->covariance->image, name); - psFree(name); - } -#endif - - pmCell *inCell = readout->parent; // Input cell - - psListAdd(cellList, PS_LIST_TAIL, inCell); - psListAdd(fpaList, PS_LIST_TAIL, inCell->parent->parent); - - cells->data[i] = psMemIncrRefCounter(inCell); - if (!filesIterateUp(config)) { - psFree(fpaList); - psFree(cellList); - psFree(covariances); - return false; - } - numGood++; - - memDump("match"); - } - psFree(sourceLists); - psFree(targetPSF); - psFree(rng); - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 2 : Generate Convolutions and Save : %f sec", psTimerClear("PPSTACK_STEPS")); - - if (numGood == 0) { - psError(PS_ERR_UNKNOWN, false, "No good images to combine."); - psFree(subKernels); - psFree(subRegions); - psFree(cells); - psFree(inputMask); - psFree(matchChi2); - psFree(fpaList); - psFree(cellList); - psFree(covariances); - return false; - } - - - // Update concepts for output - { - pmFPAview view; // View for output - view.chip = view.cell = view.readout = 0; - pmCell *outCell = pmFPAfileThisCell(config->files, &view, "PPSTACK.OUTPUT"); // Output cell - pmFPA *outFPA = outCell->parent->parent; // Output FPA - pmConceptsAverageFPAs(outFPA, fpaList); - pmConceptsAverageCells(outCell, cellList, NULL, NULL, true); - psFree(fpaList); - psFree(cellList); - - // Update the value of a concept -#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \ - psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \ - psAssert(item, "Concept should be present"); \ - psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \ - item->data.F32 = VALUE; \ - } - - UPDATE_CONCEPT(outFPA, "FPA.EXPOSURE", sumExposure); - UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", sumExposure); - UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN); - } - - - // Reject images out-of-hand on the basis of their match chi^2 - { - psVector *values = psVectorAllocEmpty(num, PS_TYPE_F32); // Values to sort - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { - continue; - } - values->data.F32[values->n++] = matchChi2->data.F32[i]; - } - assert(values->n == numGood); - if (!psVectorSortInPlace(values)) { - psError(PS_ERR_UNKNOWN, false, "Unable to sort vector."); - psFree(subKernels); - psFree(subRegions); - psFree(cells); - psFree(inputMask); - psFree(matchChi2); - psFree(values); - psFree(covariances); - return false; - } - float median = numGood % 2 ? values->data.F32[numGood / 2] : - 0.5 * (values->data.F32[numGood / 2 - 1] + values->data.F32[numGood / 2]); - - float rms = 0.74 * (values->data.F32[numGood * 3 / 4] - - values->data.F32[numGood / 4]); // Estimated RMS from interquartile range - - psFree(values); - - float rej = psMetadataLookupF32(NULL, recipe, "MATCH.REJ"); // Rejection threshold (stdevs) - if (isfinite(rej)) { - float thresh = median + rej * rms; // Threshold for rejection - psLogMsg("ppStack", PS_LOG_INFO, "chi^2 rejection threshold = %f + %f * %f = %f", - median, rej, rms, thresh); - - int numRej = 0; // Number rejected - numGood = 0; // Number of good images - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PPSTACK_MASK_ALL) { - continue; - } - if (matchChi2->data.F32[i] > thresh) { - numRej++; - inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_CHI2; - psLogMsg("ppStack", PS_LOG_INFO, "Rejecting image %d because of large matching chi^2: %f", - i, matchChi2->data.F32[i]); - } else { - psLogMsg("ppStack", PS_LOG_INFO, "Image %d has matching chi^2: %f", - i, matchChi2->data.F32[i]); - numGood++; - } - } - } - } - - if (numGood == 0) { - psError(PS_ERR_UNKNOWN, false, "No good images to combine."); - psFree(subKernels); - psFree(subRegions); - psFree(cells); - psFree(inputMask); - psFree(matchChi2); - psFree(covariances); - return false; - } - -#ifdef TESTING - psTraceSetLevel("psModules.imcombine", 7); -#endif - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 3 : Basic Image Rejection : %f sec", psTimerClear("PPSTACK_STEPS")); - - // Start threading - ppStackThreadInit(); - ppStackThreadData *stack = ppStackThreadDataSetup(cells, imageNames, maskNames, varianceNames, - covariances, config); - if (!stack) { - psError(PS_ERR_IO, false, "Unable to initialise stack threads."); - psFree(subKernels); - psFree(subRegions); - psFree(inputMask); - psFree(matchChi2); - psFree(cells); - psFree(covariances); - return false; - } - - psTimerStart("PPSTACK_INITIAL"); - - memDump("preinitial"); - - // Stack the convolved files - psTrace("ppStack", 1, "Initial stack of convolved images....\n"); - pmReadout *outRO = NULL; // Output readout - pmFPAview *view = NULL; // View to readout - psArray *inspect = NULL; // Array of arrays of pixels to inspect - { - int row0, col0; // Offset for readout - int numCols, numRows; // Size of readout - ppStackThread *thread = stack->threads->data[0]; // Representative thread - if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, thread->readouts)) { - psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(matchChi2); - psFree(cells); - psFree(covariances); - return false; - } - - pmFPAfileActivate(config->files, false, NULL); - fileActivation(config, combineFiles, true); - view = filesIterateDown(config); - if (!view) { - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(matchChi2); - psFree(cells); - psFree(covariances); - return false; - } - - pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell - outRO = pmReadoutAlloc(outCell); // Output readout - - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - if (!pmReadoutStackDefineOutput(outRO, col0, row0, numCols, numRows, true, true, maskBad)) { - psError(PS_ERR_UNKNOWN, false, "Unable to prepare output."); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(matchChi2); - psFree(view); - psFree(outRO); - psFree(cells); - psFree(covariances); - return false; - } - - psFree(cells); - - bool status; // Status of read - int numChunk; // Number of chunks - for (numChunk = 0; true; numChunk++) { - ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, overlap); - if (!status) { - // Something went wrong - psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(matchChi2); - psFree(view); - psFree(outRO); - psFree(covariances); - return false; - } - if (!thread) { - // Nothing more to read - break; - } - - // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels) - psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start - psArrayAdd(job->args, 1, thread); - psArrayAdd(job->args, 1, config); - psArrayAdd(job->args, 1, outRO); - psArrayAdd(job->args, 1, inputMask); - psArrayAdd(job->args, 1, weightings); - psArrayAdd(job->args, 1, matchChi2); - if (!psThreadJobAddPending(job)) { - psFree(job); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(weightings); - psFree(matchChi2); - psFree(view); - psFree(outRO); - psFree(covariances); - return false; - } - psFree(job); - } - - if (!psThreadPoolWait(false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do initial combination."); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(matchChi2); - psFree(view); - psFree(outRO); - psFree(covariances); - return false; - } - - // Harvest the jobs, gathering the inspection lists - inspect = psArrayAlloc(num); - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - inspect->data[i] = psArrayAllocEmpty(numChunk); - } - psThreadJob *job; // Completed job - while ((job = psThreadJobGetDone())) { - psAssert(strcmp(job->type, "PPSTACK_INITIAL_COMBINE") == 0, - "Job has incorrect type: %s", job->type); - psArray *results = job->results; // Results of job - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - inspect->data[i] = psArrayAdd(inspect->data[i], 1, results->data[i]); - } - psFree(job); - } - - memDump("initial"); - } - -#ifdef TESTING - writeImage("combined_initial.fits", NULL, outRO->image, config); - pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); -#endif - - if (stats) { - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_INITIAL", 0, - "Time to make initial stack", psTimerMark("PPSTACK_INITIAL")); - } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 4 : Make Initial Stack : %f sec", psTimerClear("PPSTACK_STEPS")); - psLogMsg("ppStack", PS_LOG_INFO, "Time to make initial stack: %f sec", psTimerClear("PPSTACK_INITIAL")); - - psTimerStart("PPSTACK_REJECT"); - - // Pixel rejection - psArray *rejected = psArrayAlloc(num); - { - int numRejected = 0; // Number of inputs rejected completely - - // Count images rejected out of hand - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - numRejected++; - } - } - - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - - psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start - psArrayAdd(job->args, 1, inspect); - PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); - if (!psThreadJobAddPending(job)) { - psFree(job); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(view); - psFree(outRO); - psFree(inspect); - psFree(rejected); - psFree(covariances); - psFree(matchChi2); - return false; - } - psFree(job); - } - - if (!psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists."); - psFree(subKernels); - psFree(subRegions); - psFree(stack); - psFree(inputMask); - psFree(view); - psFree(outRO); - psFree(inspect); - psFree(rejected); - psFree(covariances); - psFree(matchChi2); - return false; - } - - if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) { - pmStackRejectThreadsInit(); - } - - int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches - - // Reject bad pixels - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - psTimerStart("PPSTACK_REJECT"); - -#ifdef TESTING - { - psImage *mask = psPixelsToMask(NULL, inspect->data[i], - psRegionSet(0, numCols - 1, 0, numRows - 1), - 0xff); // Mask image - psString name = NULL; // Name of image - psStringAppend(&name, "inspect_%03d.fits", i); - pmStackVisualPlotTestImage(mask, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, mask, 0, NULL); - psFree(mask); - psFitsClose(fits); - } -#endif - - psPixels *reject = pmStackReject(inspect->data[i], numCols, numRows, threshold, poorFrac, stride, - subRegions->data[i], subKernels->data[i]); // Rejected pixels - -#ifdef TESTING - { - psImage *mask = psPixelsToMask(NULL, reject, psRegionSet(0, numCols - 1, 0, numRows - 1), - 0xff); // Mask image - psString name = NULL; // Name of image - psStringAppend(&name, "reject_%03d.fits", i); - pmStackVisualPlotTestImage(mask, name); - psFits *fits = psFitsOpen(name, "w"); - psFree(name); - psFitsWriteImage(fits, NULL, mask, 0, NULL); - psFree(mask); - psFitsClose(fits); - } -#endif - - psFree(inspect->data[i]); - inspect->data[i] = NULL; - - if (!reject) { - psWarning("Rejection on image %d didn't work --- reject entire image.", i); - numRejected++; - inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT; - } else { - float frac = reject->n / (float)(numCols * numRows); // Pixel fraction - psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)", - reject->n, i, frac * 100.0); - if (frac > imageRej) { - psWarning("Image %d rejected completely because rejection fraction (%.3f) " - "exceeds limit (%.3f)", i, frac, imageRej); - psFree(reject); - // reject == NULL means reject image completely - reject = NULL; - inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_BAD; - numRejected++; - } - } - - // Images without a list of rejected pixels (the list may be empty) are rejected completely - rejected->data[i] = reject; - - if (stats) { - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK, - "Time to perform rejection", psTimerMark("PPSTACK_REJECT")); - psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK, - "Number of pixels rejected", reject ? reject->n : 0); - } - psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i, - psTimerClear("PPSTACK_REJECT")); - } - - psFree(inspect); - psFree(subKernels); - psFree(subRegions); - - if (numRejected >= num - 1) { - psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed."); - psFree(stack); - psFree(rejected); - psFree(inputMask); - psFree(view); - psFree(outRO); - psFree(covariances); - psFree(matchChi2); - return false; - } - - if (stats) { - psMetadataAddS32(stats, PS_LIST_TAIL, "REJECT_IMAGES", 0, - "Number of images rejected completely", numRejected); - } - } - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 5 : Pixel Rejection : %f sec", psTimerClear("PPSTACK_STEPS")); - psTimerStart("PPSTACK_FINAL"); - - memDump("reject"); // Final combination psTrace("ppStack", 2, "Final stack of convolved images....\n"); - { - stack->lastScan = 0; // Reset read - bool status; // Status of read - for (int numChunk = 0; true; numChunk++) { - ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0); - if (!status) { - // Something went wrong - psError(PS_ERR_UNKNOWN, false, "Unable to read chunk %d", numChunk); - psFree(stack); - psFree(rejected); - psFree(inputMask); - psFree(view); - psFree(outRO); - psFree(covariances); - psFree(matchChi2); - return false; - } - if (!thread) { - // Nothing more to read - break; - } + if (!ppStackCombineFinal(stack, options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination."); + psFree(stack); + psFree(options); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("final"); - // call: ppStackReadoutFinal(config, outRO, readouts, rejected) - psThreadJob *job = psThreadJobAlloc("PPSTACK_FINAL_COMBINE"); // Job to start - psArrayAdd(job->args, 1, thread); - psArrayAdd(job->args, 1, config); - psArrayAdd(job->args, 1, outRO); - psArrayAdd(job->args, 1, inputMask); - psArrayAdd(job->args, 1, rejected); - psArrayAdd(job->args, 1, weightings); - psArrayAdd(job->args, 1, matchChi2); - if (!psThreadJobAddPending(job)) { - psFree(job); - psFree(stack); - psFree(rejected); - psFree(inputMask); - psFree(view); - psFree(outRO); - psFree(covariances); - psFree(matchChi2); - return false; - } - psFree(job); - } - if (!psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do final combination."); - psFree(stack); - psFree(inputMask); - psFree(rejected); - psFree(view); - psFree(outRO); - psFree(covariances); - psFree(matchChi2); - return false; - } + // Clean up + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + if (!ppStackCombineFinal(stack, options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to clean up."); + psFree(stack); + psFree(options); + return false; } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); - memDump("final"); - -#ifdef TESTING - pmStackVisualPlotTestImage(outRO->image, "combined_initial.fits"); - writeImage("combined_final.fits", NULL, outRO->image, config); -#endif - - // Sum covariance matrices - for (int i = 0; i < num; i++) { - if (inputMask->data.U8[i]) { - psFree(covariances->data[i]); - covariances->data[i] = NULL; - } - } - outRO->covariance = psImageCovarianceSum(covariances); - psFree(covariances); - psFree(matchChi2); - - if (stats) { - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_FINAL", 0, - "Time to make final stack", psTimerMark("PPSTACK_FINAL")); - } - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 6 : Final Stack : %f sec", psTimerClear("PPSTACK_STEPS")); - psLogMsg("ppStack", PS_LOG_INFO, "Time to make final stack: %f sec", psTimerClear("PPSTACK_FINAL")); - - psTrace("ppStack", 2, "Cleaning up after combination....\n"); - -#if 0 - // Ensure masked regions really look masked - { - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - if (!pmReadoutMaskApply(outRO, maskBad)) { - psWarning("Unable to apply mask"); - } - } -#endif - - // Close up - bool wcsDone = false; // Have we done the WCS? - for (int i = 0; i < num; i++) { - if (inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { - continue; - } - - ppStackThread *thread = stack->threads->data[0]; // Representative stack - pmReadout *inRO = thread->readouts->data[i]; // Template readout - if (inRO && !wcsDone) { - // Copy astrometry over - wcsDone = true; - pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Template HDU - pmHDU *outHDU = pmHDUFromCell(outRO->parent); // Output HDU - pmChip *outChip = outRO->parent->parent; // Output chip - pmFPA *outFPA = outChip->parent; // Output FPA - if (!outHDU || !inHDU) { - psWarning("Unable to find HDU at FPA level to copy astrometry."); - } else { - if (!pmAstromReadWCS(outFPA, outChip, inHDU->header, 1.0)) { - psErrorClear(); - psWarning("Unable to read WCS astrometry from input FPA."); - wcsDone = false; - } else { - if (!outHDU->header) { - outHDU->header = psMetadataAlloc(); - } - if (!pmAstromWriteWCS(outHDU->header, outFPA, outChip, WCS_TOLERANCE)) { - psErrorClear(); - psWarning("Unable to write WCS astrometry to output FPA."); - wcsDone = false; - } - } - } - } - - if (tempDelete) { - psString imageResolved = pmConfigConvertFilename(imageNames->data[i], config, false, false); - psString maskResolved = pmConfigConvertFilename(maskNames->data[i], config, false, false); - psString varianceResolved = pmConfigConvertFilename(varianceNames->data[i], config, false, false); - if (unlink(imageResolved) == -1 || unlink(maskResolved) == -1 || - unlink(varianceResolved) == -1) { - psWarning("Unable to delete temporary files for image %d", i); - } - psFree(imageResolved); - psFree(maskResolved); - psFree(varianceResolved); - } - } - psFree(imageNames); - psFree(maskNames); - psFree(varianceNames); - - psFree(inputMask); psFree(stack); - memDump("cleanup"); - // Generate binned JPEGs - { - int bin1 = psMetadataLookupS32(NULL, recipe, "BIN1"); // First binning level - int bin2 = psMetadataLookupS32(NULL, recipe, "BIN2"); // Second binning level - - // Target cells - pmCell *cell1 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG1"); - pmCell *cell2 = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT.JPEG2"); - psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask - - pmReadout *ro1 = pmReadoutAlloc(cell1), *ro2 = pmReadoutAlloc(cell2); // Binned readouts - if (!pmReadoutRebin(ro1, outRO, maskValue, bin1, bin1) || !pmReadoutRebin(ro2, ro1, 0, bin2, bin2)) { - psError(PS_ERR_UNKNOWN, false, "Unable to bin output."); - psFree(ro1); - psFree(ro2); - psFree(outRO); - return false; - } - psFree(ro1); - psFree(ro2); - } - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 7 : WCS & JPEGS : %f sec", psTimerClear("PPSTACK_STEPS")); - - if (psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { - psTrace("ppStack", 1, "Photometering stacked image....\n"); - - psTimerStart("PPSTACK_PHOT"); - - fileActivation(config, combineFiles, false); - fileActivation(config, photFiles, true); - pmFPAview *photView = filesIterateDown(config); - fileActivation(config, combineFiles, true); - - if (!ppStackPhotometry(config, outRO, photView)) { - psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on output."); - filesIterateUp(config); - psFree(outRO); - psFree(photView); - return false; - } - psFree(photView); - - if (stats) { - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File - pmReadout *photRO = pmFPAviewThisReadout(view, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - psMetadataAddS32(stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", sources->n); - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_PHOT", 0, - "Time to do photometry", psTimerMark("PPSTACK_PHOT")); - } - psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT")); - } - - psLogMsg("ppStack", PS_LOG_INFO, "Stage 8 : Photometry Analysis : %f sec", psTimerClear("PPSTACK_STEPS")); - - psThreadPoolFinalize(); - - memDump("phot"); - - // Statistics on output - if (stats) { - psTrace("ppStack", 1, "Gathering statistics on stacked image....\n"); - psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits for bad - psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - - ppStatsFPA(stats, outRO->parent->parent->parent, view, maskBad, config); - } - - memDump("stats"); - - // Put version information into the header - pmHDU *hdu = pmHDUFromCell(outRO->parent); - if (!hdu) { - psError(PS_ERR_UNKNOWN, false, "Unable to find HDU for output."); - psFree(outRO); - psFree(view); + // Photometry + psTrace("ppStack", 1, "Photometering stacked image....\n"); + if (!ppStackPhotometry(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry."); + psFree(options); return false; } - if (!hdu->header) { - hdu->header = psMetadataAlloc(); - } - ppStackVersionHeader(hdu->header); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("photometry"); - psFree(outRO); - psFree(view); - // Write out summary statistics - if (stats) { - psMetadataAddF32(stats, PS_LIST_TAIL, "TIME_STACK", 0, - "Time to do photometry", psTimerClear("PPSTACK_TOTAL")); - - const char *statsMDC = psMetadataConfigFormat(stats); - if (!statsMDC || strlen(statsMDC) == 0) { - psError(PS_ERR_IO, false, "Unable to get statistics MDC file.\n"); - } else { - fprintf(statsFile, "%s", statsMDC); - } - psFree((void *)statsMDC); - fclose(statsFile); - - psFree(stats); - } - - // Write out the output files - if (!filesIterateUp(config)) { + // Finish up + psTrace("ppStack", 1, "Finishing up....\n"); + if (!ppStackFinish(options, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to finish up."); + psFree(options); return false; } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 9: Final output: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("finish"); - psLogMsg("ppStack", PS_LOG_INFO, "Stage 9 : Final Output : %f sec", psTimerClear("PPSTACK_STEPS")); - - memDump("finish"); - + psFree(options); return true; } Index: /trunk/ppStack/src/ppStackLoop.h =================================================================== --- /trunk/ppStack/src/ppStackLoop.h (revision 23341) +++ /trunk/ppStack/src/ppStackLoop.h (revision 23341) @@ -0,0 +1,68 @@ +#ifndef PPSTACK_LOOP_H +#define PPSTACK_LOOP_H + +#include +#include + +#include "ppStackOptions.h" +#include "ppStackThread.h" + + +// Loop over the inputs, doing the combination +bool ppStackLoop( + pmConfig *config // Configuration + ); + +// Setup +bool ppStackSetup( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Preparation for stacking +// +// Merge input source lists, determine target PSF +bool ppStackPrepare( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Convolve inputs to match target PSF +bool ppStackConvolve( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Initial combination +bool ppStackCombineInitial( + ppStackThreadData *stack, // Stack + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Reject pixels +bool ppStackReject( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Final combination +bool ppStackCombineFinal( + ppStackThreadData *stack, // Stack + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Photometry +bool ppStackPhotometry( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +// Finish up +bool ppStackFinish( + ppStackOptions *options, // Options + pmConfig *config // Configuration + ); + +#endif Index: /trunk/ppStack/src/ppStackOptions.c =================================================================== --- /trunk/ppStack/src/ppStackOptions.c (revision 23341) +++ /trunk/ppStack/src/ppStackOptions.c (revision 23341) @@ -0,0 +1,59 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include "ppStackLoop.h" + +static void stackOptionsFree(ppStackOptions *options) +{ + psFree(options->stats); + if (options->statsFile) { + fclose(options->statsFile); + } + psFree(options->imageNames); + psFree(options->maskNames); + psFree(options->varianceNames); + psFree(options->psf); + psFree(options->inputMask); + psFree(options->sourceLists); + psFree(options->cells); + psFree(options->subKernels); + psFree(options->subRegions); + psFree(options->matchChi2); + psFree(options->weightings); + psFree(options->covariances); + psFree(options->outRO); + psFree(options->inspect); + psFree(options->rejected); + return; +} + +ppStackOptions *ppStackOptionsAlloc(void) +{ + ppStackOptions *options = psAlloc(sizeof(ppStackOptions)); // Options, to return + psMemSetDeallocator(options, (psFreeFunc)stackOptionsFree); + + options->stats = NULL; + options->statsFile = NULL; + options->imageNames = NULL; + options->maskNames = NULL; + options->varianceNames = NULL; + options->num = 0; + options->psf = NULL; + options->inputMask = NULL; + options->sourceLists = NULL; + options->cells = NULL; + options->subKernels = NULL; + options->subRegions = NULL; + options->numCols = 0; + options->numRows = 0; + options->matchChi2 = NULL; + options->weightings = NULL; + options->covariances = NULL; + options->outRO = NULL; + options->inspect = NULL; + options->rejected = NULL; + + return options; +} Index: /trunk/ppStack/src/ppStackOptions.h =================================================================== --- /trunk/ppStack/src/ppStackOptions.h (revision 23341) +++ /trunk/ppStack/src/ppStackOptions.h (revision 23341) @@ -0,0 +1,38 @@ +#ifndef PPSTACK_OPTIONS_H +#define PPSTACK_OPTIONS_H + +#include +#include + +/// Options for stacking process +typedef struct { + // Setup + psMetadata *stats; // Statistics for output + FILE *statsFile; // File to which to write statistics + psArray *imageNames, *maskNames, *varianceNames; // Filenames for the temporary convolved images + int num; // Number of inputs + // Prepare + pmPSF *psf; // Target PSF + float sumExposure; // Sum of exposure times + psVector *inputMask; // Mask for inputs + psArray *sourceLists; // Individual lists of sources for matching + // Convolve + psArray *cells; // Cells for convolved images --- a handle for reading again + psArray *subKernels; // Subtraction kernels --- required in the stacking + psArray *subRegions; // Subtraction regions --- required in the stacking + int numCols, numRows; // Size of image + psVector *matchChi2; // chi^2 for stamps from matching + psVector *weightings; // Combination weightings for images (1/noise^2) + psArray *covariances; // Covariance matrices + // Combine initial + pmReadout *outRO; // Output readout + psArray *inspect; // Array of arrays of pixels to inspect + // Rejection + psArray *rejected; // Rejected pixels +} ppStackOptions; + +/// Allocator +ppStackOptions *ppStackOptionsAlloc(void); + + +#endif Index: /trunk/ppStack/src/ppStackPhotometry.c =================================================================== --- /trunk/ppStack/src/ppStackPhotometry.c (revision 23340) +++ /trunk/ppStack/src/ppStackPhotometry.c (revision 23341) @@ -9,13 +9,10 @@ #include "ppStack.h" +#include "ppStackLoop.h" -#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \ - PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ - PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources - -bool ppStackInputPhotometry(const pmReadout *ro, const psArray *sources, const pmConfig *config) +bool ppStackPhotometry(ppStackOptions *options, pmConfig *config) { - PM_ASSERT_READOUT_NON_NULL(ro, false); - PS_ASSERT_ARRAY_NON_NULL(sources, false); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -23,98 +20,18 @@ bool mdok; // Status of MD lookup - - float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe"); - return false; - } - float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe"); - return false; - } - float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry - if (!mdok) { - psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe"); - return false; + if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { + // Nothing to do + return true; } - psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in - if (!mdok || !maskValStr) { - psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe"); - return false; - } - psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask + psTimerStart("PPSTACK_PHOT"); - psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout + ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false); + ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); + pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout + ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); - // Measure background - psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics - if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image"); - psFree(stats); - psFree(rng); - return false; - } - float bg = stats->robustMedian; // Background level - psFree(stats); - psFree(rng); - - // Photometer sources - int numCols = image->numCols, numRows = image->numRows; // Size of image - int numSources = sources->n; // Number of sources - int numGood = 0; // Number of good sources - float maxMag = -INFINITY; // Maximum magnitude - for (int j = 0; j < numSources; j++) { - pmSource *source = sources->data[j]; // Source of interest - if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) { - source->psfMag = NAN; - continue; - } - float x = source->peak->xf, y = source->peak->yf; // Coordinates of source - int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel - // Range of integration - int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius); - int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius); - - // Integrate - double sumImage = 0.0, sumKernel = 0.0; // Sums from integration - int numBadPix = 0; // Number of bad pixels - for (int v = yMin; v <= yMax; v++) { - float dy2 = PS_SQR(y - v); // Distance from centroid - for (int u = xMin; u <= xMax; u++) { - if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) { - numBadPix++; - continue; - } - float dx2 = PS_SQR(x - u); // Distance from centroid - double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value - sumImage += (image->data.F32[v][u] - bg) * kernel; - sumKernel += kernel; - } - } - - if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) { - source->psfMag = NAN; - } else { - source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius)); - if (isfinite(source->psfMag)) { - numGood++; - maxMag = PS_MAX(maxMag, source->psfMag); - } - } - } - psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag); - - return true; -} - - - -bool ppStackPhotometry(pmConfig *config, const pmReadout *readout, const pmFPAview *view) -{ - pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); - pmFPACopy(photFile->fpa, readout->parent->parent->parent); + pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry + pmFPACopy(photFile->fpa, options->outRO->parent->parent->parent); if (psMetadataLookupBool(NULL, config->arguments, "-visual")) { @@ -122,4 +39,5 @@ } +#if 0 // Need to ensure aperture residual is not calculated psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PSPHOT_RECIPE); // Recipe @@ -131,25 +49,29 @@ item->data.B = false; } +#endif // set maskValue and markValue in the psphot recipe psImageMaskType maskValue = pmConfigMaskGet("BLANK", config); // Bits to mask psImageMaskType markValue = pmConfigMaskGet("MARK.VALUE", config); // Bits to use for marking - psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "Bits to mask", maskValue); - psMetadataAddImageMask (recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, "Bits to use for mark", markValue); + psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, + "Bits to mask", maskValue); + psMetadataAddImageMask(recipe, PS_LIST_TAIL, "MARK.PSPHOT", PS_META_REPLACE, + "Bits to use for mark", markValue); - // XXX replace with psphotReadoutKnownSources (config, view) - // XXX this function requires that we construct the input source list and place it on PSPHOT.INPUT.CMF - pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); - psArray *inSources = psMetadataLookupPtr (NULL, sourcesCell->analysis, "PSPHOT.SOURCES"); + pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT"); + psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES"); if (!inSources) { psError(PS_ERR_UNKNOWN, false, "Unable to find input sources"); + psFree(photView); return false; } - if (!psphotReadoutKnownSources(config, view, inSources)) { + if (!psphotReadoutKnownSources(config, photView, inSources)) { // Clear the error, so that the output files are written. psError(PS_ERR_UNKNOWN, false, "Unable to perform photometry on stacked image."); + psFree(photView); return false; } + psFree(photView); if (!pmFPAfileDropInternal(config->files, "PSPHOT.BACKMDL") || @@ -162,4 +84,14 @@ pmFPAfileActivate(config->files, false, "PSPHOT.INPUT"); + if (options->stats) { + pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources + psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources + psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, + "Number of sources detected", sources->n); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, + "Time to do photometry", psTimerMark("PPSTACK_PHOT")); + } + psLogMsg("ppStack", PS_LOG_INFO, "Time to do photometry: %f sec", psTimerClear("PPSTACK_PHOT")); + return true; } Index: /trunk/ppStack/src/ppStackPrepare.c =================================================================== --- /trunk/ppStack/src/ppStackPrepare.c (revision 23341) +++ /trunk/ppStack/src/ppStackPrepare.c (revision 23341) @@ -0,0 +1,251 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +#define PHOT_SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_SATSTAR | PM_SOURCE_MODE_BLEND | \ + PM_SOURCE_MODE_BADPSF | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ + PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to sources + +bool ppStackInputPhotometer(const pmReadout *ro, const psArray *sources, const pmConfig *config) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PS_ASSERT_ARRAY_NON_NULL(sources, false); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + bool mdok; // Status of MD lookup + + float zpRadius = psMetadataLookupS32(&mdok, recipe, "PHOT.RADIUS"); // Radius for PHOT measurement + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.RADIUS in recipe"); + return false; + } + float zpSigma = psMetadataLookupF32(&mdok, recipe, "PHOT.SIGMA"); // Gaussian sigma for photometry + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.SIGMA in recipe"); + return false; + } + float zpFrac = psMetadataLookupF32(&mdok, recipe, "PHOT.FRAC"); // Fraction of good pixels for photometry + if (!mdok) { + psError(PS_ERR_UNKNOWN, true, "Unable to find PHOT.FRAC in recipe"); + return false; + } + + psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in + if (!mdok || !maskValStr) { + psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe"); + return false; + } + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask + + psImage *image = ro->image, *mask = ro->mask; // Image and mask from readout + + // Measure background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics + if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure background for image"); + psFree(stats); + psFree(rng); + return false; + } + float bg = stats->robustMedian; // Background level + psFree(stats); + psFree(rng); + + // Photometer sources + int numCols = image->numCols, numRows = image->numRows; // Size of image + int numSources = sources->n; // Number of sources + int numGood = 0; // Number of good sources + float maxMag = -INFINITY; // Maximum magnitude + for (int j = 0; j < numSources; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (source->mode & PHOT_SOURCE_MASK || !isfinite(source->psfMag)) { + source->psfMag = NAN; + continue; + } + float x = source->peak->xf, y = source->peak->yf; // Coordinates of source + int xCentral = x + 0.5, yCentral = y + 0.5; // Central pixel + // Range of integration + int xMin = PS_MAX(0, xCentral - zpRadius), xMax = PS_MIN(numCols - 1, xCentral + zpRadius); + int yMin = PS_MAX(0, yCentral - zpRadius), yMax = PS_MIN(numRows - 1, yCentral + zpRadius); + + // Integrate + double sumImage = 0.0, sumKernel = 0.0; // Sums from integration + int numBadPix = 0; // Number of bad pixels + for (int v = yMin; v <= yMax; v++) { + float dy2 = PS_SQR(y - v); // Distance from centroid + for (int u = xMin; u <= xMax; u++) { + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[v][u] & maskVal) { + numBadPix++; + continue; + } + float dx2 = PS_SQR(x - u); // Distance from centroid + double kernel = exp(-0.5 * (dx2 + dy2) / PS_SQR(zpSigma)); // Kernel value + sumImage += (image->data.F32[v][u] - bg) * kernel; + sumKernel += kernel; + } + } + + if (numBadPix > zpFrac * M_PI * PS_SQR(zpRadius) || !isfinite(sumImage)) { + source->psfMag = NAN; + } else { + source->psfMag = - 2.5 * log10(sumImage / sumKernel * M_PI * PS_SQR(zpRadius)); + if (isfinite(source->psfMag)) { + numGood++; + maxMag = PS_MAX(maxMag, source->psfMag); + } + } + } + psLogMsg("ppStack", PS_LOG_INFO, "Photometered %d good sources in image; max mag %f", numGood, maxMag); + + return true; +} + + +// Preparation iteration: Load the sources, and get a target PSF model +bool ppStackPrepare(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + int num = options->num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); + options->sourceLists = psArrayAlloc(num); // Individual lists of sources for matching + + options->inputMask = psVectorAlloc(num, PS_TYPE_VECTOR_MASK); // Mask for inputs + psVectorInit(options->inputMask, 0); + + if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); + pmFPAview *view = ppStackFilesIterateDown(config); + if (!view) { + return false; + } + + // Generate list of PSFs + psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, + "^PPSTACK.INPUT$"); // Iterator + psMetadataItem *fileItem; // Item from iteration + psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope + int numCols = 0, numRows = 0; // Size of image + int index = 0; // Index for file + while ((fileItem = psMetadataGetAndIncrement(fileIter))) { + assert(fileItem->type == PS_DATA_UNKNOWN); + pmFPAfile *inputFile = fileItem->data.V; // An input file + pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF + pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF + if (!psf) { + psError(PS_ERR_UNKNOWN, false, "Unable to find PSF."); + psFree(view); + psFree(fileIter); + psFree(psfs); + return false; + } + psfs->data[index] = psMemIncrRefCounter(psf); + psMetadataRemoveKey(chip->analysis, "PSPHOT.PSF"); + + pmCell *cell = pmFPAviewThisCell(view, inputFile->fpa); // Cell of interest + pmHDU *hdu = pmHDUFromCell(cell); + assert(hdu && hdu->header); + int naxis1 = psMetadataLookupS32(NULL, hdu->header, "NAXIS1"); // Number of columns + int naxis2 = psMetadataLookupS32(NULL, hdu->header, "NAXIS2"); // Number of rows + if (naxis1 <= 0 || naxis2 <= 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine size of image from PSF."); + psFree(view); + psFree(fileIter); + psFree(psfs); + return false; + } + if (numCols == 0 && numRows == 0) { + numCols = naxis1; + numRows = naxis2; + } + + pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources + psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources + if (!sources) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout."); + return NULL; + } + options->sourceLists->data[index] = psMemIncrRefCounter(sources); + + // Re-do photometry if we don't trust the source lists + if (psMetadataLookupBool(NULL, recipe, "PHOT")) { + psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num); + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File + pmFPAview *photView = ppStackFilesIterateDown(config); + if (!photView) { + psFree(view); + return false; + } + + pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest + + if (!ppStackInputPhotometer(ro, sources, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); + psFree(view); + psFree(photView); + return false; + } + + psFree(photView); + if (!ppStackFilesIterateUp(config)) { + psFree(view); + return false; + } + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); + } + + index++; + } + psFree(fileIter); + + // Zero point calibration + options->sumExposure = ppStackSourcesTransparency(options->sourceLists, options->inputMask, + view, config); + if (!isfinite(options->sumExposure) || options->sumExposure <= 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences"); + psFree(view); + return false; + } + + // Generate target PSF + options->psf = ppStackPSF(config, numCols, numRows, psfs, options->inputMask); + psFree(psfs); + if (!options->psf) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); + psFree(view); + return false; + } + psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, + "Target PSF for stack", options->psf); + + pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip + psMetadataAddPtr(outChip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, + "Target PSF", options->psf); + outChip->data_exists = true; + + psFree(view); + + if (!ppStackFilesIterateUp(config)) { + return false; + } + } + + return true; +} Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 23340) +++ /trunk/ppStack/src/ppStackReadout.c (revision 23341) @@ -8,4 +8,5 @@ #include +#include "ppStackThread.h" #include "ppStack.h" @@ -16,9 +17,11 @@ psArray *args = job->args; // Arguments ppStackThread *thread = args->data[0]; // Thread - pmConfig *config = args->data[1]; // Configuration - pmReadout *outRO = args->data[2]; // Output readout - psVector *mask = args->data[3]; // Mask for inputs - psVector *weightings = args->data[4]; // Weightings (1/noise^2) for each image - psVector *addVariance = args->data[5]; // Additional variance when rejecting + ppStackOptions *options = args->data[1]; // Options + pmConfig *config = args->data[2]; // Configuration + + pmReadout *outRO = options->outRO; // Output readout + psVector *mask = options->inputMask; // Mask for inputs + psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image + psVector *addVariance = options->matchChi2; // Additional variance when rejecting psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask, @@ -37,10 +40,12 @@ psArray *args = job->args; // Arguments ppStackThread *thread = args->data[0]; // Thread - pmConfig *config = args->data[1]; // Configuration - pmReadout *outRO = args->data[2]; // Output readout - psVector *mask = args->data[3]; // Mask for inputs - psArray *rejected = args->data[4]; // Rejected pixels - psVector *weightings = args->data[5]; // Weighting (1/noise^2) for each image - psVector *addVariance = args->data[6]; // Weighting (1/noise^2) for each image + ppStackOptions *options = args->data[1]; // Options + pmConfig *config = args->data[2]; // Configuration + + pmReadout *outRO = options->outRO; // Output readout + psVector *mask = options->inputMask; // Mask for inputs + psArray *rejected = options->rejected; // Rejected pixels + psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image + psVector *addVariance = options->matchChi2; // Additional variance when rejecting bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected, Index: /trunk/ppStack/src/ppStackReject.c =================================================================== --- /trunk/ppStack/src/ppStackReject.c (revision 23341) +++ /trunk/ppStack/src/ppStackReject.c (revision 23341) @@ -0,0 +1,157 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +bool ppStackReject(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + int num = options->num; // Number of inputs + options->rejected = psArrayAlloc(num); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(ppsub, "We've thrown an error on this before."); + + float threshold = psMetadataLookupF32(NULL, recipe, "THRESHOLD.MASK"); // Threshold for mask deconvolution + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + float imageRej = psMetadataLookupF32(NULL, recipe, "IMAGE.REJ"); // Maximum fraction of image to reject + // before rejecting entire image + int stride = psMetadataLookupS32(NULL, ppsub, "STRIDE"); // Size of convolution patches + + // Count images rejected out of hand + int numRejected = 0; // Number of inputs rejected completely + for (int i = 0; i < num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + numRejected++; + } + } + + // Concatenate inspection lists + for (int i = 0; i < num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + + psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start + psArrayAdd(job->args, 1, options->inspect); + PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); + if (!psThreadJobAddPending(job)) { + psFree(job); + return false; + } + psFree(job); + } + if (!psThreadPoolWait(true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to concatenate inspection lists."); + return false; + } + + + // Reject bad pixels + if (psMetadataLookupS32(NULL, config->arguments, "-threads") > 0) { + pmStackRejectThreadsInit(); + } + for (int i = 0; i < num; i++) { + if (options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + continue; + } + psTimerStart("PPSTACK_REJECT"); + +#ifdef TESTING + { + psImage *mask = psPixelsToMask(NULL, options->inspect->data[i], + psRegionSet(0, options->numCols - 1, 0, options->numRows - 1), + 0xff); // Mask image + psString name = NULL; // Name of image + psStringAppend(&name, "inspect_%03d.fits", i); + pmStackVisualPlotTestImage(mask, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, mask, 0, NULL); + psFree(mask); + psFitsClose(fits); + } +#endif + + psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, + threshold, poorFrac, stride, options->subRegions->data[i], + options->subKernels->data[i]); // Rejected pixels + +#ifdef TESTING + { + psImage *mask = psPixelsToMask(NULL, reject, + psRegionSet(0, options->numCols - 1, 0, options->numRows - 1), + 0xff); // Mask image + psString name = NULL; // Name of image + psStringAppend(&name, "reject_%03d.fits", i); + pmStackVisualPlotTestImage(mask, name); + psFits *fits = psFitsOpen(name, "w"); + psFree(name); + psFitsWriteImage(fits, NULL, mask, 0, NULL); + psFree(mask); + psFitsClose(fits); + } +#endif + + psFree(options->inspect->data[i]); + options->inspect->data[i] = NULL; + + if (!reject) { + psWarning("Rejection on image %d didn't work --- reject entire image.", i); + numRejected++; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_REJECT; + } else { + float frac = reject->n / (float)(options->numCols * options->numRows); // Pixel fraction + psLogMsg("ppStack", PS_LOG_INFO, "%ld pixels rejected from image %d (%.1f%%)", + reject->n, i, frac * 100.0); + if (frac > imageRej) { + psWarning("Image %d rejected completely because rejection fraction (%.3f) " + "exceeds limit (%.3f)", i, frac, imageRej); + psFree(reject); + // reject == NULL means reject image completely + reject = NULL; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD; + numRejected++; + } + } + + // Images without a list of rejected pixels (the list may be empty) are rejected completely + options->rejected->data[i] = reject; + + if (options->stats) { + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_REJECT", PS_META_DUPLICATE_OK, + "Time to perform rejection", psTimerMark("PPSTACK_REJECT")); + psMetadataAddS32(options->stats, PS_LIST_TAIL, "REJECT_PIXELS", PS_META_DUPLICATE_OK, + "Number of pixels rejected", reject ? reject->n : 0); + } + psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i, + psTimerClear("PPSTACK_REJECT")); + } + + psFree(options->inspect); options->inspect = NULL; + psFree(options->subKernels); options->subKernels = NULL; + psFree(options->subRegions); options->subRegions = NULL; + + if (numRejected >= num - 1) { + psError(PS_ERR_UNKNOWN, true, "All inputs completely rejected; unable to proceed."); + return false; + } + + if (options->stats) { + psMetadataAddS32(options->stats, PS_LIST_TAIL, "REJECT_IMAGES", 0, + "Number of images rejected completely", numRejected); + } + + return true; +} Index: /trunk/ppStack/src/ppStackSetup.c =================================================================== --- /trunk/ppStack/src/ppStackSetup.c (revision 23341) +++ /trunk/ppStack/src/ppStackSetup.c (revision 23341) @@ -0,0 +1,84 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppStack.h" +#include "ppStackLoop.h" + +bool ppStackSetup(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + + int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + options->num = num; + + bool mdok; // Status of MD lookup + const char *statsName = psMetadataLookupStr(&mdok, config->arguments, "STATS"); // Filename for statistics + if (statsName && strlen(statsName) > 0) { + psString resolved = pmConfigConvertFilename(statsName, config, true, true); // Resolved filename + options->statsFile = fopen(resolved, "w"); + if (!options->statsFile) { + psError(PS_ERR_IO, true, "Unable to open statistics file %s for writing.\n", resolved); + psFree(resolved); + return false; + } + psFree(resolved); + options->stats = psMetadataAlloc(); + } + + const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images + if (!tempDir) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe"); + return false; + } + + psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, + "OUTPUT")); // Name for temporary files + const char *tempName = basename(outputName); + if (!tempName) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files."); + psFree(outputName); + return false; + } + + const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images + const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks + const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp var maps + if (!tempImage || !tempMask || !tempVariance) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe"); + psFree(outputName); + return false; + } + + options->imageNames = psArrayAlloc(num); + options->maskNames = psArrayAlloc(num); + options->varianceNames = psArrayAlloc(num); + for (int i = 0; i < num; i++) { + psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images + psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage); + psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask); + psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance); + psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName); + options->imageNames->data[i] = imageName; + options->maskNames->data[i] = maskName; + options->varianceNames->data[i] = varianceName; + } + psFree(outputName); + + if (!pmConfigMaskSetBits(NULL, NULL, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to determine mask value."); + return false; + } + + return true; +} Index: /trunk/ppStack/src/ppStackThread.c =================================================================== --- /trunk/ppStack/src/ppStackThread.c (revision 23340) +++ /trunk/ppStack/src/ppStackThread.c (revision 23341) @@ -9,4 +9,7 @@ #include "ppStack.h" +#include "ppStackOptions.h" +#include "ppStackThread.h" + #define THREAD_WAIT 10000 // Microseconds to wait if thread is not available @@ -53,8 +56,15 @@ } -ppStackThreadData *ppStackThreadDataSetup(const psArray *cells, const psArray *imageNames, - const psArray *maskNames, const psArray *varianceNames, - const psArray *covariances, const pmConfig *config) -{ +ppStackThreadData *ppStackThreadDataSetup(const ppStackOptions *options, const pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + const psArray *cells = options->cells; // Array of input cells + const psArray *imageNames = options->imageNames; // Names of images to read + const psArray *maskNames = options->maskNames; // Names of masks to read + const psArray *varianceNames = options->varianceNames; // Names of variance maps to read + const psArray *covariances = options->covariances; // Covariance matrices (already read) + PS_ASSERT_ARRAY_NON_NULL(cells, NULL); PS_ASSERT_ARRAYS_SIZE_EQUAL(cells, imageNames, NULL); @@ -245,5 +255,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 6); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INITIAL_COMBINE", 3); task->function = &ppStackReadoutInitialThread; psThreadTaskAdd(task); @@ -259,5 +269,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 3); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task); Index: /trunk/ppStack/src/ppStackThread.h =================================================================== --- /trunk/ppStack/src/ppStackThread.h (revision 23341) +++ /trunk/ppStack/src/ppStackThread.h (revision 23341) @@ -0,0 +1,52 @@ +#ifndef PPSTACK_THREAD_H +#define PPSTACK_THREAD_H + +#include +#include + +#include "ppStackOptions.h" + +// Thread for stacking chunks +// +// Each input file contributes a readout, into which is read a chunk from that file +typedef struct { + psArray *readouts; // Input readouts to read and stack + bool read; // Has the scan been read? + bool busy; // Is the scan being processed? + int firstScan; // First row of the chunk to be read for this group + int lastScan; // Last row of the chunk to be read for this group +} ppStackThread; + +// Allocator +ppStackThread *ppStackThreadAlloc( + psArray *readouts // Inputs readouts to read and stack + ); + +// Data for threads +typedef struct { + psArray *threads; // Threads doing stacking + int lastScan; // Last row that's been read + psArray *imageFits; // FITS file pointers for images + psArray *maskFits; // FITS file pointers for masks + psArray *varianceFits; // FITS file pointers for variances +} ppStackThreadData; + +// Set up thread data +ppStackThreadData *ppStackThreadDataSetup( + const ppStackOptions *options, // Options + const pmConfig *config // Configuration + ); + +// Read chunk into the first available file thread +ppStackThread *ppStackThreadRead(bool *status, // Status of read + ppStackThreadData *stack, // Stacks available for reading + pmConfig *config, // Configuration + int numChunk, // Chunk number (only for interest) + int overlap // Overlap between subsequent scans + ); + +// Initialise the threads +void ppStackThreadInit(void); + + +#endif