Index: /trunk/psModules/src/detrend/pmMaskBadPixels.c =================================================================== --- /trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 9980) +++ /trunk/psModules/src/detrend/pmMaskBadPixels.c (revision 9981) @@ -4,4 +4,5 @@ #include +#include #include @@ -67,2 +68,110 @@ return true; } + + +psImage *pmMaskFlagSuspectPixels(psImage *out, const pmReadout *readout, float rej, + psMaskType maskVal, float frac, psRandom *rng) +{ + PS_ASSERT_PTR_NON_NULL(readout, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(frac, 0.0, 1.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); + 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); + } + + psImage *image = readout->image; // Image of interest + psImage *mask = readout->mask; // Corresponding mask + + if (rng) { + psMemIncrRefCounter(rng); + } else { + rng = psRandomAlloc(PS_RANDOM_TAUS, 0); + } + + psStats *stats = psImageBackground(image, mask, maskVal, 0.25, 0.75, + frac * image->numCols * image->numRows, rng); // Image statistics + if (!stats || !isfinite(stats->robustMedian) || !isfinite(stats->robustUQ) || + !isfinite(stats->robustLQ)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure image statistics.\n"); + psFree(stats); + psFree(rng); + return NULL; + } + psFree(rng); + + float median = stats->robustMedian; // Median value + float stdev = 0.74*(stats->robustUQ - stats->robustLQ); // Estimate of the standard deviation + psFree(stats); + + if (!out) { + out = psImageAlloc(image->numCols, image->numRows, PS_TYPE_S32); + psImageInit(out, 0); + } + + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + 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, float thresh, psMaskType maskVal) +{ + PS_ASSERT_IMAGE_NON_NULL(suspects, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(suspects, NULL); + PS_ASSERT_IMAGE_TYPE(suspects, PS_TYPE_S32, NULL); + + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics + if (!psImageStats(stats, suspects, NULL, 0) || !isfinite(stats->sampleMean) || + !isfinite(stats->sampleStdev)) { + psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics.\n"); + psFree(stats); + return NULL; + } + + float mean = stats->sampleMean; // Mean value + float stdev = stats->sampleStdev; // Standard deviation + + psImage *badpix = psImageAlloc(suspects->numCols, suspects->numRows, PS_TYPE_MASK); // Bad pixel mask + psImageInit(badpix, 0); + + if (psTraceGetLevel("psModules.detrend") > 9) { + psStats *stats = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics + psImageStats(stats, suspects, NULL, 0); + psHistogram *histo = psHistogramAlloc(-0.5, stats->max + 0.5, stats->max + 1); + psImageHistogram(histo, suspects, NULL, 0); + for (int i = 0; i < histo->nums->n; i++) { + printf("%f --> %f : %f\n", histo->bounds->data.F32[i], histo->bounds->data.F32[i + 1], + histo->nums->data.F32[i]); + } + psFree(stats); + psFree(histo); + printf("Threshold: %f\n", mean + thresh * stdev); + } + + for (int y = 0; y < suspects->numRows; y++) { + for (int x = 0; x < suspects->numCols; x++) { + if (suspects->data.S32[y][x] - mean >= thresh * stdev) { + badpix->data.PS_TYPE_MASK_DATA[y][x] = maskVal; + } + } + } + + return badpix; +} Index: /trunk/psModules/src/detrend/pmMaskBadPixels.h =================================================================== --- /trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 9980) +++ /trunk/psModules/src/detrend/pmMaskBadPixels.h (revision 9981) @@ -8,6 +8,6 @@ /// @author Eugene Magnier, IfA /// -/// @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ -/// @date $Date: 2006-10-17 22:17:38 $ +/// @version $Revision: 1.9 $ $Name: not supported by cvs2svn $ +/// @date $Date: 2006-11-15 00:27:52 $ /// /// Copyright 2004 Institute for Astronomy, University of Hawaii @@ -33,3 +33,30 @@ ); +/// Find pixels outlying from the background, flagging suspect pixels +/// +/// Pixels more than "rej" standard deviations from the background level (in flat-fielded images) have the +/// corresponding pixel in the "suspect pixels" image incremented. After accumulating over a suitable sample +/// of images, bad pixels should have a high value in the suspect pixels image, allowing them to be +/// identified. The suspect pixels image is of type S32. +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 + float frac, ///< Fraction of pixels to consider + psRandom *rng ///< Random number generator + ); + +/// 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). +psImage *pmMaskIdentifyBadPixels(const psImage *suspects, ///< Accumulated suspect pixels image + float thresh, ///< Threshold for bad pixel (standard deviations) + psMaskType maskVal ///< Value to set for bad pixels + ); + + + + #endif