Index: branches/tap_branches/ppStack/src/ppStackReadout.c =================================================================== --- branches/tap_branches/ppStack/src/ppStackReadout.c (revision 25900) +++ branches/tap_branches/ppStack/src/ppStackReadout.c (revision 27838) @@ -23,13 +23,12 @@ 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 - psArray *inspect = ppStackReadoutInitial(config, outRO, thread->readouts, mask, - weightings, addVariance); - - job->results = inspect; + job->results = ppStackReadoutInitial(config, outRO, thread->readouts, mask, + weightings, exposures, addVariance); thread->busy = false; - return inspect ? true : false; + return job->results ? true : false; } @@ -40,18 +39,23 @@ 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 + psArray *reject = args->data[1]; // Rejected pixels for each image + ppStackOptions *options = args->data[2]; // Options + pmConfig *config = args->data[3]; // Configuration + bool safety = PS_SCALAR_VALUE(args->data[4], U8); // Safety switch on? + bool normalise = PS_SCALAR_VALUE(args->data[5], U8); // Normalise images? + psVector *mask = options->inputMask; // Mask for inputs - psArray *rejected = options->rejected; // Rejected pixels psVector *weightings = options->weightings; // Weightings (1/noise^2) for each image + psVector *exposures = options->exposures; // Exposure times for each image psVector *addVariance = options->matchChi2; // Additional variance when rejecting - - bool status = ppStackReadoutFinal(config, outRO, thread->readouts, mask, rejected, - weightings, addVariance); // Status of operation + psVector *norm = normalise ? options->norm : NULL; // Normalisations to apply to images + + bool status = ppStackReadoutFinal(config, options->outRO, options->expRO, thread->readouts, mask, reject, + weightings, exposures, addVariance, safety, norm); // Status of operation thread->busy = false; + psAssert(status, "Stacking failed."); + return status; } @@ -63,24 +67,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; @@ -89,5 +104,6 @@ psArray *ppStackReadoutInitial(const pmConfig *config, pmReadout *outRO, const psArray *readouts, - const psVector *mask, const psVector *weightings, const psVector *addVariance) + const psVector *mask, const psVector *weightings, const psVector *exposures, + const psVector *addVariance) { assert(config); @@ -104,5 +120,5 @@ bool mdok; // Status of MD lookup - int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + float iter = psMetadataLookupF32(NULL, recipe, "COMBINE.ITER"); // Rejection iterations float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error @@ -116,4 +132,6 @@ psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels @@ -129,15 +147,4 @@ } -#if 0 - // This doesn't seem to work, so getting the weightings directly from a vector - float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight - if (!mdok || !isfinite(weighting)) { - psWarning("No weighting supplied for image %d --- set to unity.", i); - weighting = 1.0; - } else { - psLogMsg("ppStack", PS_LOG_INFO, "Weighting for image %d is %f", i, weighting); - } -#endif - // Ensure there is a mask, or pmStackCombine will complain if (!ro->mask) { @@ -146,9 +153,10 @@ } - stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], addVariance->data.F32[i]); - } - - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, kernelSize, iter, - combineRej, combineSys, combineDiscard, true, useVariance, safe, false)) { + stack->data[i] = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], + addVariance->data.F32[i]); + } + + if (!pmStackCombine(outRO, NULL, stack, maskVal | maskBad, maskSuspect, maskBad, kernelSize, iter, + combineRej, combineSys, combineDiscard, useVariance, safe, false)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts with rejection."); psFree(stack); @@ -156,6 +164,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 @@ -168,4 +177,5 @@ } inspect->data[i] = psMemIncrRefCounter(data->inspect); + reject->data[i] = psMemIncrRefCounter(data->reject); } psFree(stack); @@ -173,15 +183,21 @@ sectionNum++; - return inspect; -} - - - -bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, const psArray *readouts, + psArray *results = psArrayAlloc(2); // Array of results + results->data[0] = inspect; + results->data[1] = reject; + + return results; +} + + + +bool ppStackReadoutFinal(const pmConfig *config, pmReadout *outRO, pmReadout *expRO, const psArray *readouts, const psVector *mask, const psArray *rejected, const psVector *weightings, - const psVector *addVariance) + const psVector *exposures, const psVector *addVariance, bool safety, + const psVector *norm) { assert(config); assert(outRO); + assert(expRO); assert(readouts); assert(!rejected || readouts->n == rejected->n); @@ -196,8 +212,4 @@ bool mdok; // Status of MD lookup - int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations - float combineRej = psMetadataLookupF32(NULL, recipe, "COMBINE.REJ"); // Combination threshold - float combineSys = psMetadataLookupF32(NULL, recipe, "COMBINE.SYS"); // Combination systematic error - float combineDiscard = psMetadataLookupF32(NULL, recipe, "COMBINE.DISCARD"); // Olympic discard fraction bool useVariance = psMetadataLookupBool(&mdok, recipe, "VARIANCE"); // Use variance for rejection? bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels @@ -205,4 +217,6 @@ psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch + psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits + psImageMaskType maskSuspect = pmConfigMaskGet(maskSuspectStr, config); // Suspect bits psString maskBadStr = psMetadataLookupStr(NULL, recipe, "MASK.BAD"); // Name of bits to mask for bad psImageMaskType maskBad = pmConfigMaskGet(maskBadStr, config); // Bits to mask for bad pixels @@ -211,34 +225,17 @@ psArray *stack = psArrayAlloc(num); // Array for stacking - bool entire = true; // Combine entire image? - if (rejected) { - // We have rejection from a previous combination: combine without flagging pixels to inspect - entire = false; - safe = false; - iter = 0; - combineRej = NAN; - combineSys = NAN; - } + // We have rejection from a previous combination: combine without flagging pixels to inspect + safe &= safety; + int iter = 0; + float combineRej = NAN; + float combineSys = NAN; + float combineDiscard = NAN; 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 - entire = true; - continue; - } else if (mask->data.U8[i]) { - // Image completely rejected before original combination - continue; - } - -#if 0 - // This doesn't seem to work, so getting the weightings directly from a vector - bool mdok; // Status of MD lookup - float weighting = psMetadataLookupF32(&mdok, ro->analysis, "PPSTACK.WEIGHTING"); // Relative weight - if (!mdok || !isfinite(weighting)) { - psWarning("No WEIGHTING supplied for image %d --- set to unity.", i); - weighting = 1.0; - } -#endif + if (mask->data.U8[i]) { + // Image completely rejected + continue; + } // Ensure there is a mask, or pmStackCombine will complain @@ -248,13 +245,18 @@ } - pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], + pmStackData *data = pmStackDataAlloc(ro, weightings->data.F32[i], exposures->data.F32[i], addVariance ? addVariance->data.F32[i] : NAN); data->reject = rejected ? psMemIncrRefCounter(rejected->data[i]) : NULL; stack->data[i] = data; - } - - if (!pmStackCombine(outRO, stack, maskVal | maskBad, maskBad, 0, - iter, combineRej, combineSys, combineDiscard, - entire, useVariance, safe, !rejected)) { + + if (norm) { + float normalise = powf(10.0, -0.4 * norm->data.F32[i]); // Normalisation + psBinaryOp(ro->image, ro->image, "*", psScalarAlloc(normalise, PS_TYPE_F32)); + psBinaryOp(ro->variance, ro->variance, "*", psScalarAlloc(PS_SQR(normalise), PS_TYPE_F32)); + } + } + + if (!pmStackCombine(outRO, expRO, stack, maskVal | maskBad, maskSuspect, maskBad, 0, iter, combineRej, + combineSys, combineDiscard, useVariance, safe, rejected)) { psError(PS_ERR_UNKNOWN, false, "Unable to combine input readouts."); psFree(stack); @@ -264,9 +266,14 @@ pmCell *outCell = outRO->parent; // Output cell pmChip *outChip = outCell->parent; // Output chip - outRO->data_exists = true; outCell->data_exists = true; outChip->data_exists = true; + pmCell *expCell = expRO->parent; // Exposure cell + pmChip *expChip = expCell->parent; // Exposure chip + expRO->data_exists = true; + expCell->data_exists = true; + expChip->data_exists = true; + psFree(stack);