Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 41743) +++ trunk/psModules/src/detrend/pmPattern.c (revision 41892) @@ -9,4 +9,107 @@ #define PATTERN_ROW_BKG_FIX 1 +/* some in-line notes: + + patternMaskRow sets the data value to NAN, inconsistent with new plan + + here is the outline of pmPatternRow + + * at this point, we have already done overscan subtraction, right? + + * measure stats on the full cell (MEDIAN, STDEV) + ** subsample? + ** if it fails, it masks the entire cell and set the value to NAN + + * calculate an upper and lower threshold (median +/- T * sigma) + * define a normalized x-coordinate ('index') : + ** see note below on chebys + + * each row is treated independently + * pixels are masked for the fit if they are out-of-range + or if they are already masked + + ** note that the clipping threshold will be larger if there + are pixels which have astronomical structures + + a possible better option would be to set the threshold based on the median + and a sigma calculated from Poisson stats (do we know the gain?) + + ** fit is allowed to proceed if even N+1 pixels exist, which is clearly too low + + ** Remaining pixels are fitted with clip-fit + + ** solution is subtracted from the data + (this is implemented with psPolynomial1DEvalVector) + perhaps faster if we fixed the order to 2 and hardwired the result + + * after each row is fitted, the intercept (A value) is fitted + as a function of the y-coordinate and the result is subtracted + + * the slope value is also fitted as a function of the column + and added back in -- I'm not sure I understand this step. + + ***************** + + ** what we calculate are related to chebychevs (domain is -1 : +1) + *** T0(x) = 1 + *** T1(x) = x + *** T2(x) = 2x^2 - 1 + + *** we calculate y = A + Bx + Cx^2 + + a_0 + a_1 x + a_2 (2x^2 - 1) = A + B + Cx^2 + + a_1 = B + a_0 - a_2 = A + 2 a_2 = C + + a_0 = A + C/2 + a_1 = B + a_2 = C/2 + + ***************** + + I have 3 goals in re-working the code: + + 1) improve overall speed + 2) improve reliability of the fit + 3) skip fit if we can + + Let's assume the signal in the cell is light + bias drift + + The bias drift has an amplitude of ~5 - 10 DN + + That makes a detectable source with ~N * a few counts (multiple pixels in a row) + + So, the effective flux is ~10 * 5 = 50 DN + for which sky level is this value - 3 sigma? + + 50 / sqrt(sky sigma^2 * effective area) + + area ~ 5pixels, sky sigma^2 = sky + + 10 * Npix / sqrt(sky * Npix) = 3 + + Npix = sky * (S/N)^2 / (peak^2) + + sky = Npix * peak^2 / SN^2 + + if (sky < Npix * peak^2 / SN^2), we should skip: + + Npix ~ 5 + peak ~ 10 + SN ~ 3 (or even less) + + sky < 5 * 100 / 9 = 55 or so + + + To address these in order: + + 1) speed: + * the analysis threaded is not threaded: thread across cells + + 2) + + */ // Mask a row as bad @@ -33,9 +136,47 @@ } +// Comparison and swap functions for sorting values directly +#define SORT_COMPARE(A,B) (sampleArray[A] < sampleArray[B]) +#define SORT_SWAP(TYPE,A,B) { \ + if (A != B) { \ + TYPE temp = sampleArray[A]; \ + sampleArray[A] = sampleArray[B]; \ + sampleArray[B] = temp; \ + } \ +} + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Measurement and application ////////////////////////////////////////////////////////////////////////////////////////////////////////////// +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); + + +// XXX allow user choice of binned vs unbinned analysis? bool pmPatternRow(pmReadout *ro, int order, int iter, float rej, float thresh, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) { + + bool status = false; + if (true) { + status = pmPatternRowBinned(ro, order, iter, rej, thresh, clipMean, clipStdev, maskVal, maskBad); + } else { + 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 pmPatternRowUnbinned(pmReadout *ro, int order, int iter, float rej, float thresh, psStatsOptions clipMean, psStatsOptions clipStdev, psImageMaskType maskVal, psImageMaskType maskBad) @@ -48,4 +189,11 @@ PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + psImage *image = ro->image; // Image to correct psImage *mask = ro->mask; // Mask for image @@ -55,20 +203,57 @@ psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { - psWarning("Unable to calculate statistics on readout; masking entire readout."); + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); psErrorClear(); psFree(stats); psFree(rng); - psImageInit(image, NAN); - if (mask) { - psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); - } - if (ro->variance) { - psImageInit(image, NAN); - } return true; } + +# if (USE_BACKGROUND_STDEV) float lower = stats->robustMedian - thresh * stats->robustStdev; // Lower bound for data float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data float background = stats->robustMedian; +# else + // the signal we are looking for is a small variation on top of the background. if + // the background is uniform with only read noise + sky noise, then the pixel-to-pixel + // stdev should only be due to known noise sources and predictable. If the + // pixel-to-pixel variations are from other features, then those variations will + // probably dominate the row-by-row bias variations. + + // instead of using the image pixel statistics to measure the stdev, lets assume only + // dark noise plus poisson sky noise. we are not carrying in the read noise, but it is + // fairly modest for GPC1 (~10 DN) + + // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have + // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17 + + // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2) + // so smaller by a factor of 1.4 than what we predict, which is not very large + + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + + // XXX retrieve noise and gain from 'concepts' if possible +# define READNOISE 10 /* arbitrary number */ + float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data + float background = stats->robustMedian; + + // if the noise from the background is too large + float significance = nominalAmplitude / sigma; + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + return true; + } + fprintf (stderr, "correcting pattern row background %s: %f - %f - %f : %f %f %f\n", cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); +# endif + psFree(stats); psFree(rng); @@ -109,18 +294,35 @@ psVectorInit(yaxisMask, 0); #endif + + // we really need more than order + 1 points (= 4). + // this should be tunable, but let's try 5 - 10% + int validNmin = numCols * 0.1; + for (int y = 0; y < numRows; y++) { psVectorInit(clipMask, 0); data = psImageRow(data, image, y); int num = 0; // Number of good pixels + + // if the unmasked pixels only span a small range in x then we cannot fit the + // 2nd order polynomial variations very well. Require a minimum fractional range + float validXmin = +1; + float validXmax = -1; + + // XXX can we do just as well fitting 1/3 of the pixels? (NOT REALLY) + // (x % 3) || for (int x = 0; x < numCols; x++) { - if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || - data->data.F32[x] < lower || data->data.F32[x] > upper) { - clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; + if ((mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) || + data->data.F32[x] < lower || data->data.F32[x] > upper) { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0xFF; } else { clipMask->data.PS_TYPE_VECTOR_MASK_DATA[x] = 0; num++; + validXmin = PS_MIN(indices->data.F32[x], validXmin); + validXmax = PS_MAX(indices->data.F32[x], validXmax); } } - if (num < order + 1) { + + // XXX how much time is spent in the fitting + if (num < validNmin) { // Not enough points to fit patternMaskRow(ro, y, maskBad); @@ -131,4 +333,5 @@ continue; } + // 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)) { psWarning("Unable to fit polynomial to row %d", y); @@ -215,5 +418,5 @@ for (int x = 0; x < numCols; x++) { image->data.F32[y][x] += solution->data.F32[y]; - psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]); + psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]); } corr->data.F64[y][0] -= solution->data.F32[y]; @@ -241,5 +444,6 @@ for (int x = 0; x < numCols; x++) { image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x]; - psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],indices->data.F32[x]); + // XXX EAM : this was [x] which is wrong + psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]); } corr->data.F64[y][1] -= solution->data.F32[y] ; @@ -262,4 +466,339 @@ psFree(indices); + psFree(clip); + psFree(clipMask); + psFree(poly); + psFree(data); + + return true; +} + +# define NPIX 15 + +// 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, + psStatsOptions clipMean, psStatsOptions clipStdev, + psImageMaskType maskVal, psImageMaskType maskBad) +{ + 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); + PS_ASSERT_FLOAT_LARGER_THAN(thresh, 0.0, false); + + bool mdok; // Status of MD lookup + + pmCell *cell = ro->parent; + pmChip *chip = cell->parent; + const char *chipName = psMetadataLookupStr(&mdok, chip->concepts, "CHIP.NAME"); // Name of chip + const char *cellName = psMetadataLookupStr(&mdok, cell->concepts, "CELL.NAME"); // Name of cell + + psImage *image = ro->image; // Image to correct + psImage *mask = ro->mask; // Mask for image + int numCols = image->numCols, numRows = image->numRows; // Size of image + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + if (!psImageBackground(stats, NULL, ro->image, ro->mask, maskVal, rng)) { + psWarning("Unable to calculate statistics on readout; skipping pattern correction for %s, %s.", chipName, cellName); + psErrorClear(); + psFree(stats); + psFree(rng); + + // EAM 20211011 : we used to mask cells which fail the above, but this seems excessive + // psImageInit(image, NAN); + // if (mask) { + // psBinaryOp(mask, mask, "|", psScalarAlloc(maskBad, PS_TYPE_IMAGE_MASK)); + // } + // if (ro->variance) { + // psImageInit(image, NAN); + // } + return true; + } + + // 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 + float upper = stats->robustMedian + thresh * stats->robustStdev; // Upper bound for data + float background = stats->robustMedian; +# else + // the signal we are looking for is a small variation on top of the background. if + // the background is uniform with only read noise + sky noise, then the pixel-to-pixel + // stdev should only be due to known noise sources and predictable. If the + // pixel-to-pixel variations are from other features, then those variations will + // probably dominate the row-by-row bias variations. + + // instead of using the image pixel statistics to measure the stdev, lets assume only + // dark noise plus poisson sky noise. we are not carrying in the read noise, but it is + // fairly modest for GPC1 (~10 DN) + + // if we assume a gain of 1 and the read noise of 10 DN, then a sky of 200 would have + // a noise of N = sqrt (1 * 200 + 10^2) = sqrt (300) ~ 17 + + // if the gain were as much as 2, then the noise in DN would be N = sqrt(2 * (200 + 100)) / 2 = sqrt(300) / sqrt(2) + // so smaller by a factor of 1.4 than what we predict, which is not very large + + // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file) + float nominalAmplitude = psMetadataLookupF32 (&mdok, cell->analysis, "PTN.ROW.AMP"); + if (!mdok) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + if (!isfinite(nominalAmplitude)) nominalAmplitude = 20; // XXX EAM : somewhat arbitrary number + // If we cannot determine the nominal amplitude, we fall-back on the worst case + + // retrieve noise and gain from 'concepts' if possible +# define READNOISE 10 /* arbitrary number */ + float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); + float lower = stats->robustMedian - delta; // Lower bound for data + float upper = stats->robustMedian + delta; // Upper bound for data + float background = stats->robustMedian; + + // if the noise from the background is too large + float significance = nominalAmplitude / sigma; + + // XXX EAM : arbitrary number + if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psFree(stats); + psFree(rng); + return true; + } + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); +# endif + + psFree(stats); + psFree(rng); + + // 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(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; + 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 + clip->clipIter = iter; + clip->clipSigma = rej; + 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 + + psImage *corr = psImageAlloc(order + 1, numRows, PS_TYPE_F64); // Corrections applied + psImageInit(corr, NAN); + + // CZW: 2011-11-30 + // Define the vectors to hold the "x" and "y" slope trends. + // Briefly, the slope trend in the y-axis is a due to variations in the 0-th order term + // of the PATTERN.ROW fit between individual rows across the cell. Similarly, the 1-st + // order term of the PATTERN.ROW fit defines the trend in the x-axis (as that's what we + // are fitting with PATTERN.ROW in the first place). However, the thing we're trying to + // fix with PATTERN.ROW is the detector level bias wiggles. These should be overlaid on + // the true sky level. Therefore, simply applying the PATTERN.ROW correction will + // introduce cell-to-cell sky variations as these two trends are removed. To avoid this, + // We store the 0th and 1st order values used for each row, and then fit a polynomial to + // these results. By re-adding these systematic trends back, we can remove the row-to-row + // variations without improperly removing the real sky trend. + psVector *yaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the constant term + psVector *yaxisMask = psVectorAlloc(numRows, PS_TYPE_VECTOR_MASK); // Mask for rows with no fit + psVector *xaxisData = psVectorAlloc(numRows, PS_TYPE_F32); // Data to fit to the linear term + psVectorInit(yaxisMask, 0); + + // validNmin is the minimum number of samples needed to measure the trend. + // this should be tunable, but let's try 5 - 10% + int validNmin = PS_MAX (nSamples * 0.1, order + 2); + + for (int y = 0; y < numRows; y++) { + psVectorInit(clipMask, 0); + data = psImageRow(data, image, y); + int num = 0; // Number of good pixels + + // if the unmasked pixels only span a small range in x then we cannot fit the + // 2nd order polynomial variations very well. Require a minimum fractional range + float validXmin = +1; + float validXmax = -1; + + for (int sample = 0; sample < nSamples; sample ++) { + + // 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; // real pixel elements in x-dir + psAssert (pix >= 0, "invalid pix value"); + psAssert (pix < numCols, "invalid pix value"); + 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]; // store the value to be sorted + seq ++; + } + if (seq < 1) { + clipMask->data.PS_TYPE_VECTOR_MASK_DATA[sample] = 0xFF; + 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; + if (seq % 2 == 1) { psAssert (midPt >= 0, "invalid midPt"); } + if (seq % 2 == 0) { psAssert (midPt >= 1, "invalid midPt"); } + psAssert (midPt < NPIX, "invalid midPt"); + + float medValue = (seq % 2) ? sampleArray[midPt] : 0.5*(sampleArray[midPt] + sampleArray[midPt-1]); + yFit->data.F32[sample] = medValue; + } + + // If not enough points are valid, skip + if (num < validNmin) { + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + // 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, yFit, NULL, xFit)) { + psWarning("Unable to fit polynomial to row %d", y); + psErrorClear(); + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + // Store the results we found for this row. + yaxisData->data.F32[y] = poly->coeff[0]; + xaxisData->data.F32[y] = poly->coeff[1]; + psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]); + + memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); + psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector + if (!solution) { + psWarning("Unable to evaluate polynomial for row %d", y); + psErrorClear(); + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; + continue; + } + + psAssert (solution->n == numCols, "oops"); + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] -= solution->data.F32[x]; + psTrace("pattern",5,"A: %d %d %g\n",x,y,solution->data.F32[x]); + } + psFree(solution); + } + + // Put the global trends back that were removed by the PATTERN.ROW correction. + // Set up the indices for the polynomial + psVector *yaxisIndices = psVectorAlloc(numRows, PS_TYPE_F32); + norm = 2.0 / (float)numRows; + for (int y = 0; y < numRows; y++) { + yaxisIndices->data.F32[y] = y * norm - 1.0; + psTrace("psModules.detrend.pattern",10,"%d %f %f\n",y,yaxisIndices->data.F32[y],yaxisData->data.F32[y]); + } + + // Fit the trend of the constant term, producing the y-axis global trend + psStatsInit(clip); + psPolynomial1D *yaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(yaxisPoly,clip,yaxisMask,0xFF,yaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to y-axis trend"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + } + corr->data.F64[y][0] -= background; + } + } else { + psVector *solution = psPolynomial1DEvalVector(yaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + // If we've failed, we need to do something, so add back in the background level, and + // expect that the final image will have background mismatches. + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += background; + } + corr->data.F64[y][0] -= background; + } + } else { + psAssert (solution->n == numRows, "oops"); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y]; + // XXX EAM : this was [x], which is wrong + psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[y]); + } + corr->data.F64[y][0] -= solution->data.F32[y]; + } + } + psFree(solution); + } + + // Fit the trend of the linear term, producing the x-axis global trend + // We can use the same mask vector, as the same rows failed the row-fit earlier. + psStatsInit(clip); + psPolynomial1D *xaxisPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); // Polynomial to fit. + if (!psVectorClipFitPolynomial1D(xaxisPoly,clip,yaxisMask,0xFF,xaxisData, NULL, yaxisIndices)) { + psWarning("Unable to fit polynomial to x-axis trend"); + psErrorClear(); + } else { + psVector *solution = psPolynomial1DEvalVector(xaxisPoly,yaxisIndices); + if (!solution) { + psWarning("Unable to evaluate polynomial"); + psErrorClear(); + } else { + psAssert (solution->n == numRows, "oops"); + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += solution->data.F32[y] * indices->data.F32[x]; + // XXX EAM : this was set to [x] which is wrong (numCols > numRows) + psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[y],indices->data.F32[x]); + } + corr->data.F64[y][1] -= solution->data.F32[y] ; + } + } + psFree(solution); + } + psFree(yaxisPoly); + psFree(xaxisPoly); + psFree(yaxisIndices); + psFree(yaxisMask); + psFree(yaxisData); + psFree(xaxisData); + + psMetadataAddImage(ro->analysis, PS_LIST_TAIL, PM_PATTERN_ROW_CORRECTION, PS_META_REPLACE, + "Pattern row correction", corr); + psFree(corr); + + psFree(indices); + psFree(xFit); + psFree(yFit); psFree(clip); psFree(clipMask); @@ -1316,2 +1855,3 @@ return true; } +