Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 27676) +++ trunk/psModules/src/detrend/pmPattern.c (revision 32970) @@ -6,4 +6,6 @@ #include "pmPattern.h" + +#define PATTERN_ROW_BKG_FIX 1 @@ -89,4 +91,22 @@ psImageInit(corr, NAN); +#ifdef PATTERN_ROW_BKG_FIX + // 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); +#endif for (int y = 0; y < numRows; y++) { psVectorInit(clipMask, 0); @@ -105,4 +125,8 @@ // Not enough points to fit patternMaskRow(ro, y, maskBad); +#ifdef PATTERN_ROW_BKG_FIX + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif continue; } @@ -111,8 +135,16 @@ psErrorClear(); patternMaskRow(ro, y, maskBad); - continue; - } - - poly->coeff[0] -= background; +#ifdef PATTERN_ROW_BKG_FIX + // Ignore this row in our subsequent fits, because the fit failed. + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif + continue; + } +#ifndef PATTERN_ROW_BKG_FIX + poly->coeff[0] -= background; +#else + // Store the results we found for this row. + yaxisData->data.F32[y] = poly->coeff[0]; +#endif memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64)); psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector @@ -121,4 +153,7 @@ psErrorClear(); patternMaskRow(ro, y, maskBad); +#ifdef PATTERN_ROW_BKG_FIX + yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF; +#endif continue; } @@ -130,4 +165,87 @@ } +#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 + 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, 2); // 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 { + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + image->data.F32[y][x] += 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, 2); // 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 { + 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]; + corr->data.F64[y][0] -= solution->data.F32[y] * indices->data.F32[x]; + } + } + } + psFree(solution); + } + psFree(yaxisPoly); + psFree(xaxisPoly); + psFree(yaxisIndices); + psFree(yaxisMask); + 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, "Pattern row correction", corr);