Index: trunk/ppStack/src/ppStack.c =================================================================== --- trunk/ppStack/src/ppStack.c (revision 27318) +++ trunk/ppStack/src/ppStack.c (revision 27319) @@ -70,5 +70,6 @@ ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); ppStackFileActivation(config, PPSTACK_FILES_CONVOLVE, true); - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + ppStackFileActivation(config, PPSTACK_FILES_STACK, true); + ppStackFileActivation(config, PPSTACK_FILES_UNCONV, true); ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); if (!ppStackFilesIterateUp(config)) { Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 27318) +++ trunk/ppStack/src/ppStack.h (revision 27319) @@ -27,5 +27,6 @@ PPSTACK_FILES_PREPARE, // Files for preparation PPSTACK_FILES_CONVOLVE, // Files for convolution - PPSTACK_FILES_COMBINE, // Files for combination + PPSTACK_FILES_STACK, // Stack files + PPSTACK_FILES_UNCONV, // Unconvolved stack files PPSTACK_FILES_PHOT // Files for photometry } ppStackFileList; Index: trunk/ppStack/src/ppStackCleanup.c =================================================================== --- trunk/ppStack/src/ppStackCleanup.c (revision 27318) +++ trunk/ppStack/src/ppStackCleanup.c (revision 27319) @@ -6,4 +6,5 @@ #include #include +#include #include "ppStack.h" @@ -59,4 +60,18 @@ } + // 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); + } + return true; } Index: trunk/ppStack/src/ppStackCombineFinal.c =================================================================== --- trunk/ppStack/src/ppStackCombineFinal.c (revision 27318) +++ trunk/ppStack/src/ppStackCombineFinal.c (revision 27319) @@ -15,9 +15,31 @@ bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances, - ppStackOptions *options, pmConfig *config, bool safe, bool normalise) + ppStackOptions *options, pmConfig *config, bool safe, bool normalise, bool grow) { psAssert(stack, "Require stack"); psAssert(options, "Require options"); psAssert(config, "Require configuration"); + + int numCols = target->image->numCols, numRows = target->image->numRows; // Size of image + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + + // Grow the list of rejected pixels, if desired + psArray *reject = psArrayAlloc(options->num); // Pixels rejected for each image + for (int i = 0; i < options->num; i++) { + if (grow) { + reject->data[i] = pmStackRejectGrow(options->rejected->data[i], numCols, numRows, poorFrac, + options->regions->data[i], options->kernels->data[i]); + if (!reject->data[i]) { + psError(psErrorCodeLast(), false, "Unable to grow rejected pixels for image %d", i); + psFree(reject); + return false; + } + } else { + reject->data[i] = psMemIncrRefCounter(options->rejected->data[i]); + } + } if (!target->mask) { @@ -32,4 +54,5 @@ // Something went wrong psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk); + psFree(reject); return false; } @@ -43,4 +66,5 @@ psArrayAdd(job->args, 1, target); psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, reject); psArrayAdd(job->args, 1, options); psArrayAdd(job->args, 1, config); @@ -49,4 +73,5 @@ if (!psThreadJobAddPending(job)) { psFree(job); + psFree(reject); return false; } @@ -56,6 +81,9 @@ if (!psThreadPoolWait(true)) { psError(psErrorCodeLast(), false, "Unable to do final combination."); + psFree(reject); return false; } + + psFree(reject); // Sum covariance matrices Index: trunk/ppStack/src/ppStackCombinePrepare.c =================================================================== --- trunk/ppStack/src/ppStackCombinePrepare.c (revision 27318) +++ trunk/ppStack/src/ppStackCombinePrepare.c (revision 27319) @@ -10,5 +10,6 @@ #include "ppStackLoop.h" -bool ppStackCombinePrepare(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +bool ppStackCombinePrepare(pmReadout **ro, const char *name, ppStackFileList files, ppStackThreadData *stack, + ppStackOptions *options, pmConfig *config) { psAssert(stack, "Require stack"); @@ -25,5 +26,5 @@ pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + ppStackFileActivation(config, files, true); pmFPAview *view = ppStackFilesIterateDown(config); // View to readout if (!view) { @@ -31,9 +32,6 @@ } - pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); // Output cell - options->outRO = pmReadoutAlloc(outCell); // Output readout - - pmCell *unconvCell = pmFPAfileThisCell(config->files, view, "PPSTACK.UNCONV"); // Unconvolved cell - options->unconvRO = pmReadoutAlloc(unconvCell); // Unconvolved readout + pmCell *cell = pmFPAfileThisCell(config->files, view, name); // Output cell + *ro = pmReadoutAlloc(cell); // Output readout psFree(view); @@ -44,16 +42,9 @@ psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels - if (!pmReadoutStackDefineOutput(options->outRO, col0, row0, numCols, numRows, true, true, maskBad)) { + if (!pmReadoutStackDefineOutput(*ro, col0, row0, numCols, numRows, true, true, maskBad)) { psError(PPSTACK_ERR_ARGUMENTS, false, "Unable to prepare output."); return false; } - options->unconvRO->image = psImageCopy(NULL, options->outRO->image, PS_TYPE_F32); -// options->unconvRO->mask = psImageCopy(NULL, options->outRO->mask, PS_TYPE_IMAGE_MASK); - options->unconvRO->col0 = options->outRO->col0; - options->unconvRO->row0 = options->outRO->row0; - - - return true; } Index: trunk/ppStack/src/ppStackFiles.c =================================================================== --- trunk/ppStack/src/ppStackFiles.c (revision 27318) +++ trunk/ppStack/src/ppStackFiles.c (revision 27319) @@ -21,8 +21,9 @@ // "PPSTACK.CONV.KERNEL", NULL }; -/// Output files for the combination -static char *filesCombine[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", - "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", - "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; +/// Regular (convolved) stack files +static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", + "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", NULL }; +/// Unconvolved stack files +static char *filesUnconv[] = { "PPSTACK.UNCONV", "PPSTACK.UNCONV.MASK", "PPSTACK.UNCONV.VARIANCE", NULL }; /// Files for photometry @@ -37,5 +38,6 @@ case PPSTACK_FILES_PREPARE: return filesPrepare; case PPSTACK_FILES_CONVOLVE: return filesConvolve; - case PPSTACK_FILES_COMBINE: return filesCombine; + case PPSTACK_FILES_STACK: return filesStack; + case PPSTACK_FILES_UNCONV: return filesUnconv; case PPSTACK_FILES_PHOT: return filesPhot; default: Index: trunk/ppStack/src/ppStackFinish.c =================================================================== --- trunk/ppStack/src/ppStackFinish.c (revision 27318) +++ trunk/ppStack/src/ppStackFinish.c (revision 27319) @@ -7,5 +7,4 @@ #include #include -#include #include @@ -49,22 +48,4 @@ } - - // 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; return true; Index: trunk/ppStack/src/ppStackLoop.c =================================================================== --- trunk/ppStack/src/ppStackLoop.c (revision 27318) +++ trunk/ppStack/src/ppStackLoop.c (revision 27319) @@ -111,5 +111,6 @@ // Prepare for combination - if (!ppStackCombinePrepare(stack, options, config)) { + if (!ppStackCombinePrepare(&options->outRO, "PPSTACK.OUTPUT", PPSTACK_FILES_STACK, + stack, options, config)) { psError(psErrorCodeLast(), false, "Unable to prepare for combination."); psFree(stack); @@ -148,5 +149,6 @@ // Final combination psTrace("ppStack", 2, "Final stack of convolved images....\n"); - if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, false, false)) { + if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, + false, false, true)) { psError(psErrorCodeLast(), false, "Unable to perform final combination."); psFree(stack); @@ -163,5 +165,27 @@ return false; } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + + // Photometry + psTrace("ppStack", 1, "Photometering stacked image....\n"); + if (!ppStackPhotometry(options, config)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("photometry"); + + if (!ppStackFilesIterateUp(config)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + return false; + } + ppStackFileActivation(config, PPSTACK_FILES_STACK, false); + ppStackFileActivation(config, PPSTACK_FILES_PHOT, false); + options->outRO->data_exists = false; + options->outRO->parent->data_exists = false; + options->outRO->parent->parent->data_exists = false; + psFree(options->outRO); + options->outRO = NULL; + + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("cleanup"); @@ -177,12 +201,34 @@ return false; } + + // Prepare for combination + if (!ppStackCombinePrepare(&options->unconvRO, "PPSTACK.UNCONV", PPSTACK_FILES_UNCONV, + stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to prepare for combination."); + psFree(stack); + return false; + } + psTrace("ppStack", 2, "Stack of unconvolved images....\n"); - if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, false, true)) { + if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, + false, true, false)) { psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); psFree(stack); return false; } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("unconv"); + + if (!ppStackFilesIterateUp(config)) { + psError(psErrorCodeLast(), false, "Unable to close files."); + psFree(stack); + return false; + } + ppStackFileActivation(config, PPSTACK_FILES_UNCONV, false); + options->unconvRO->data_exists = false; + options->unconvRO->parent->data_exists = false; + options->unconvRO->parent->parent->data_exists = false; + psFree(options->unconvRO); + options->unconvRO = NULL; psFree(stack); @@ -190,13 +236,4 @@ psFree(options->cells); options->cells = NULL; #endif - - // Photometry - psTrace("ppStack", 1, "Photometering stacked image....\n"); - if (!ppStackPhotometry(options, config)) { - psError(psErrorCodeLast(), false, "Unable to perform photometry."); - return false; - } - psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); - ppStackMemDump("photometry"); // Finish up Index: trunk/ppStack/src/ppStackLoop.h =================================================================== --- trunk/ppStack/src/ppStackLoop.h (revision 27318) +++ trunk/ppStack/src/ppStackLoop.h (revision 27319) @@ -7,4 +7,5 @@ #include "ppStackOptions.h" #include "ppStackThread.h" +#include "ppStack.h" @@ -37,4 +38,7 @@ // Prepare for combination bool ppStackCombinePrepare( + pmReadout **readout, // Readout to set + const char *name, // Name of file + ppStackFileList files, // Files of interest ppStackThreadData *stack, // Stack ppStackOptions *options, // Options @@ -63,5 +67,6 @@ pmConfig *config, // Configuration bool safe, // Allow safe combination? - bool norm // Normalise images? + bool norm, // Normalise images? + bool grow // Grow rejection masks? ); Index: trunk/ppStack/src/ppStackPhotometry.c =================================================================== --- trunk/ppStack/src/ppStackPhotometry.c (revision 27318) +++ trunk/ppStack/src/ppStackPhotometry.c (revision 27319) @@ -26,8 +26,8 @@ psTimerStart("PPSTACK_PHOT"); - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, false); + pmFPAfileActivate(config->files, false, NULL); ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout - ppStackFileActivation(config, PPSTACK_FILES_COMBINE, true); + ppStackFileActivation(config, PPSTACK_FILES_STACK, true); pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 27318) +++ trunk/ppStack/src/ppStackReadout.c (revision 27319) @@ -39,16 +39,16 @@ pmReadout *target = args->data[0]; // Output readout ppStackThread *thread = args->data[1]; // Thread - ppStackOptions *options = args->data[2]; // Options - pmConfig *config = args->data[3]; // Configuration - bool safety = PS_SCALAR_VALUE(args->data[4], U8); // Safety switch on? - bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images? + psArray *reject = args->data[2]; // Rejected pixels for each image + ppStackOptions *options = args->data[3]; // Options + pmConfig *config = args->data[4]; // Configuration + bool safety = PS_SCALAR_VALUE(args->data[5], U8); // Safety switch on? + bool normalise = PS_SCALAR_VALUE(args->data[6], U8); // Normalise images? 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 psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images - bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected, + bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, reject, weightings, addVariance, safety, norm); // Status of operation Index: trunk/ppStack/src/ppStackReject.c =================================================================== --- trunk/ppStack/src/ppStackReject.c (revision 27318) +++ trunk/ppStack/src/ppStackReject.c (revision 27319) @@ -31,5 +31,4 @@ 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 @@ -109,5 +108,5 @@ psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, - threshold, poorFrac, stride, options->regions->data[i], + threshold, stride, options->regions->data[i], options->kernels->data[i]); // Rejected pixels @@ -130,11 +129,9 @@ options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_BAD; numRejected++; + } else { + // Add to list of pixels already rejected + reject = psPixelsConcatenate(reject, options->rejected->data[i]); + options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject); } - } - - if (reject) { - // Add to list of pixels already rejected - reject = psPixelsConcatenate(reject, options->rejected->data[i]); - options->rejected->data[i] = psPixelsDuplicates(options->rejected->data[i], reject); } @@ -170,6 +167,4 @@ psFree(options->inspect); options->inspect = NULL; - psFree(options->kernels); options->kernels = NULL; - psFree(options->regions); options->regions = NULL; if (numRejected >= num) { Index: trunk/ppStack/src/ppStackThread.c =================================================================== --- trunk/ppStack/src/ppStackThread.c (revision 27318) +++ trunk/ppStack/src/ppStackThread.c (revision 27319) @@ -284,5 +284,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task); Index: trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.c (revision 27318) +++ trunk/psModules/src/imcombine/pmStackReject.c (revision 27319) @@ -10,4 +10,6 @@ #include "pmSubtractionThreads.h" #include "pmSubtractionKernels.h" + +#include "pmStackReject.h" #define PIXEL_LIST_BUFFER 100 // Number of pixels to add to list at a time @@ -115,6 +117,6 @@ -psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, float poorFrac, - int stride, const psArray *subRegions, const psArray *subKernels) +psPixels *pmStackReject(const psPixels *in, int numCols, int numRows, float threshold, int stride, + const psArray *subRegions, const psArray *subKernels) { PS_ASSERT_PIXELS_NON_NULL(in, NULL); @@ -223,6 +225,17 @@ psTrace("psModules.imcombine", 7, "Found %ld bad pixels", bad->n); - // Now, grow the mask to include everything that touches a bad pixel in the convolution - psImage *source = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), + return bad; +} + + +psPixels *pmStackRejectGrow(const psPixels *in, int numCols, int numRows, float poorFrac, + const psArray *subRegions, const psArray *subKernels) +{ + PS_ASSERT_PIXELS_NON_NULL(in, NULL); + PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL); + PS_ASSERT_ARRAY_NON_NULL(subKernels, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, subKernels, NULL); + + psImage *source = psPixelsToMask(NULL, in, psRegionSet(0, numCols - 1, 0, numRows - 1), PM_STACK_MASK_BAD); // Mask image to grow @@ -243,5 +256,5 @@ bool oldThreads = psImageConvolveSetThreads(false); // Old value of threading for psImageColvolve - psImage *target = psImageRecycle(convolved, numCols, numRows, PS_TYPE_IMAGE_MASK); // Grown image + psImage *target = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); // Grown image psImageInit(target, 0); for (int i = 0; i < subRegions->n; i++) { @@ -326,5 +339,5 @@ psFree(source); - bad = psPixelsFromMask(bad, target, PM_STACK_MASK_ALL); + psPixels *bad = psPixelsFromMask(NULL, target, PM_STACK_MASK_ALL); // All bad pixels psFree(target); psTrace("psModules.imcombine", 7, "Total %ld bad pixels", bad->n); Index: trunk/psModules/src/imcombine/pmStackReject.h =================================================================== --- trunk/psModules/src/imcombine/pmStackReject.h (revision 27318) +++ trunk/psModules/src/imcombine/pmStackReject.h (revision 27319) @@ -12,8 +12,16 @@ int numCols, int numRows, ///< Size of image of interest float threshold, ///< Threshold on convolved image, 0..1 - float poorFrac, ///< Fraction for "poor" int stride, ///< Size of convolution patches const psArray *regions, ///< Array of image regions for image const psArray *kernels ///< Array of kernel parameters for each region + ); + +/// Given a list of pixels from the convolved image, we grow them by convolution to get the list of all pixels +/// which should be rejected. +psPixels *pmStackRejectGrow(const psPixels *in, ///< List of pixels in the convolved image + int numCols, int numRows, ///< Size of image of interest + float poorFrac, ///< Fraction for "poor" + const psArray *regions, ///< Array of image regions for image + const psArray *kernels ///< Array of kernel parameters for each region );