Index: trunk/ppStack/src/ppStackReadout.c =================================================================== --- trunk/ppStack/src/ppStackReadout.c (revision 36886) +++ 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; +} + +