Index: /trunk/psModules/src/detrend/pmDark.c =================================================================== --- /trunk/psModules/src/detrend/pmDark.c (revision 18858) +++ /trunk/psModules/src/detrend/pmDark.c (revision 18859) @@ -98,48 +98,43 @@ } - -bool pmDarkCombine(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept, - int iter, float rej, psMaskType maskVal) -{ - PS_ASSERT_PTR_NON_NULL(output, false); - PS_ASSERT_ARRAY_NON_NULL(inputs, false); - PS_ASSERT_ARRAY_NON_NULL(ordinates, false); - PS_ASSERT_INT_NONNEGATIVE(iter, false); - PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); - - bool inRange = true; - - // Extract fitting orders - int numOrdinates = ordinates->n; // Number of ordinates - int numInputs = inputs->n; // Number of inputs +// this creates and saves: values, roMask, norm, poly, counts, sigma, and saves the on output->analysis +bool pmDarkCombinePrepare(pmCell *output, const psArray *inputs, psArray *ordinates, const char *normConcept) +{ + psArray *values = psArrayAlloc(inputs->n); + psVector *roMask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for bad readouts + psVector *norm = normConcept ? psVectorAlloc(inputs->n, PS_TYPE_F32) : NULL; // Normalizations for each + psVector *orders = psVectorAlloc(ordinates->n, PS_TYPE_U8); // Orders for each concept + + psVectorInit(roMask, 0); + + bool inRange = false; int numBadInputs = 0; // Number of bad inputs - psArray *values = psArrayAlloc(numInputs); - psVector *roMask = psVectorAlloc(numInputs, PS_TYPE_U8); // Mask for bad readouts - psVectorInit(roMask, 0); - psVector *norm = normConcept ? psVectorAlloc(numInputs, PS_TYPE_F32) : NULL; // Normalisations for each - for (int i = 0; i < numInputs; i++) { - values->data[i] = psVectorAlloc(numOrdinates, PS_TYPE_F32); - if (norm) { - pmReadout *ro = inputs->data[i]; // Readout of interest - float normValue; // Normalisation value - if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, ro)) { - psWarning("Unable to find value of %s for readout %d", normConcept, i); - roMask->data.U8[i] = 0xff; - norm->data.F32[i] = NAN; - numBadInputs++; - continue; - } - if (normValue == 0.0) { - psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); - roMask->data.U8[i] = 0xff; - norm->data.F32[i] = NAN; - numBadInputs++; - continue; - } - norm->data.F32[i] = 1.0 / normValue; - } - } - psVector *orders = psVectorAlloc(numOrdinates, PS_TYPE_U8); // Orders for each concept - for (int i = 0; i < numOrdinates; i++) { + + // build the 'norm' vector and the 'values' vectors, count the number of bad inputs + for (int i = 0; i < inputs->n; i++) { + values->data[i] = psVectorAlloc(ordinates->n, PS_TYPE_F32); + if (!norm) continue; + + pmReadout *readout = inputs->data[i]; // Readout of interest + float normValue; // Normalisation value + if (!ordinateLookup(&normValue, &inRange, normConcept, false, NAN, NAN, readout)) { + psWarning("Unable to find value of %s for readout %d", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + if (normValue == 0.0) { + psWarning("Normalisation value (%s) for readout %d is zero", normConcept, i); + roMask->data.U8[i] = 0xff; + norm->data.F32[i] = NAN; + numBadInputs++; + continue; + } + norm->data.F32[i] = 1.0 / normValue; + } + + // build the 'orders' vector and set the array of 'values' + for (int i = 0; i < ordinates->n; i++) { pmDarkOrdinate *ord = ordinates->data[i]; // Ordinate information if (ord->order <= 0) { @@ -148,17 +143,10 @@ psFree(roMask); psFree(orders); + psFree(norm); return false; } orders->data.U8[i] = ord->order; - // Mask the readout and move on - #define MASK_READOUT_VALUE { \ - roMask->data.U8[j] = 0xff; \ - val->data.F32[i] = NAN; \ - numBadInputs++; \ - continue; \ - } - - for (int j = 0; j < numInputs; j++) { + for (int j = 0; j < inputs->n; j++) { psVector *val = values->data[j]; // Value vector for readout if (roMask->data.U8[j]) { @@ -167,7 +155,7 @@ } - pmReadout *ro = inputs->data[j]; // Readout of interest + pmReadout *readout = inputs->data[j]; // Readout of interest float value = NAN; // Value of ordinate - if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, ro)) { + if (!ordinateLookup(&value, &inRange, ord->name, ord->scale, ord->min, ord->max, readout)) { roMask->data.U8[j] = 0xff; val->data.F32[i] = NAN; @@ -186,19 +174,19 @@ if (psTraceGetLevel("psModules.detrend") > 9) { - for (int i = 0; i < numInputs; i++) { + for (int i = 0; i < inputs->n; i++) { psVector *val = values->data[i]; (void) val; // avoid unused variable message when tracing is compiled out - for (int j = 0; j < numOrdinates; j++) { + for (int j = 0; j < ordinates->n; j++) { psTrace("psModules.detrend", 9, "Image %d, ordinate %d: %f\n", i, j, val->data.F32[j]); } } } - psPolynomialMD *poly = psPolynomialMDAlloc(orders); // Polynomial for fitting psFree(orders); + int numTerms = poly->coeff->n; // Number of terms in polynomial - if (numTerms > numInputs - numBadInputs) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%d) to fit polynomial terms (%d).", - numInputs - numBadInputs, numTerms); + if (numTerms > inputs->n - numBadInputs) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Insufficient inputs (%ld) to fit polynomial terms (%d).", + inputs->n - numBadInputs, numTerms); psFree(values); psFree(roMask); @@ -207,5 +195,95 @@ } - // Set up output + // determine the output image size based on the input images + int row0, col0, numCols, numRows; + if (!pmReadoutStackSetOutputSize(&col0, &row0, &numCols, &numRows, inputs)) { + psError(PS_ERR_UNKNOWN, false, "problem setting output readout size."); + return false; + } + + // the output is potentially a cube (depending on the dimensionality of the fit) + if (output->readouts->n != numTerms) { + output->readouts = psArrayRealloc(output->readouts, numTerms); + } + + // generate the required output images based on the specified sizes + for (int i = 0; i < numTerms; i++) { + pmReadout *readout = output->readouts->data[i]; // Readout to update + if (!readout) { + readout = output->readouts->data[i] = pmReadoutAlloc(output); + } + + pmReadoutStackDefineOutput(readout, col0, row0, numCols, numRows, false, false, 0); + psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", col0, row0); + } + + // these calls allocate and save the requested images on the output analysis metadata + psImage *counts = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT, numCols, numRows, PS_TYPE_U16, 0); + if (!counts) { + return false; + } + psImage *sigma = pmReadoutSetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA, numCols, numRows, PS_TYPE_F32, NAN); + if (!sigma) { + return false; + } + + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); + psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, PS_META_REPLACE, "Dark normalisation", normConcept); + + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.VALUES", PS_DATA_ARRAY | PS_META_REPLACE, "Dark values", values); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.RO.MASK", PS_DATA_VECTOR | PS_META_REPLACE, "Dark Readout Mask", roMask); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.NORM", PS_DATA_VECTOR | PS_META_REPLACE, "Dark norm", norm); + psMetadataAddPtr(output->analysis, PS_LIST_TAIL, "DARK.POLY", PS_DATA_UNKNOWN | PS_META_REPLACE, "Dark poly", poly); + + for (int i = 0; i < numTerms; i++) { + pmReadout *readout = output->readouts->data[i]; // Readout to update + readout->data_exists = true; + } + output->data_exists = true; + output->parent->data_exists = true; + + psFree(norm); + psFree(roMask); + psFree(poly); + psFree(values); + + return true; +} + +// do the combine work for this portion of the output (range is set by input data) +bool pmDarkCombine(pmCell *output, const psArray *inputs, int iter, float rej, psMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(output, false); + PS_ASSERT_PTR_NON_NULL(output->readouts, false); + PS_ASSERT_INT_NONNEGATIVE(output->readouts->n, false); + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + + pthread_t id = pthread_self(); + char name[64]; + sprintf (name, "%x", (unsigned int) id); + psTimerStart (name); + + bool mdok = false; + + // retrieve the required parameter vectors + psArray *values = psMetadataLookupPtr(&mdok, output->analysis, "DARK.VALUES"); psAssert (values, "values not supplied"); + psVector *roMask = psMetadataLookupPtr(&mdok, output->analysis, "DARK.RO.MASK"); psAssert (roMask, "roMask not supplied"); + psPolynomialMD *poly = psMetadataLookupPtr(&mdok, output->analysis, "DARK.POLY"); psAssert (poly, "orders not supplied"); + + // retrieve the norm vector, if supplied + psVector *norm = psMetadataLookupPtr(&mdok, output->analysis, "DARK.NORM"); + + // retrieve the 'counts' and 'sigma' images + psImage *counts = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT); + if (!counts) { + return false; + } + psImage *sigma = pmReadoutGetAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA); + if (!sigma) { + return false; + } + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine int xSize, ySize; // Size of the output image @@ -213,55 +291,24 @@ inputs)) { psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); - psFree(values); - psFree(roMask); - psFree(norm); - return false; - } - if (output->readouts->n != numTerms) { - output->readouts = psArrayRealloc(output->readouts, numTerms); - } - int outRow0 = 0, outCol0 = 0; // Output row0 and col0 - for (int i = 0; i < numTerms; i++) { - pmReadout *ro = output->readouts->data[i]; // Readout to update - if (!ro) { - ro = output->readouts->data[i] = pmReadoutAlloc(output); - } - - pmReadoutUpdateSize(ro, minInputCols, minInputRows, xSize, ySize, false, false, 0); - if (i == 0) { - outRow0 = ro->row0; - outCol0 = ro->col0; - } else { - assert(ro->row0 == outRow0); - assert(ro->col0 == outCol0); - } - } - - psImage *counts = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_COUNT, - xSize, ySize, PS_TYPE_U16, 0); - if (!counts) { - return false; - } - psImage *sigma = pmReadoutAnalysisImage(output->readouts->data[0], PM_READOUT_STACK_ANALYSIS_SIGMA, - xSize, ySize, PS_TYPE_F32, NAN); - if (!sigma) { - psFree(counts); - return false; - } + return false; + } + + pmReadout *outReadout = output->readouts->data[0]; // Iterate over pixels, fitting polynomial - psVector *pixels = psVectorAlloc(numInputs, PS_TYPE_F32); // Stack of pixels - psVector *mask = psVectorAlloc(numInputs, PS_TYPE_MASK); // Mask for stack + psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels + psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_MASK); // Mask for stack for (int i = minInputRows; i < maxInputRows; i++) { - int yOut = i - outRow0; // y position on output readout -#ifdef SHOW_BUSY + int yOut = i - outReadout->row0; // y position on output readout + +# ifdef SHOW_BUSY if (psTraceGetLevel("psModules.detrend") > 9) { printf("Processing row %d\r", i); fflush(stdout); } -#endif +# endif for (int j = minInputCols; j < maxInputCols; j++) { - int xOut = j - outCol0; // x position on output readout + int xOut = j - outReadout->col0; // x position on output readout psVectorInit(mask, 0); @@ -303,5 +350,5 @@ psVectorInit(poly->coeff, NAN); } - for (int k = 0; k < numTerms; k++) { + for (int k = 0; k < poly->coeff->n; k++) { pmReadout *ro = output->readouts->data[k]; // Readout of interest ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; @@ -312,27 +359,10 @@ } - psFree(norm); - psFree(roMask); - psFree(poly); psFree(pixels); psFree(mask); - psFree(values); - - psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_ORDINATES, - PS_DATA_ARRAY | PS_META_REPLACE, "Dark ordinates", ordinates); - psMetadataAddStr(output->analysis, PS_LIST_TAIL, PM_DARK_ANALYSIS_NORM, - PS_META_REPLACE, "Dark normalisation", normConcept); - - for (int i = 0; i < numTerms; i++) { - pmReadout *ro = output->readouts->data[i]; // Readout to update - ro->data_exists = true; - } - output->data_exists = true; - output->parent->data_exists = true; - + + fprintf (stderr, "done with combine %x : %f sec\n", (unsigned int) id, psTimerMark (name)); return true; } - - bool pmDarkApply(pmReadout *readout, const pmCell *dark, psMaskType bad) Index: /trunk/psModules/src/detrend/pmDark.h =================================================================== --- /trunk/psModules/src/detrend/pmDark.h (revision 18858) +++ /trunk/psModules/src/detrend/pmDark.h (revision 18859) @@ -25,9 +25,14 @@ -// Combine darks +// Combine darks -- preparation step +bool pmDarkCombinePrepare(pmCell *output, // Output cell; readouts will be attached + const psArray *inputs, // Input readouts for combination + psArray *ordinates, // Ordinates for fitting + const char *normConcept // Concept name to use to divide input pixel values + ); + +// Combine darks -- do the actual work bool pmDarkCombine(pmCell *output, // Output cell; readouts will be attached const psArray *inputs, // Input readouts for combination - psArray *ordinates, // Ordinates for fitting - const char *normConcept, // Concept name to use to divide input pixel values int iter, // Number of rejection iterations float rej, // Rejection threshold (standard deviations)