Index: /branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c
===================================================================
--- /branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c	(revision 42652)
+++ /branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c	(revision 42653)
@@ -222,4 +222,6 @@
 }
 
+static int saveVectors = 1;
+
 bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) {
 
@@ -339,6 +341,6 @@
     // We are performing a row-by-row overscan subtraction
     int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction
-    if ((cellreaddir != 1) && (cellreaddir != 2)) {
-	psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n");
+    if ((cellreaddir != -1) && (cellreaddir != 1) && (cellreaddir != 2)) {
+	psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols) or 0 (2D)\n");
 	return false;
     }
@@ -369,4 +371,5 @@
 		values = psVectorRealloc(values, values->n + overscan->numCols);
 		values->n += overscan->numCols;
+		// XXX double-check the range of values->n here
 		memcpy(&values->data.F32[index], overscan->data.F32[j],
 		       overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
@@ -545,4 +548,176 @@
     }
 
+    // 2D bias subtraction with x-dir readout direction: the
+    // bias is constructed by combining a 1D pattern in the
+    // readout direction from the top overscan region and a second
+    // 1D pattern in the cross direction from the overscan
+    if (cellreaddir == -1) {
+	// we require 2 overscan regions, and they must match these directions:
+	if (overscans->n != 2) {
+	    psLogMsg (__func__, PS_LOG_ERROR, "CELL.READDIR = -1, requires 2 overscan regions but %d supplied", (int) overscans->n);
+	    psLogMsg (__func__, PS_LOG_ERROR, "e.g.: CELL.BIASSEC STR [591:624,1:608],[1:590,598:608]");
+	    psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to generate overscan vector.\n");
+	    psFree(stats);
+	    return false;
+	}
+
+	// the serial (fast readout) direction is columns (x-dir)
+	psImage *yscan = psListGet (overscans, 0); // overscan region spanning all rows
+	psImage *xscan = psListGet (overscans, 1); // overscan region spanning all columns
+
+	// Extract the y-dir overscan vector.  The overscan and image might not be aligned:
+	// diff represents the offset between the rows in the image data and the overscan.
+	// pixels->data represents the image row pixels. For example, the image region may be
+	// inset in the y-direction but the overscan could cover the entire y-range
+
+	// The read direction is rows
+	psArray *yscanPixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
+	for (int i = 0; i < yscanPixels->n; i++) {
+	    yscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32);
+	}
+
+	// XXX this code allows multiple yscans to be appended, but this does not
+	// match the concept of how they are assigned above: biassec[0] = yscan
+	int yDiff = yscan->row0 - image->row0; // Offset between the two regions
+	for (int i = PS_MAX(0,yDiff); i < PS_MIN(image->numRows, yscan->numRows + yDiff); i++) {
+	    int j = i - yDiff;
+	    // i is row on image, j is row on yscan
+	    psVector *values = yscanPixels->data[i];
+	    int index = values->n; // Index in the vector
+	    values = psVectorRealloc(values, values->n + yscan->numCols);
+	    values->n += yscan->numCols;
+	    // XXX double-check the range of values->n here
+	    memcpy(&values->data.F32[index], yscan->data.F32[j],
+		   yscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+	    index += yscan->numCols;
+	    yscanPixels->data[i] = values; // Update the pointer in case it's moved
+	}
+
+	// Extract the x-dir overscan vector.  The overscan and image might not be aligned:
+	// diff represents the offset between the rows in the image data and the overscan.
+	// pixels->data represents the image row pixels. For example, the image region may be
+	// inset in the x-direction but the overscan could cover the entire y-range
+
+	// Extract the top region as a vector of the columns
+	psArray *xscanPixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
+	for (int i = 0; i < xscanPixels->n; i++) {
+	    xscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32);
+	}
+
+	int xDiff = xscan->col0 - image->col0; // Offset between the two regions
+	for (int ix = PS_MAX(0,xDiff); ix < PS_MIN(image->numCols, xscan->numCols + xDiff); ix++) {
+	    int jx = ix - xDiff;
+	    // ix is row on image, jx is column on xscan
+	    psVector *values = xscanPixels->data[ix];
+	    values = psVectorRealloc(values, xscan->numRows);
+	    values->n = xscan->numRows;
+	    for (int iy = 0; iy < xscan->numRows; iy++) {
+		values->data.F32[iy] = xscan->data.F32[iy][jx];
+	    }
+	    xscanPixels->data[ix] = values; // Update the pointer in case it's moved
+	}
+
+	// Reduce the overscans
+	// XXX need to save 2 different chi-square values
+	psVector *xReduced = pmOverscanVector(&chi2, overscanOpts, xscanPixels, stats);
+	psFree(xscanPixels);
+	if (!xReduced) {
+	    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n");
+	    psFree(stats);
+	    return false;
+	}
+
+	// Reduce the overscans
+	// XXX need to save 2 different chi-square values
+	psVector *yReduced = pmOverscanVector(&chi2, overscanOpts, yscanPixels, stats);
+	psFree(yscanPixels);
+	if (!yReduced) {
+	    psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n");
+	    psFree(stats);
+	    return false;
+	}
+
+
+	if (saveVectors) {
+	    FILE *foutX = fopen ("xdir.dat", "w");
+	    int fdX = fileno (foutX);
+	    p_psVectorPrint (fdX, xReduced, "xscan");
+	    fclose (foutX);
+	    FILE *foutY = fopen ("ydir.dat", "w");
+	    int fdY = fileno (foutY);
+	    p_psVectorPrint (fdY, yReduced, "yscan");
+	    fclose (foutY);
+	    saveVectors = 0;
+	}
+
+	// write info about stats of y-dir overscan vector in header
+	{ 
+	    psString comment = NULL;    // Comment to add
+	    psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+	    if (!psVectorStats (vectorStats, yReduced, NULL, NULL, 0)) {
+		psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+		return false;
+	    }
+	    psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean);
+	    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+	    psFree(comment);
+
+	    // write metadata header value
+	    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean);
+	    psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev);
+
+	    // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
+	    // declare the readout dead and mask
+	  
+	    if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) {
+		fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean);
+		psImage *mask = input->mask;
+		for (int y = 0; y < mask->numRows; y++) {
+		    for (int x = 0; x < mask->numCols; x++) {
+			mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal;
+		    }
+		}
+	    }
+	    psFree (vectorStats);
+	}
+
+	// XXX apply the 2D bias correction here
+	for (int i = 0; i < image->numRows; i++) {
+	    for (int j = 0; j < image->numCols; j++) {
+		image->data.F32[i][j] -= yReduced->data.F32[i];
+	    }
+	}
+
+	// subtract the y-dir vector from the y-dir overscan regions (why?)
+	{
+	    // the overscan and image might not be aligned.
+	    int diff = yscan->row0 - image->row0; // Offset between the two regions
+	    for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, yscan->numRows + diff); i++) {
+		int j = i - diff;
+		// i is row on image (aligned with yReduced)
+		// j is row on yscan
+		for (int k = 0; k < yscan->numCols; k++) {
+		    yscan->data.F32[j][k] -= yReduced->data.F32[i];
+		}
+	    }
+	}
+
+	// subtract the x-dir vector from the x-dir overscan regions (why?)
+	{
+	    // the overscan and image might not be aligned.
+	    int diff = xscan->col0 - image->col0; // Offset between the two regions
+	    for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, xscan->numCols + diff); i++) {
+		int j = i - diff;
+		// i is column on image (aligned with xReduced)
+		// j is column on xscan
+		for (int k = 0; k < xscan->numRows; k++) {
+		    xscan->data.F32[k][j] -= xReduced->data.F32[i];
+		}
+	    }
+	}
+	psFree(xReduced);
+	psFree(yReduced);
+    } 
+
     pmOverscanUpdateHeader (hdu, overscanOpts, chi2);
     psFree(stats);