Index: /trunk/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9335) +++ /trunk/psModules/src/detrend/pmShutterCorrection.c (revision 9336) @@ -60,4 +60,5 @@ #include +#include "pmFPA.h" #include "pmShutterCorrection.h" #include "psVectorBracket.h" @@ -193,5 +194,9 @@ const psVector *counts, const psVector *cntError, - float offref + const psVector *mask, + float offref, + int nIter, + float rej, + psMaskType maskVal ) { @@ -232,6 +237,11 @@ line->coeff[1][1] = 0; - if (!psVectorFitPolynomial2D(line, NULL, 0, counts, cntError, x, y)) { + psStats *stats = psStatsAlloc(0); + stats->clipSigma = rej; + stats->clipIter = nIter; + + if (!psVectorClipFitPolynomial2D(line, stats, mask, maskVal, counts, cntError, x, y)) { psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to fit shutter correction.\n"); + psFree(stats); psFree(x); psFree(y); @@ -239,4 +249,5 @@ return NULL; } + psFree(stats); pmShutterCorrection *corr = pmShutterCorrectionAlloc(); @@ -355,4 +366,6 @@ psStatsOptions meanStat, // Statistic to use for mean psStatsOptions stdevStat, // Statistic to use for stdev + int nIter, // Number of iterations per pixel + float rej, // Rejection limit psMaskType maskVal // Mask value ) @@ -382,4 +395,8 @@ if (!image) { continue; + } + if (image->type.type != PS_TYPE_F32) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for image: %x\n", image->type.type); + goto MEASURE_ERROR; } if (numRows == 0 || numCols == 0) { @@ -401,4 +418,19 @@ goto MEASURE_ERROR; } + if (masks) { + psImage *mask = masks->data[i]; + if (mask) { + if (mask->type.type != PS_TYPE_U8) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Bad type for mask: %x\n", mask->type.type); + goto MEASURE_ERROR; + } + if (mask->numRows != numRows || mask->numCols != numCols) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Mask sizes don't match: %dx%d vs %dx%d\n", mask->numCols, mask->numRows, + numCols, numRows); + goto MEASURE_ERROR; + } + } + } // Measure statistics @@ -407,4 +439,5 @@ } refs->data.F32[i] = psStatsGetValue(stats, meanStat); + psTrace("psModules.detrend", 3, "Reference value for image %ld = %f\n", i, refs->data.F32[i]); if (refs->data.F32[i] <= 0.0) { psError(PS_ERR_UNKNOWN, true, "Measured non-positive reference value.\n"); @@ -424,4 +457,6 @@ samplesMean->data.F32[j][i] = psStatsGetValue(stats, meanStat) * refs->data.F32[i]; samplesStdev->data.F32[j][i] = psStatsGetValue(stats, stdevStat) * refs->data.F32[i]; + psTrace("psModules.detrend", 5, "Image %ld, sample %d: %f +/- %f\n", i, j, + samplesMean->data.F32[j][i], samplesStdev->data.F32[j][i]); } } @@ -437,10 +472,11 @@ errors = psImageRow(errors, samplesStdev, i); pmShutterCorrection *guess = pmShutterCorrectionGuess(exptimes, counts); // Guess at correction - pmShutterCorrection *params = pmShutterCorrectionFullFit(exptimes, counts, errors, guess); // Correc'n - if (isfinite(params->offref)) { - meanRef += params->offref; + pmShutterCorrection *corr = pmShutterCorrectionFullFit(exptimes, counts, errors, guess); // Correct'n + psTrace("psModules.detrend", 5, "Sample reference value: %f\n", corr->offref); + if (isfinite(corr->offref)) { + meanRef += corr->offref; numGood++; } - psFree(params); + psFree(corr); psFree(guess); } @@ -455,16 +491,26 @@ } meanRef /= (float)numGood; + psTrace("psModules.detrend", 3, "Mean reference value: %f\n", meanRef); psImage *shutter = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Shutter correction image //psImage *pattern = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Illumination pattern + psVector *mask = psVectorAlloc(num, PS_TYPE_U8); // Mask for each image + mask->n = num; + psVectorInit(mask, 0); + psTrace("psModules.detrend", 2, "Performing linear fit on individual pixels...\n"); for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { for (int i = 0; i < num; i++) { psImage *image = images->data[i]; // Image of interest + psImage *maskImage; // Mask image + if (masks && (maskImage = masks->data[i])) { + mask->data.U8[i] = maskImage->data.U8[y][x]; + } counts->data.F32[i] = image->data.F32[y][x] * refs->data.F32[i]; errors->data.F32[i] = sqrtf(image->data.F32[y][x]) * refs->data.F32[i]; } - pmShutterCorrection *corr = pmShutterCorrectionLinFit(exptimes, counts, errors, meanRef); + pmShutterCorrection *corr = pmShutterCorrectionLinFit(exptimes, counts, errors, mask, meanRef, + nIter, rej, maskVal); shutter->data.F32[y][x] = corr->offset; //pattern->data.F32[y][x] = corr->scale; @@ -488,2 +534,50 @@ return NULL; } + + +bool pmShutterCorrectionApply(pmReadout *readout, // Readout to which to apply shutter correction + const pmReadout *shutter // Shutter correction readout + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_PTR_NON_NULL(shutter, false); + PS_ASSERT_IMAGE_NON_NULL(readout->image, false); + PS_ASSERT_IMAGE_NON_NULL(shutter->image, false); + PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(shutter->image, PS_TYPE_F32, false); + + psRegion region = psRegionSet(readout->col0, readout->col0 + readout->image->numCols, + readout->row0, readout->row0 + readout->image->numRows); // Detector region + + pmCell *cell = readout->parent; // Parent cell + if (!cell) { + psError(PS_ERR_BAD_PARAMETER_NULL, true, + "Parent cell is NULL --- unable to determine exposure time.\n"); + return false; + } + float exptime = psMetadataLookupF32(NULL, cell->concepts, "CELL.EXPOSURE"); // Exposure time + if (!isfinite(exptime) || exptime < 0.0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Bad exposure time: %f.\n", exptime); + return false; + } + + psImage *shutterImage = psImageSubset(shutter->image, region); // Subimage with shutter + if (!shutterImage) { + psString regionString = psRegionToString(region); + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Size mismatch: %s vs %dx%d\n", + regionString, shutter->image->numCols, shutter->image->numRows); + psFree(regionString); + psFree(shutterImage); + return false; + } + psImage *image = readout->image; // Image to correct + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + image->data.F32[y][x] *= exptime / (exptime + shutterImage->data.F32[y][x]); + } + } + + psFree(shutterImage); + + return true; +} Index: /trunk/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9335) +++ /trunk/psModules/src/detrend/pmShutterCorrection.h (revision 9336) @@ -48,6 +48,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-10-05 23:36:26 $ + * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-10-06 03:06:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -82,5 +82,9 @@ const psVector *counts, // Counts for each exposure const psVector *cntError, // Error in the counts, for each exp. - float offref // Reference time offset + const psVector *mask, // Mask for each exposure + float offref, // Reference time offset + int nIter, // Number of iterations + float rej, // Rejection threshold (sigma) + psMaskType maskVal // Mask value ); @@ -99,6 +103,13 @@ psStatsOptions meanStat, // Statistic to use for mean psStatsOptions stdevStat, // Statistic to use for stdev + int nIter, // Number of iterations + float rej, // Rejection threshold (sigma) psMaskType maskVal // Mask value ); +// Apply a shutter correction +bool pmShutterCorrectionApply(pmReadout *readout, // Readout to which to apply shutter correction + const pmReadout *shutter // Shutter correction readout + ); + #endif