Index: trunk/psModules/src/detrend/pmPattern.c
===================================================================
--- trunk/psModules/src/detrend/pmPattern.c	(revision 32970)
+++ trunk/psModules/src/detrend/pmPattern.c	(revision 33243)
@@ -146,4 +146,10 @@
 	// 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));
@@ -161,4 +167,5 @@
         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);
@@ -177,5 +184,5 @@
     // 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.
+    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");
@@ -209,4 +216,5 @@
 	    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]);
 	  }
 	}
@@ -218,5 +226,5 @@
     // 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.
+    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");
@@ -233,5 +241,6 @@
 	  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];
+	    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]);
 	  }
 	}
@@ -500,2 +509,346 @@
 
 
+
+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
+    
+    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); }
+	  }
+	  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.F32[y][x],B->data.F32[y][x],C->data.F32[y][x],D->data.F32[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;
+      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
+	psTrace("psModules.detrend.cont",5,"CmD %d %d %d %d %g %g", // diagnostic
+		i,x,y,j,
+		C->data.F64[y][x],
+		D->data.F64[y][x+1]
+		);
+	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:
+	  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;
+	psTrace("psModules.detrend.cont",5,"DmC %d %d %d %d %g %g",
+		i,x,y,j,
+		D->data.F64[y][x],
+		C->data.F64[y][x-1]
+		);
+
+	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]))) {
+	  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 (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
+	    (isfinite(A->data.F64[y][x]))&&(isfinite(B->data.F64[y+1][x]))) {
+	  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 (!(meanMask->data.PS_TYPE_VECTOR_MASK_DATA[j] & PM_PATTERN_IGNORE)&&
+	    (isfinite(B->data.F64[y][x]))&&(isfinite(A->data.F64[y-1][x]))) {
+	  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;
+	}
+      }
+    }
+    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 (0)
+    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);
+
+    /* 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;
+}
+
+