Index: /branches/pap/ppStack/src/ppStack.h =================================================================== --- /branches/pap/ppStack/src/ppStack.h (revision 25963) +++ /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 25963) +++ /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 25963) +++ /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 25963) +++ /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 25963) +++ /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 25963) +++ /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 25963) +++ /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 25963) +++ /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); Index: /branches/pap/psModules/src/imcombine/pmStack.c =================================================================== --- /branches/pap/psModules/src/imcombine/pmStack.c (revision 25963) +++ /branches/pap/psModules/src/imcombine/pmStack.c (revision 25964) @@ -290,4 +290,5 @@ // Mark a pixel for inspection +// Value in pixel doesn't seem to agree with the stack, so need to look closer static inline void combineInspect(const psArray *inputs, // Stack data int x, int y, // Pixel @@ -301,4 +302,21 @@ } data->inspect = psPixelsAdd(data->inspect, data->inspect->nalloc, x, y); + return; +} + +// Mark a pixel for rejection +// Cannot possibly inspect this pixel and confirm that it's good. +// e.g., Only a single input +static inline void combineReject(const psArray *inputs, // Stack data + int x, int y, // Pixel + int source // Source image index + ) +{ + pmStackData *data = inputs->data[source]; // Stack data of interest + if (!data) { + psWarning("Can't find input data for source %d", source); + return; + } + data->reject = psPixelsAdd(data->reject, data->reject->nalloc, x, y); return; } @@ -319,8 +337,8 @@ float rej, // Number of standard deviations at which to reject float sys, // Relative systematic error - float discard,// Fraction of values to discard (Olympic weighted mean) + float olympic,// Fraction of values to discard (Olympic weighted mean) bool useVariance, // Use variance for rejection when combining? bool safe, // Combine safely? - bool rejectInspect, // Reject values marked for inspection from combination? + bool rejection, // Reject values marked for inspection from combination? combineBuffer *buffer // Buffer for combination; to avoid multiple allocations ) @@ -419,13 +437,15 @@ } maskValue = 0; - } -#ifdef TESTING - else { + } else { + if (!rejection) { + combineReject(inputs, x, y, pixelSources->data.U16[0]); + } +#ifdef TESTING numRejected = 1; if (x == TEST_X && y == TEST_Y) { fprintf(stderr, "Single input to combine, safety on, pixel is bad.\n"); } +#endif } -#endif break; } @@ -459,5 +479,5 @@ if (PS_SQR(diff) > PS_SQR(rej) * sigma2) { // Not consistent: mark both for inspection - if (rejectInspect) { + if (rejection) { imageValue = NAN; varianceValue = NAN; @@ -550,5 +570,5 @@ } #endif - median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, discard, sort); + median = combinationWeightedOlympic(pixelData, pixelWeights, pixelMasks, olympic, sort); } @@ -563,5 +583,5 @@ #define MASK_PIXEL_FOR_INSPECTION() \ pixelMasks->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; \ - if (!rejectInspect) { \ + if (!rejection) { \ combineInspect(inputs, x, y, pixelSources->data.U16[j]); \ } \ @@ -605,5 +625,5 @@ } - if (rejectInspect && totalClipped > 0) { + if (rejection && totalClipped > 0) { // Get rid of the masked values // The alternative to this is to make combinationMeanVariance() accept a mask @@ -658,8 +678,8 @@ // Ensure the input array of pmStackData is valid, and get some details out of it static bool validateInputData(bool *haveVariances, // Do we have variance maps in the sky images? - bool *haveRejects, // Do we have lists of rejected pixels? int *num, // Number of inputs int *numCols, int *numRows, // Size of (sky) images - psArray *input // Input array of pmStackData to validate + const psArray *input, // Input array of pmStackData to validate + const pmReadout *output // Output readout ) { @@ -687,5 +707,5 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } - *haveRejects = (data->reject != NULL); + bool haveRejects = (data->reject != NULL); // Do we have rejected pixels? // Make sure the rest correspond with the first @@ -700,5 +720,5 @@ return false; } - if ((*haveRejects && !data->reject) || (data->reject && !*haveRejects)) { + if ((haveRejects && !data->reject) || (data->reject && !haveRejects)) { psError(PS_ERR_UNEXPECTED_NULL, true, "The rejected pixels are specified in some but not all inputs."); @@ -716,4 +736,13 @@ PS_ASSERT_IMAGE_TYPE(data->readout->variance, PS_TYPE_F32, false); } + } + + PM_ASSERT_READOUT_NON_NULL(output, false); + if (output->image) { + PS_ASSERT_IMAGE_NON_NULL(output->image, false); + PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(output->mask, false); + PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_IMAGE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(output->image, output->mask, false); } @@ -802,13 +831,11 @@ /// Stack input images bool pmStackCombine(pmReadout *combined, psArray *input, psImageMaskType maskVal, psImageMaskType bad, - int kernelSize, int numIter, float rej, float sys, float discard, - bool entire, bool useVariance, bool safe, bool rejectInspect) -{ - PS_ASSERT_PTR_NON_NULL(combined, false); + int kernelSize, int numIter, float rej, float sys, float olympic, + bool useVariance, bool safe, bool rejection) +{ bool haveVariances; // Do we have the variance maps? - bool haveRejects; // Do we have lists of rejected pixels? int num; // Number of inputs int numCols, numRows; // Size of (sky) images - if (!validateInputData(&haveVariances, &haveRejects, &num, &numCols, &numRows, input)) { + if (!validateInputData(&haveVariances, &num, &numCols, &numRows, input, combined)) { return false; } @@ -821,12 +848,4 @@ PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); } - if (haveRejects) { - // This is a subsequent combination, so expect that the image and mask already exist - PS_ASSERT_IMAGE_NON_NULL(combined->image, false); - PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false); - PS_ASSERT_IMAGE_NON_NULL(combined->mask, false); - PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_IMAGE_MASK, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); - } if (useVariance && !haveVariances) { psWarning("Unable to use variance in rejection if no variance maps supplied --- option turned off"); @@ -852,6 +871,12 @@ } #endif - if (!haveRejects && !data->inspect) { - data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + if (!rejection) { + // Ensure pixels can be put on the appropriate list + if (!data->inspect) { + data->inspect = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } + if (!data->reject) { + data->reject = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } } } @@ -882,91 +907,56 @@ combineBuffer *buffer = combineBufferAlloc(num); - if (haveRejects) { - psImage *combinedImage = combined->image; // Combined image - psImage *combinedMask = combined->mask; // Combined mask - psImage *combinedVariance = combined->variance; // Combined variance map - - psArray *pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, - minInputRows, maxInputRows); // Map of pixels to source - psPixels *pixels = NULL; // Total list of pixels, with no duplicates + // Pull the products out, allocate if necessary + psImage *combinedImage = combined->image; // Combined image + if (!combinedImage) { + combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedImage = combined->image; + } + psImage *combinedMask = combined->mask; // Combined mask + if (!combinedMask) { + combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + combinedMask = combined->mask; + } + + psImage *combinedVariance = combined->variance; // Combined variance map + if (haveVariances && !combinedVariance) { + combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedVariance = combined->variance; + } + + // Set up rejection list + psArray *pixelMap = NULL; // Map of pixels to source + if (rejection) { + pixelMap = pixelMapGenerate(input, minInputCols, maxInputCols, minInputRows, maxInputRows); + } + + // Combine each pixel + for (int y = minInputRows; y < maxInputRows; y++) { + for (int x = minInputCols; x < maxInputCols; x++) { + psVector *reject = NULL; // Images to reject for this pixel + if (rejection) { + reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, x, y); + } + combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, + addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, olympic, + useVariance, safe, rejection, buffer); + } + } + + psFree(pixelMap); + psFree(weights); + psFree(buffer); + +#ifndef PS_NO_TRACE + if (!rejection && psTraceGetLevel("psModules.imcombine") >= 5) { for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels - if (!data) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data || !data->inspect) { continue; } - pixels = psPixelsConcatenate(pixels, data->reject); - } - pixels = psPixelsDuplicates(pixels, pixels); - - if (entire) { - // Combine entire image - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - } else { - // Only combine previously rejected pixels - for (int i = 0; i < pixels->n; i++) { - // Pixel coordinates are in the frame of the original image - int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest - if (x < minInputCols || x >= maxInputCols || y < minInputRows || y >= maxInputRows) { - continue; - } - psVector *reject = pixelMapQuery(pixelMap, minInputCols, minInputRows, - x, y); // Reject these images - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, reject, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - psFree(pixels); - psFree(pixelMap); - } else { - // Pull the products out, allocate if necessary - psImage *combinedImage = combined->image; // Combined image - if (!combinedImage) { - combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedImage = combined->image; - } - psImage *combinedMask = combined->mask; // Combined mask - if (!combinedMask) { - combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); - combinedMask = combined->mask; - } - - psImage *combinedVariance = combined->variance; // Combined variance map - if (haveVariances && !combinedVariance) { - combined->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); - combinedVariance = combined->variance; - } - - for (int y = minInputRows; y < maxInputRows; y++) { - for (int x = minInputCols; x < maxInputCols; x++) { - combinePixels(combinedImage, combinedMask, combinedVariance, input, weights, - addVariance, NULL, x, y, maskVal, bad, numIter, rej, sys, discard, - useVariance, safe, rejectInspect, buffer); - } - } - -#ifndef PS_NO_TRACE - if (psTraceGetLevel("psModules.imcombine") >= 5) { - for (int i = 0; i < num; i++) { - pmStackData *data = input->data[i]; // Stack data for this input - if (!data || !data->inspect) { - continue; - } - psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); - } - } -#endif - } - - psFree(weights); - psFree(buffer); + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->inspect->n); + } + } +#endif return true; Index: /branches/pap/psModules/src/imcombine/pmStack.h =================================================================== --- /branches/pap/psModules/src/imcombine/pmStack.h (revision 25963) +++ /branches/pap/psModules/src/imcombine/pmStack.h (revision 25964) @@ -30,4 +30,5 @@ psPixels *reject; ///< Pixels to reject psPixels *inspect; ///< Pixels to inspect + psPixels *discard; ///< Pixels to discard float weight; ///< Relative weighting for image float addVariance; ///< Additional variance when rejecting @@ -49,9 +50,8 @@ float rej, ///< Rejection limit (standard deviations) float sys, ///< Relative systematic error - float discard, ///< Fraction of values to discard for Olympic weighted mean - bool entire, ///< Combine entire image even if rejection lists provided? + float olympic, ///< Fraction of values to discard for Olympic weighted mean bool useVariance, ///< Use variance values for rejection? bool safe, ///< Play safe with small numbers of input pixels (mask if N <= 2)? - bool rejectInspect ///< Reject pixels instead of marking them for inspection? + bool rejection ///< Reject pixels instead of marking them for inspection? );