Index: trunk/psModules/src/detrend/pmPattern.c
===================================================================
--- trunk/psModules/src/detrend/pmPattern.c	(revision 35038)
+++ trunk/psModules/src/detrend/pmPattern.c	(revision 35081)
@@ -938,3 +938,366 @@
 }
 
-
+// This should reuse code more efficiently
+bool pmPatternContinuityBackground(pmFPAfile *in, pmFPAfile *out,  psStatsOptions bgStat, psStatsOptions cellStat,
+				   psImageMaskType maskVal, psImageMaskType maskBad, int edgeWidth)
+{
+    PS_ASSERT_PTR_NON_NULL(in, false);
+    PS_ASSERT_PTR_NON_NULL(out, false);
+
+    int numChips = out->fpa->chips->n;            // Number of cells
+
+    psVector *meanMask = psVectorAlloc(numChips, 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
+    };
+
+    // Measure mean of each cell edge, and use that to determine the cell offsets.
+    psStatsOptions stat = cellStat;     // Define which statistic to use.
+    
+    psStats *bgStats = psStatsAlloc(stat); // 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);
+
+    // Corners don't exist.
+    A->data.F64[0][0] = NAN;
+    B->data.F64[0][0] = NAN;
+    C->data.F64[0][0] = NAN;
+    D->data.F64[0][0] = NAN;
+    A->data.F64[7][0] = NAN;
+    B->data.F64[7][0] = NAN;
+    C->data.F64[7][0] = NAN;
+    D->data.F64[7][0] = NAN;
+    A->data.F64[0][7] = NAN;
+    B->data.F64[0][7] = NAN;
+    C->data.F64[0][7] = NAN;
+    D->data.F64[0][7] = NAN;
+    A->data.F64[7][7] = NAN;
+    B->data.F64[7][7] = NAN;
+    C->data.F64[7][7] = NAN;
+    D->data.F64[7][7] = NAN;
+    
+    for (int i = 0; i < numChips; i++) {
+      pmChip *chip = out->fpa->chips->data[i];
+      pmCell *cell = chip->cells->data[0]; // 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
+      // This is wrong for chips
+      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
+      int xParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.XPARITY");
+      int yParity          = psMetadataLookupS16(NULL,chip->concepts, "CHIP.YPARITY");
+      int x = chipName[2] - '0';
+      int y = chipName[3] - '0';
+
+      
+      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);
+
+	if (!psImageBackground(bgStats, NULL, subset, NULL, 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);
+	  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,stat) < 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 (xParity == -1) {
+	    if (j == 2)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	    else if (j == 3) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	  }
+	  if (xParity == 1) {
+	    if (j == 3)      { C->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	    else if (j == 2) { D->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	  }
+	  if (yParity == 1) {
+	    if (j == 0)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	    else if (j == 1) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	  }
+	  if (yParity == -1) {
+	    if (j == 1)      { B->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	    else if (j == 0) { A->data.F64[y][x] = psStatsGetValue(bgStats,stat); }
+	  }
+	  
+	}
+#if (0)	
+	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],0);
+	  }
+	}
+#endif	
+	psFree(subset);
+	
+      }
+      psTrace("psModules.detrend.cont",5, "OTA: %d (%d %d) (%d %d) A: %f B: %f C: %f D: %f",
+	      i,x,y,xParity,yParity,
+	      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 < numChips; 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.
+      pmChip *chip = out->fpa->chips->data[i];
+      const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
+      int x = chipName[2] - '0';
+      int y = chipName[3] - '0';
+
+      int j;
+      int k = 8 * x + y;
+      double critical_value = 0.0;
+
+      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 (// 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[k] += C->data.F64[y][x] - D->data.F64[y][x+1];     // Take the difference
+	  XX->data.F64[k][k] += 1;                                              // increment our relation with ourself
+	  XX->data.F64[k][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 (
+	    (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[k] += D->data.F64[y][x] - C->data.F64[y][x-1];
+	  XX->data.F64[k][k] += 1;
+	  XX->data.F64[k][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 (
+	    (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[k] += A->data.F64[y][x] - B->data.F64[y+1][x];
+	  XX->data.F64[k][k] += 1;
+	  XX->data.F64[k][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 (
+	    (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[k] += B->data.F64[y][x] - A->data.F64[y-1][x];
+	  XX->data.F64[k][k] += 1;
+	  XX->data.F64[k][j] += -1;
+	}
+      }
+    }
+    double max_XX = 0;
+    double solution_V = 0;
+    int i_peak = -1;
+    for (int i = 0; i < numChips + 4; 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);
+    
+#if (1)
+    for (int i = 0; i < numChips + 4; 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 < numChips + 4; 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 < numChips + 4; i++) { // print vector b
+      psTrace("psModules.detrend.cont",5,"x: %d %f",
+	      i,
+	      solution->data.F64[i]
+	      );
+    }
+#endif
+    
+    // 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 < numChips; i++) {
+	pmChip *chip = out->fpa->chips->data[i];
+        pmCell *cell = chip->cells->data[0]; // Cell of interest
+        pmReadout *ro = cell->readouts->data[0]; // Readout of interest
+	const char *chipName = psMetadataLookupStr(NULL,chip->concepts, "CHIP.NAME");
+	int x = chipName[2] - '0';
+	int y = chipName[3] - '0';
+	int k = 8 * x + y;
+        float correction = solution->data.F64[k];
+        psLogMsg("psModules.detrend", PS_LOG_DETAIL, "Correcting background of chip %s by %f",
+                 chipName, 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 chip correction solution", correction);
+    }
+
+    psFree(solution);
+    psFree(meanMask);
+
+    return true;
+}
Index: trunk/psModules/src/detrend/pmPattern.h
===================================================================
--- trunk/psModules/src/detrend/pmPattern.h	(revision 35038)
+++ trunk/psModules/src/detrend/pmPattern.h	(revision 35081)
@@ -14,4 +14,10 @@
 #include <pslib.h>
 #include <pmFPA.h>
+#include <pmFPAview.h>
+#include <pmFPAfile.h>
+#include <pmFPAfileFitsIO.h>
+#include <pmFPAHeader.h>
+#include <pmHDU.h>
+#include <pmHDUUtils.h>
 
 /// @addtogroup detrend Detrend Creation and Application
@@ -65,4 +71,12 @@
     );
 
+bool pmPatternContinuityBackground(pmFPAfile *in,
+				   pmFPAfile *out,
+				   psStatsOptions bgStat,
+				   psStatsOptions cellStat,
+				   psImageMaskType maskVal,
+				   psImageMaskType maskBad,
+				   int edgeWidth);
+
 /// Apply previously measured cell pattern correction
 bool pmPatternContinuityApply(pmReadout *ro,          ///< Readout to correct
Index: trunk/psModules/src/imcombine/pmStack.c
===================================================================
--- trunk/psModules/src/imcombine/pmStack.c	(revision 35038)
+++ trunk/psModules/src/imcombine/pmStack.c	(revision 35081)
@@ -1385,4 +1385,5 @@
 
   psVector *pixelData = psVectorAlloc(input->n,PS_TYPE_F32);
+  psVector *pixelMask = psVectorAlloc(input->n,PS_TYPE_VECTOR_MASK);
   psStats  *stats     = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
 
@@ -1392,10 +1393,25 @@
 	pmReadout *ro  = stack->data[i];
 	psImage *image = ro->image;
+	pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
 	pixelData->data.F32[i] = image->data.F32[y][x];
+	if (isfinite(image->data.F32[y][x])&&
+	    (fabs(image->data.F32[y][x]) < 1e5)) {
+	  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+	}
+	else {
+	  pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
+	}
+#if (0)
+      if ((x == 59)&&(y > 40)&&(y < 50)) {
+	  fprintf(stderr,"%d %d %d %d %g\n",
+		  x,y,i,pixelMask->data.PS_TYPE_VECTOR_MASK_DATA[i],pixelData->data.F32[i]);
+	}
       }
-      if (!psVectorStats(stats,pixelData,NULL,NULL,0)) {
+#endif
+      if (!psVectorStats(stats,pixelData,NULL,pixelMask,1)) {
 	psError(PS_ERR_UNKNOWN, false, "Unable to calculate median");
 	psFree(stats);
 	psFree(pixelData);
+	psFree(pixelMask);
 	psFree(stack);
 	return(false);
@@ -1403,4 +1419,11 @@
       combined->image->data.F32[y][x] = stats->robustMedian;
       combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0;
+#if (0)
+      if ((x == 59)&&(y > 40)&&(y < 50)) {
+	fprintf(stderr,"%d %d %d %d %g\n",
+		x,y,-1,combined->mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x],
+		combined->image->data.F32[y][x]);
+      }
+#endif
     }
   }
@@ -1408,4 +1431,5 @@
   psFree(stats);
   psFree(pixelData);
+  psFree(pixelMask);
   psFree(stack);
   return (true);