Changeset 42889 for trunk/psModules

Timestamp:

Jun 5, 2025, 3:33:01 PM (14 months ago)

Author:

tdeboer

Message:

updates to detrend functions for 2D overscan correction

Location:

trunk/psModules/src/detrend

Files:

• pmOverscan.c (modified) (1 diff)
• pmOverscan.h (modified) (2 diffs)
• pmPattern.c (modified) (5 diffs)

trunk/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, 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
+                // XXX HG: investigate whether the y-overscan is being used or the x-overscan
+                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, false);
+                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, false);
+                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");
+
+                //parse 2D overscan BISASEC region information
+                psRegion yRegion = overscanOpts->secondary->biassecslow;
+                psRegion xRegion = overscanOpts->secondary->biassecfast;
+                            
+                // Cut out relevant region image. Need to use the parent uncut image
+                psImage *parent = (psImage *) image->parent;
+
+                // the serial (fast readout) direction is columns (x-dir)
+                psImage *yscan = psImageSubset(parent, yRegion);  // overscan region spanning all rows
+                psImage *xscan = psImageSubset(parent, xRegion);  // 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, true);
+                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, true);
+                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 how to normalize the xReduced vector
+                        // slice the part overlapping withe yoverscan from xReduced
+                        long startIndex = yscan->col0 - xscan->col0;
+                        long normSize = xReduced->n - startIndex; // Calculate the size of the new vector
+                        // Allocate the new vector
+                        psVector *normVector = psVectorAlloc(normSize, PS_TYPE_F32);
+                        // Copy the data from xReduced starting from startIndex to the end
+                        for (long i = startIndex; i < xReduced->n; i++)
+                        {
+                                normVector->data.F32[i - startIndex] = xReduced->data.F32[i];
+                        }
+                        // Now normVector holds the sliced data; then, compute its robust mean
+                        psStatsOptions statistic = PS_STAT_ROBUST_MEDIAN;
+                        psStats *stats = psStatsAlloc(statistic);
+                        if (!psVectorStats(stats, normVector, NULL, NULL, 0))
+                        {
+                                psError(PS_ERR_UNKNOWN, false, "failure to measure robust median as the normalization of xReduced");
+                                psFree(stats);
+                                psFree(normVector);
+                                return NULL;
+                        }
+                        float xReducedNormalized = psStatsGetValue(stats, statistic);
+                        psFree(stats);
+                        psFree(normVector);
+                        // 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] - xReducedNormalized + 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;
+
+} // End of overscan subtraction
 
 static void pmOverscanOptionsFree(pmOverscanOptions *options)
 {
-    psFree(options->stat);
-    psFree(options->poly);
-    psFree(options->spline);
+        psFree(options->primary);
+        psFree(options->secondary);
 }
 
-pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat,
-                                          int boxcar, float gauss)
+static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options)
 {
-    pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions));
-    psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree);
-
-    // Inputs
-    opts->single = single;
-    opts->constant = false;
-    opts->fitType = fitType;
-    opts->order = order;
-    opts->stat = psMemIncrRefCounter(stat);
-
-    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
-
-    // Smoothing
-    opts->boxcar = boxcar;
-    opts->gauss = gauss;
-
-    // Outputs
-    opts->poly = NULL;
-    opts->spline = NULL;
-
-    return opts;
+        psFree(options->stat);
+        psFree(options->poly);
+        psFree(options->spline);
 }
 
+// Globally defined
+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;
+}
+
+// Globally defined
+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;
+}
+
 // Produce an overscan vector from an array of pixels
-psVector *pmOverscanVector(float *chi2, // chi^2 from fit
-                           pmOverscanOptions *overscanOpts, // Overscan options
-                           const psArray *pixels, // Array of vectors containing the pixel values
-                           psStats *myStats // Statistic to use in reducing the overscan
-    )
+psVector *pmOverscanVector(float *chi2, // chi^2 from fit
+                           pmOverscanStatOptions *overscanOpts, // Overscan statistic options
+                           const psArray *pixels,                               // Array of vectors containing the pixel values
+                           bool robust)                                                 // Use robust statistics for boxcar smoothing
 {
-    assert(overscanOpts);
-    assert(pixels);
-    assert(myStats);
-
-    psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use
-    assert(statistic != 0);
-
-    // 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(myStats, values, NULL, NULL, 0)) {
+        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, robust); // 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;
+        }
+
+        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;
+}
+
+// XXX EAM 2024.04.13 : this function seems poorly conceived.  it is replacing the stats from
+// 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;
+}
+
+// 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;
-            }
-            reduced->data.F32[i] = psStatsGetValue(myStats, 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");
-            return NULL;
-        } 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);
-        return NULL;
-        break;
-    }
-
-    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);
-
-    return fitted;
-}
-
-bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *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 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));
-
-        return true;
-    }
-
-    // we are performing a statitical analysis of the overscan region
-
-    // Check for an unallowable pmFit.
-    if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD &&
-        overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) {
-        psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d).  Returning original image.\n",
-                overscanOpts->fitType);
-        return false;
-    }
-
-    psList *overscans = input->bias; // List of the overscan images
-
-    psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use
-    if (statistic == 0) {
-        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
-
-    psString comment = NULL;    // Comment to add
-    psStringAppend(&comment, "Subtracting overscan (stat %x; type %x; order %d)",
-                   statistic, overscanOpts->fitType, overscanOpts->order);
-    psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK,
-                     comment, "");
-    psTrace ("psModules.detrend", 4, "%s\n", comment);
-    psFree(comment);
-
-    // Reduce all overscan pixels to a single value
-    if (overscanOpts->single) {
-        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);
-
-        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);
+        }
+        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 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);
+        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 ((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);
+        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 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 (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;
-                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, pixels, stats);
-        psFree(pixels);
-        if (! reduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
-            psFree(stats);
-            return false;
-        }
-
-        // generate stats of overscan vector for header
-        { 
-            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);
-        }
-
-        // 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 (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, pixels, stats);
-        psFree(pixels);
-        if (! reduced) {
-            psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n");
-            psFree(stats);
-            return false;
-        }
-
-        // generate stats of overscan vector for header
-        { 
-          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 (3) %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);
-        }
-
-        // 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);
-    }
-
-    pmOverscanUpdateHeader (hdu, overscanOpts, chi2);
-    psFree(stats);
-
-    return true;
-
-} // End of overscan subtraction
-    
+}

trunk/psModules/src/detrend/pmOverscan.h

@@ -18 +18 @@
 
 /// Type of fit to perform
-typedef enum {
-    PM_FIT_NONE,                        ///< No fit
-    PM_FIT_POLY_ORD,                    ///< Fit ordinary polynomial
-    PM_FIT_POLY_CHEBY,                  ///< Fit Chebyshev polynomial
-    PM_FIT_SPLINE                       ///< Fit cubic splines
+typedef enum
+{
+    PM_FIT_NONE,       ///< No fit
+    PM_FIT_POLY_ORD,   ///< Fit ordinary polynomial
+    PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial
+    PM_FIT_SPLINE      ///< Fit cubic splines
 } pmFit;
 
@@ -35 +36 @@
 {
     // Inputs
-    bool single;                        ///< Reduce all overscan regions to a single value?
-    bool constant;                      ///< use a supplied constant value (do not measure region)
-    float value;                        ///< supplied value if needed (per above)
-    pmFit fitType;                      ///< Type of fit to overscan
-    unsigned int order;                 ///< Order of polynomial, or number of spline pieces
-    psStats *stat;                      ///< Statistic to use when reducing the minor direction
-    int boxcar;                         ///< Boxcar smoothing radius
-    float gauss;                        ///< Gaussian smoothing sigma
-    float minValid;                     ///< if overscan is too low, readout is dead : mask
-    float maxValid;                     ///< if overscan is too high, readout is dead : mask
-    psImageMaskType maskVal;            ///< Mask value to give dead readouts
+    pmFit fitType;      ///< Type of fit to overscan
+    unsigned int order; ///< Order of polynomial, or number of spline pieces
+    psStats *stat;      ///< Statistic to use when reducing the minor direction
+    int boxcar;         ///< Boxcar smoothing radius
+    float gauss;        ///< Gaussian smoothing sigma
+    psRegion biassecslow; ///< The 2D slow bias section
+    psRegion biassecfast; ///< The 2D fast bias section
 
-    // Outputs
-    psPolynomial1D *poly;               ///< Result of polynomial fit
-    psSpline1D *spline;                 ///< Result of spline fit
-}
-pmOverscanOptions;
+    
+    // These are used to carry the fitted functions, but are only used to generate
+    // an updated reduced vector
+    psPolynomial1D *poly; ///< Result of polynomial fit
+    psSpline1D *spline;   ///< Result of spline fit
+} pmOverscanStatOptions;
+
+/// Options for overscan subtraction
+///
+/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if
+/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the
+/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without
+/// fitting (if the structure defies characterisation).  In any case, statistics are required to reduce
+/// multiple values to a single value (either for the scan, or for the entire overscan regions).
+typedef struct
+{
+    // Inputs
+    bool single;             ///< Reduce all overscan regions to a single value?
+    bool constant;           ///< use a supplied constant value (do not measure region)
+    bool TwoD;               ///< use a supplied constant value (do not measure region)
+    float value;             ///< supplied value if needed (per above)
+    float minValid;          ///< if overscan is too low, readout is dead : mask
+    float maxValid;          ///< if overscan is too high, readout is dead : mask
+    psImageMaskType maskVal; ///< Mask value to give dead readouts
+
+    pmOverscanStatOptions *primary;
+    pmOverscanStatOptions *secondary;
+} pmOverscanOptions;
 
 /// Allocator for overscan options
-pmOverscanOptions *pmOverscanOptionsAlloc(bool single, ///< Reduce all overscan regions to a single value?
-                                          pmFit fitType, ///< Type of fit to overscan
-                                          unsigned int order, ///< Order of polynomial, or number of splines
-                                          psStats *stat, ///< Statistic to use
-                                          int boxcar, ///< Boxcar smoothing radius
-                                          float gauss ///< Gaussian smoothing sigma
-                                         );
+pmOverscanOptions *pmOverscanOptionsAlloc(void);
 
-psVector *pmOverscanVector(float *chi2, // chi^2 from fit
-                           pmOverscanOptions *overscanOpts, // Overscan options
-                           const psArray *pixels, // Array of vectors containing the pixel values
-                           psStats *myStats // Statistic to use in reducing the overscan
-    );
+pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void);
 
-bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2);
+psVector *pmOverscanVector(float *chi2,                         // chi^2 from fit
+                           pmOverscanStatOptions *overscanOpts, // Overscan options
+                           const psArray *pixels,               // Array of vectors containing the pixel values
+                           bool robust                          // use robust statistics for boxcar smoothing
+);
 
-bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts);
+// bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2);
+
+bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan);
 
 /// @}
 #endif
-

trunk/psModules/src/detrend/pmPattern.c

@@ -239 +239 @@
 # define READNOISE 10 /* arbitrary number */
     float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    if(isnan(sigma)) sigma=20.;
     float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); 
     float lower = stats->robustMedian - delta; // Lower bound for data
@@ -246 +247 @@
     // if the noise from the background is too large 
     float significance = nominalAmplitude / sigma;
+    significance = abs(significance);
    
     // XXX EAM : arbitrary number
     if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      return true; 
+    }
+    //Add some safeties in case of bad fits
+    if (isnan(significance) || (background < 0.0)) {
       psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
       return true; 
@@ -554 +561 @@
 # define READNOISE 10 /* arbitrary number */
     float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE));
+    if(isnan(sigma)) sigma=20.; // XXX EAM : somewhat arbitrary number
     float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); 
     float lower = stats->robustMedian - delta; // Lower bound for data
@@ -561 +569 @@
     // if the noise from the background is too large 
     float significance = nominalAmplitude / sigma;
-   
+    significance = abs(significance);
+      
     // XXX EAM : arbitrary number
     if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) {
@@ -569 +578 @@
       return true; 
     }
+    //Add some safeties in case of bad fits
+    if (isnan(significance) || (background < 0.0)) {
+      psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction due to bad fit for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
+      psFree(stats);
+      psFree(rng);
+      return true; 
+    }
+    
     psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance);
 # endif