Changeset 33415 for branches/meh_branches/ppstack_test/psModules/src/detrend/pmPattern.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psModules (modified) (1 prop)
• psModules/src/detrend/pmPattern.c (modified) (7 diffs)

branches/meh_branches/ppstack_test/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];
+        xaxisData->data.F32[y] = poly->coeff[1];
+        psTrace("pattern",1,"%d %g %g\n",y,poly->coeff[0],poly->coeff[1]);
+        
+        //      yaxisData->data.F32[y] = 0.0;
+/*      xaxisData->data.F32[y] = 0.0; */
+        
+#endif
         memcpy(corr->data.F64[y], poly->coeff, (order + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
         psVector *solution = psPolynomial1DEvalVector(poly, indices); // Solution vector
@@ -121 +159 @@
             psErrorClear();
             patternMaskRow(ro, y, maskBad);
+#ifdef PATTERN_ROW_BKG_FIX
+            yaxisMask->data.PS_TYPE_VECTOR_MASK_DATA[y] = 0xFF;
+#endif
             continue;
         }
@@ -126 +167 @@
         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);
     }
 
+#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, 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 {
+        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];
+            psTrace("pattern",5,"B: %d %d %g\n",x,y,solution->data.F32[x]);
+          }
+        }
+      }
+      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 {
+        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][1]  -= solution->data.F32[y] ;
+            psTrace("pattern",5,"C: %d %d %g %g\n",x,y,solution->data.F32[x],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);
@@ -382 +509 @@
 
 
+
+bool pmPatternContinuity(pmChip *chip, const psVector *tweak, psStatsOptions bgStat, psStatsOptions cellStat,
+                         psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth)
+{
+    PS_ASSERT_PTR_NON_NULL(chip, false);
+    PS_ASSERT_VECTOR_NON_NULL(tweak, false);
+    PS_ASSERT_VECTOR_SIZE(tweak, chip->cells->n, false);
+    PS_ASSERT_VECTOR_TYPE(tweak, PS_TYPE_U8, false);
+
+    int numCells = tweak->n;            // Number of cells
+
+    psVector *meanMask = psVectorAlloc(numCells, PS_TYPE_VECTOR_MASK); // Mask for means
+    psVectorInit(meanMask, 0);
+
+    // Mask bits
+    enum {
+        PM_PATTERN_IGNORE = 0x01,       // Ignore this cell
+        PM_PATTERN_TWEAK  = 0x02,       // Tweak this cell
+        PM_PATTERN_ERROR  = 0x04,       // Error in calculating background
+        PM_PATTERN_ALL    = 0xFF,       // All causes
+    };
+
+    // Count number of cells to tweak
+    int numTweak = 0;                   // Number of cells to tweak
+    int numIgnore = 0;                  // Number of cells to ignore
+    for (int i = 0; i < numCells; i++) {
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+        if (!cell || !cell->data_exists || !cell->process ||
+            cell->readouts->n == 0 || cell->readouts->n > 1 || !cell->readouts->data[0]) {
+            numIgnore++;
+            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_IGNORE;
+            continue;
+        }
+        if (tweak->data.U8[i]) {
+            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PM_PATTERN_TWEAK;
+            numTweak++;
+        }
+    }
+    if (numTweak == 0) {
+        // Nothing to do
+        psFree(meanMask);
+        return true;
+    }
+
+    // Measure mean of each cell edge, and use that to determine the cell offsets.
+
+    psStats *bgStats = psStatsAlloc(bgStat); // Statistics on background
+    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator
+
+    psRegion region = {0,0,0,0};
+
+    /* These images hold the edge data for the OTA structure.  */
+    psImage *A = psImageAlloc(8,8,PS_TYPE_F64); // Top edge
+    psImage *B = psImageAlloc(8,8,PS_TYPE_F64); // Bottom edge
+    psImage *C = psImageAlloc(8,8,PS_TYPE_F64); // Right edge
+    psImage *D = psImageAlloc(8,8,PS_TYPE_F64); // Left edge
+    psImageInit(A,0.0);
+    psImageInit(B,0.0);
+    psImageInit(C,0.0);
+    psImageInit(D,0.0);
+    
+    for (int i = 0; i < numCells; i++) {
+        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
+            continue;
+        }
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+
+        psStatsInit(bgStats);
+
+        // Convert cell iterator i into an xy coordinate on the grid of cells
+        int y = (i % 8);
+        int x = (i - y) / 8;
+        
+        for (int j = 0; j < 4; j++) {
+          if (j == 0) {  // Region B
+            region = psRegionSet(0,ro->image->numCols,
+                                 0,edgeWidth);
+          }
+          else if (j == 1) { // Region A
+            region = psRegionSet(0,ro->image->numCols,
+                                 ro->image->numRows - edgeWidth,ro->image->numRows);
+          }
+          else if (j == 2) { // Region D
+            region = psRegionSet(0,edgeWidth,
+                                 0,ro->image->numRows);
+          }
+          else if (j == 3) { // Region C
+            region = psRegionSet(ro->image->numCols - edgeWidth,ro->image->numCols,
+                                 0,ro->image->numRows);
+          }
+          psImage *subset  = psImageSubset(ro->image,region);
+          psImage *submask = psImageSubset(ro->mask,region);
+
+          if (!psImageBackground(bgStats, NULL, subset, submask, maskVal, rng)) {
+            psWarning("Unable to measure background for cell %d on edge %d\n", i, j);
+            psErrorClear();
+            meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PM_PATTERN_ERROR;
+            if (j == 0)      { B->data.F64[y][x] = NAN; }
+            else if (j == 1) { A->data.F64[y][x] = NAN; }
+            else if (j == 2) { C->data.F64[y][x] = NAN; }
+            else if (j == 3) { D->data.F64[y][x] = NAN; }
+            psFree(subset);
+            psFree(submask);
+            continue; // Move on to next edge, as only part of this cell may be a problem
+          }
+ 
+          // If the returned value is zero, assume something is wrong.  Do I still need this?
+          if (psStatsGetValue(bgStats,bgStat) < 1e-6) {
+            if (j == 0)      { B->data.F64[y][x] = NAN; }
+            else if (j == 1) { A->data.F64[y][x] = NAN; }
+            else if (j == 2) { C->data.F64[y][x] = NAN; }
+            else if (j == 3) { D->data.F64[y][x] = NAN; }
+          }
+          // If we have an error for this cell/edge, make sure we mask the value
+          if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_ERROR) {
+            if (j == 0)      { B->data.F64[y][x] = NAN; }
+            else if (j == 1) { A->data.F64[y][x] = NAN; }
+            else if (j == 2) { C->data.F64[y][x] = NAN; }
+            else if (j == 3) { D->data.F64[y][x] = NAN; }
+          }
+          else { // Set the value to match what we got from the edge box.
+            if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
+            else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
+            else if (j == 2) { C->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
+            else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,bgStat); }
+          }
+
+          for (int u = 0; u < subset->numCols; u++) {
+            for (int v = 0; v < subset->numRows; v++) {
+              psTrace("psModules.detrend.cont",10,"BOX: %d %d (%d %d) (%d %d) %f %d",
+                      i,j,x,y,u,v,subset->data.F32[v][u],submask->data.PS_TYPE_IMAGE_MASK_DATA[v][u]);
+            }
+          }       
+          
+          psFree(subset);
+          psFree(submask);
+
+        }
+        psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) A: %f B: %f C: %f D: %f",
+                i,x,y,
+                A->data.F64[y][x],B->data.F64[y][x],C->data.F64[y][x],D->data.F64[y][x]);               
+    }
+    psFree(bgStats);
+    psFree(rng);
+
+    // We've now allocated all the edge values, so we can now minimize the offsets.
+    // This involves solving the equation A x = b, where
+    // A is the (64x64 for GPC1) matrix containing the edges that match for each cell
+    // x is the solution vector
+    // b is the combination of offsets across each cell boundary for each cell.
+    // Below "XX" is used as the matrix A, and "solution" is used as both b and x
+    //   (due to the way psMatrixLUSolve operates).
+    psVector *solution = psVectorAlloc(64,PS_TYPE_F64);
+    psImage  *XX       = psImageAlloc(64,64,PS_TYPE_F64);
+    psVectorInit(solution,0.0);
+    psImageInit(XX,0.0);
+    
+    for (int i = 0; i < numCells; i++) {
+      // Accumulate all the possible edge differences we can for this cell.
+      // As we do so, make a note of the correlations by incrementing the element of the matrix.
+      int y = (i % 8);
+      int x = (i - y) / 8;
+      int j;
+      double critical_value = 0.0;
+      if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
+        continue;
+      }
+      if (x + 1 < 8) {  // We have a neighbor adjacent in the +x direction
+        j = 8 * (x + 1) + y; // Determine that neighbor's index
+        if (fabs(C->data.F64[y][x]) > fabs(D->data.F64[y][x+1])) {
+          critical_value = 2.0 * fabs(D->data.F64[y][x+1]);
+        }
+        else {
+          critical_value = 2.0 * fabs(C->data.F64[y][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g %g", // diagnostic
+                i,x,y,j,
+                C->data.F64[y][x],
+                D->data.F64[y][x+1],
+                critical_value
+                );
+        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&  // If there are no errors with the neighbor,
+            (isfinite(C->data.F64[y][x]))&&(isfinite(D->data.F64[y][x+1]))&&     // and all edges have valid values,
+            (fabs(C->data.F64[y][x] - D->data.F64[y][x+1]) < critical_value)     // and there are no large discontinuities,
+            ) {    
+          solution->data.F64[i] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
+          XX->data.F64[i][i] += 1;                                              // increment our relation with ourself
+          XX->data.F64[i][j] += -1;                                             // decrement our relation with the neighbor
+        }
+      }
+      if (x - 1 > -1) { // etc.
+        j = 8 * (x - 1) + y;
+        if (fabs(C->data.F64[y][x-1]) > fabs(D->data.F64[y][x])) {
+          critical_value = 2.0 * fabs(D->data.F64[y][x]);
+        }
+        else {
+          critical_value = 2.0 * fabs(C->data.F64[y][x-1]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g %g",
+                i,x,y,j,
+                D->data.F64[y][x],
+                C->data.F64[y][x-1],
+                critical_value
+                );
+
+        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
+            (isfinite(D->data.F64[y][x]))&&(isfinite(C->data.F64[y][x-1]))&&
+            (fabs(D->data.F64[y][x] - C->data.F64[y][x-1]) < critical_value)
+            ) {
+          solution->data.F64[i] += D->data.F64[y][x] - C->data.F64[y][x-1];
+          XX->data.F64[i][i] += 1;
+          XX->data.F64[i][j] += -1;
+        }
+      }
+      if (y + 1 < 8) {
+        j = 8 * x + (y + 1);
+        psTrace("psModules.detrend.cont",5,"AmB %d %d %d %d %g %g",
+                i,x,y,j,
+                A->data.F64[y][x],
+                B->data.F64[y+1][x]
+                );
+        if (fabs(A->data.F64[y][x]) > fabs(B->data.F64[y+1][x])) {
+          critical_value = 2.0 * fabs(B->data.F64[y+1][x]);
+        }
+        else {
+          critical_value = 2.0 * fabs(A->data.F64[y][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
+            (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))&&
+            (fabs(A->data.F64[y][x] - B->data.F64[y+1][x]) < critical_value)
+            ) {
+          solution->data.F64[i] += A->data.F64[y][x] - B->data.F64[y+1][x];
+          XX->data.F64[i][i] += 1;
+          XX->data.F64[i][j] += -1;
+        }
+      }
+      if (y - 1 > -1) {
+        j = 8 * x +  (y - 1);
+        psTrace("psModules.detrend.cont",5,"BmA %d %d %d %d %g %g",
+                i,x,y,j,
+                B->data.F64[y][x],
+                A->data.F64[y-1][x]
+                );
+        if (fabs(A->data.F64[y-1][x]) > fabs(B->data.F64[y][x])) {
+          critical_value = 2.0 * fabs(B->data.F64[y][x]);
+        }
+        else {
+          critical_value = 2.0 * fabs(A->data.F64[y-1][x]);
+        }
+        if (critical_value < 25) { critical_value = 25; }
+        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
+            (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))&&
+            (fabs(B->data.F64[y][x] - A->data.F64[y-1][x]) < critical_value)
+            ) {
+          solution->data.F64[i] += B->data.F64[y][x] - A->data.F64[y-1][x];
+          XX->data.F64[i][i] += 1;
+          XX->data.F64[i][j] += -1;
+        }
+      }
+    }
+    double max_XX = 0;
+    double solution_V = 0;
+    int i_peak = -1;
+    for (int i = 0; i < numCells; i++) { // If any cells have no value of themself, set the matrix to 1.0.
+      if (XX->data.F64[i][i] == 0.0) {
+        XX->data.F64[i][i] = 1.0;
+      }
+      if (XX->data.F64[i][i] > max_XX) {
+        max_XX = XX->data.F64[i][i];
+        solution_V = solution->data.F64[i];
+        i_peak = i;
+      }
+    }
+    psTrace("psModules.detrend.cont",5,"fixed point: %d %g\n",
+            i_peak,solution_V);
+
+    for (int i = 0; i < numCells; i++) {
+/*        if (!((XX->data.F64[i][i] == 1.0)&& */
+/*          (solution->data.F64[i] == 0.0))) { */
+        solution->data.F64[i] -= solution_V;
+        if (i != i_peak) {
+          for (int j = 0; j < numCells; j++) {
+            XX->data.F64[i][j] -= XX->data.F64[i_peak][j];
+          }
+        }
+/*        } */
+    }
+    for (int i = 0; i < numCells; i++) {
+      XX->data.F64[i_peak][i] = 0.0;
+    }
+    XX->data.F64[i_peak][i_peak] = 1.0;
+    
+    
+#if (1)
+    for (int i = 0; i < numCells; i++) { // print matrix A
+      psTrace("psModules.detrend.cont",5,"A: %3d % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f % 2.0f",
+              i,
+              XX->data.F64[i][0],             XX->data.F64[i][1],             XX->data.F64[i][2],             XX->data.F64[i][3],
+              XX->data.F64[i][4],             XX->data.F64[i][5],             XX->data.F64[i][6],             XX->data.F64[i][7],
+              XX->data.F64[i][8],             XX->data.F64[i][9],             XX->data.F64[i][10],            XX->data.F64[i][11],
+              XX->data.F64[i][12],            XX->data.F64[i][13],            XX->data.F64[i][14],            XX->data.F64[i][15],
+              XX->data.F64[i][16],            XX->data.F64[i][17],            XX->data.F64[i][18],            XX->data.F64[i][19],
+              XX->data.F64[i][20],            XX->data.F64[i][21],            XX->data.F64[i][22],            XX->data.F64[i][23],
+              XX->data.F64[i][24],            XX->data.F64[i][25],            XX->data.F64[i][26],            XX->data.F64[i][27],
+              XX->data.F64[i][28],            XX->data.F64[i][29],            XX->data.F64[i][30],            XX->data.F64[i][31],
+              XX->data.F64[i][32],            XX->data.F64[i][33],            XX->data.F64[i][34],            XX->data.F64[i][35],
+              XX->data.F64[i][36],            XX->data.F64[i][37],            XX->data.F64[i][38],            XX->data.F64[i][39],
+              XX->data.F64[i][40],            XX->data.F64[i][41],            XX->data.F64[i][42],            XX->data.F64[i][43],
+              XX->data.F64[i][44],            XX->data.F64[i][45],            XX->data.F64[i][46],            XX->data.F64[i][47],
+              XX->data.F64[i][48],            XX->data.F64[i][49],            XX->data.F64[i][50],            XX->data.F64[i][51],
+              XX->data.F64[i][52],            XX->data.F64[i][53],            XX->data.F64[i][54],            XX->data.F64[i][55],
+              XX->data.F64[i][56],            XX->data.F64[i][57],            XX->data.F64[i][58],            XX->data.F64[i][59],
+              XX->data.F64[i][60],            XX->data.F64[i][61],            XX->data.F64[i][62],            XX->data.F64[i][63]
+              );
+    }
+
+    for (int i = 0; i < numCells; i++) { // print vector b
+      psTrace("psModules.detrend.cont",5,"b: %d %f",
+              i,
+              solution->data.F64[i]
+              );
+    }
+#endif    
+    
+    // Solve the Ax=b equation
+    //    psMatrixLUSolve(XX,solution);
+    psMatrixGJSolve(XX,solution);
+#if (1)
+    for (int i = 0; i < numCells; i++) { // print vector b
+      psTrace("psModules.detrend.cont",5,"x: %d %f",
+              i,
+              solution->data.F64[i]
+              );
+    }
+#endif
+    
+    /* old code to remove the minimum solution value from the set, to give a "minimal set of offsets." Mathematically unnecessary. */
+/*     double min = 99e99; */
+/*     for (int i = 0; i < numCells; i++) { */
+/*       if (solution->data.F64[i] < min) { */
+/*      min = solution->data.F64[i]; */
+/*       } */
+/*       psTrace("psModules.detrend.cont",5,"x: %d %f %f ", */
+/*            i, */
+/*            solution->data.F64[i],min */
+/*            ); */
+/*     } */
+/*     for (int i = 0; i < numCells; i++) { */
+/*      if (solution->data.F64[i] != 0.0) { */
+/*        solution->data.F64[i] -= min; */
+/*      } */
+/*     } */
+
+    // Cleanup
+    psFree(XX);
+    psFree(A);
+    psFree(B);
+    psFree(C);
+    psFree(D);
+
+    // Correct cells based on the offsets calculated, and store the result in the analysis metadata.
+    for (int i = 0; i < numCells; i++) {
+        if (meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_IGNORE) {
+            continue;
+        }
+        if (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PM_PATTERN_TWEAK)) {
+            continue;
+        }
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+
+        float correction = solution->data.F64[i];
+        const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME"); // Name of cell
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of cell %s by %f",
+                 cellName, correction);
+        psBinaryOp(ro->image, ro->image, "-", psScalarAlloc(correction, PS_TYPE_F32));
+        psMetadataAddF32(ro->analysis, PS_LIST_TAIL, PM_PATTERN_CELL_CORRECTION, PS_META_REPLACE,
+                         "Pattern cell correction solution", correction);
+    }
+
+    psFree(solution);
+    psFree(meanMask);
+
+    return true;
+}
+
+bool pmPatternContinuityApply(pmReadout *ro, psImageMaskType maskBad)
+{
+    PM_ASSERT_READOUT_NON_NULL(ro, false);
+    PM_ASSERT_READOUT_IMAGE(ro, false);
+
+    bool mdok;                          // Status of MD lookup
+    float corr = psMetadataLookupF32(&mdok, ro->analysis, PM_PATTERN_CELL_CORRECTION); // Correction to apply
+    if (!mdok) {
+        // No correction to apply
+        return true;
+    }
+
+    psImage *image = ro->image, *mask = ro->mask; // Image and mask of interest
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    if (!isfinite(corr)) {
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                image->data.F32[y][x] = NAN;
+            }
+        }
+        if (mask) {
+            for (int y = 0; y < numRows; y++) {
+                for (int x = 0; x < numCols; x++) {
+                    mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskBad;
+                }
+            }
+        }
+    } else {
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                image->data.F32[y][x] += corr;
+            }
+        }
+    }
+
+    return true;
+}
+
+