Index: branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c (revision 41814) +++ branches/eam_branches/ipp-dev-20210817/psModules/src/detrend/pmPattern.c (revision 41817) @@ -137,10 +137,10 @@ // Comparison and swap functions for sorting values directly -#define SORT_COMPARE(A,B) (value[A] < value[B]) -#define SORT_SWAP(A,B) { \ +#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B]) +#define SORT_SWAP(TYPE,A,B) { \ if (A != B) { \ - float temp = value[A]; \ - value[A] = value[B]; \ - value[B] = temp; \ + TYPE temp = sampleArray[A]; \ + sampleArray[A] = sampleArray[B]; \ + sampleArray[B] = temp; \ } \ } @@ -150,9 +150,28 @@ ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad); + + +bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad); + + +bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) { + +// bool status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + bool status = pmPatternRowUnbinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + return status; +} + // USE_BACKGROUND_STDEV: if true, the analysis will use the measured robust stdev to clip the out-of-range pixles // if false, the stdev will be estimated based on the Poisson statistics of the median (sky level). # define USE_BACKGROUND_STDEV 0 -bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, +bool pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh, psStatsOptions clipMean, psStatsOptions clipStdev, psImageMaskType maskVal, psImageMaskType maskBad) @@ -233,6 +252,5 @@ psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics - // XXX clip->clipIter = iter; - clip->clipIter = 1; // XXX skip iteration for a test + clip->clipIter = iter; clip->clipSigma = rej; psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping @@ -440,5 +458,7 @@ } -// bin by Npix in the x-direction to reduce the number of calculations needed to measure +# define NPIX 31 + +// bin by NPIX in the x-direction to reduce the number of calculations needed to measure // the pattern bool pmPatternRowBinned(pmReadout *ro, int order, int iter, float rej, float thresh, @@ -474,4 +494,9 @@ } + // if USE_BACKGROUND_STDEV is TRUE, the observed standard deviation is used to set the + // thresholds. this is going to be an overestimate if there is any structure in the + // image. If FALSE, the thresholds are set based on poisson stats for the background + // level. We assume the gain is 1, so this is an overestimate if the gain is > 1 + # if (USE_BACKGROUND_STDEV) float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data @@ -497,6 +522,7 @@ # define READNOISE 10 float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); - float lower = stats->robustMedian - thresh * sigma; // Lower bound for data - float upper = stats->robustMedian + thresh * sigma; // Upper bound for data + float delta = PS_MIN (thresh * sigma, 40); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data float background = stats->robustMedian; # endif @@ -513,24 +539,28 @@ psFree(rng); -# define NPIX 15 - - // Indices are distributed [-1:1] [-1 = 0, +1 = numCols = indices[nSamples] + // the vector 'indices' maps the x-coordinate to a range [-1:1]. the element number (i) of indices + // related to the x-coordinate (column number) by x = (i + 0.5) * NPIX + int nSamples = numCols / NPIX; - psVector *indices = psVectorAlloc(nSamples, PS_TYPE_F32); // Indices for fitting - psVector *fitData = psVectorAlloc(nSAmples, PS_TYPE_F32); // Data to fit - - // indices elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX + + psVector *indices = psVectorAlloc(numCols, PS_TYPE_F32); // Indices for fit solutions + psVector *xFit = psVectorAlloc(nSamples, PS_TYPE_F32); // x-coordinate for fitting + psVector *yFit = psVectorAlloc(nSamples, PS_TYPE_F32); // flux values for fitting + + // xFit elements run from 0 - nSamples, element 'sample' corresponds to the middle of the bin sample*NPIX + 0.5*NPIX float norm = 2.0 / (float)numCols; // Normalisation for indices for (int sample = 0; sample < nSamples; sample ++) { int x = (sample + 0.5)*NPIX; - indices->data.F32[sample] = x * norm - 1.0; + xFit->data.F32[sample] = x * norm - 1.0; + } + for (int x = 0; x < numCols; x ++) { + indices->data.F32[x] = x * norm - 1.0; } psStats *clip = psStatsAlloc(clipMean | clipStdev); // Clipping statistics - // XXX clip->clipIter = iter; - clip->clipIter = 1; // XXX skip iteration for a test + clip->clipIter = iter; clip->clipSigma = rej; - psVector *clipMask = psVectorAlloc(numCols, PS_TYPE_VECTOR_MASK); // Mask for clipping + psVector *clipMask = psVectorAlloc(nSamples, PS_TYPE_VECTOR_MASK); // Mask for clipping psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, order); // Polynomial to fit psVector *data = psVectorAlloc(numCols, PS_TYPE_F32); // Data to fit @@ -539,5 +569,4 @@ psImageInit(corr, NAN); -#ifdef PATTERN_ROW_BKG_FIX // CZW: 2011-11-30 // Define the vectors to hold the "x" and "y" slope trends. @@ -556,9 +585,8 @@ psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term psVectorInit(yaxisMask, 0); -#endif - - // we really need more than order + 1 points (= 4). + + // validNmin is the minimum number of samples needed to measure the trend. // this should be tunable, but let's try 5 - 10% - int validNmin = numCols * 0.1; + int validNmin = PS_MAX (nSamples * 0.1, order + 2); for (int y = 0; y < numRows; y++) { @@ -572,33 +600,39 @@ float validXmax = -1; - // XXX bookkeeping might be easier if this loop is over elements of 'indices' - // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY) - // (x % 3) || for (int sample = 0; sample < nSamples; sample ++) { - // store valid samples in the array + // store valid samples in the array to be sorted float sampleArray[NPIX]; int seq = 0; for (int j = 0; j < NPIX; j++) { - int pix = sample * NPIX + j; + int pix = sample * NPIX + j; // real pixel elements in x-dir if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][pix] & maskVal)) continue; if (data->data.F32[pix] < lower || data->data.F32[pix] > upper) continue; - sampleArray[seq] = data->data.F32[pix]; + sampleArray[seq] = data->data.F32[pix]; // store the value to be sorted seq ++; } if (seq < 1) { clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF; - } else { - clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0; - num++; - validXmin = PS_MIN(indices->data.F32[sample], validXmin); - validXmax = PS_MAX(indices->data.F32[sample], validXmax); - } - - PSSORT (seq, sampleArray - - } - - // XXX how much time is spent in the fitting + yFit->data.F32[sample] = NAN; + continue; + } + // note that we are treating the x-coordinate as the center + // of the binned pixel group, even if some or most pixels have + // been masked. compared to the amplitude of the slope, this + // error is small + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0; + validXmin = PS_MIN(xFit->data.F32[sample], validXmin); + validXmax = PS_MAX(xFit->data.F32[sample], validXmax); + num++; + + // PSSORT operates on sampleArray (see define of macro SORT_SWAP above) + PSSORT (seq, SORT_COMPARE, SORT_SWAP, float); + + int midPt = 0.5 * seq; + float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]); + yFit->data.F32[sample] = medValue; + } + + // XXX how much time is spent in the fitting? to test: make this always true if (num < validNmin) { // Not enough points to fit @@ -609,5 +643,5 @@ } // XXX does this need to be a clipped fit if we are clipping based on the median poisson noise? - if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, data, NULL, indices)) { + if (!psVectorClipFitPolynomial1D(poly, clip, clipMask, 0xFF, yFit, NULL, xFit)) { psWarning("Unable to fit polynomial to row %d", y); psErrorClear(); @@ -628,7 +662,5 @@ psErrorClear(); patternMaskRow(ro, y, maskBad); -#ifdef PATTERN_ROW_BKG_FIX yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; -#endif continue; } @@ -641,5 +673,4 @@ } -#ifdef PATTERN_ROW_BKG_FIX // Put the global trends back that were removed by the PATTERN.ROW correction. // Set up the indices for the polynomial @@ -665,6 +696,5 @@ corr->data.F64[y][0] -= background; } - } - else { + } else { psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices); if (!solution) { @@ -679,6 +709,5 @@ corr->data.F64[y][0] -= background; } - } - else { + } else { for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { @@ -699,12 +728,10 @@ psWarning("Unable to fit polynomial to x-axis trend"); psErrorClear(); - } - else { + } else { psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices); if (!solution) { psWarning("Unable to evaluate polynomial"); psErrorClear(); - } - else { + } else { for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { @@ -723,6 +750,4 @@ psFree(yaxisData); psFree(xaxisData); - // End PATTERN_ROW_BKG_FIX global trend replacement -#endif psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, @@ -731,4 +756,6 @@ psFree(indices); + psFree(xFit); + psFree(yFit); psFree(clip); psFree(clipMask);