Index: trunk/ippconfig/recipes/ppImage.config =================================================================== --- trunk/ippconfig/recipes/ppImage.config (revision 24890) +++ trunk/ippconfig/recipes/ppImage.config (revision 24891) @@ -25,4 +25,12 @@ USE.DEBURNED.IMAGE BOOL FALSE # use burntool-repaired image? TILTYSTREAK.APPLY BOOL FALSE # apply the 'tiltystreak' tool + +PATTERN BOOL FALSE # Fit and remove pattern noise? +PATTERN.ORDER S32 1 # Polynomial order +PATTERN.ITER S32 3 # Rejection iterations +PATTERN.REJ F32 2.0 # Rejection threshold +PATTERN.MEAN STR ROBUST_MEDIAN # Statistic for mean +PATTERN.STDEV STR ROBUST_STDEV # Statistic for standard deviation + TILTYSTREAK.BY.CLASS METADATA # apply the 'tiltystreak' tool Index: trunk/ppImage/src/ppImageDetrendReadout.c =================================================================== --- trunk/ppImage/src/ppImageDetrendReadout.c (revision 24890) +++ trunk/ppImage/src/ppImageDetrendReadout.c (revision 24891) @@ -64,27 +64,27 @@ if (options->doVarianceBuild) { // create the target mask and variance images - psImage *noiseImage = NULL; - if (options->doNoiseMap) { - // XXX convert the noiseMap image to a binned image - pmReadout *noiseMap = NULL; - noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP"); - noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32); - psImageInit (noiseImage, 0.0); - - // XXX this works, but is not really quite right: the model shoud include the - // offset information, we are not really getting exactly the right mapping from the - // original file. - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = noiseMap->image->numCols; - binning->nYruff = noiseMap->image->numRows; - binning->nXfine = input->image->numCols; - binning->nYfine = input->image->numRows; - psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT); - - psImageUnbin (noiseImage, noiseMap->image, binning); - psFree (binning); - } + psImage *noiseImage = NULL; + if (options->doNoiseMap) { + // XXX convert the noiseMap image to a binned image + pmReadout *noiseMap = NULL; + noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP"); + noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32); + psImageInit (noiseImage, 0.0); + + // XXX this works, but is not really quite right: the model shoud include the + // offset information, we are not really getting exactly the right mapping from the + // original file. + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = noiseMap->image->numCols; + binning->nYruff = noiseMap->image->numRows; + binning->nXfine = input->image->numCols; + binning->nYfine = input->image->numRows; + psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT); + + psImageUnbin (noiseImage, noiseMap->image, binning); + psFree (binning); + } pmReadoutGenerateVariance(input, noiseImage, true); - psFree (noiseImage); + psFree (noiseImage); } @@ -123,4 +123,15 @@ } } + + // Pattern noise correction + if (options->doPattern) { + if (!pmPatternRow(input, options->patternOrder, +options->patternIter, options->patternRej, options->patternMean, +options->patternStdev)) { + psFree(detview); + return false; + } + } + // Normalization by a single (known) constant Index: trunk/psModules/src/detrend/Makefile.am =================================================================== --- trunk/psModules/src/detrend/Makefile.am (revision 24890) +++ trunk/psModules/src/detrend/Makefile.am (revision 24891) @@ -16,5 +16,6 @@ pmShifts.c \ pmDark.c \ - pmRemnance.c + pmRemnance.c \ + pmPattern.c # pmSkySubtract.c @@ -33,5 +34,6 @@ pmShifts.h \ pmDark.h \ - pmRemnance.h + pmRemnance.h \ + pmPattern.h # pmSkySubtract.h Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 24891) +++ trunk/psModules/src/detrend/pmPattern.c (revision 24891) @@ -0,0 +1,84 @@ +#include +#include + +#include "pmFPA.h" + +#include "pmPattern.h" + + +bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, + psStatsOptions clipMean, psStatsOptions clipStdev) +{ + PM_ASSERT_READOUT_NON_NULL(ro, false); + PM_ASSERT_READOUT_IMAGE(ro, false); + PS_ASSERT_INT_NONNEGATIVE(order, false); + PS_ASSERT_INT_NONNEGATIVE(iter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + + psImage *image = ro->image; // Image to correct + int numCols = image->numCols, numRows = image->numRows; // Size of image + + // Indices are distributed [-1:1) + psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fitting + float norm = 2.0 / (float)numCols; // Normalisation for indices + for (int x = 0; x < numCols; x++) { + indices->data.F32[x] = x * norm - 1.0; + } + + psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics + clip->clipIter = iter; + clip->clipSigma = rej; + psVector *mask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for fitting + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit + psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit + + for (int y = 0; y < numRows; y++) { + psVectorInit(mask, 0); + data = psImageRow(data, image, y); + int num = 0; // Number of good pixels + for (int x = 0; x < numCols; x++) { + if (ro->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] > 0) { + mask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; + } else { + mask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0; + num++; + } + } + if (num < order + 1) { + // Not enough points to fit + continue; + } + if (!psVectorClipFitPolynomial1D(poly, clip, mask, 0xFF, data, NULL, indices)) { + psError(PS_ERR_UNKNOWN, false, "Unable to fit polynomial to row %d", y); + psFree(indices); + psFree(clip); + psFree(mask); + psFree(poly); + psFree(data); + return false; + } + psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "Unable to evaluate polynomial for row %d", y); + psFree(indices); + psFree(clip); + psFree(mask); + psFree(poly); + psFree(data); + return false; + } + + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] -= solution->data.F32[x]; + } + psFree(solution); + } + + psFree(indices); + psFree(clip); + psFree(mask); + psFree(poly); + psFree(data); + + return true; +} Index: trunk/psModules/src/detrend/pmPattern.h =================================================================== --- trunk/psModules/src/detrend/pmPattern.h (revision 24891) +++ trunk/psModules/src/detrend/pmPattern.h (revision 24891) @@ -0,0 +1,34 @@ +/* @file pmPattern.h + * @brief Fit and remove pattern noise + * + * @author Paul Price, IfA + * + * @version $Revision: 1.16 $ + * @date $Date: 2009-02-12 19:25:52 $ + * Copyright 2004-2009 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_PATTERN_H +#define PM_PATTERN_H + +#include +#include + +/// @addtogroup detrend Detrend Creation and Application +/// @{ + +/// Fit and remove pattern noise over rows +bool pmPatternRow( + pmReadout *ro, ///< Readout to correct + int order, ///< Polynomial order + int iter, ///< Number of clipping iterations + float rej, ///< Rejection threshold for clipping + psStatsOptions clipMean, ///< Statistic to use for mean + psStatsOptions clipStdev ///< Statistic to use for standard deviation + ); + +/// @} +#endif + + + Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 24890) +++ trunk/psModules/src/psmodules.h (revision 24891) @@ -77,4 +77,5 @@ #include #include +#include // the following headers are from psModule:astrom