Index: /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.c =================================================================== --- /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.c (revision 17154) +++ /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.c (revision 17155) @@ -3,8 +3,40 @@ #endif +#include +#include #include #include "pmReadoutStack.h" +psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image + const char *name, // Name of image in analysis metadata + int numCols, int numRows, // Expected size of image + psElemType type, // Expected type of image + double init // Initial value + ) +{ + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_STRING_NON_EMPTY(name, false); + + bool mdok; // Status of MD lookup + psImage *image = psMetadataLookupPtr(&mdok, readout->analysis, name); + if (!image) { + image = psImageAlloc(numCols, numRows, type); + psMetadataAddImage(readout->analysis, PS_LIST_TAIL, name, 0, "Analysis image from " __FILE__, image); + psImageInit(image, init); + return image; + } + if (image->numCols != numCols || image->numRows != numRows) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Analysis image %s has incorrect size (%dx%d vs %dx%d)", + name, image->numCols, image->numRows, numCols, numRows); + return NULL; + } + if (image->type.type != type) { + psError(PS_ERR_BAD_PARAMETER_TYPE, true, "Analysis image %s has incorrect type (%x vs %x)", + name, image->type.type, type); + return NULL; + } + return psMemIncrRefCounter(image); +} bool pmReadoutUpdateSize(pmReadout *readout, int minCols, int minRows, Index: /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.h =================================================================== --- /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.h (revision 17154) +++ /branches/pap_branch_080320/psModules/src/camera/pmReadoutStack.h (revision 17155) @@ -4,4 +4,7 @@ #include "pmHDU.h" #include "pmFPA.h" + +#define PM_READOUT_STACK_ANALYSIS_COUNT "STACK.COUNT" // Name for count image in analysis metadata +#define PM_READOUT_STACK_ANALYSIS_SIGMA "STACK.SIGMA" // Name for sigma image in analysis metadata /// Update an output readout (for a stack) with the correct col0,row0 and the image size @@ -21,4 +24,11 @@ ); +/// Return an image from analysis metadata, produced while stacking +psImage *pmReadoutAnalysisImage(pmReadout *readout, // Readout containing image + const char *name, // Name of image in analysis metadata + int numCols, int numRows, // Expected size of image + psElemType type, // Expected type of image + double init // Initial value + ); #endif Index: /branches/pap_branch_080320/psModules/src/detrend/pmDark.c =================================================================== --- /branches/pap_branch_080320/psModules/src/detrend/pmDark.c (revision 17154) +++ /branches/pap_branch_080320/psModules/src/detrend/pmDark.c (revision 17155) @@ -227,4 +227,16 @@ } + 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; + } + // Iterate over pixels, fitting polynomial psVector *pixels = psVectorAlloc(numInputs, PS_TYPE_F32); // Stack of pixels @@ -269,4 +281,6 @@ ro->image->data.F32[yOut][xOut] = poly->coeff->data.F64[k]; } + counts->data.U16[yOut][xOut] = poly->numFit; + sigma->data.F32[yOut][xOut] = poly->stdevFit; } } Index: /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.c (revision 17154) +++ /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.c (revision 17155) @@ -81,50 +81,45 @@ -psImage *pmMaskFlagSuspectPixels(psImage *out, const pmReadout *readout, float rej, - psMaskType maskVal) +bool pmMaskFlagSuspectPixels(pmReadout *output, const pmReadout *readout, float median, float stdev, + float rej, psMaskType maskVal) { - PS_ASSERT_PTR_NON_NULL(readout, NULL); - PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, NULL); - PS_ASSERT_IMAGE_NON_NULL(readout->image, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(readout->image, NULL); - PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, NULL); + PS_ASSERT_PTR_NON_NULL(readout, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + PS_ASSERT_IMAGE_NON_NULL(readout->image, false); + PS_ASSERT_IMAGE_NON_EMPTY(readout->image, false); + PS_ASSERT_IMAGE_TYPE(readout->image, PS_TYPE_F32, false); if (readout->mask) { - PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, NULL); - PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, NULL); - } - if (out) { - PS_ASSERT_IMAGE_NON_EMPTY(out, NULL); - PS_ASSERT_IMAGE_TYPE(out, PS_TYPE_S32, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, out, NULL); - } - - bool status; + PS_ASSERT_IMAGE_NON_EMPTY(readout->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, readout->mask, false); + PS_ASSERT_IMAGE_TYPE(readout->mask, PS_TYPE_MASK, false); + } + PS_ASSERT_PTR_NON_NULL(output, false); + + bool mdok; // Status of MD lookup + psImage *suspect = psMetadataLookupPtr(&mdok, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img + if (suspect) { + PS_ASSERT_IMAGE_NON_EMPTY(output->image, false); + PS_ASSERT_IMAGE_TYPE(output->image, PS_TYPE_S32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(readout->image, output->image, false); + psMemIncrRefCounter(suspect); + } else { + suspect = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32); + psImageInit(suspect, 0); + psMetadataAddImage(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_SUSPECT, 0, + "Suspect pixels", suspect); + psMetadataAddS32(output->analysis, PS_LIST_TAIL, PM_MASK_ANALYSIS_NUM, 0, + "Number of input images", 0); + } + + if (!isfinite(median) || !isfinite(stdev)) { + // If we get down here and the statistics are missing, then we should go and mask the entire image + psWarning("Missing statistics --- flagging entire image as suspect."); + return (psImage*)psBinaryOp(suspect, suspect, "+", psScalarAlloc(1.0, PS_TYPE_S32)); + } + psImage *image = readout->image; // Image of interest psImage *mask = readout->mask; // Corresponding mask - if (!out) { - out = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S32); - psImageInit(out, 0); - } - - bool whole = false; // Mask whole image? - float median = psMetadataLookupF32 (&status, readout->analysis, PM_MASK_ANALYSIS_MEAN); - if (!status || !isfinite(median)) { - whole = true; - } - float stdev = psMetadataLookupF32 (&status, readout->analysis, PM_MASK_ANALYSIS_STDEV); - if (!status || !isfinite(stdev)) { - whole = true; - } - - psTrace ("psModules.detrend", 3, "suspect: %f +/- %f\n", median, stdev); - - if (whole) { - // If we get down here and the statistics are missing, then we should go and mask the entire image - psWarning("Missing statistics --- flagging entire image as suspect."); - return (psImage*)psBinaryOp(out, out, "+", psScalarAlloc(1.0, PS_TYPE_S32)); - } for (int y = 0; y < image->numRows; y++) { @@ -132,20 +127,38 @@ if (fabs((image->data.F32[y][x] - median) / stdev) >= rej && (!mask || !(mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal))) { - out->data.S32[y][x]++; - } - } - } - - return out; -} - -psImage *pmMaskIdentifyBadPixels(const psImage *suspects, psMaskType maskVal, int nTotal, float thresh, pmMaskIdentifyMode mode) + suspect->data.S32[y][x]++; + } + } + } + psFree(suspect); // Drop reference + + psMetadataItem *numItem = psMetadataLookup(output->analysis, PM_MASK_ANALYSIS_NUM); // Item with number + assert(numItem); + numItem->data.S32++; + + return true; +} + +bool pmMaskIdentifyBadPixels(pmReadout *output, psMaskType maskVal, float thresh, pmMaskIdentifyMode mode) { - PS_ASSERT_IMAGE_NON_NULL(suspects, NULL); - PS_ASSERT_IMAGE_NON_EMPTY(suspects, NULL); - PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, NULL); + PS_ASSERT_PTR_NON_NULL(output, false); + psImage *suspects = psMetadataLookupPtr(NULL, output->analysis, PM_MASK_ANALYSIS_SUSPECT); // Suspect img + if (!suspects) { + psError(PS_ERR_UNKNOWN, false, "Unable to find image with suspected bad pixels."); + return false; + } + PS_ASSERT_IMAGE_NON_EMPTY(suspects, false); + PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, false); + if (output->mask) { + PS_ASSERT_IMAGE_NON_EMPTY(output->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(output->mask, suspects, false); + PS_ASSERT_IMAGE_TYPE(output->mask, PS_TYPE_MASK, false); + } else { + output->mask = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); + } + int num = psMetadataLookupS32(NULL, output->analysis, PM_MASK_ANALYSIS_NUM); // Number of inputs + PS_ASSERT_INT_POSITIVE(num, false); float limit = NAN; // Limit for masking - switch (mode) { case PM_MASK_ID_VALUE: @@ -154,5 +167,5 @@ case PM_MASK_ID_FRACTION: - limit = thresh * nTotal; + limit = thresh * num; break; @@ -200,5 +213,6 @@ limit = max + 1.0 - thresh * sqrtf((float)max + 1.0); - psTrace ("psModules.detrend", 3, "bad: mode: %d, stdev: %f, limit: %f\n", max, sqrtf((float)max + 1.0), limit); + psTrace ("psModules.detrend", 3, "bad: mode: %d, stdev: %f, limit: %f\n", + max, sqrtf((float)max + 1.0), limit); break; } @@ -207,7 +221,4 @@ return NULL; } - - psImage *badpix = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); // Bad pixel mask - psImageInit(badpix, 0); if (psTraceGetLevel("psModules.detrend") > 9) { @@ -227,4 +238,7 @@ psTrace ("psModules.detrend", 3, "bad pixel threshold: %f", limit); + psImage *badpix = output->mask; // Bad pixel mask + psImageInit(badpix, 0); + for (int y = 0; y < suspects->numRows; y++) { for (int x = 0; x < suspects->numCols; x++) { @@ -235,19 +249,19 @@ } - return badpix; + return true; } pmMaskIdentifyMode pmMaskIdentifyModeFromString (const char *string) { - if (!strcasecmp (string, "VALUE")) { + if (!strcasecmp(string, "VALUE")) { return PM_MASK_ID_VALUE; } - if (!strcasecmp (string, "FRACTION")) { + if (!strcasecmp(string, "FRACTION")) { return PM_MASK_ID_FRACTION; } - if (!strcasecmp (string, "SIGMA")) { + if (!strcasecmp(string, "SIGMA")) { return PM_MASK_ID_SIGMA; } - if (!strcasecmp (string, "POISSON")) { + if (!strcasecmp(string, "POISSON")) { return PM_MASK_ID_POISSON; } Index: /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.h (revision 17154) +++ /branches/pap_branch_080320/psModules/src/detrend/pmMaskBadPixels.h (revision 17155) @@ -5,6 +5,6 @@ * @author Eugene Magnier, IfA * - * @version $Revision: 1.15.6.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2008-03-21 03:24:32 $ + * @version $Revision: 1.15.6.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-03-27 23:28:38 $ * Copyright 2004 Institute for Astronomy, University of Hawaii */ @@ -16,6 +16,6 @@ /// @{ -#define PM_MASK_ANALYSIS_MEAN "MASK.MEAN" // Readout analysis metadata keyword with mean value -#define PM_MASK_ANALYSIS_STDEV "MASK.STDEV" // Readout analysis metadata keyword with stdev value +#define PM_MASK_ANALYSIS_SUSPECT "MASK.SUSPECT" // Readout analysis metadata keyword for suspect image +#define PM_MASK_ANALYSIS_NUM "MASK.NUM" // Readout analysis metadata keyword for number of inputs @@ -51,21 +51,21 @@ /// image is of type S32. The relevant median and standard deviation must be supplied in the /// readout->analysis metadata as READOUT.MEDIAN, READOUT.STDEVe -psImage *pmMaskFlagSuspectPixels(psImage *out, ///< Suspected bad pixels image, or NULL - const pmReadout *readout, ///< Readout to inspect - float rej, ///< Rejection threshold (standard deviations) - psMaskType maskVal ///< Mask value for statistics - ); +bool pmMaskFlagSuspectPixels(pmReadout *output, ///< Output readout, optionally with suspect pixels image + const pmReadout *readout, ///< Readout to inspect + float median, ///< Image median + float stdev, ///< Image standard deviation + float rej, ///< Rejection threshold (standard deviations) + psMaskType maskVal ///< Mask value for statistics + ); /// Identify bad pixels from the suspect pixels image /// -/// Bad pixels are identified from the suspect pixels image (accumulated over a large number of images). -/// Pixels marked as suspect in more than "thresh" standard deviations from the mean are identified as bad -/// pixels (output image). If "thresh" is negative, a Poisson is assumed. -psImage *pmMaskIdentifyBadPixels(const psImage *suspects, ///< Accumulated suspect pixels image - psMaskType maskVal, ///< Value to set for bad pixels - int nTotal, - float thresh, ///< Threshold for bad pixel (standard deviations) - pmMaskIdentifyMode mode - ); +/// Bad pixels are identified from the suspect pixels image (accumulated over a large number of images), +/// according to the chosen mode. +bool pmMaskIdentifyBadPixels(pmReadout *output, ///< Output readout, with suspect pixels imageOut + psMaskType maskVal, ///< Value to set for bad pixels + float thresh, ///< Threshold for bad pixel + pmMaskIdentifyMode mode ///< Mode for identifying bad pixels + ); /// @} #endif Index: /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.c =================================================================== --- /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.c (revision 17154) +++ /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.c (revision 17155) @@ -78,4 +78,6 @@ corr->offset = 0.0; corr->offref = 0.0; + corr->num = 0; + corr->stdev = NAN; return corr; @@ -190,13 +192,7 @@ // linear fit to the counts and exptime, given a value for offref -pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, - const psVector *counts, - const psVector *cntError, - const psVector *mask, - float offref, - int nIter, - float rej, - psMaskType maskVal - ) +pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, const psVector *counts, + const psVector *cntError, const psVector *mask, float offref, + int nIter, float rej, psMaskType maskVal) { PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); @@ -249,5 +245,4 @@ return NULL; } - psFree(stats); pmShutterCorrection *corr = pmShutterCorrectionAlloc(); @@ -255,5 +250,8 @@ corr->scale = line->coeff[1][0]; corr->offset = line->coeff[0][1] / line->coeff[1][0]; - + corr->num = stats->clippedNvalues; + corr->stdev = stats->clippedStdev; + + psFree(stats); psFree(x); psFree(y); @@ -263,7 +261,5 @@ } -static psF32 pmShutterCorrectionModel(psVector *deriv, - const psVector *params, - const psVector *x) +static psF32 pmShutterCorrectionModel(psVector *deriv, const psVector *params, const psVector *x) { // This is in a tight loop, so we won't assert here. @@ -283,9 +279,6 @@ // non-linear fit to the counts and exptime, given a guess for the three parameters -pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, - const psVector *counts, - const psVector *cntError, - const pmShutterCorrection *guess - ) +pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, const psVector *counts, + const psVector *cntError, const pmShutterCorrection *guess) { PS_ASSERT_VECTOR_NON_NULL(exptime, NULL); @@ -919,4 +912,16 @@ } + psImage *nums = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize, + PS_TYPE_U16, 0); // Image with number fitted per pixel + if (!nums) { + return false; + } + psImage *sigma = pmReadoutAnalysisImage(shutter, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize, + PS_TYPE_F32, NAN); // Image with stdev per pixel + if (!sigma) { + psFree(nums); + return false; + } + psImage *shutterImage = shutter->image; // Shutter correction image psImage *patternImage; // Illumination pattern @@ -962,8 +967,12 @@ shutterImage->data.F32[yOut][xOut] = NAN; patternImage->data.F32[yOut][xOut] = NAN; + nums->data.U16[yOut][xOut] = 0; + sigma->data.F32[yOut][xOut] = NAN; continue; } shutterImage->data.F32[yOut][xOut] = corr->offset; patternImage->data.F32[yOut][xOut] = corr->scale; + nums->data.U16[yOut][xOut] = corr->num; + sigma->data.F32[yOut][xOut] = corr->stdev; psFree(corr); } @@ -972,4 +981,6 @@ psFree(errors); psFree(counts); + psFree(nums); + psFree(sigma); // Update the "concepts" Index: /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.h =================================================================== --- /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.h (revision 17154) +++ /branches/pap_branch_080320/psModules/src/detrend/pmShutterCorrection.h (revision 17155) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-19 03:40:48 $ + * @version $Revision: 1.14.22.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-03-27 23:28:38 $ * Copyright 2006 Institute for Astronomy, University of Hawaii */ @@ -61,11 +61,11 @@ /// Shutter correction parameters, applicable for a single pixel -typedef struct -{ +typedef struct { double scale; ///< The normalisation for an exposure, A(k) double offset; ///< The time offset, dTk double offref; ///< The reference time offset, dTo -} -pmShutterCorrection; + int num; ///< Number of points used + float stdev; ///< Standard deviation +} pmShutterCorrection; /// Allocator for shutter correction parameters @@ -76,7 +76,8 @@ /// This function is used before doing the full non-linear fit, to get parameters close to the true. Assumes /// exptime vector is sorted (ascending order; longest is last) prior to input. -pmShutterCorrection *pmShutterCorrectionGuess(const psVector *exptime, ///< Exposure times for each exposure - const psVector *counts ///< Counts for each exposure - ); +pmShutterCorrection *pmShutterCorrectionGuess( + const psVector *exptime, ///< Exposure times for each exposure + const psVector *counts ///< Counts for each exposure + ); /// Generate shutter correction based on a linear fit @@ -84,13 +85,14 @@ /// Performs a linear fit to counts as a function of exposure time, with the reference time offset fixed (so /// that the system is linear). Performs iterative clipping, if nIter > 1. -pmShutterCorrection *pmShutterCorrectionLinFit(const psVector *exptime, ///< Exposure times for each exposure - const psVector *counts, ///< Counts for each exposure - const psVector *cntError, ///< Error in the counts - 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 - ); +pmShutterCorrection *pmShutterCorrectionLinFit( + const psVector *exptime, ///< Exposure times for each exposure + const psVector *counts, ///< Counts for each exposure + const psVector *cntError, ///< Error in the counts + 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 + ); /// Generate shutter correction based on a full non-linear fit @@ -99,9 +101,10 @@ /// the reference time offset, so that future fits may be performed using linear fitting with the reference /// time offset fixed. -pmShutterCorrection *pmShutterCorrectionFullFit(const psVector *exptime, ///< Exposure times for each exposure - const psVector *counts, ///< Counts for each exposure - const psVector *cntError, ///< Error in the counts - const pmShutterCorrection *guess ///< Initial guess - ); +pmShutterCorrection *pmShutterCorrectionFullFit( + const psVector *exptime, ///< Exposure times for each exposure + const psVector *counts, ///< Counts for each exposure + const psVector *cntError, ///< Error in the counts + const pmShutterCorrection *guess ///< Initial guess + ); /// Measure a shutter correction image from an array of images @@ -111,12 +114,13 @@ /// measuring the reference time offset using the full non-linear fit for a small number of representative /// regions (middle and corners), and then using that to perform a linear fit to each pixel. -bool pmShutterCorrectionMeasure(pmReadout *output, ///< Output readout - const psArray *readouts, ///< Array of readouts - int size, ///< Size of samples for statistics for non-linear fit - 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 +bool pmShutterCorrectionMeasure( + pmReadout *output, ///< Output readout + const psArray *readouts, ///< Array of readouts + int size, ///< Size of samples for statistics for non-linear fit + 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 ); @@ -124,7 +128,8 @@ /// /// Given a shutter correction (with dT for each pixel), applies this correction to an input image. -bool pmShutterCorrectionApply(pmReadout *readout, ///< Readout to which to apply shutter correction - const pmReadout *shutter ///< Shutter correction readout, with dT for each pixel - ); +bool pmShutterCorrectionApply( + pmReadout *readout, ///< Readout to which to apply shutter correction + const pmReadout *shutter ///< Shutter correction readout, with dT for each pixel + ); ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -156,14 +161,16 @@ /// /// Performs statistics on the readout, recording the data -bool pmShutterCorrectionAddReadout(pmShutterCorrectionData *data, ///< Correction data - const pmReadout *readout, ///< Readout to add - psStatsOptions meanStat, ///< Statistic to use for mean - psStatsOptions stdevStat, ///< Statistic to use for stdev - psMaskType maskVal, ///< Mask value - psRandom *rng ///< Random number generator +bool pmShutterCorrectionAddReadout( + pmShutterCorrectionData *data, ///< Correction data + const pmReadout *readout, ///< Readout to add + psStatsOptions meanStat, ///< Statistic to use for mean + psStatsOptions stdevStat, ///< Statistic to use for stdev + psMaskType maskVal, ///< Mask value + psRandom *rng ///< Random number generator ); /// Calculate the reference shutter time from the correction data -float pmShutterCorrectionReference(const pmShutterCorrectionData *data ///< Correction data +float pmShutterCorrectionReference( + const pmShutterCorrectionData *data ///< Correction data ); @@ -171,12 +178,13 @@ /// /// Performs the linear fit to each pixel in the stack. -bool pmShutterCorrectionGenerate(pmReadout *shutter, ///< Shutter correction - pmReadout *pattern, ///< Background pattern (or NULL) - const psArray *inputs, ///< Stack of input pmReadouts - float reference, ///< Reference shutter time (from pmShutterCorrectionRef) - const pmShutterCorrectionData *data, ///< Correction data - int nIter, ///< Number of iterations - float rej, ///< Rejection threshold (sigma) - psMaskType maskVal ///< Mask value +bool pmShutterCorrectionGenerate( + pmReadout *shutter, ///< Shutter correction + pmReadout *pattern, ///< Background pattern (or NULL) + const psArray *inputs, ///< Stack of input pmReadouts + float reference, ///< Reference shutter time (from pmShutterCorrectionRef) + const pmShutterCorrectionData *data, ///< Correction data + int nIter, ///< Number of iterations + float rej, ///< Rejection threshold (sigma) + psMaskType maskVal ///< Mask value ); Index: /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c (revision 17154) +++ /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.c (revision 17155) @@ -18,4 +18,5 @@ //#define SHOW_BUSY 1 // Show that the function is busy + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// @@ -94,5 +95,24 @@ } - psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination + psStatsOptions combineStdev = 0; // Statistics option for weights + switch (params->combine) { + case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEDIAN: + combineStdev = PS_STAT_SAMPLE_STDEV; + break; + case PS_STAT_ROBUST_MEDIAN: + combineStdev = PS_STAT_ROBUST_STDEV; + break; + case PS_STAT_FITTED_MEAN: + combineStdev = PS_STAT_FITTED_STDEV; + break; + case PS_STAT_CLIPPED_MEAN: + combineStdev = PS_STAT_CLIPPED_STDEV; + break; + default: + psAbort("Should never get here --- checked params->combine before.\n"); + } + + psStats *stats = psStatsAlloc(params->combine | combineStdev); // The statistics to use in the combination if (params->combine == PS_STAT_CLIPPED_MEAN) { stats->clipSigma = params->rej; @@ -107,4 +127,8 @@ } } + if (params->weights && first) { + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + "Using input weights to combine images", ""); + } int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine @@ -120,31 +144,17 @@ psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); - psStatsOptions combineStdev = 0; // Statistics option for weights - if (params->weights) { - if (first) { - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - "Using input weights to combine images", ""); - } - - // Get the correct statistics option for weights - switch (params->combine) { - case PS_STAT_SAMPLE_MEAN: - case PS_STAT_SAMPLE_MEDIAN: - combineStdev = PS_STAT_SAMPLE_STDEV; - break; - case PS_STAT_ROBUST_MEDIAN: - combineStdev = PS_STAT_ROBUST_STDEV; - break; - case PS_STAT_FITTED_MEAN: - combineStdev = PS_STAT_FITTED_STDEV; - break; - case PS_STAT_CLIPPED_MEAN: - combineStdev = PS_STAT_CLIPPED_STDEV; - break; - default: - psAbort("Should never get here --- checked params->combine before.\n"); - } - stats->options |= combineStdev; - } + psImage *counts = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_COUNT, xSize, ySize, + PS_TYPE_U16, 0); + if (!counts) { + return false; + } + psImage *sigma = pmReadoutAnalysisImage(output, PM_READOUT_STACK_ANALYSIS_SIGMA, xSize, ySize, + PS_TYPE_F32, NAN); + if (!sigma) { + psFree(counts); + return false; + } + + stats->options |= combineStdev; // We loop through each pixel in the output image. We loop through each input readout. We determine if @@ -273,4 +283,6 @@ outputMask[yOut][xOut] = params->blank; outputImage[yOut][xOut] = NAN; + counts->data.U16[yOut][xOut] = 0; + sigma->data.F32[yOut][xOut] = NAN; continue; } @@ -288,4 +300,5 @@ maskData[indexData[k]] = 1; numMasked++; + numValid--; } } @@ -293,10 +306,12 @@ for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) { // Don't count the ones that are already masked - if (! maskData[indexData[k]]) { + if (!maskData[indexData[k]]) { maskData[indexData[k]] = 1; numMasked++; + numValid--; } } } + counts->data.U16[yOut][xOut] = numValid; // XXXXX this step probably is very expensive : convert errors to variance everywhere? @@ -314,4 +329,5 @@ outputWeight[yOut][xOut] = NAN; } + sigma->data.F32[yOut][xOut] = NAN; } else { outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); @@ -323,4 +339,5 @@ // also, the weighted mean is not obviously the correct thing here } + sigma->data.F32[yOut][xOut] = psStatsGetValue(stats, combineStdev); } } @@ -338,4 +355,6 @@ psFree(stats); psFree(invScale); + psFree(counts); + psFree(sigma); // Update the "concepts" Index: /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.h (revision 17154) +++ /branches/pap_branch_080320/psModules/src/imcombine/pmReadoutCombine.h (revision 17155) @@ -5,6 +5,6 @@ * @author Paul Price, IfA * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-06-02 03:51:03 $ + * @version $Revision: 1.12.24.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-03-27 23:28:38 $ * Copyright 2004-2006 Institute for Astronomy, University of Hawaii */ @@ -20,6 +20,5 @@ /// These values define how the combination is performed, and should not vary by detector, so that it can be /// re-used for multiple combinations. -typedef struct -{ +typedef struct { psStatsOptions combine; ///< Statistic to use when performing the combination psMaskType maskVal; ///< Mask value @@ -31,6 +30,5 @@ float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) bool weights; ///< Use the supplied weights (instead of calculated stdev)? -} -pmCombineParams; +} pmCombineParams; // Allocator for pmCombineParams