Index: /trunk/ppStack/src/Makefile.am =================================================================== --- /trunk/ppStack/src/Makefile.am (revision 42093) +++ /trunk/ppStack/src/Makefile.am (revision 42094) @@ -29,4 +29,6 @@ ppStackPrepare.c \ ppStackConvolve.c \ + ppStackLoopByPercent.c \ + ppStackCombinePercent.c \ ppStackCombinePrepare.c \ ppStackCombineInitial.c \ Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 42093) +++ /trunk/ppStack/src/ppStack.h (revision 42094) @@ -115,4 +115,17 @@ bool ppStackReadoutFinalThread(psThreadJob *job // Job to process ); + +// Thread entry point for ppStackReadoutPercent +bool ppStackReadoutPercentThread(psThreadJob *job); + +bool ppStackReadoutPercent(const pmConfig *config, + pmReadout *outRO, + pmReadout *expRO, + const psArray *readouts, + const psVector *mask, + const psVector *weightings, + const psVector *exposures, + const psVector *addVariance); + // Perform median stacking for background bool ppStackCombineBackground(ppStackThreadData *stack, ppStackOptions *options, pmConfig *config); Index: /trunk/ppStack/src/ppStackArguments.c =================================================================== --- /trunk/ppStack/src/ppStackArguments.c (revision 42093) +++ /trunk/ppStack/src/ppStackArguments.c (revision 42094) @@ -197,5 +197,5 @@ } - // stack-type : used to define the stack for PSPS + // stack-type : used to define the stack for PSPS (added as metadata to header key STK_TYPE) if ((argNum = psArgumentGet (argc, argv, "-stack-type"))) { if (argc <= argNum+1) { Index: /trunk/ppStack/src/ppStackCombineFinal.c =================================================================== --- /trunk/ppStack/src/ppStackCombineFinal.c (revision 42093) +++ /trunk/ppStack/src/ppStackCombineFinal.c (revision 42094) @@ -2,5 +2,5 @@ // This is the doomsday switch. -// #define TESTING // Enable test output +// # define TESTING // Enable test output //MEH -- adhoc addition to blank mask border of final stack since rejection different/none on order of KERNEL.SIZE with overlap in CombineInitial Index: /trunk/ppStack/src/ppStackCombineInitial.c =================================================================== --- /trunk/ppStack/src/ppStackCombineInitial.c (revision 42093) +++ /trunk/ppStack/src/ppStackCombineInitial.c (revision 42094) @@ -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: /trunk/ppStack/src/ppStackCombinePercent.c =================================================================== --- /trunk/ppStack/src/ppStackCombinePercent.c (revision 42094) +++ /trunk/ppStack/src/ppStackCombinePercent.c (revision 42094) @@ -0,0 +1,62 @@ +#include "ppStack.h" + +# define USE_THREADS 1 +// #define TESTING + +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 (USE_THREADS) { + if (!psThreadJobAddPending(job)) { + return false; + } + } else { + ppStackReadoutPercentThread(job); + } + } + + if (USE_THREADS && !psThreadPoolWait(true, true)) { + psError(psErrorCodeLast(), false, "Unable to do initial combination."); + return false; + } + + 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: /trunk/ppStack/src/ppStackCombinePrepare.c =================================================================== --- /trunk/ppStack/src/ppStackCombinePrepare.c (revision 42093) +++ /trunk/ppStack/src/ppStackCombinePrepare.c (revision 42094) @@ -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: /trunk/ppStack/src/ppStackFiles.c =================================================================== --- /trunk/ppStack/src/ppStackFiles.c (revision 42093) +++ /trunk/ppStack/src/ppStackFiles.c (revision 42094) @@ -1,3 +1,4 @@ #include "ppStack.h" +# define DEBUG 0 // Here follows lists of files for activation/deactivation at various stages. Each must be NULL-terminated. @@ -28,5 +29,4 @@ static char *filesStack[] = { "PPSTACK.OUTPUT", "PPSTACK.OUTPUT.MASK", "PPSTACK.OUTPUT.VARIANCE", "PPSTACK.OUTPUT.EXP", "PPSTACK.OUTPUT.EXPNUM", "PPSTACK.OUTPUT.EXPWT", - "PPSTACK.OUTPUT.JPEG1", "PPSTACK.OUTPUT.JPEG2", "PPSTACK.OUTPUT.BKGMODEL", NULL }; @@ -45,14 +45,14 @@ { switch (list) { - case PPSTACK_FILES_NONE: return filesNOP; - case PPSTACK_FILES_PREPARE: return filesPrepare; - case PPSTACK_FILES_TARGET: return filesTarget; - case PPSTACK_FILES_CONVOLVE: return filesConvolve; - case PPSTACK_FILES_STACK: return filesStack; - case PPSTACK_FILES_UNCONV: return filesUnconv; - case PPSTACK_FILES_PHOT: return filesPhot; - case PPSTACK_FILES_BKG: return filesBkg; - case PPSTACK_FILES_MEDIAN_IN: return filesMedianIn; - case PPSTACK_FILES_MEDIAN_OUT: return filesMedianOut; + case PPSTACK_FILES_NONE: return filesNOP; + case PPSTACK_FILES_PREPARE: return filesPrepare; + case PPSTACK_FILES_TARGET: return filesTarget; + case PPSTACK_FILES_CONVOLVE: return filesConvolve; + case PPSTACK_FILES_STACK: return filesStack; + case PPSTACK_FILES_UNCONV: return filesUnconv; + case PPSTACK_FILES_PHOT: return filesPhot; + case PPSTACK_FILES_BKG: return filesBkg; + case PPSTACK_FILES_MEDIAN_IN: return filesMedianIn; + case PPSTACK_FILES_MEDIAN_OUT: return filesMedianOut; default: psAbort("Unrecognised file list: %x", list); @@ -99,4 +99,11 @@ char **files = stackFiles(list); // Files to turn on/off for (int i = 0; files[i] != NULL; i++) { + if (DEBUG) { + if (state) { + fprintf (stderr, "activate %s\n", files[i]); + } else { + fprintf (stderr, "deativate %s\n", files[i]); + } + } pmFPAfileActivate(config->files, state, files[i]); } Index: /trunk/ppStack/src/ppStackJPEGs.c =================================================================== --- /trunk/ppStack/src/ppStackJPEGs.c (revision 42093) +++ /trunk/ppStack/src/ppStackJPEGs.c (revision 42094) @@ -33,4 +33,7 @@ psFree(ro2); + pmFPAfileActivate(config->files, true, "PPSTACK.OUTPUT.JPEG1"); + pmFPAfileActivate(config->files, true, "PPSTACK.OUTPUT.JPEG2"); + return true; } Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 42093) +++ /trunk/ppStack/src/ppStackLoop.c (revision 42094) @@ -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,9 +132,11 @@ 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 - if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, true, true, true)) { + // MEH -- see below for comment on ppStackCombineFinal and bscaleOffset + // last 4 parameters are: safe, normalize, grow, bscaleoffset + // EAM : looks like we do NOT need to normalize here eithr because ppStackConvolve is still called and + // the inputs are normalized + if (!ppStackCombineFinal(stack, options->convCovars, options, config, false, false, true, true)) { psError(psErrorCodeLast(), false, "Unable to perform final combination."); psFree(stack); @@ -154,5 +162,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 +193,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: /trunk/ppStack/src/ppStackLoop.h =================================================================== --- /trunk/ppStack/src/ppStackLoop.h (revision 42093) +++ /trunk/ppStack/src/ppStackLoop.h (revision 42094) @@ -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: /trunk/ppStack/src/ppStackLoopByPercent.c =================================================================== --- /trunk/ppStack/src/ppStackLoopByPercent.c (revision 42094) +++ /trunk/ppStack/src/ppStackLoopByPercent.c (revision 42094) @@ -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 because + // it performs the normalization of the inputs). + 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 + + // Define threading elements + ppStackThreadData *stack = ppStackThreadDataSetup(options, config, true); + if (!stack) { + psError(psErrorCodeLast(), false, "Unable to initialise stack threads."); + return false; + } + + // Prepare for combination (define outputs, generate output cells) + // NOTE: expmaps will be skipped if "PPSTACK.OUTPUT.EXP" -> NULL + ppStackFileList stackFiles = options->convolve ? PPSTACK_FILES_STACK : PPSTACK_FILES_UNCONV; + char *outName = options->convolve ? "PPSTACK.OUTPUT" : "PPSTACK.UNCONV"; + char *expName = options->convolve ? "PPSTACK.OUTPUT.EXP" : "PPSTACK.UNCONV.EXP"; + char *bckName = options->doBackground ? "PPSTACK.OUTPUT.BKGMODEL" : NULL; + + if (!ppStackCombinePrepare(outName, expName, bckName, stackFiles, 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 + for (int i = 0; i < options->num; i++) { + pmCellFreeData(options->cells->data[i]); + } + + 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, stackFiles, PPSTACK_FILES_PHOT, true)) { + psError(psErrorCodeLast(), false, "Unable to clean up."); + psFree(stack); + return false; + } + psFree(stack); + psLogMsg("ppStack", PS_LOG_INFO, "Stage 7: Cleanup, WCS & JPEGS: %f sec", psTimerClear("PPSTACK_STEPS")); + ppStackMemDump("cleanup"); + + // Unconvolved stack --- it's cheap to calculate, compared to everything else! + // XXX unconvolved stack is currently using the convolved mask! oops! (EAM: is this true?) + 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: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 42093) +++ /trunk/ppStack/src/ppStackMatch.c (revision 42094) @@ -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: /trunk/ppStack/src/ppStackOptions.c =================================================================== --- /trunk/ppStack/src/ppStackOptions.c (revision 42093) +++ /trunk/ppStack/src/ppStackOptions.c (revision 42094) @@ -56,4 +56,7 @@ options->num = 0; options->quality = 0; + + options->clipPercent = false; + options->psf = NULL; options->sumExposure = NAN; Index: /trunk/ppStack/src/ppStackOptions.h =================================================================== --- /trunk/ppStack/src/ppStackOptions.h (revision 42093) +++ /trunk/ppStack/src/ppStackOptions.h (revision 42094) @@ -16,4 +16,7 @@ int num; // Number of inputs int quality; // Bad data quality flag + + bool clipPercent; // use percentile range to clip? + // Prepare pmPSF *psf; // Target PSF Index: /trunk/ppStack/src/ppStackPhotometry.c =================================================================== --- /trunk/ppStack/src/ppStackPhotometry.c (revision 42093) +++ /trunk/ppStack/src/ppStackPhotometry.c (revision 42094) @@ -19,5 +19,7 @@ ppStackFileActivation(config, PPSTACK_FILES_PHOT, true); pmFPAview *photView = ppStackFilesIterateDown(config); // View to readout - ppStackFileActivation(config, PPSTACK_FILES_STACK, true); + + ppStackFileList stackFiles = options->convolve ? PPSTACK_FILES_STACK : PPSTACK_FILES_UNCONV; + ppStackFileActivation(config, stackFiles, true); pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); // File for photometry @@ -57,5 +59,5 @@ // psphotSetNFrames pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); - if (photRO && options->expRO->mask) { + if (photRO && options->expRO && options->expRO->mask) { psMetadataAddImage(photRO->analysis, PS_LIST_TAIL, "EXPNUM", PS_META_REPLACE, "EXPNUM image", options->expRO->mask); } Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 42093) +++ /trunk/ppStack/src/ppStackReadout.c (revision 42094) @@ -19,4 +19,8 @@ weightings, exposures, addVariance); thread->busy = false; + + thread->status = job->results ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE; + + psAssert(job->results, "Stacking failed."); return job->results ? true : false; @@ -48,4 +52,6 @@ thread->busy = false; + thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE; + psAssert(status, "Stacking failed."); @@ -53,4 +59,36 @@ } +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 + pmReadout *expRO = options->expRO; // Expmap 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 + + // fprintf (stderr, "start thread %d : scan range %d - %d\n", thread->entry, thread->firstScan, thread->lastScan); + + bool status = ppStackReadoutPercent(config, outRO, expRO, thread->readouts, mask, weightings, exposures, addVariance); + + thread->busy = false; + + thread->status = status ? PPSTACK_THREAD_SUCCESS : PPSTACK_THREAD_FAILURE; + + // fprintf (stderr, "finished thread %d, status %d : scan range %d - %d\n", thread->entry, thread->status, thread->firstScan, thread->lastScan); + + psAssert(status, "Stacking failed."); + + return status; +} + +//////////////////////////////////// bool ppStackInspect(psThreadJob *job) @@ -145,5 +183,5 @@ int num = readouts->n; // Number of inputs - psArray *stack = psArrayAlloc(num); // Array for stacking + psArray *stackData = psArrayAlloc(num); // Array of data to be stacked for (int i = 0; i < num; i++) { @@ -160,12 +198,19 @@ } - stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], + // stackData is an array of pmStackData structures + stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], addVariance->data.F32[i]); } - if (!pmStackCombine(outRO, NULL, stack, maskBad, maskSuspect, maskBlank, kernelSize, iter, + if (!pmStackCombine(outRO, NULL, stackData, maskBad, maskSuspect, maskBlank, kernelSize, iter, combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); - psFree(stack); + psFree(stackData); + + // XXX EAM : if pmStackCombine fails above this return will cause the thread to block. + // The failure status results in a psThread.fault. the psThreadLauncher function, which + // each background thread runs as a loop, will sleep until the fault is cleared externally + // by the handler. psThreadPoolWait() catches, counts, and clears these faults, but + // ppStackThreadRead does not. return false; } @@ -175,5 +220,5 @@ psArray *reject = psArrayAlloc(num); // List of pixels rejected for (int i = 0; i < num; i++) { - pmStackData *data = stack->data[i]; // Data for this image + pmStackData *data = stackData->data[i]; // Data for this image if (!data) { continue; @@ -186,5 +231,5 @@ reject->data[i] = psMemIncrRefCounter(data->reject); } - psFree(stack); + psFree(stackData); //MEH change to trace @@ -247,5 +292,5 @@ int num = readouts->n; // Number of inputs - psArray *stack = psArrayAlloc(num); // Array for stacking + psArray *stackData = psArrayAlloc(num); // Array of data to be stacked // We have rejection from a previous combination: combine without flagging pixels to inspect @@ -272,5 +317,5 @@ addVariance ? addVariance->data.F32[i] : NAN); data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL; - stack->data[i] = data; + stackData->data[i] = data; //MEH -- apply bscale offset before norm @@ -288,8 +333,8 @@ } - if (!pmStackCombine(outRO, expRO, stack, maskBad, maskSuspect, maskBlank, 0, iter, combineRej, + if (!pmStackCombine(outRO, expRO, stackData, maskBad, maskSuspect, maskBlank, 0, iter, combineRej, combineSys, combineDiscard, useVariance, safe, nminpix, rejected)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); - psFree(stack); + psFree(stackData); return false; } @@ -307,5 +352,5 @@ expChip->data_exists = true; - psFree(stack); + psFree(stackData); //MEH change to trace @@ -316,2 +361,86 @@ return true; } + +// NOTE: expRO is allowed to be NULL +bool ppStackReadoutPercent(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, 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 + + // fprintf (stderr, "starting ReadoutPercent, %d\n", sectionNum); + + // 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 + psF64 rejectFraction = psMetadataLookupF32(&mdok, recipe, "COMBINE.REJECT.FRACTION"); // fraction of outliers to reject + + char defaultBlankStr[16] = "BLANK"; + + psString maskBadStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask for bad + psString maskSuspectStr = psMetadataLookupStr(&mdok, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psString maskBlankStr = psMetadataLookupStr(&mdok, recipe, "MASK.BLANK"); // Name of bits to set for empty pixels + if (!maskBlankStr) { + maskBlankStr = psMetadataLookupStr(&mdok, recipe, "MASK.BAD"); // Old name for MASK.BLANK + } + if (!maskBlankStr) { + maskBlankStr = defaultBlankStr; // this is statically allocated above + } + + psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits + psImageMaskType maskBlank = pmConfigMaskGet(maskBlankStr, config); // Bits to mask for bad pixels + + int num = readouts->n; // Number of inputs + psArray *stackData = psArrayAlloc(num); // Array of data to be stacked + + for (int i = 0; i < num; i++) { + stackData->data[i] = NULL; + + 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); + } + + stackData->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], + addVariance->data.F32[i]); + } + + if (!pmStackCombineByPercentile(outRO, expRO, stackData, rejectFraction, maskBad, maskSuspect, maskBlank)) { + psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); + psFree(stackData); + return false; + } + + outRO->data_exists = true; // output readout + outRO->parent->data_exists = true; // output cell + outRO->parent->parent->data_exists = true; // output chip + + if (expRO) { + expRO->data_exists = true; // expmap readout + expRO->parent->data_exists = true; // expmap cell + expRO->parent->parent->data_exists = true; // expmap chip + } + + psFree(stackData); + sectionNum++; + + return true; +} + + Index: /trunk/ppStack/src/ppStackSetup.c =================================================================== --- /trunk/ppStack/src/ppStackSetup.c (revision 42093) +++ /trunk/ppStack/src/ppStackSetup.c (revision 42094) @@ -18,4 +18,6 @@ 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? if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { Index: /trunk/ppStack/src/ppStackSources.c =================================================================== --- /trunk/ppStack/src/ppStackSources.c (revision 42093) +++ /trunk/ppStack/src/ppStackSources.c (revision 42094) @@ -341,7 +341,9 @@ options->norm = psVectorAlloc(num, PS_TYPE_F32); for (int i = 0; i < num; i++) { - if (!isfinite(trans->data.F32[i])) { - continue; - } + if (!isfinite(trans->data.F32[i])) { + psLogMsg("ppStack", PS_LOG_INFO, "Non-finite transparency, skipping correction for image %d: %f mag (%f) (with inputMask value %d)\n", + i, NAN, trans->data.F32[i], inputMask->data.U8[i]); + continue; + } psArray *sources = sourceLists->data[i]; // Sources of interest float magCorr = zp->data.F32[i] - trans->data.F32[i] - 2.5*log10(options->sumExposure) - airmassTerm * airmassTarget; Index: /trunk/ppStack/src/ppStackThread.c =================================================================== --- /trunk/ppStack/src/ppStackThread.c (revision 42093) +++ /trunk/ppStack/src/ppStackThread.c (revision 42094) @@ -19,4 +19,6 @@ thread->firstScan = 0; thread->lastScan = 0; + thread->entry = -1; // this gets assigned after allocation + thread->status = PPSTACK_THREAD_NEW; return thread; } @@ -124,5 +126,7 @@ readouts->data[j] = ro; } - threads->data[i] = ppStackThreadAlloc(readouts); + ppStackThread *thread = ppStackThreadAlloc(readouts); + thread->entry = i; + threads->data[i] = thread; psFree(readouts); // Drop reference } @@ -157,6 +161,8 @@ while (true) { // check for any groups which can read data + // NOTE: threads->n is set (in ppStackThreadDataSetup) to be 1 more than the argument to -threads for (int j = 0; j < threads->n; j++) { ppStackThread *thread = threads->data[j]; + // fprintf (stderr, "grab thread %d, scan: %d - %d : busy: %d read: %d status: %d\n", thread->entry, thread->firstScan, thread->lastScan, thread->busy, thread->read, thread->status); if (thread->read) { continue; @@ -166,4 +172,5 @@ thread->lastScan = stack->lastScan + rows; stack->lastScan = thread->lastScan; + thread->status = PPSTACK_THREAD_READ; psArray *readouts = thread->readouts; @@ -172,4 +179,5 @@ psTrace("ppStack", 2, "Reading data for chunk %d into group %d....\n", numChunk, j); + // fprintf (stderr, "Reading data for chunk %d into group %d (thread %d, scan: %d - %d)....\n", numChunk, j, thread->entry, thread->firstScan, thread->lastScan); for (int i = 0; i < readouts->n; i++) { pmReadout *ro = readouts->data[i]; // Input readout @@ -178,15 +186,17 @@ } + // fprintf (stderr, "read data for readout %d\n", i); + // override the recorded last scan - ro->thisImageScan = thread->firstScan; + ro->thisImageScan = thread->firstScan; ro->thisVarianceScan = thread->firstScan; - ro->thisMaskScan = thread->firstScan; - ro->lastImageScan = thread->lastScan; - ro->lastMaskScan = thread->lastScan; + ro->thisMaskScan = thread->firstScan; + ro->lastImageScan = thread->lastScan; + ro->lastMaskScan = thread->lastScan; ro->lastVarianceScan = thread->lastScan; - ro->forceScan = true; - - psFits *imageFits = stack->imageFits->data[i]; // FITS file for image - psFits *maskFits = stack->maskFits->data[i]; // FITS file for mask + ro->forceScan = true; + + psFits *imageFits = stack->imageFits->data[i]; // FITS file for image + psFits *maskFits = stack->maskFits->data[i]; // FITS file for mask psFits *varianceFits = stack->varianceFits->data[i]; // FITS file for variance @@ -232,4 +242,5 @@ thread->read = thread->busy = true; + thread->status = PPSTACK_THREAD_RUN; return thread; } @@ -239,7 +250,10 @@ for (int j = 0; j < threads->n; j++) { ppStackThread *thread = threads->data[j]; + // fprintf (stderr, "test thread %d, scan: %d - %d : busy: %d, read: %d, status: %d\n", thread->entry, thread->firstScan, thread->lastScan, thread->busy, thread->read, thread->status); + if (thread->busy) { continue; } + thread->status = PPSTACK_THREAD_NEW; thread->read = false; wait = false; @@ -283,4 +297,10 @@ } + { + psThreadTask *task = psThreadTaskAlloc("PPSTACK_PERCENT_COMBINE", 3); + task->function = &ppStackReadoutPercentThread; + psThreadTaskAdd(task); + psFree(task); + } return; } Index: /trunk/ppStack/src/ppStackThread.h =================================================================== --- /trunk/ppStack/src/ppStackThread.h (revision 42093) +++ /trunk/ppStack/src/ppStackThread.h (revision 42094) @@ -11,4 +11,6 @@ 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 + int entry; // thread number for tracking progress + int status; // thread number for tracking progress } ppStackThread; @@ -47,3 +49,6 @@ +enum {PPSTACK_THREAD_NEW, PPSTACK_THREAD_READ, PPSTACK_THREAD_RUN, PPSTACK_THREAD_SUCCESS, PPSTACK_THREAD_FAILURE}; + + #endif