Changeset 32970 for trunk/psModules/src/detrend/pmPattern.c

Timestamp:

Dec 19, 2011, 12:50:30 PM (14 years ago)

Author:

watersc1

Message:

Fix for pattern correction that I forgot to commit previously.

File:

• trunk/psModules/src/detrend/pmPattern.c (modified) (6 diffs)

trunk/psModules/src/detrend/pmPattern.c

@@ -6 +6 @@
 
 #include "pmPattern.h"
+
+#define PATTERN_ROW_BKG_FIX 1
 
 
@@ -89 +91 @@
     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 +125 @@
             // 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 +135 @@
             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 +153 @@
             psErrorClear();
             patternMaskRow(ro, y, maskBad);
+#ifdef PATTERN_ROW_BKG_FIX
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+#endif
             continue;
         }
@@ -130 +165 @@
     }
 
+#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);