Changeset 42679 for branches/2dbias/psModules/src/detrend/pmOverscan.c

Timestamp:

May 21, 2024, 8:50:38 PM (2 years ago)

Author:

hgao

Message:

implement 2D overscan subtraction per cell as configured in OVERSCAN.2D.SUBSET

File:

• branches/2dbias/psModules/src/detrend/pmOverscan.c (modified) (5 diffs)

branches/2dbias/psModules/src/detrend/pmOverscan.c

@@ -17 +17 @@
 #include "pmOverscan.h"
 
-#define SMOOTH_NSIGMA 4.0               // Number of Gaussian sigma the smoothing kernel extends
+#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends
 
 static void pmOverscanOptionsFree(pmOverscanOptions *options);
 static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options);
-bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced);
-
-bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) {
-
-    assert (input);
-
-    if (overscanOpts == NULL) return true; // no overscan subtraction requested
-
-    pmHDU *hdu = pmHDUFromReadout(input);  // HDU of interest
-    psImage *image = input->image;
-
-    // check for 'soft bias' (simple, fixed offset to be subtracted)
-    if (overscanOpts->constant) {
-        // write metadata header value
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value);
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
-
-        // NOTE psBinaryOp frees arg2 if it is a scalar
-        (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32));
-
+bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced);
+
+bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan)
+{
+
+        assert(input);
+
+        if (overscanOpts == NULL)
+                return true; // no overscan subtraction requested
+
+        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
+        psImage *image = input->image;
+
+        // check for 'soft bias' (simple, fixed offset to be subtracted)
+        if (overscanOpts->constant)
+        {
+                // write metadata header value
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value);
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
+
+                // NOTE psBinaryOp frees arg2 if it is a scalar
+                (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32));
+
+                return true;
+        }
+
+        // we are performing a statitical analysis of the overscan region
+
+        // Check for an unallowable pmFit.
+        if (overscanOpts->primary)
+        {
+                if (overscanOpts->primary->fitType != PM_FIT_NONE &&
+                        overscanOpts->primary->fitType != PM_FIT_POLY_ORD &&
+                        overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY &&
+                        overscanOpts->primary->fitType != PM_FIT_SPLINE)
+                {
+                        psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n",
+                                        overscanOpts->primary->fitType);
+                        return false;
+                }
+        }
+        if (overscanOpts->secondary)
+        {
+                if (overscanOpts->secondary->fitType != PM_FIT_NONE &&
+                        overscanOpts->secondary->fitType != PM_FIT_POLY_ORD &&
+                        overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY &&
+                        overscanOpts->secondary->fitType != PM_FIT_SPLINE)
+                {
+                        psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d).  Returning original image.\n",
+                                        overscanOpts->secondary->fitType);
+                        return false;
+                }
+        }
+
+        psList *overscans = input->bias; // List of the overscan images
+
+        // Reduce all overscan pixels to a single value
+        if (overscanOpts->single)
+        {
+
+                // extract overscan pixels to a single vector
+                psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL;                                                                                                       // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter)))
+                {
+                        int index = pixels->n; // Index
+                        pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
+                        pixels->n += overscan->numRows * overscan->numCols;
+                        for (int i = 0; i < overscan->numRows; i++)
+                        {
+                                memcpy(&pixels->data.F32[index], overscan->data.F32[i],
+                                           overscan->numCols * sizeof(psF32));
+                                index += overscan->numCols;
+                        }
+                }
+                psFree(iter);
+
+                // statistic to be calculated
+                psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use
+                if (!statistic)
+                {
+                        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n",
+                                        overscanOpts->primary->stat);
+                        return false;
+                }
+                psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
+
+                if (!psVectorStats(stats, pixels, NULL, NULL, 0))
+                {
+                        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+                        return false;
+                }
+                psFree(pixels);
+                double reduced = psStatsGetValue(stats, statistic); // Result of statistics
+
+                psString comment = NULL; // Comment to add
+                psStringAppend(&comment, "Overscan value: %f", reduced);
+                psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+                psFree(comment);
+
+                // write metadata header value
+                // XXX EAM : this could / should write the stdev of the overscan region
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced);
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
+
+                psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32);
+                psBinaryOp(image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
+
+                // subtract the measured value from each overscan region as well
+                iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                overscan = NULL;                                                                                        // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter)))
+                {
+                        psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
+                }
+                psFree(iter);
+                psFree(reducedScalar);
+
+                // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
+                // declare the readout dead and mask
+
+                if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid))
+                {
+                        fprintf(stderr, "bad overscan (1) %f, masking readout\n", reduced);
+                        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(stats);
+                return true;
+        }
+
+        bool mdok = false;
+
+        // 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");
+                return false;
+        }
+
+        float chi2 = NAN; // chi^2 from fit
+
+        // adjust operation depending on the read direction : need to re-org pixels for columns
+        if (!doTwoDOverscan && (cellreaddir == 1))
+        {
+                // The read direction is rows
+                psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++)
+                {
+                        pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32);
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL;                                                                                                       // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter)))
+                {
+                        // the overscan and image might not be aligned.  pixels->data represents
+                        // the image row pixels.
+                        int diff = overscan->row0 - image->row0; // Offset between the two regions
+                        for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++)
+                        {
+                                int j = i - diff;
+                                // i is row on image
+                                // j is row on overscan
+                                psVector *values = pixels->data[i];
+                                int index = values->n; // Index in the vector
+                                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));
+                                index += overscan->numCols;
+                                pixels->data[i] = values; // Update the pointer in case it's moved
+                        }
+                }
+                psFree(iter);
+
+                // Reduce the overscans
+                psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
+                psFree(pixels);
+                if (!reduced)
+                {
+                        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
+                        return false;
+                }
+
+                if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced))
+                {
+                        psError(PS_ERR_UNKNOWN, false, "failure to update header");
+                        return false;
+                }
+
+                // Subtract row by row
+                for (int i = 0; i < image->numRows; i++)
+                {
+                        for (int j = 0; j < image->numCols; j++)
+                        {
+                                image->data.F32[i][j] -= reduced->data.F32[i];
+                        }
+                }
+
+                // subtract from the overscan regions
+                {
+                        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                        psImage *overscan = NULL;                                                                                                       // Overscan image from iterator
+                        while ((overscan = psListGetAndIncrement(iter)))
+                        {
+                                // the overscan and image might not be aligned.
+                                int diff = overscan->row0 - image->row0; // Offset between the two regions
+                                for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++)
+                                {
+                                        int j = i - diff;
+                                        // i is row on image
+                                        // j is row on overscan
+                                        for (int k = 0; k < overscan->numCols; k++)
+                                        {
+                                                overscan->data.F32[j][k] -= reduced->data.F32[j];
+                                        }
+                                }
+                        }
+                        psFree(iter);
+                }
+                psFree(reduced);
+        }
+
+        if (!doTwoDOverscan && (cellreaddir == 2))
+        {
+                // The read direction is columns
+                psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
+                for (int i = 0; i < pixels->n; i++)
+                {
+                        psVector *values = psVectorAlloc(0, PS_TYPE_F32);
+                        pixels->data[i] = values;
+                }
+
+                // Pull the pixels out into the vectors
+                psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                psImage *overscan = NULL;                                                                                                       // Overscan image from iterator
+                while ((overscan = psListGetAndIncrement(iter)))
+                {
+                        // the overscan and image might not be aligned.  pixels->data represents
+                        // the image row pixels.
+                        int diff = overscan->col0 - image->col0; // Offset between the two regions
+                        for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++)
+                        {
+                                int iFixed = i - diff;
+                                // i is column on image
+                                // iFixed is column on overscan
+                                psVector *values = pixels->data[i];
+                                int index = values->n; // Index in the vector
+                                values = psVectorRealloc(values, values->n + overscan->numRows);
+                                for (int j = 0; j < overscan->numRows; j++)
+                                {
+                                        values->data.F32[index++] = overscan->data.F32[j][iFixed];
+                                }
+                                values->n += overscan->numRows;
+                                pixels->data[i] = values; // Update the pointer in case it's moved
+                        }
+                }
+                psFree(iter);
+
+                // Reduce the overscans
+                psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
+                psFree(pixels);
+                if (!reduced)
+                {
+                        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
+                        return false;
+                }
+
+                if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced))
+                {
+                        psError(PS_ERR_UNKNOWN, false, "failure to update header");
+                        return false;
+                }
+
+                // Subtract column by column
+                for (int j = 0; j < image->numRows; j++)
+                {
+                        for (int i = 0; i < image->numCols; i++)
+                        {
+                                image->data.F32[j][i] -= reduced->data.F32[i];
+                        }
+                }
+
+                // subtract from the overscan regions
+                {
+                        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
+                        psImage *overscan = NULL;                                                                                                       // Overscan image from iterator
+                        while ((overscan = psListGetAndIncrement(iter)))
+                        {
+                                // the overscan and image might not be aligned.
+                                int diff = overscan->col0 - image->col0; // Offset between the two regions
+                                for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++)
+                                {
+                                        int j = i - diff;
+                                        // i is col on image
+                                        // j is col on overscan
+                                        for (int k = 0; i < overscan->numRows; k++)
+                                        {
+                                                overscan->data.F32[k][j] -= reduced->data.F32[j];
+                                        }
+                                }
+                        }
+                        psFree(iter);
+                }
+
+                psFree(reduced);
+        }
+
+        // 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 (doTwoDOverscan && (cellreaddir == 1))
+        {
+                psAssert(overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters");
+                // we require 2 overscan regions, and they must match these directions:
+                if (overscans->n != 2)
+                {
+                        psLogMsg(__func__, PS_LOG_ERROR, "OVERSCAN.2D 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:624,599:608]");
+                        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to generate overscan vector.\n");
+                        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(yscan->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 = 0; i < yscanPixels->n; i++)
+                {
+                        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[i],
+                                   yscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+                        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(xscan->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 = 0; ix < xscanPixels->n; ix++)
+                {
+                        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][ix];
+                        }
+                        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 *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels);
+                psFree(yscanPixels);
+                if (!yReduced)
+                {
+                        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n");
+                        return false;
+                }
+                if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, yReduced))
+                {
+                        psError(PS_ERR_UNKNOWN, false, "failure to update header");
+                        return false;
+                }
+
+                // Reduce the overscans
+                // XXX need to save 2 different chi-square values
+                psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels);
+                psFree(xscanPixels);
+                if (!xReduced)
+                {
+                        psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n");
+                        return false;
+                }
+
+                // subtract the 2D bias from the image
+                if (yscan->col0 >= xscan->col0 && yscan->col0 + yscan->numCols <= xscan->col0 + xscan->numCols)
+                {
+                        // define where to normalize the xReduced vector
+                        int j0 = yscan->col0 - xscan->col0 + (int)(yscan->numCols / 2) - 1;
+                        // XXX apply the 2D bias correction here
+                        int yDiff = yscan->row0 - image->row0; // y offset between the science and the yoverscan region
+                        int xDiff = xscan->col0 - image->col0; // x offset between the science and the xoverscan region
+                        for (int i = 0; i < image->numRows; i++)
+                        {
+                                for (int j = 0; j < image->numCols; j++)
+                                {
+                                        int iy = i + yDiff;
+                                        int jx = j + xDiff;
+                                        image->data.F32[i][j] -= yReduced->data.F32[iy] - xReduced->data.F32[j0] + xReduced->data.F32[jx];
+                                }
+                        }
+                }
+                else
+                {
+                        psError(PS_ERR_UNKNOWN, true, "x dimension of yscan is not fully contained by xscan\n");
+                        return false;
+                }
+
+                // 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 = 0; i < yscan->numRows; i++)
+                        {
+                                for (int j = 0; j < yscan->numCols; j++)
+                                {
+                                        yscan->data.F32[i][j] -= 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 = 0; i < xscan->numCols; i++)
+                        {
+                                // int j = i - diff;
+                                // i is column on image (aligned with xReduced)
+                                // j is column on xscan
+                                for (int j = 0; j < xscan->numRows; j++)
+                                {
+                                        xscan->data.F32[j][i] -= xReduced->data.F32[i];
+                                }
+                        }
+                }
+                psFree(xReduced);
+                psFree(yReduced);
+        }
+        // pmOverscanUpdateHeader (hdu, overscanOpts, chi2);
         return true;
-    }
-
-    // we are performing a statitical analysis of the overscan region
-
-    // Check for an unallowable pmFit.
-    if (overscanOpts->primary) {
-        if (overscanOpts->primary->fitType != PM_FIT_NONE &&
-            overscanOpts->primary->fitType != PM_FIT_POLY_ORD &&
-            overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY &&
-            overscanOpts->primary->fitType != PM_FIT_SPLINE) {
-            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n",
-                    overscanOpts->primary->fitType);
-            return false;
-        }
-    }
-    if (overscanOpts->secondary) {
-        if (overscanOpts->secondary->fitType != PM_FIT_NONE &&
-            overscanOpts->secondary->fitType != PM_FIT_POLY_ORD &&
-            overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY &&
-            overscanOpts->secondary->fitType != PM_FIT_SPLINE) {
-            psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d).  Returning original image.\n",
-                    overscanOpts->secondary->fitType);
-            return false;
-        }
-    }
-
-    psList *overscans = input->bias; // List of the overscan images
-
-    // Reduce all overscan pixels to a single value
-    if (overscanOpts->single) {
-
-        // extract overscan pixels to a single vector
-        psVector *pixels = psVectorAlloc(0, PS_TYPE_F32);
-        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-        psImage *overscan = NULL;   // Overscan image from iterator
-        while ((overscan = psListGetAndIncrement(iter))) {
-            int index = pixels->n;  // Index
-            pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols);
-            pixels->n += overscan->numRows * overscan->numCols;
-            for (int i = 0; i < overscan->numRows; i++) {
-                memcpy(&pixels->data.F32[index], overscan->data.F32[i],
-                       overscan->numCols * sizeof(psF32));
-                index += overscan->numCols;
-            }
-        }
-        psFree(iter);
-
-        // statistic to be calculated
-        psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use
-        if (!statistic) {
-            psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n",
-                    overscanOpts->primary->stat);
-            return false;
-        }
-        psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
-
-        if (!psVectorStats(stats, pixels, NULL, NULL, 0)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-            return false;
-        }
-        psFree(pixels);
-        double reduced = psStatsGetValue(stats, statistic); // Result of statistics
-
-        psString comment = NULL;    // Comment to add
-        psStringAppend(&comment, "Overscan value: %f", reduced);
-        psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-        psFree(comment);
-
-        // write metadata header value
-        // XXX EAM : this could / should write the stdev of the overscan region
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced);
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN);
-
-        psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32);
-        psBinaryOp (image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
-
-        // subtract the measured value from each overscan region as well
-        iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-        overscan = NULL;   // Overscan image from iterator
-        while ((overscan = psListGetAndIncrement(iter))) {
-            psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use
-        }
-        psFree(iter);
-        psFree(reducedScalar);
-
-        // EAM 2022.03.29 : if the calculated overscan value is below the threshold,
-        // declare the readout dead and mask
-
-        if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid)) {
-            fprintf (stderr, "bad overscan (1) %f, masking readout\n", reduced);
-            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(stats);
-        return true;
-    } 
-
-    bool mdok = false;
-
-    // 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");
-        return false;
-    }
-
-    float chi2 = NAN;           // chi^2 from fit
-
-    // adjust operation depending on the read direction : need to re-org pixels for columns
-    if (!overscanOpts->TwoD && (cellreaddir == 1)) {
-        // The read direction is rows
-        psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels
-        for (int i = 0; i < pixels->n; i++) {
-            pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32);
-        }
-
-        // Pull the pixels out into the vectors
-        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-        psImage *overscan = NULL; // Overscan image from iterator
-        while ((overscan = psListGetAndIncrement(iter))) {
-            // the overscan and image might not be aligned.  pixels->data represents
-            // the image row pixels.
-            int diff = overscan->row0 - image->row0; // Offset between the two regions
-            for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) {
-                int j = i - diff;
-                // i is row on image
-                // j is row on overscan
-                psVector *values = pixels->data[i];
-                int index = values->n; // Index in the vector
-                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));
-                index += overscan->numCols;
-                pixels->data[i] = values; // Update the pointer in case it's moved
-            }
-        }
-        psFree(iter);
-
-        // Reduce the overscans
-        psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
-        psFree(pixels);
-        if (! reduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
-            return false;
-        }
-
-        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to update header");
-            return false;
-        }
-
-        // Subtract row by row
-        for (int i = 0; i < image->numRows; i++) {
-            for (int j = 0; j < image->numCols; j++) {
-                image->data.F32[i][j] -= reduced->data.F32[i];
-            }
-        }
-
-        // subtract from the overscan regions
-        {
-            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-            psImage *overscan = NULL; // Overscan image from iterator
-            while ((overscan = psListGetAndIncrement(iter))) {
-                // the overscan and image might not be aligned.
-                int diff = overscan->row0 - image->row0; // Offset between the two regions
-                for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) {
-                    int j = i - diff;
-                    // i is row on image
-                    // j is row on overscan
-                    for (int k = 0; k < overscan->numCols; k++) {
-                        overscan->data.F32[j][k] -= reduced->data.F32[j];
-                    }
-                }
-            }
-            psFree(iter);
-        }
-        psFree(reduced);
-    } 
-
-    if (!overscanOpts->TwoD && (cellreaddir == 2)) {
-        // The read direction is columns
-        psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels
-        for (int i = 0; i < pixels->n; i++) {
-            psVector *values = psVectorAlloc(0, PS_TYPE_F32);
-            pixels->data[i] = values;
-        }
-
-        // Pull the pixels out into the vectors
-        psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-        psImage *overscan = NULL; // Overscan image from iterator
-        while ((overscan = psListGetAndIncrement(iter))) {
-            // the overscan and image might not be aligned.  pixels->data represents
-            // the image row pixels.
-            int diff = overscan->col0 - image->col0; // Offset between the two regions
-            for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) {
-                int iFixed = i - diff;
-                // i is column on image
-                // iFixed is column on overscan
-                psVector *values = pixels->data[i];
-                int index = values->n; // Index in the vector
-                values = psVectorRealloc(values, values->n + overscan->numRows);
-                for (int j = 0; j < overscan->numRows; j++) {
-                    values->data.F32[index++] = overscan->data.F32[j][iFixed];
-                }
-                values->n += overscan->numRows;
-                pixels->data[i] = values; // Update the pointer in case it's moved
-            }
-        }
-        psFree(iter);
-
-        // Reduce the overscans
-        psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels);
-        psFree(pixels);
-        if (! reduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
-            return false;
-        }
-
-        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, reduced)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to update header");
-            return false;
-        }
-
-        // Subtract column by column 
-        for (int j = 0; j < image->numRows; j++) {
-            for (int i = 0; i < image->numCols; i++) {
-                image->data.F32[j][i] -= reduced->data.F32[i];
-            }
-        }
-
-        // subtract from the overscan regions
-        {
-            psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator
-            psImage *overscan = NULL; // Overscan image from iterator
-            while ((overscan = psListGetAndIncrement(iter))) {
-                // the overscan and image might not be aligned.
-                int diff = overscan->col0 - image->col0; // Offset between the two regions
-                for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) {
-                    int j = i - diff;
-                    // i is col on image
-                    // j is col on overscan
-                    for (int k = 0; i < overscan->numRows; k++) {
-                        overscan->data.F32[k][j] -= reduced->data.F32[j];
-                    }
-                }
-            }
-            psFree(iter);
-        }
-
-        psFree(reduced);
-    }
-
-    // 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 (overscanOpts->TwoD && (cellreaddir == 1)) {
-        psAssert (overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters");
-        // we require 2 overscan regions, and they must match these directions:
-        if (overscans->n != 2) {
-            psLogMsg (__func__, PS_LOG_ERROR, "OVERSCAN.2D 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");
-            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 *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels);
-        psFree(yscanPixels);
-        if (!yReduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n");
-            return false;
-        }
-        if (!pmOverscanUpdateHeaderVector (input, hdu, overscanOpts, yReduced)) {
-            psError(PS_ERR_UNKNOWN, false, "failure to update header");
-            return false;
-        }
-
-        // Reduce the overscans
-        // XXX need to save 2 different chi-square values
-        psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels);
-        psFree(xscanPixels);
-        if (!xReduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n");
-            return false;
-        }
-
-        // 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);
-    return true;
 
 } // End of overscan subtraction
@@ -437 +500 @@
 static void pmOverscanOptionsFree(pmOverscanOptions *options)
 {
-    psFree(options->primary);
-    psFree(options->secondary);
+        psFree(options->primary);
+        psFree(options->secondary);
 }
 
 static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options)
 {
-    psFree(options->stat);
-    psFree(options->poly);
-    psFree(options->spline);
+        psFree(options->stat);
+        psFree(options->poly);
+        psFree(options->spline);
 }
 
@@ -451 +514 @@
 pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void)
 {
-    pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions));
-    psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree);
-
-    // Inputs
-    opts->fitType = PM_FIT_NONE;
-    opts->order   = 0;
-    opts->stat    = NULL;
-
-    // Smoothing
-    opts->boxcar = 0;
-    opts->gauss  = 0.0;
-
-    // Outputs
-    opts->poly = NULL;
-    opts->spline = NULL;
-
-    return opts;
+        pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions));
+        psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree);
+
+        // Inputs
+        opts->fitType = PM_FIT_NONE;
+        opts->order = 0;
+        opts->stat = NULL;
+
+        // Smoothing
+        opts->boxcar = 0;
+        opts->gauss = 0.0;
+
+        // Outputs
+        opts->poly = NULL;
+        opts->spline = NULL;
+
+        return opts;
 }
 
@@ -473 +536 @@
 pmOverscanOptions *pmOverscanOptionsAlloc(void)
 {
-    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
-    psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree);
-
-    // Inputs
-    opts->single   = false;
-    opts->constant = false;
-    opts->TwoD     = false;
-
-    opts->value    = 0.0;
-
-    opts->minValid = 0.0; // default value if not defined
-    opts->maxValid = (float) 0x10000; // default value if not defined
-    opts->maskVal  = 0x0001; // default value if not defined
-
-    // stat options
-    opts->primary = NULL;
-    opts->secondary = NULL;
-
-    return opts;
+        pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
+        psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree);
+
+        // Inputs
+        opts->single = false;
+        opts->constant = false;
+        opts->TwoD = false;
+
+        opts->value = 0.0;
+
+        opts->minValid = 0.0;                    // default value if not defined
+        opts->maxValid = (float)0x10000; // default value if not defined
+        opts->maskVal = 0x0001;                  // default value if not defined
+
+        // stat options
+        opts->primary = NULL;
+        opts->secondary = NULL;
+
+        return opts;
 }
 
 // Produce an overscan vector from an array of pixels
-psVector *pmOverscanVector(float *chi2, // chi^2 from fit
-                           pmOverscanStatOptions *overscanOpts, // Overscan statistic options
-                           const psArray *pixels // Array of vectors containing the pixel values
-    )
+psVector *pmOverscanVector(float *chi2,                                                 // chi^2 from fit
+                                                   pmOverscanStatOptions *overscanOpts, // Overscan statistic options
+                                                   const psArray *pixels                                // Array of vectors containing the pixel values
+)
 {
-    assert(overscanOpts);
-    assert(pixels);
-    
-    // statisctic to be calculated
-    psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use
-    if (!statistic) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat);
-        return false;
-    }
-    psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
-
-    // Reduce the overscans
-    psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row
-    psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate
-    psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting
-
-    for (int i = 0; i < pixels->n; i++) {
-        psVector *values = pixels->data[i]; // Vector with overscan values
-        if (values->n > 0) {
-            mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
-            ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1]
-            if (!psVectorStats(stats, values, NULL, NULL, 0)) {
-                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+        assert(overscanOpts);
+        assert(pixels);
+
+        // statisctic to be calculated
+        psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use
+        if (!statistic)
+        {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat);
+                return false;
+        }
+        psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original
+
+        // Reduce the overscans
+        psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32);              // Overscan for each row
+        psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32);             // Ordinate
+        psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting
+
+        for (int i = 0; i < pixels->n; i++)
+        {
+                psVector *values = pixels->data[i]; // Vector with overscan values
+                if (values->n > 0)
+                {
+                        mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0;
+                        ordinate->data.F32[i] = 2.0 * (float)i / (float)pixels->n - 1.0; // Scale to [-1,1]
+                        if (!psVectorStats(stats, values, NULL, NULL, 0))
+                        {
+                                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
+                                goto escape;
+                        }
+                        reduced->data.F32[i] = psStatsGetValue(stats, statistic);
+                }
+                else
+                {
+                        if (overscanOpts->fitType == PM_FIT_NONE)
+                        {
+                                psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n");
+                                goto escape;
+                        }
+                        else
+                        {
+                                // We'll fit this one out
+                                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
+                        }
+                }
+        }
+        // Smooth the reduced vector
+        if (overscanOpts->boxcar > 0)
+        {
+                psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector
+                psFree(reduced);
+                reduced = smoothed;
+        }
+        if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0)
+        {
+                if (overscanOpts->boxcar > 0)
+                {
+                        psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it.");
+                }
+                psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA);
+                psFree(reduced);
+                reduced = smoothed;
+        }
+
+        // Fit the overscan, if required
+        psVector *fitted = NULL; // Fitted overscan values
+        switch (overscanOpts->fitType)
+        {
+        case PM_FIT_NONE:
+                // No fitting --- that's easy.
+                fitted = psMemIncrRefCounter(reduced);
+                break;
+        case PM_FIT_POLY_ORD:
+                if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
+                                                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD))
+                {
+                        psFree(overscanOpts->poly);
+                        overscanOpts->poly = NULL;
+                }
+                if (!overscanOpts->poly)
+                {
+                        overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
+                }
+                psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
+                fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+                break;
+        case PM_FIT_POLY_CHEBY:
+                if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
+                                                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB))
+                {
+                        psFree(overscanOpts->poly);
+                        overscanOpts->poly = NULL;
+                }
+                if (!overscanOpts->poly)
+                {
+                        overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
+                }
+                psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
+                fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
+                break;
+        case PM_FIT_SPLINE:
+
+                // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes
+                // a knot for every input point.  it needs an argument like 'number of knots' for the
+                // output spline.  EAM: still true 2023.01.22
+
+                // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
+                // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate);
+                psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n");
+                break;
+        default:
+                psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
                 goto escape;
-            }
-            reduced->data.F32[i] = psStatsGetValue(stats, statistic);
-        } else {
-            if (overscanOpts->fitType == PM_FIT_NONE) {
-                psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n");
-                goto escape;
-            } else {
-                // We'll fit this one out
-                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1;
-            }
-        }
-    }
-    // Smooth the reduced vector
-    if (overscanOpts->boxcar > 0) {
-        psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector
+        }
+
+        if (chi2)
+        {
+                *chi2 = 0.0; // chi^2 (sort of)
+                for (int i = 0; i < reduced->n; i++)
+                {
+                        *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]);
+                }
+        }
+
         psFree(reduced);
-        reduced = smoothed;
-    }
-    if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) {
-        if (overscanOpts->boxcar > 0) {
-            psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it.");
-        }
-        psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA);
+        psFree(ordinate);
+        psFree(mask);
+        psFree(stats);
+        return fitted;
+
+escape:
         psFree(reduced);
-        reduced = smoothed;
-    }
-
-    // Fit the overscan, if required
-    psVector *fitted = NULL;                   // Fitted overscan values
-    switch (overscanOpts->fitType) {
-      case PM_FIT_NONE:
-        // No fitting --- that's easy.
-        fitted = psMemIncrRefCounter(reduced);
-        break;
-      case PM_FIT_POLY_ORD:
-        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
-                                   overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) {
-            psFree(overscanOpts->poly);
-            overscanOpts->poly = NULL;
-        }
-        if (! overscanOpts->poly) {
-            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order);
-        }
-        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
-        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
-        break;
-      case PM_FIT_POLY_CHEBY:
-        if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order ||
-                                   overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) {
-            psFree(overscanOpts->poly);
-            overscanOpts->poly = NULL;
-        }
-        if (! overscanOpts->poly) {
-            overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order);
-        }
-        psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate);
-        fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate);
-        break;
-      case PM_FIT_SPLINE:
-
-        // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes
-        // a knot for every input point.  it needs an argument like 'number of knots' for the
-        // output spline.  EAM: still true 2023.01.22
-
-        // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate);
-        // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate);
-        psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n");
-        break;
-      default:
-        psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType);
-        goto escape;
-    }
-
-    if (chi2) {
-        *chi2 = 0.0;                    // chi^2 (sort of)
-        for (int i = 0; i < reduced->n; i++) {
-            *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]);
-        }
-    }
-
-    psFree(reduced);
-    psFree(ordinate);
-    psFree(mask);
-    psFree(stats);
-    return fitted;
-
-escape:
-    psFree(reduced);
-    psFree(ordinate);
-    psFree(mask);
-    psFree(stats);
-    return NULL;
+        psFree(ordinate);
+        psFree(mask);
+        psFree(stats);
+        return NULL;
 }
 
@@ -621 +703 @@
 // the reduced overscan vector with the 0-order element from the polynomial fit.  Not clear is
 // adding any useful information
-bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2) {
-
-    psString comment = NULL;    // Comment to add
-
-    switch (overscanOpts->fitType) {
-      case PM_FIT_POLY_ORD:
-      case PM_FIT_POLY_CHEBY: {
-          psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2);
-          psPolynomial1D *poly = overscanOpts->poly; // The polynomial
-          for (int i = 0; i < poly->nX; i++) {
-              psStringAppend(&comment, "%.1f ", poly->coeff[i]);
-          }
-          psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
-          psFree(comment);
-          comment = NULL;
-
-          // write metadata header value
-          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]);
-          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]);
-          break;
-      }
-      case PM_FIT_SPLINE: {
-          /*
-            psSpline1D *spline = overscanOpts->spline; // The spline
-            for (int i = 0; i < spline->n; i++) {
-            psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i);
-            psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
-            for (int j = 0; j < poly->nX; j++) {
-            psStringAppend(&comment, "%.1f ", poly->coeff[i]);
-            }
-            psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-            comment, "");
-            psFree(comment);
-            comment = NULL;
-            }
-          */
-          // write metadata header value
-          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE,
-                           "Overscan value", NAN);
-          psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE,
-                           "Overscan stdev", NAN);
-          break;
-      }
-      case PM_FIT_NONE:
-        break;
-      default:
-        psAbort("Should never get here!!!\n");
-    }
-    return true;
+bool pmOverscanUpdateHeader(pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2)
+{
+
+        psString comment = NULL; // Comment to add
+
+        switch (overscanOpts->fitType)
+        {
+        case PM_FIT_POLY_ORD:
+        case PM_FIT_POLY_CHEBY:
+        {
+                psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2);
+                psPolynomial1D *poly = overscanOpts->poly; // The polynomial
+                for (int i = 0; i < poly->nX; i++)
+                {
+                        psStringAppend(&comment, "%.1f ", poly->coeff[i]);
+                }
+                psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
+                psFree(comment);
+                comment = NULL;
+
+                // write metadata header value
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]);
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]);
+                break;
+        }
+        case PM_FIT_SPLINE:
+        {
+                /*
+                  psSpline1D *spline = overscanOpts->spline; // The spline
+                  for (int i = 0; i < spline->n; i++) {
+                  psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i);
+                  psPolynomial1D *poly = spline->spline[i]; // i-th polynomial
+                  for (int j = 0; j < poly->nX; j++) {
+                  psStringAppend(&comment, "%.1f ", poly->coeff[i]);
+                  }
+                  psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
+                  comment, "");
+                  psFree(comment);
+                  comment = NULL;
+                  }
+                */
+                // write metadata header value
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE,
+                                                 "Overscan value", NAN);
+                psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE,
+                                                 "Overscan stdev", NAN);
+                break;
+        }
+        case PM_FIT_NONE:
+                break;
+        default:
+                psAbort("Should never get here!!!\n");
+        }
+        return true;
 }
 
 // generate stats of overscan vector for header
 // reduced: 1D vector with overscan stats
-bool pmOverscanUpdateHeaderVector (pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced) {
-
-    psString comment = NULL;    // Comment to add
-    psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
-    if (!psVectorStats (vectorStats, reduced, 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);
-    return true;
+bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced)
+{
+
+        psString comment = NULL; // Comment to add
+        psStats *vectorStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
+        if (!psVectorStats(vectorStats, reduced, 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);
+        return true;
 }
-