Index: branches/eam_branches/ppStack.20211015/src/ppStackCombineInitial.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackCombineInitial.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackCombineInitial.c (revision 41841) @@ -9,10 +9,4 @@ psAssert(options, "Require options"); psAssert(config, "Require configuration"); - -/* if (!options->convolve) { */ -/* // No need to do initial combination when we haven't convolved */ -/* // XXX either allocate inspect and rejected here, or do not require them downstream */ -/* return true; */ -/* } */ psTimerStart("PPSTACK_INITIAL"); Index: branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c (revision 41841) +++ branches/eam_branches/ppStack.20211015/src/ppStackCombinePercent.c (revision 41841) @@ -0,0 +1,75 @@ +#include "ppStack.h" + +// This is the doomsday switch. +// #define TESTING // Enable test output + +bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config) +{ + psAssert(stack, "Require stack"); + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psTimerStart("PPSTACK_PERCENT"); + + bool status; // Status of read + int numChunk; // Number of chunks + for (numChunk = 0; true; numChunk++) { + ppStackThread *thread = ppStackThreadRead(&status, stack, config, numChunk, 0); + if (!status) { + // Something went wrong + psError(psErrorCodeLast(), false, "Unable to read chunk %d", numChunk); + return false; + } + if (!thread) { + // Nothing more to read + break; + } + + psThreadJob *job = psThreadJobAlloc("PPSTACK_PERCENT_COMBINE"); // Job to start + psArrayAdd(job->args, 1, thread); + psArrayAdd(job->args, 1, options); + psArrayAdd(job->args, 1, config); + if (!psThreadJobAddPending(job)) { + return false; + } + } + + if (!psThreadPoolWait(false, true)) { + psError(psErrorCodeLast(), false, "Unable to do initial combination."); + return false; + } + + // Harvest the jobs, gathering the inspection lists + options->inspect = psArrayAlloc(options->num); + options->rejected = 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); + options->rejected->data[i] = psArrayAllocEmpty(numChunk); + } + + psThreadJob *job; // Completed job + while ((job = psThreadJobGetDone())) { + psAssert(strcmp(job->type, "PPSTACK_PERCENT_COMBINE") == 0, + "Job has incorrect type: %s", job->type); + psFree(job); + } + + ppStackMemDump("initial"); + +#ifdef TESTING + ppStackWriteImage("combined_image_initial.fits", NULL, options->outRO->image, config); + ppStackWriteImage("combined_mask_initial.fits", NULL, options->outRO->mask, config); + ppStackWriteImage("combined_variance_initial.fits", NULL, options->outRO->variance, config); + + pmStackVisualPlotTestImage(options->outRO->image, "combined_image_initial.fits"); +#endif + + if (options->stats) { + // psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_INITIAL", PS_META_REPLACE, "Time to make initial stack", psTimerMark("PPSTACK_INITIAL")); + } + + return true; +} Index: branches/eam_branches/ppStack.20211015/src/ppStackCombinePrepare.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackCombinePrepare.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackCombinePrepare.c (revision 41841) @@ -31,7 +31,5 @@ options->expRO = pmReadoutAlloc(expCell); //Output readout } -/* else { */ -/* options->expRO = NULL; */ -/* } */ + int bkg_r0 = 0,bkg_c0 = 0; int bkg_nC = 0,bkg_nR = 0; @@ -53,6 +51,6 @@ options->bkgRO = pmReadoutAlloc(bkgCell); // BKG readout - // if (!pmHDUGenerateForFPA(options->bkgRO->parent->parent->parent)) { options->bkgRO->parent->parent->parent->hdu = pmHDUAlloc(NULL); + if (!options->bkgRO->parent->parent->parent->hdu) { fprintf(stderr,"failed to generate a HDU for this thing.\n"); @@ -68,6 +66,5 @@ ro->parent->parent->parent->concepts); - } - else { + } else { options->bkgRO = NULL; } @@ -97,7 +94,4 @@ return false; } - - - } Index: branches/eam_branches/ppStack.20211015/src/ppStackLoop.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackLoop.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackLoop.c (revision 41841) @@ -31,4 +31,10 @@ ppStackMemDump("prepare"); if (options->quality) return true; // Can't do anything else + + // if we want to skip the rejection steps, and instead use a percentile clipping: + if (options->clipPercent) { + bool result = ppStackLoopByPercent (config, options); + return result; + } // Convolve inputs @@ -88,5 +94,5 @@ pmCellFreeData(options->cells->data[i]); } - //MEH -- must uncomment back out -- unclear is should be moved after pixel rejection + // MEH -- must uncomment back out -- unclear is should be moved after pixel rejection psFree(stack); @@ -126,8 +132,7 @@ return false; } - } - else { + } else { // Since we haven't convolved, I believe we do need to normalize here. - //MEH -- see below for comment on ppStackCombineFinal and bscaleOffset + // MEH -- see below for comment on ppStackCombineFinal and bscaleOffset if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) { psError(psErrorCodeLast(), false, "Unable to perform final combination."); @@ -154,5 +159,6 @@ return false; } - ppStackFileActivation(config,PPSTACK_FILES_BKG ,false); + ppStackFileActivation(config, PPSTACK_FILES_BKG, false); + // Photometry psTrace("ppStack", 1, "Photometering stacked image....\n"); @@ -184,5 +190,5 @@ return false; } - //MEH -- also must uncomment back out.. + // MEH -- also must uncomment back out.. psFree(stack); psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); Index: branches/eam_branches/ppStack.20211015/src/ppStackLoop.h =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackLoop.h (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackLoop.h (revision 41841) @@ -7,4 +7,5 @@ ppStackOptions *options // Options for stacking ); + // Median only loop. bool ppStackMedianLoop( @@ -12,4 +13,10 @@ ppStackOptions *options // Options for stacking ); + +// stack by percentile range +bool ppStackLoopByPercent( + pmConfig *config, + ppStackOptions *options + ); // Setup @@ -92,4 +99,6 @@ +bool ppStackCombinePercent(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); + bool ppStackUpdateHeader(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); bool ppStackJPEGs(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); Index: branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c (revision 41841) +++ branches/eam_branches/ppStack.20211015/src/ppStackLoopByPercent.c (revision 41841) @@ -0,0 +1,159 @@ +#include "ppStack.h" + +bool ppStackLoopByPercent(pmConfig *config, ppStackOptions *options) +{ + assert(config); + + // Convolve inputs (this must be called even if the convolution is turned off) + psTrace("ppStack", 1, "Convolving inputs to target PSF....\n"); + if (!ppStackConvolve(options, config)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 2: Generate Convolutions and Save: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("convolve"); + if (options->quality) return true; // Can't do anything else + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + + // Define threading elements + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Prepare for combination + if (!ppStackCombinePrepare("PPSTACK.OUTPUT", "PPSTACK.OUTPUT.EXP", options->doBackground ? "PPSTACK.OUTPUT.BKGMODEL" : NULL, PPSTACK_FILES_STACK, stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to prepare for combination."); + psFree(stack); + return false; + } + + // Initial combination + psTrace("ppStack", 1, "Initial stack of convolved images....\n"); + if (!ppStackCombinePercent(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to perform initial combination."); + psFree(stack); + return false; + } + ppStackMemDump("percent"); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 3: Make Initial Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + + // Done with stack inputs for now + // XXX is this where we are leaking?? + for (int i = 0; i < options->num; i++) { + pmCellFreeData(options->cells->data[i]); + } + // MEH -- must uncomment back out -- unclear is should be moved after pixel rejection + psFree(stack); + + psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("final"); + + // Update Header + // Do this before performing photometry so that the cmf header gets all of the information. + if (!ppStackUpdateHeader(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to update header."); + psFree(stack); + return false; + } + + // Generate median background stack here. + if (options->doBackground && !ppStackCombineBackground(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to generate median of background images."); + psFree(stack); + return false; + } + ppStackFileActivation(config, PPSTACK_FILES_BKG, false); + + // Photometry + psTrace("ppStack", 1, "Photometering stacked image....\n"); + if (!ppStackPhotometry(options, config)) { + psError(psErrorCodeLast(), false, "Unable to perform photometry."); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 6: Photometry Analysis: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("photometry"); + + // Create JPEGS + if (!ppStackJPEGs(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to make jpegs."); + psFree(stack); + return false; + } + // Assemble Stats + if (!ppStackStats(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to assemble statistics."); + psFree(stack); + return false; + } + + // Clean up + psTrace("ppStack", 2, "Cleaning up after combination....\n"); + if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_STACK, PPSTACK_FILES_PHOT, true)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + // MEH -- also must uncomment back out.. + psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); + + // Unconvolved stack --- it's cheap to calculate, compared to everything else! + // XXX unconvolved stack is currently using the convolved mask! oops! + if (options->convolve) { + // Start threading + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, false); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Prepare for combination + if (!ppStackCombinePrepare("PPSTACK.UNCONV", "PPSTACK.UNCONV.EXP", NULL, PPSTACK_FILES_UNCONV, + stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to prepare for combination."); + psFree(stack); + return false; + } + + // generate the unconvolved stack. NOTE: this one must be normalized since the inputs have not been + psTrace("ppStack", 2, "Stack of unconvolved images....\n"); + if (!ppStackCombinePercent(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to perform unconvolved combination."); + psFree(stack); + return false; + } + psLogMsg("ppStack", PS_LOG_INFO, "Stage 8: Unconvolved Stack: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("unconv"); + + // Update Header + if (!ppStackUpdateHeader(stack, options, config)) { + psError(psErrorCodeLast(), false, "Unable to update header."); + psFree(stack); + return false; + } + // Clean up unconvolved stack + psTrace("ppStack", 2, "Cleaning up after unconvolved stack....\n"); + if (!ppStackCleanupFiles(stack, options, config, PPSTACK_FILES_UNCONV, PPSTACK_FILES_NONE, false)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + psFree(stack); + } + psFree(options->cells); options->cells = NULL; + + // Finish up + psTrace("ppStack", 1, "Finishing up....\n"); + if (!ppStackFinish(options, config)) { + psError(psErrorCodeLast(), false, "Unable to finish up."); + return false; + } + ppStackMemDump("finish"); + + return true; +} + Index: branches/eam_branches/ppStack.20211015/src/ppStackMatch.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackMatch.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackMatch.c (revision 41841) @@ -524,4 +524,5 @@ psFree(conv); } else { + // this branch is done if options->convolve is FALSE // Match the normalisation float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation Index: branches/eam_branches/ppStack.20211015/src/ppStackOptions.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackOptions.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackOptions.c (revision 41841) @@ -56,4 +56,9 @@ options->num = 0; options->quality = 0; + + options->clipPercent = false; + options->clipMinRange = 0.16; // 1-sigma default + options->clipMaxRange = 0.84; // 1-sigma default + options->psf = NULL; options->sumExposure = NAN; Index: branches/eam_branches/ppStack.20211015/src/ppStackOptions.h =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackOptions.h (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackOptions.h (revision 41841) @@ -16,4 +16,9 @@ int num; // Number of inputs int quality; // Bad data quality flag + + bool clipPercent; // use percentile range to clip? + psF32 clipMinRange; // min percentile + psF32 clipMaxRange; // max percentile + // Prepare pmPSF *psf; // Target PSF Index: branches/eam_branches/ppStack.20211015/src/ppStackReadout.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackReadout.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackReadout.c (revision 41841) @@ -53,4 +53,27 @@ } +bool ppStackReadoutPercentThread(psThreadJob *job) +{ + PS_ASSERT_THREAD_JOB_NON_NULL(job, false); + + psArray *args = job->args; // Arguments + ppStackThread *thread = args->data[0]; // Thread + 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 *exposures = options->exposures; // Exposure times for each image + psVector *addVariance = options->matchChi2; // Additional variance when rejecting + + job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask, + weightings, exposures, addVariance); + thread->busy = false; + + return job->results ? true : false; +} + +//////////////////////////////////// bool ppStackInspect(psThreadJob *job) @@ -316,2 +339,83 @@ return true; } + +psArray *ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + const psVector *mask, const psVector *weightings, const psVector *exposures, + const psVector *addVariance) +{ + assert(config); + assert(outRO); + assert(readouts); + assert(mask && mask->n == readouts->n && mask->type.type == PS_TYPE_VECTOR_MASK); + assert(weightings && weightings->n == readouts->n && weightings->type.type == PS_TYPE_F32); + assert(addVariance && addVariance->n == readouts->n && addVariance->type.type == PS_TYPE_F32); + static int sectionNum = 0; // Section number; for debugging outputs + + // Get the recipe values + 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 iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations + float minRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MIN.RANGE"); // Combination threshold + float maxRange = psMetadataLookupF32(NULL, recipe, "COMBINE.MAX.RANGE"); // Combination threshold + + psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size + + psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask for bad + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + + psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits + + bool status = false; + psImageMaskType maskBlank; + psString maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels + if (maskBlankStr) { + maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels + } else { + maskBlankStr = psMetadataLookupStr(&status, recipe, "MASK.BAD"); // Old name for MASK.BLANK + if (maskBlankStr) { + maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels + } else { + maskBlank = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + } + } + + int num = readouts->n; // Number of inputs + psArray *stack = psArrayAlloc(num); // Array for stacking + + for (int i = 0; i < num; i++) { + pmReadout *ro = readouts->data[i]; + if (!ro || mask->data.PS_TYPE_VECTOR_MASK_DATA[i]) { + // Bad image + continue; + } + + // Ensure there is a mask, or pmStackCombine will complain + if (!ro->mask) { + ro->mask = psImageAlloc(ro->image->numCols, ro->image->numRows, PS_TYPE_IMAGE_MASK); + psImageInit(ro->mask, 0); + } + + stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], + addVariance->data.F32[i]); + } + + if (!pmStackCombineByPercentile(outRO, stack, minRange, maxRange, maskBad, maskSuspect, maskBlank)) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); + psFree(stack); + return false; + } + + sectionNum++; + + psArray *results = psArrayAlloc(2); // Array of results + results->data[0] = inspect; + results->data[1] = reject; + + return results; +} + + Index: branches/eam_branches/ppStack.20211015/src/ppStackSetup.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackSetup.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackSetup.c (revision 41841) @@ -18,4 +18,9 @@ options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images? options->doBackground = psMetadataLookupBool(NULL, recipe, "BACKGROUND.MODEL"); + + options->clipPercent = psMetadataLookupBool(NULL, recipe, "COMBINE.PERCENT"); // use percentile range to clip? + options->clipMinRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MIN.RANGE"); // min percentile + options->clipMaxRange = psMetadataLookupF32 (NULL, recipe, "COMBINE.MAX.RANGE"); // max percentile + if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { Index: branches/eam_branches/ppStack.20211015/src/ppStackThread.c =================================================================== --- branches/eam_branches/ppStack.20211015/src/ppStackThread.c (revision 41840) +++ branches/eam_branches/ppStack.20211015/src/ppStackThread.c (revision 41841) @@ -283,4 +283,10 @@ } + { + psThreadTask *task = psThreadTaskAlloc("PPSTACK_PERCENT_COMBINE", 3); + task->function = &ppStackReadoutPercentThread; + psThreadTaskAdd(task); + psFree(task); + } return; }