Index: branches/pap/ppStack/src/ppStack.h =================================================================== --- branches/pap/ppStack/src/ppStack.h (revision 25959) +++ branches/pap/ppStack/src/ppStack.h (revision 25964) @@ -59,5 +59,5 @@ // Perform stacking on a readout // -// Returns an array of pixels to inspect for each input image +// Returns two arrays: pixels to inspect for each input image, and pixels to reject for each input image. psArray *ppStackReadoutInitial(const pmConfig *config, // Configuration pmReadout *outRO, // Output readout @@ -84,5 +84,4 @@ const psVector *weightings, // Weighting factors for each image const psVector *addVariance, // Additional variance for rejection - bool full, // Combine full image? bool safety, // Enable safety switch? const psVector *norm // Normalisations to apply Index: branches/pap/ppStack/src/ppStackCombineFinal.c =================================================================== --- branches/pap/ppStack/src/ppStackCombineFinal.c (revision 25959) +++ branches/pap/ppStack/src/ppStackCombineFinal.c (revision 25964) @@ -13,5 +13,5 @@ bool ppStackCombineFinal(pmReadout *target, ppStackThreadData *stack, psArray *covariances, - ppStackOptions *options, pmConfig *config, bool full, bool safe, bool normalise) + ppStackOptions *options, pmConfig *config, bool safe, bool normalise) { psAssert(stack, "Require stack"); @@ -43,5 +43,4 @@ psArrayAdd(job->args, 1, options); psArrayAdd(job->args, 1, config); - PS_ARRAY_ADD_SCALAR(job->args, full, PS_TYPE_U8); PS_ARRAY_ADD_SCALAR(job->args, safe, PS_TYPE_U8); PS_ARRAY_ADD_SCALAR(job->args, normalise, PS_TYPE_U8); Index: branches/pap/ppStack/src/ppStackCombineInitial.c =================================================================== --- branches/pap/ppStack/src/ppStackCombineInitial.c (revision 25959) +++ branches/pap/ppStack/src/ppStackCombineInitial.c (revision 25964) @@ -62,4 +62,5 @@ // 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]) { @@ -67,4 +68,5 @@ } options->inspect->data[i] = psArrayAllocEmpty(numChunk); + options->rejected->data[i] = psArrayAllocEmpty(numChunk); } psThreadJob *job; // Completed job @@ -73,9 +75,13 @@ "Job has incorrect type: %s", job->type); psArray *results = job->results; // Results of job + psAssert(results->n == 2, "Results array has wrong size!"); + psArray *inspect = results->data[0]; // Pixels to inspect + psArray *reject = results->data[1]; // Pixels to reject 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]); + options->inspect->data[i] = psArrayAdd(options->inspect->data[i], 1, inspect->data[i]); + options->rejected->data[i] = psArrayAdd(options->rejected->data[i], 1, reject->data[i]); } psFree(job); Index: branches/pap/ppStack/src/ppStackLoop.c =================================================================== --- branches/pap/ppStack/src/ppStackLoop.c (revision 25959) +++ branches/pap/ppStack/src/ppStackLoop.c (revision 25964) @@ -97,6 +97,5 @@ // Final combination psTrace("ppStack", 2, "Final stack of convolved images....\n"); - if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, - true, false, false)) { + if (!ppStackCombineFinal(options->outRO, stack, options->convCovars, options, config, false, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform final combination."); psFree(stack); @@ -106,5 +105,4 @@ psLogMsg("ppStack", PS_LOG_INFO, "Stage 5: Final Stack: %f sec", psTimerClear("PPSTACK_STEPS")); ppStackMemDump("final"); - // Clean up @@ -121,5 +119,4 @@ psFree(stack); - #if 1 // Unconvolved stack --- it's cheap to calculate, compared to everything else! @@ -133,6 +130,5 @@ } psTrace("ppStack", 2, "Stack of unconvolved images....\n"); - if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, - true, true, true)) { + if (!ppStackCombineFinal(options->unconvRO, stack, options->origCovars, options, config, true, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to perform unconvolved combination."); psFree(stack); @@ -158,5 +154,4 @@ ppStackMemDump("photometry"); - // Finish up psTrace("ppStack", 1, "Finishing up....\n"); @@ -170,4 +165,5 @@ psFree(options); + return true; } Index: branches/pap/ppStack/src/ppStackLoop.h =================================================================== --- branches/pap/ppStack/src/ppStackLoop.h (revision 25959) +++ branches/pap/ppStack/src/ppStackLoop.h (revision 25964) @@ -61,5 +61,4 @@ ppStackOptions *options, // Options pmConfig *config, // Configuration - bool full, // Combine full image? bool safe, // Allow safe combination? bool norm // Normalise images? Index: branches/pap/ppStack/src/ppStackReadout.c =================================================================== --- branches/pap/ppStack/src/ppStackReadout.c (revision 25959) +++ branches/pap/ppStack/src/ppStackReadout.c (revision 25964) @@ -25,11 +25,9 @@ psVector *addVariance = options->matchChi2; // Additional variance when rejecting - psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask, - weightings, addVariance); - - job->results = inspect; + job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask, + weightings, addVariance); thread->busy = false; - return inspect ? true : false; + return job->results ? true : false; } @@ -43,7 +41,6 @@ ppStackOptions *options = args->data[2]; // Options pmConfig *config = args->data[3]; // Configuration - bool full = PS_SCALAR_VALUE(args->data[4], U8); // Combine full image? - bool safety = PS_SCALAR_VALUE(args->data[5], U8); // Safety switch on? - bool normalise = PS_SCALAR_VALUE(args->data[6], U8); // Normalise images? + bool safety = PS_SCALAR_VALUE(args->data[4], U8); // Safety switch on? + bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images? psVector *mask = options->inputMask; // Mask for inputs @@ -54,5 +51,5 @@ bool status = ppStackReadoutFinal(config, target, thread->readouts, mask, rejected, - weightings, addVariance, full, safety, norm); // Status of operation + weightings, addVariance, safety, norm); // Status of operation thread->busy = false; @@ -67,24 +64,35 @@ psArray *args = job->args; // Input arguments - psArray *inspect = args->data[0]; // Array of pixel arrays - int index = PS_SCALAR_VALUE(args->data[1], S32); // Index of interest - - psArray *inputs = inspect->data[index]; // Array of interest - psPixels *output = NULL; // Output pixel list - for (int i = 0; i < inputs->n; i++) { - psPixels *input = inputs->data[i]; // Input pixel list - if (!input || input->n == 0) { - continue; - } - output = psPixelsConcatenate(output, input); - } - - if (!output) { - // If there are no pixels to inspect, then just fake it - output = psPixelsAllocEmpty(0); - } - - psFree(inputs); - inspect->data[index] = output; + psArray *inspects = args->data[0]; // Array of pixel arrays + psArray *rejects = args->data[1]; // Array of pixel arrays + int index = PS_SCALAR_VALUE(args->data[2], S32); // Index of interest + + psArray *inInspects = inspects->data[index]; // Array of interest + psArray *inRejects = rejects->data[index]; // Array of interest + psAssert(inInspects->n == inRejects->n, "Size should be the same"); + psPixels *outInspect = NULL, *outReject = NULL; // Output pixel lists + for (int i = 0; i < inInspects->n; i++) { + psPixels *inInspect = inInspects->data[i]; // Input pixel list + if (inInspect && inInspect->n > 0) { + outInspect = psPixelsConcatenate(outInspect, inInspect); + } + psPixels *inReject = inRejects->data[i]; // Input pixel list + if (inReject && inReject->n > 0) { + outReject = psPixelsConcatenate(outReject, inReject); + } + } + + // If there are no pixels to inspect, then just fake it + if (!outInspect) { + outInspect = psPixelsAllocEmpty(0); + } + if (!outReject) { + outReject = psPixelsAllocEmpty(0); + } + + psFree(inspects->data[index]); + inspects->data[index] = outInspect; + psFree(rejects->data[index]); + rejects->data[index] = outReject; return true; @@ -154,5 +162,5 @@ if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, - combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) { + combineRej, combineSys, combineDiscard, useVariance, safe, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); psFree(stack); @@ -160,6 +168,7 @@ } - // Save list of pixels to inspect + // Save lists of pixels psArray *inspect = psArrayAlloc(num); // List of pixels to inspect + psArray *reject = psArrayAlloc(num); // List of pixels rejected for (int i = 0; i < num; i++) { pmStackData *data = stack->data[i]; // Data for this image @@ -172,4 +181,5 @@ } inspect->data[i] = psMemIncrRefCounter(data->inspect); + reject->data[i] = psMemIncrRefCounter(data->reject); } psFree(stack); @@ -177,5 +187,9 @@ sectionNum++; - return inspect; + psArray *results = psArrayAlloc(2); // Array of results + results->data[0] = inspect; + results->data[1] = reject; + + return results; } @@ -184,5 +198,5 @@ bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, - const psVector *addVariance, bool full, bool safety, const psVector *norm) + const psVector *addVariance, bool safety, const psVector *norm) { assert(config); @@ -225,10 +239,6 @@ for (int i = 0; i < num; i++) { pmReadout *ro = readouts->data[i]; - if (mask->data.U8[i] & (PPSTACK_MASK_REJECT | PPSTACK_MASK_BAD)) { - // Image completely rejected since previous combination - full = true; - continue; - } else if (mask->data.U8[i]) { - // Image completely rejected before original combination + if (mask->data.U8[i]) { + // Image completely rejected continue; } @@ -262,7 +272,6 @@ } - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, - iter, combineRej, combineSys, combineDiscard, - full, useVariance, safe, !rejected)) { + if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, iter, combineRej, + combineSys, combineDiscard, useVariance, safe, !rejected)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); psFree(stack); Index: branches/pap/ppStack/src/ppStackReject.c =================================================================== --- branches/pap/ppStack/src/ppStackReject.c (revision 25959) +++ branches/pap/ppStack/src/ppStackReject.c (revision 25964) @@ -23,5 +23,4 @@ int num = options->num; // Number of inputs - options->rejected = psArrayAlloc(num); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -53,4 +52,5 @@ psThreadJob *job = psThreadJobAlloc("PPSTACK_INSPECT"); // Job to start psArrayAdd(job->args, 1, options->inspect); + psArrayAdd(job->args, 1, options->rejected); PS_ARRAY_ADD_SCALAR(job->args, i, PS_TYPE_S32); if (!psThreadJobAddPending(job)) { @@ -92,15 +92,11 @@ #endif - psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, - threshold, poorFrac, stride, options->regions->data[i], - options->kernels->data[i]); // Rejected pixels - #ifdef TESTING { - psImage *mask = psPixelsToMask(NULL, reject, + psImage *mask = psPixelsToMask(NULL, options->rejected->data[i], psRegionSet(0, options->numCols - 1, 0, options->numRows - 1), 0xff); // Mask image psString name = NULL; // Name of image - psStringAppend(&name, "reject_%03d.fits", i); + psStringAppend(&name, "pre_reject_%03d.fits", i); pmStackVisualPlotTestImage(mask, name); psFits *fits = psFitsOpen(name, "w"); @@ -111,4 +107,8 @@ } #endif + + psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, + threshold, poorFrac, stride, options->regions->data[i], + options->kernels->data[i]); // Rejected pixels psFree(options->inspect->data[i]); @@ -127,5 +127,4 @@ "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; @@ -134,6 +133,25 @@ } - // Images without a list of rejected pixels (the list may be empty) are rejected completely - 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); + } + +#ifdef TESTING + { + psImage *mask = psPixelsToMask(NULL, options->rejected->data[i], + 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 if (options->stats) { @@ -143,4 +161,6 @@ "Number of pixels rejected", reject ? reject->n : 0); } + + psFree(reject); psLogMsg("ppStack", PS_LOG_INFO, "Time to perform rejection on image %d: %f sec", i, psTimerClear("PPSTACK_REJECT")); Index: branches/pap/ppStack/src/ppStackThread.c =================================================================== --- branches/pap/ppStack/src/ppStackThread.c (revision 25959) +++ branches/pap/ppStack/src/ppStackThread.c (revision 25964) @@ -275,5 +275,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 2); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_INSPECT", 3); task->function = &ppStackInspect; psThreadTaskAdd(task); @@ -282,5 +282,5 @@ { - psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 7); + psThreadTask *task = psThreadTaskAlloc("PPSTACK_FINAL_COMBINE", 6); task->function = &ppStackReadoutFinalThread; psThreadTaskAdd(task);