Changeset 42679

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

Location:

branches/2dbias

Files:

• ippconfig/gpc1/format_20100723.config (modified) (2 diffs)
• ippconfig/gpc1/ppImage.config (modified) (1 diff)
• ippconfig/recipes/ppImage.config (modified) (1 diff)
• ppImage/src/ppImage.h (modified) (1 diff)
• ppImage/src/ppImageDetrendPattern.c (modified) (16 diffs)
• ppImage/src/ppImageDetrendReadout.c (modified) (12 diffs)
• psModules/src/detrend/pmOverscan.c (modified) (5 diffs)
• psModules/src/detrend/pmOverscan.h (modified) (1 diff)

branches/2dbias/ippconfig/gpc1/format_20100723.config

@@ -166 +166 @@
 # 2D bias subtraction requires ppImage.config OVERSCAN.2D to be TRUE
                CELL.BIASSEC.SOURCE     STR     VALUE
-               CELL.BIASSEC            STR     [591:624,1:598],[1:590,598:608]
+               CELL.BIASSEC            STR     [591:624,1:598],[1:624,599:608]
 #
 # a single region will be used for constant and 1D bias subtraction
@@ -174 +174 @@
 END
 
+OVERSCAN.2D.SUBSET METADATA
+      XY00    BOOL    TRUE
+      XY01    BOOL    TRUE
+      XY02    BOOL    TRUE
+      XY03    BOOL    TRUE
+      XY04    BOOL    TRUE
+      XY05    BOOL    TRUE
+      XY06    BOOL    TRUE
+      XY07    BOOL    TRUE
+      XY10    BOOL    TRUE
+      XY11    BOOL    TRUE
+      XY12    BOOL    TRUE
+      XY13    BOOL    TRUE
+      XY14    BOOL    TRUE
+      XY15    BOOL    TRUE
+      XY16    BOOL    TRUE
+      XY17    BOOL    TRUE
+      XY20    BOOL    TRUE
+      XY21    BOOL    TRUE
+      XY22    BOOL    TRUE
+      XY23    BOOL    TRUE
+      XY24    BOOL    TRUE
+      XY25    BOOL    TRUE
+      XY26    BOOL    TRUE
+      XY27    BOOL    TRUE
+      XY30    STR     1111100011111000111100001111000011110000111110001111000011110000
+      XY31    BOOL    TRUE
+      XY32    BOOL    TRUE
+      XY33    BOOL    TRUE
+      XY34    BOOL    TRUE
+      XY35    BOOL    TRUE
+      XY36    BOOL    TRUE
+      XY37    BOOL    TRUE
+      XY40    BOOL    TRUE
+      XY41    BOOL    TRUE
+      XY42    BOOL    TRUE
+      XY43    BOOL    TRUE
+      XY44    BOOL    TRUE
+      XY45    BOOL    TRUE
+      XY46    BOOL    TRUE
+      XY47    BOOL    TRUE
+      XY50    BOOL    TRUE
+      XY51    BOOL    TRUE
+      XY52    BOOL    TRUE
+      XY53    BOOL    TRUE
+      XY54    BOOL    TRUE
+      XY55    BOOL    TRUE
+      XY56    BOOL    TRUE
+      XY57    BOOL    TRUE
+      XY60    BOOL    TRUE
+      XY61    BOOL    TRUE
+      XY62    BOOL    TRUE
+      XY63    BOOL    TRUE
+      XY64    BOOL    TRUE
+      XY65    BOOL    TRUE
+      XY66    BOOL    TRUE
+      XY67    BOOL    TRUE
+      XY70    BOOL    TRUE
+      XY71    BOOL    TRUE
+      XY72    BOOL    TRUE
+      XY73    BOOL    TRUE
+      XY74    BOOL    TRUE
+      XY75    BOOL    TRUE
+      XY76    BOOL    TRUE
+      XY77    BOOL    TRUE
+END
 
 # How to translate PS concepts into FITS headers

branches/2dbias/ippconfig/gpc1/ppImage.config

@@ -9 +9 @@
 OVERSCAN.STAT           STR     MEDIAN
 OVERSCAN.BOXCAR         S32     3
+
+OVERSCAN.2D.SUBSET METADATA
+      XY00    BOOL    TRUE
+      XY01    BOOL    TRUE
+      XY02    BOOL    TRUE
+      XY03    BOOL    TRUE
+      XY04    BOOL    TRUE
+      XY05    BOOL    TRUE
+      XY06    BOOL    TRUE
+      XY07    BOOL    TRUE
+      XY10    BOOL    TRUE
+      XY11    BOOL    TRUE
+      XY12    BOOL    TRUE
+      XY13    BOOL    TRUE
+      XY14    BOOL    TRUE
+      XY15    BOOL    TRUE
+      XY16    BOOL    TRUE
+      XY17    BOOL    TRUE
+      XY20    BOOL    TRUE
+      XY21    BOOL    TRUE
+      XY22    BOOL    TRUE
+      XY23    BOOL    TRUE
+      XY24    BOOL    TRUE
+      XY25    BOOL    TRUE
+      XY26    BOOL    TRUE
+      XY27    BOOL    TRUE
+      XY30    BOOL    TRUE
+      XY31    BOOL    TRUE
+      XY32    BOOL    TRUE
+      XY33    BOOL    TRUE
+      XY34    BOOL    TRUE
+      XY35    BOOL    TRUE
+      XY36    BOOL    TRUE
+      XY37    BOOL    TRUE
+      XY40    BOOL    TRUE
+      XY41    BOOL    TRUE
+      XY42    BOOL    TRUE
+      XY43    BOOL    TRUE
+      XY44    BOOL    TRUE
+      XY45    BOOL    TRUE
+      XY46    BOOL    TRUE
+      XY47    BOOL    TRUE
+      XY50    BOOL    TRUE
+      XY51    BOOL    TRUE
+      XY52    BOOL    TRUE
+      XY53    BOOL    TRUE
+      XY54    BOOL    TRUE
+      XY55    BOOL    TRUE
+      XY56    BOOL    TRUE
+      XY57    BOOL    TRUE
+      XY60    BOOL    TRUE
+      XY61    BOOL    TRUE
+      XY62    BOOL    TRUE
+      XY63    BOOL    TRUE
+      XY64    BOOL    TRUE
+      XY65    BOOL    TRUE
+      XY66    BOOL    TRUE
+      XY67    BOOL    TRUE
+      XY70    BOOL    TRUE
+      XY71    BOOL    TRUE
+      XY72    BOOL    TRUE
+      XY73    BOOL    TRUE
+      XY74    BOOL    TRUE
+      XY75    BOOL    TRUE
+      XY76    BOOL    TRUE
+      XY77    BOOL    TRUE
+END
 
 # Normalization class

branches/2dbias/ippconfig/recipes/ppImage.config

@@ -110 +110 @@
 
 OVERSCAN.2D             BOOL    FALSE           # use a 2D model for the overscan subtraction?
+OVERSCAN.2D.SUBSET METADATA
+END
 
 # statistics to use to convert 2D top-overscan region to a vector

branches/2dbias/ppImage/src/ppImage.h

@@ -188 +188 @@
 
 bool ppImageDetrendFree(pmConfig *config, pmFPAview *view);
+bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue);
 bool ppImageFringeFree(pmConfig *config, pmFPAview *view);
 

branches/2dbias/ppImage/src/ppImageDetrendPattern.c

@@ -5 +5 @@
 #include "ppImage.h"
 
-#define ESCAPE(STATUS,...) {                            \
-        psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__);   \
-        psFree(view);                                   \
-        return false;                                   \
-    }
-
-static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue);
+#define ESCAPE(STATUS, ...)                           \
+    {                                                 \
+        psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \
+        psFree(view);                                 \
+        return false;                                 \
+    }
+
+bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipename, const char *recipevalue);
 
 bool ppImageDetrendPatternRowApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options);
@@ -26 +27 @@
     assert(inputView->readout == -1);
 
-    if (!ppImageDetrendPatternRowApply (config, chip, inputView, options)) {
-        return false;
-    }
-    if (!ppImageDetrendPatternContinuityApply (config, chip, inputView, options)) {
-        return false;
-    }
-    if (!ppImageDetrendPatternCellApply (config, chip, inputView, options)) {
-        return false;
-    }
-    if (!ppImageDetrendPatternDeadCellsApply (config, chip, inputView, options)) {
-        return false;
-    }
-    return(true);
+    if (!ppImageDetrendPatternRowApply(config, chip, inputView, options))
+    {
+        return false;
+    }
+    if (!ppImageDetrendPatternContinuityApply(config, chip, inputView, options))
+    {
+        return false;
+    }
+    if (!ppImageDetrendPatternCellApply(config, chip, inputView, options))
+    {
+        return false;
+    }
+    if (!ppImageDetrendPatternDeadCellsApply(config, chip, inputView, options))
+    {
+        return false;
+    }
+    return (true);
 }
 
@@ -45 +50 @@
     bool status;
 
-    if (!options->doPatternRow) return true;
-    
+    if (!options->doPatternRow)
+        return true;
+
     assert(inputView->chip != -1);
     assert(inputView->cell == -1);
@@ -73 +79 @@
 
     pmHDU *hdu = pmHDUFromChip(chip);
-    if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW")) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done.");
-        return true;
+    if (psMetadataLookupBool(&status, hdu->header, "PTRN_ROW"))
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Not performing row pattern correction as it has already been done.");
+        return true;
     }
 
@@ -85 +92 @@
 
     // grab the PATTERN.ROW.AMP file
-    pmFPAfile *PRAfile = psMetadataLookupPtr (&status, config->files, "PPIMAGE.PATTERN.ROW.AMP");
-    if (!PRAfile) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) ");
-    }
-        
+    pmFPAfile *PRAfile = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.ROW.AMP");
+    if (!PRAfile)
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Row Amplitude file not found, applying to all (with limits from ROW.SUBSET) ");
+    }
+
     pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
     *view = *inputView;
 
     pmCell *cell = NULL;
-    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-        // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
-        // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
-        // psWarning ("Looping through %s, %s\n", chipName, cellName);
-
-        if (!cell->process || !cell->file_exists) {
-            continue;
-        }
-
-        // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?)
-        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
-            ESCAPE(false, "load failure for Cell");
-        }
-            
-        if (!cell->data_exists) {
-            continue;
-        }
-
-        if (cell->readouts->n > 1) {
-            psWarning ("Skipping Video Cell for ppImageDetrendPatternApply");
-            continue;
-        }
-
-        // grab the corresponding cell
-        if (PRAfile) {
-            pmCell *PRAcell = pmFPAviewThisCell (view, PRAfile->fpa);
-            psAssert (PRAcell, "found Pattern Row Amplitude file, but not cell?");
-
-            // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
-            float amplitude = psMetadataLookupF32 (&status, PRAcell->analysis, "PTN.ROW.AMP");
-            if (!status) amplitude = NAN;
-            
-            // put the value on the science cell
-            psMetadataAddF32 (cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
-        }
-
-        bool doPattern = false;
-        if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET")) {
-            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
-        }
-        if (!doPattern) continue;
-
-        // switch to test threaded version
-        if (true) {
-            // I need to allocate a view here to be freed by the
-            // called function below.  
-            pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing
-            *myView = *view;
-
-            // allocate a job, construct the arguments for this job
-            psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL");
-            psArrayAdd(job->args, 1, config);
-            psArrayAdd(job->args, 1, input->fpa);
-            psArrayAdd(job->args, 1, chip);
-            psArrayAdd(job->args, 1, cell);
-            psArrayAdd(job->args, 1, myView);
-            psArrayAdd(job->args, 1, options);
-            if (!psThreadJobAddPending(job)) {
-                return false;
-            }
-        } else {
-            // bump the counter since it must be freed by the function below.  
-            psMemIncrRefCounter (view); // View for local processing
-            if (!ppImageDetrendPatternApplyCell (config, input->fpa, chip, cell, view, options)) {
-                return false;
-            }
-        }
+    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL)
+    {
+        // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME");
+        // const char *cellName = psMetadataLookupStr(NULL, cell->concepts, "CELL.NAME");
+        // psWarning ("Looping through %s, %s\n", chipName, cellName);
+
+        if (!cell->process || !cell->file_exists)
+        {
+            continue;
+        }
+
+        // this forces pmFPAfileRead of the PATTERN.ROW.AMP file (XXX but is this needed?)
+        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE))
+        {
+            ESCAPE(false, "load failure for Cell");
+        }
+
+        if (!cell->data_exists)
+        {
+            continue;
+        }
+
+        if (cell->readouts->n > 1)
+        {
+            psWarning("Skipping Video Cell for ppImageDetrendPatternApply");
+            continue;
+        }
+
+        // grab the corresponding cell
+        if (PRAfile)
+        {
+            pmCell *PRAcell = pmFPAviewThisCell(view, PRAfile->fpa);
+            psAssert(PRAcell, "found Pattern Row Amplitude file, but not cell?");
+
+            // find the nominal signal amplitude (check the ghost and/or crosstalk recipe file)
+            float amplitude = psMetadataLookupF32(&status, PRAcell->analysis, "PTN.ROW.AMP");
+            if (!status)
+                amplitude = NAN;
+
+            // put the value on the science cell
+            psMetadataAddF32(cell->analysis, PS_LIST_TAIL, "PTN.ROW.AMP", PS_META_REPLACE, "", amplitude);
+        }
+
+        bool doPattern = false;
+        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.ROW.SUBSET"))
+        {
+            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
+        }
+        if (!doPattern)
+            continue;
+
+        // switch to test threaded version
+        if (true)
+        {
+            // I need to allocate a view here to be freed by the
+            // called function below.
+            pmFPAview *myView = pmFPAviewAlloc(0); // View for local processing
+            *myView = *view;
+
+            // allocate a job, construct the arguments for this job
+            psThreadJob *job = psThreadJobAlloc("PPIMAGE_PATTERN_ROW_CELL");
+            psArrayAdd(job->args, 1, config);
+            psArrayAdd(job->args, 1, input->fpa);
+            psArrayAdd(job->args, 1, chip);
+            psArrayAdd(job->args, 1, cell);
+            psArrayAdd(job->args, 1, myView);
+            psArrayAdd(job->args, 1, options);
+            if (!psThreadJobAddPending(job))
+            {
+                return false;
+            }
+        }
+        else
+        {
+            // bump the counter since it must be freed by the function below.
+            psMemIncrRefCounter(view); // View for local processing
+            if (!ppImageDetrendPatternApplyCell(config, input->fpa, chip, cell, view, options))
+            {
+                return false;
+            }
+        }
     }
 
     // wait here for the threaded jobs to finish
-    // if no threads are allocated, this 
-    if (!psThreadPoolWait(true, true)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");
-        return false;
-    }
-
-    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW",PS_META_REPLACE,"PATTERN.ROW correction applied",true);
+    // if no threads are allocated, this
+    if (!psThreadPoolWait(true, true))
+    {
+        psError(PS_ERR_UNKNOWN, false, "Unable to apply bias correction.");
+        return false;
+    }
+
+    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_ROW", PS_META_REPLACE, "PATTERN.ROW correction applied", true);
     psFree(view);
-   
+
     return true;
 }
@@ -180 +203 @@
     bool status;
 
-    if (!options->doPatternContinuity) return true;
+    if (!options->doPatternContinuity)
+        return true;
 
     assert(inputView->chip != -1);
@@ -188 +212 @@
     // see the comment for PATTERN.ROW; the same rules apply for PATTERN.CELL
 
-    int numCells = chip->cells->n;       // Number of cells
+    int numCells = chip->cells->n;                         // Number of cells
     psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell?
-    pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
+    pmFPAview *view = pmFPAviewAlloc(0);                   // View for local processing
     *view = *inputView;
 
     pmHDU *hdu = pmHDUFromChip(chip);
-    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON")) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done.");
-        return true;
+    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CON"))
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction as it has already been done.");
+        return true;
     }
 
@@ -202 +227 @@
     psLogMsg("ppImage", PS_LOG_INFO, "Performing pattern continuity correction for %s\n", chipName);
 
-    for (int i = 0; i < chip->cells->n; i++) {
-        view->cell = i;
-
-        pmCell *cell = chip->cells->data[i]; // Cell of interest
-
-        if (cell->readouts->n > 1) {
-            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell.");
-            continue;
-        }
-
-        bool doPattern = false;
-        if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET")) {
-            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
-        }
-
-        if (doPattern) {
-            tweak->data.U8[i] = 0xFF;
-        }
+    for (int i = 0; i < chip->cells->n; i++)
+    {
+        view->cell = i;
+
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+
+        if (cell->readouts->n > 1)
+        {
+            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell continuity correction on video cell.");
+            continue;
+        }
+
+        bool doPattern = false;
+        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CONTINUITY.SUBSET"))
+        {
+            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
+        }
+
+        if (doPattern)
+        {
+            tweak->data.U8[i] = 0xFF;
+        }
     }
 
     // Tweak the cells
     if (!pmPatternContinuity(chip, tweak, options->patternCellBG, options->patternCellMean,
-                             options->maskValue, options->darkMask,options->patternContinuityEdgeWidth)) {
-        psFree(tweak);
-        psFree(view);
-        return false;
+                             options->maskValue, options->darkMask, options->patternContinuityEdgeWidth))
+    {
+        psFree(tweak);
+        psFree(view);
+        return false;
     }
     psFree(tweak);
     psFree(view);
 
-    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON",PS_META_REPLACE,"PATTERN.CONTINUITY correction applied",true);
+    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CON", PS_META_REPLACE, "PATTERN.CONTINUITY correction applied", true);
     return true;
 }
@@ -240 +270 @@
     bool status;
 
-    if (!options->doPatternCell) return true;
+    if (!options->doPatternCell)
+        return true;
 
     assert(inputView->chip != -1);
@@ -246 +277 @@
     assert(inputView->readout == -1);
 
-    int numCells = chip->cells->n;       // Number of cells
+    int numCells = chip->cells->n;                         // Number of cells
     psVector *tweak = psVectorAlloc(numCells, PS_TYPE_U8); // Tweak cell?
-    pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
+    pmFPAview *view = pmFPAviewAlloc(0);                   // View for local processing
     *view = *inputView;
 
     pmHDU *hdu = pmHDUFromChip(chip);
-    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL")) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done.");
-        return true;
+    if (psMetadataLookupBool(&status, hdu->header, "PTRN_CEL"))
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction as it has already been done.");
+        return true;
     }
 
@@ -260 +292 @@
     psLogMsg("ppImage", PS_LOG_INFO, "Performing cell pattern correction for %s\n", chipName);
 
-    for (int i = 0; i < chip->cells->n; i++) {
-        view->cell = i;
-
-        pmCell *cell = chip->cells->data[i]; // Cell of interest
-
-        if (cell->readouts->n > 1) {
-            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell.");
-            continue;
-        }
-
-        bool doPattern = false;
-        if (!doPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET")) {
-            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
-        }
-        if (doPattern) {
-            tweak->data.U8[i] = 0xFF;
-        }
+    for (int i = 0; i < chip->cells->n; i++)
+    {
+        view->cell = i;
+
+        pmCell *cell = chip->cells->data[i]; // Cell of interest
+
+        if (cell->readouts->n > 1)
+        {
+            psLogMsg("ppImage", PS_LOG_INFO, "Not performing cell pattern correction on video cell.");
+            continue;
+        }
+
+        bool doPattern = false;
+        if (!ppImageDoPatternForView(&doPattern, config, chip, view, RECIPE_NAME, "PATTERN.CELL.SUBSET"))
+        {
+            ESCAPE(false, "Unable to determine whether row pattern matching should be applied.");
+        }
+        if (doPattern)
+        {
+            tweak->data.U8[i] = 0xFF;
+        }
     }
 
     // Tweak the cells
     if (!pmPatternCell(chip, tweak, options->patternCellBG, options->patternCellMean,
-                       options->maskValue, options->darkMask)) {
-        psFree(tweak);
-        psFree(view);
-        return false;
+                       options->maskValue, options->darkMask))
+    {
+        psFree(tweak);
+        psFree(view);
+        return false;
     }
     psFree(tweak);
     psFree(view);
 
-    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL",PS_META_REPLACE,"PATTERN.CELL correction applied",true);
+    psMetadataAddBool(hdu->header, PS_LIST_TAIL, "PTRN_CEL", PS_META_REPLACE, "PATTERN.CELL correction applied", true);
     return true;
 }
 
-bool ppImageDetrendPatternDeadCellsMask (pmChip *chip, psImageMaskType maskVal) {
-
-    int numCells = chip->cells->n;       // Number of cells
+bool ppImageDetrendPatternDeadCellsMask(pmChip *chip, psImageMaskType maskVal)
+{
+
+    int numCells = chip->cells->n; // Number of cells
 
     // now mask bad cells
-    for (int i = 0; i < numCells; i++) {
-        pmCell *cell = chip->cells->data[i]; // Cell of interest
+    for (int i = 0; i < numCells; i++)
+    {
+        pmCell *cell = chip->cells->data[i];     // Cell of interest
         pmReadout *ro = cell->readouts->data[0]; // Readout of interest
 
-        psImage *mask = ro->mask; // mask of interest
-        int numCols = mask->numCols, numRows = mask->numRows; // Size of image
-        
-        for (int y = 0; y < numRows; y++) {
-            for (int x = 0; x < numCols; x++) {
-                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
-            }
-        }
+        psImage *mask = ro->mask;                             // mask of interest
+        int numCols = mask->numCols, numRows = mask->numRows; // Size of image
+
+        for (int y = 0; y < numRows; y++)
+        {
+            for (int x = 0; x < numCols; x++)
+            {
+                mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal;
+            }
+        }
     }
     return true;
 }
 
-bool ppImageDetrendPatternDeadCellsXY42 (pmChip *chip, psImageMaskType maskVal) {
+bool ppImageDetrendPatternDeadCellsXY42(pmChip *chip, psImageMaskType maskVal)
+{
 
     bool status;
 
     // extract the MEDIAN_CELL_BACKGROUND values to check for problems
-    int numCells = chip->cells->n;       // Number of cells
-    psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32);
-
-    for (int i = 0; i < numCells; i++) {
-        pmCell *cell = chip->cells->data[i];;
-
-        if (!cell->process || !cell->file_exists || !cell->data_exists) continue;
-        
-        // select the HDU for this cell
-        pmHDU *cellHDU = pmHDUFromCell(cell);  // HDU of interest
-        
-        psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
-        if (!status) continue;
-
-        psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR");
-        if (!status) continue;
-
-        if (BackValue == 0.0) continue;
-        if (BackError == 0.0) continue;
-
-        psVectorAppend (cellBackground, BackValue / BackError);
-    }
-
-    if (cellBackground->n < 0.375*numCells) {
-        // Chip is bad, mask the whole thing
-        ppImageDetrendPatternDeadCellsMask (chip, maskVal);
-        psFree (cellBackground);
-        return true;
+    int numCells = chip->cells->n; // Number of cells
+    psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32);
+
+    for (int i = 0; i < numCells; i++)
+    {
+        pmCell *cell = chip->cells->data[i];
+        ;
+
+        if (!cell->process || !cell->file_exists || !cell->data_exists)
+            continue;
+
+        // select the HDU for this cell
+        pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest
+
+        psF32 BackValue = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
+        if (!status)
+            continue;
+
+        psF32 BackError = psMetadataLookupF32(&status, cellHDU->header, "BACK_ERR");
+        if (!status)
+            continue;
+
+        if (BackValue == 0.0)
+            continue;
+        if (BackError == 0.0)
+            continue;
+
+        psVectorAppend(cellBackground, BackValue / BackError);
+    }
+
+    if (cellBackground->n < 0.375 * numCells)
+    {
+        // Chip is bad, mask the whole thing
+        ppImageDetrendPatternDeadCellsMask(chip, maskVal);
+        psFree(cellBackground);
+        return true;
     }
 
     // Second, calculate the median
-    psVectorSortInPlace (cellBackground);
+    psVectorSortInPlace(cellBackground);
     int midPt = cellBackground->n / 2.0;
-    float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5*(cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt-1]);
-    psFree (cellBackground);
+    float median = cellBackground->n % 2 ? cellBackground->data.F32[midPt] : 0.5 * (cellBackground->data.F32[midPt] + cellBackground->data.F32[midPt - 1]);
+    psFree(cellBackground);
 
     // XXX hardwired number
-    if (median < 3.0) {
-        // Chip is bad, mask the whole thing
-        ppImageDetrendPatternDeadCellsMask (chip, maskVal);
-    }
-    psFree (cellBackground);
+    if (median < 3.0)
+    {
+        // Chip is bad, mask the whole thing
+        ppImageDetrendPatternDeadCellsMask(chip, maskVal);
+    }
+    psFree(cellBackground);
     return true;
 }
 
 // to apply the dead cell pattern, we need to transfer the pattern for this chip from the
-// pmFPAfile for the pattern to the one for the image being processed.  
+// pmFPAfile for the pattern to the one for the image being processed.
 bool ppImageDetrendPatternDeadCellsApply(pmConfig *config, pmChip *chip, const pmFPAview *inputView, ppImageOptions *options)
 {
     bool status;
 
-    if (!options->doPatternDeadCells) return true;
+    if (!options->doPatternDeadCells)
+        return true;
 
     assert(inputView->chip != -1);
@@ -380 +432 @@
 
     pmHDU *hdu = pmHDUFromChip(chip);
-    if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED")) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done.");
-        return true;
-    }
-
-    if (!strcmp (chipName, "XY40") || !strcmp (chipName, "XY42")) {
-        // special case : check for BACK_VAL / BACK_ERR > 3 or < 3
-        psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42");
-        ppImageDetrendPatternDeadCellsXY42 (chip, options->blankMask);
-        return true;
+    if (psMetadataLookupBool(&status, hdu->header, "PTRN_DED"))
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Not performing dead cell pattern correction as it has already been done.");
+        return true;
+    }
+
+    if (!strcmp(chipName, "XY40") || !strcmp(chipName, "XY42"))
+    {
+        // special case : check for BACK_VAL / BACK_ERR > 3 or < 3
+        psLogMsg("ppImage", PS_LOG_INFO, "Using special case for XY40 and XY42");
+        ppImageDetrendPatternDeadCellsXY42(chip, options->blankMask);
+        return true;
     }
 
     pmFPAfile *pattern = psMetadataLookupPtr(&status, config->files, "PPIMAGE.PATTERN.DEAD.CELLS");
-    if (!pattern) {
-        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping");
-        return true;
+    if (!pattern)
+    {
+        psLogMsg("ppImage", PS_LOG_INFO, "Pattern Dead Cells file not found, skipping");
+        return true;
     }
 
@@ -401 +456 @@
     *view = *inputView;
 
-    pmChip *patternChip = pmFPAviewThisChip (view, pattern->fpa);
+    pmChip *patternChip = pmFPAviewThisChip(view, pattern->fpa);
 
     // grab the pattern from the input pattern file chip
-    psImage *deadCellPattern = (psImage *) psMetadataLookupPtr (&status, patternChip->analysis, "PTN.DEAD.CELL");
-    if (!deadCellPattern) {
+    psImage *deadCellPattern = (psImage *)psMetadataLookupPtr(&status, patternChip->analysis, "PTN.DEAD.CELL");
+    if (!deadCellPattern)
+    {
         psLogMsg("psModules.detrend", PS_LOG_DETAIL, "No DEAD CELL pattern for chip, skipping\n");
-        psFree (view);
-        return true;
+        psFree(view);
+        return true;
     }
 
     // copy the pattern pointer to the chip of the image being processed
-    psMetadataAddImage (chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern);
+    psMetadataAddImage(chip->analysis, PS_LIST_TAIL, "PTN.DEAD.CELL", PS_META_REPLACE, "", deadCellPattern);
 
     // extract the MEDIAN_CELL_BACKGROUND values to check for problems
-    int numCells = chip->cells->n;       // Number of cells
-    psVector *cellBackground = psVectorAllocEmpty (numCells, PS_DATA_F32);
+    int numCells = chip->cells->n; // Number of cells
+    psVector *cellBackground = psVectorAllocEmpty(numCells, PS_DATA_F32);
 
     pmCell *cell = NULL;
-    while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL) {
-        if (!cell->process || !cell->file_exists || !cell->data_exists) {
-            psVectorAppend (cellBackground, NAN);
-            continue;
-        }
-        
-        // select the HDU for this cell
-        pmHDU *cellHDU = pmHDUFromCell(cell);  // HDU of interest
-        
-        psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
-        if (!status) {
-            psVectorAppend (cellBackground, NAN);
-            continue;
-        }
-        psVectorAppend (cellBackground, value);
+    while ((cell = pmFPAviewNextCell(view, chip->parent, 1)) != NULL)
+    {
+        if (!cell->process || !cell->file_exists || !cell->data_exists)
+        {
+            psVectorAppend(cellBackground, NAN);
+            continue;
+        }
+
+        // select the HDU for this cell
+        pmHDU *cellHDU = pmHDUFromCell(cell); // HDU of interest
+
+        psF32 value = psMetadataLookupF32(&status, cellHDU->header, "BACK_VAL");
+        if (!status)
+        {
+            psVectorAppend(cellBackground, NAN);
+            continue;
+        }
+        psVectorAppend(cellBackground, value);
     }
 
     // match cellBackground pattern to registered patterns and mask as needed
-    if (!pmPatternDeadCells(chip, cellBackground, options->blankMask)) {
-        psFree(cellBackground);
-        psFree(view);
-        return false;
+    if (!pmPatternDeadCells(chip, cellBackground, options->blankMask))
+    {
+        psFree(cellBackground);
+        psFree(view);
+        return false;
     }
     psFree(cellBackground);
@@ -449 +509 @@
 }
 
-static bool doPatternForView (bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue) {
+bool ppImageDoPatternForView(bool *doit, const pmConfig *config, const pmChip *chip, const pmFPAview *view, const char *recipeName, const char *recipeValue)
+{
 
     *doit = false;
@@ -456 +517 @@
 
     doPattern = pmConfigRecipeValueByView(config, recipeName, recipeValue, chip->parent, view);
-    
-    if (!doPattern) {
+
+    if (!doPattern)
+    {
         psError(PS_ERR_UNKNOWN, false, "Unable to determine whether row pattern matching should be applied.");
         return false;
     }
-    if (doPattern->type == PS_DATA_BOOL) {
+    if (doPattern->type == PS_DATA_BOOL)
+    {
         *doit = doPattern->data.B;
         return true;
     }
-    if (doPattern->type == PS_DATA_STRING) {
+    if (doPattern->type == PS_DATA_STRING)
+    {
         // expect a string of the form "000110001001" with at least view->cell entries
         char *string = doPattern->data.str;
-        if (strlen(string) < view->cell) {
-            psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int) strlen(string));
+        if (strlen(string) < view->cell)
+        {
+            psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string (too few elements %d)", (int)strlen(string));
             return false;
         }
-        switch (string[view->cell]) {
-          case '0':
-          case 'f':
-          case 'F':
-          case 'n':
-          case 'N':
+        switch (string[view->cell])
+        {
+        case '0':
+        case 'f':
+        case 'F':
+        case 'n':
+        case 'N':
             *doit = false;
             return true;
-          case '1':
-          case 't':
-          case 'T':
-          case 'y':
-          case 'Y':
+        case '1':
+        case 't':
+        case 'T':
+        case 'y':
+        case 'Y':
             *doit = true;
             return true;
-          default:
+        default:
             psError(PS_ERR_UNKNOWN, true, "error in PATTERN.ROW.SUBSET chip string %s (unknown value %c))", string, string[view->cell]);
             return false;
@@ -497 +563 @@
 }
 
-// thread safety : 
-bool ppImageDetrendPatternApplyCell (pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options) {
-  
+// thread safety :
+bool ppImageDetrendPatternApplyCell(pmConfig *config, pmFPA *fpa, pmChip *chip, pmCell *cell, pmFPAview *view, ppImageOptions *options)
+{
+
     // process each of the readouts
-    pmReadout *readout;         // Readout from cell
-    while ((readout = pmFPAviewNextReadout (view, fpa, 1)) != NULL) {
-        if (!readout->data_exists) {
-            continue;
-        }
-
-        // Perform pattern correction
-        if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter,
-                          options->patternRowRej, options->patternRowThresh, options->patternRowMean,
-                          options->patternRowStdev, options->maskValue, options->darkMask)) {
-            psFree (view);
-            return false;
-        }
-    }
-    psFree (view); // supplied view must be freeable
+    pmReadout *readout; // Readout from cell
+    while ((readout = pmFPAviewNextReadout(view, fpa, 1)) != NULL)
+    {
+        if (!readout->data_exists)
+        {
+            continue;
+        }
+
+        // Perform pattern correction
+        if (!pmPatternRow(readout, options->patternRowOrder, options->patternRowIter,
+                          options->patternRowRej, options->patternRowThresh, options->patternRowMean,
+                          options->patternRowStdev, options->maskValue, options->darkMask))
+        {
+            psFree(view);
+            return false;
+        }
+    }
+    psFree(view); // supplied view must be freeable
     return true;
 }
 
-psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView) {
-    
+psVector *ppImageDetrendPatternCellFailures(pmFPAfile *input, const pmFPAview *inputView)
+{
+
     bool status = false;
     pmCell *cell = NULL;
-    
+
     pmFPAview *view = pmFPAviewAlloc(0); // View for local processing
     *view = *inputView;
 
     // extract the MEDIAN_CELL_BACKGROUND values to check for problems
-    psVector *cellBackground = psVectorAllocEmpty (64, PS_DATA_F32);
-    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL) {
-        if (!cell->process || !cell->file_exists || !cell->data_exists) {
-            psVectorAppend (cellBackground, NAN);
-            continue;
-        }
-        
-        // select the HDU for this cell
-        pmHDU *hdu = pmHDUFromCell(cell);  // HDU of interest
-        
-        psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL");
-        psVectorAppend (cellBackground, value);
+    psVector *cellBackground = psVectorAllocEmpty(64, PS_DATA_F32);
+    while ((cell = pmFPAviewNextCell(view, input->fpa, 1)) != NULL)
+    {
+        if (!cell->process || !cell->file_exists || !cell->data_exists)
+        {
+            psVectorAppend(cellBackground, NAN);
+            continue;
+        }
+
+        // select the HDU for this cell
+        pmHDU *hdu = pmHDUFromCell(cell); // HDU of interest
+
+        psF32 value = psMetadataLookupF32(&status, hdu->header, "BACK_VAL");
+        psVectorAppend(cellBackground, value);
     }
     return cellBackground;

branches/2dbias/ppImage/src/ppImageDetrendReadout.c

@@ -5 +5 @@
 #include "ppImage.h"
 
+#define ESCAPE(STATUS, ...)                           \
+    {                                                 \
+        psError(PS_ERR_UNKNOWN, STATUS, __VA_ARGS__); \
+        psFree(view);                                 \
+        return false;                                 \
+    }
+
 bool ppImageDetrendReadout(pmConfig *config, ppImageOptions *options, pmFPAview *view)
 {
-  // psTimerStart("detrend.readout");
+    // psTimerStart("detrend.readout");
 
     // construct a view for the detrend images (which have only one readout)
@@ -18 +25 @@
 
     // Check that the gain is set (this is used by both pmReadoutGenerateMask and pmReadoutGenerateVariance)
-    { 
-      float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell
-      if (!isfinite(gain)) {
-        psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe
-        psAssert(recipe, "Should be there!");
-        bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain?
-        if (override) {
-          psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout);
-          psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item
-          psAssert(item, "Should be there!");
-          item->data.F32 = 1.0;
-
-          // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated
-          psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE");
-        } else {
-          psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout);
-        }
-      }
+    {
+        float gain = psMetadataLookupF32(NULL, input->parent->concepts, "CELL.GAIN"); // Gain for cell
+        if (!isfinite(gain))
+        {
+            psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, RECIPE_NAME); // Recipe
+            psAssert(recipe, "Should be there!");
+            bool override = psMetadataLookupBool(NULL, recipe, "GAIN.OVERRIDE"); // Override the bad gain?
+            if (override)
+            {
+                psWarning("CELL.GAIN is not set for readout (%d,%d,%d) --- setting to unity.", view->chip, view->cell, view->readout);
+                psMetadataItem *item = psMetadataLookup(input->parent->concepts, "CELL.GAIN"); // Gain item
+                psAssert(item, "Should be there!");
+                item->data.F32 = 1.0;
+
+                // for unity gain, there is no modification for the readnoise, note that it has (effectively) been updated
+                psMetadataRemoveKey(input->parent->concepts, "CELL.READNOISE.UPDATE");
+            }
+            else
+            {
+                psWarning("CELL.GAIN is NAN for readout (%d,%d,%d), image will be masked.", view->chip, view->cell, view->readout);
+            }
+        }
     }
     // psLogMsg ("ppImage", 6, "check gain: %f sec\n", psTimerMark ("detrend.readout"));
@@ -43 +54 @@
     bool hasVideo = false;
     {
-      // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config?
-      bool status = false;
-      if (!input) goto done_video_check;
-      if (!input->parent) goto done_video_check;
-      if (!input->parent->parent) goto done_video_check;
-      if (!input->parent->parent->hdu) goto done_video_check;
-      if (!input->parent->parent->hdu->header) goto done_video_check;
-      char *ptr = psMetadataLookupStr(&status,input->parent->parent->hdu->header,"CELLMODE");
-      if (status) {
-        psLogMsg ("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
-        char *Vptr = strchr(ptr,'V');
-        if (Vptr) {
-          hasVideo = options->hasVideo = true;
-          psLogMsg ("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n",(int) options->hasVideo,(int) options->useVideoDark, (int) options->useVideoMask);
-        }
-      }
-    }
-    done_video_check:
+        // XXX test for GPC1? or CAN_HAVE_VIDEO in camera config?
+        bool status = false;
+        if (!input)
+            goto done_video_check;
+        if (!input->parent)
+            goto done_video_check;
+        if (!input->parent->parent)
+            goto done_video_check;
+        if (!input->parent->parent->hdu)
+            goto done_video_check;
+        if (!input->parent->parent->hdu->header)
+            goto done_video_check;
+        char *ptr = psMetadataLookupStr(&status, input->parent->parent->hdu->header, "CELLMODE");
+        if (status)
+        {
+            psLogMsg("ppImage.detrend", PS_LOG_DETAIL, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask);
+            char *Vptr = strchr(ptr, 'V');
+            if (Vptr)
+            {
+                hasVideo = options->hasVideo = true;
+                psLogMsg("ppImage.detrend", PS_LOG_INFO, "VIDEO: %d %d %d\n", (int)options->hasVideo, (int)options->useVideoDark, (int)options->useVideoMask);
+            }
+        }
+    }
+done_video_check:
     // psLogMsg ("ppImage", 6, "check video: %f sec\n", psTimerMark ("detrend.readout"));
 
@@ -74 +92 @@
 
     // Masking on the basis of pixel value needs to be done before anything else, so the values are pristine.
-    if (options->doMaskBuild) {
+    if (options->doMaskBuild)
+    {
         psImageMaskType satMask = options->doMaskSat ? options->satMask : 0;
         psImageMaskType lowMask = options->doMaskLow ? options->lowMask : 0;
-        if (!pmReadoutGenerateMask(input, satMask, lowMask)) {
-          psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask.");
-          psFree(detview);
-          return false;
-        }
-        // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout"));
+        if (!pmReadoutGenerateMask(input, satMask, lowMask))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate a mask.");
+            psFree(detview);
+            return false;
+        }
+        // psLogMsg ("ppImage", 6, "generate mask: %f sec\n", psTimerMark ("detrend.readout"));
     }
     // apply the externally supplied mask to the input->mask pixels
-    if (options->doMask) {
-      pmReadout *mask;
-      if ((options->useVideoMask)&&(hasVideo)) {
-        mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK");
-      }
-      else {
-        mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
-      }
-      pmMaskBadPixels(input, mask, options->maskValue);
-      // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout"));
+    if (options->doMask)
+    {
+        pmReadout *mask;
+        if ((options->useVideoMask) && (hasVideo))
+        {
+            mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEOMASK");
+        }
+        else
+        {
+            mask = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.MASK");
+        }
+        pmMaskBadPixels(input, mask, options->maskValue);
+        // psLogMsg ("ppImage", 6, "apply mask: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
@@ -101 +124 @@
         // extern bool ppImageBurntoolApply(pmConfig *, ppImageOptions *, pmFPAview *, pmReadout *);
         ppImageBurntoolApply(config, options, view, input);
-        // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-    
-    if (options->doMaskBurntool) {
-        if (options->doApplyBurntool) {
+        // psLogMsg ("ppImage", 6, "apply burntool: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
+    if (options->doMaskBurntool)
+    {
+        if (options->doApplyBurntool)
+        {
             // build burntool mask from data input burntool table
-            ppImageBurntoolMaskFromTable(config,options,view,input);
-        } else {
+            ppImageBurntoolMaskFromTable(config, options, view, input);
+        }
+        else
+        {
             // build burntool mask from data in the input image's fits extension
-            ppImageBurntoolMask(config,options,view,input);
-        }
-        // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout"));
+            ppImageBurntoolMask(config, options, view, input);
+        }
+        // psLogMsg ("ppImage", 6, "apply burntool mask: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
     // Subtract the overscan
-    if (options->doOverscan) {
-      if (!pmOverscanSubtract (input, options->overscan)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
-        psFree(detview);
-        return false;
-      }
-      // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout"));
+    if (options->doOverscan)
+    {
+        bool doTwoDOverscan = false;
+        if (options->overscan->TwoD)
+        {
+            pmChip *chip = pmFPAfileThisChip(config->files, view, "PPIMAGE.INPUT");
+            if (!ppImageDoPatternForView(&doTwoDOverscan, config, chip, view, RECIPE_NAME, "OVERSCAN.2D.SUBSET"))
+            {
+                ESCAPE(false, "Unable to determine whether 2D Overscan subtraction should be applied.");
+            }
+        }
+        if (!pmOverscanSubtract(input, options->overscan, doTwoDOverscan))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Unable to subtract overscan.");
+            psFree(detview);
+            return false;
+        }
+        // psLogMsg ("ppImage", 6, "subtract overscan: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
     // measure the overscan-subtracted readoutBackground here (or subtract the overscan value?)
     // XXX this is the measurements and should be independent of the pattern masking
-    if (options->doPatternDeadCells) {
-        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
-        psImageBackground (stats, NULL, input->image, NULL, 0xffff, rng);
-
-        // save this value somewhere
-        pmHDU *hdu = pmHDUFromReadout(input);  // HDU of interest
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian);
-        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev);
-        psFree (stats);
-        psFree (rng);
+    if (options->doPatternDeadCells)
+    {
+        psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
+        psImageBackground(stats, NULL, input->image, NULL, 0xffff, rng);
+
+        // save this value somewhere
+        pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_VAL", PS_META_REPLACE, "Median cell background", stats->robustMedian);
+        psMetadataAddF32(hdu->header, PS_LIST_TAIL, "BACK_ERR", PS_META_REPLACE, "Stdev of cell background", stats->robustStdev);
+        psFree(stats);
+        psFree(rng);
     }
 
     // Non-linearity correction
-    if (options->doNonLin) {
-      if (!ppImageDetrendNonLinear(input,detview,config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
-        psFree(detview);
-        return(false);
-      }
-      // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
+    if (options->doNonLin)
+    {
+        if (!ppImageDetrendNonLinear(input, detview, config))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Unable to correct NonLinearity");
+            psFree(detview);
+            return (false);
+        }
+        // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
     }
     // New Non-linearity correction (exclusive of the above)
-    if (options->doNewNonLin) {
-      if (!ppImageDetrendNewNonLinear(input, detview, config)) {
-        psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)");
-        psFree(detview);
-        return(false);
-      }
-      // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
+    if (options->doNewNonLin)
+    {
+        if (!ppImageDetrendNewNonLinear(input, detview, config))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Unable to correct Non-Linearity with new version (2023)");
+            psFree(detview);
+            return (false);
+        }
+        // psLogMsg ("ppImage", 6, "nonlinear correction: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
     // set up the dark and bias
-    pmCell *dark = NULL;                // Multi-dark
-    pmReadout *oldDark = NULL;          // Old-fashioned dark
-    pmReadout  *bias = NULL;
-    if (options->doBias) {
+    pmCell *dark = NULL;       // Multi-dark
+    pmReadout *oldDark = NULL; // Old-fashioned dark
+    pmReadout *bias = NULL;
+    if (options->doBias)
+    {
         bias = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.BIAS");
     }
-    if (options->doDark) {
-        bool mdok;                      // Status of MD lookup
-        psMetadata *recipe = psMetadataLookupPtr (&mdok, config->recipes, RECIPE_NAME);
+    if (options->doDark)
+    {
+        bool mdok; // Status of MD lookup
+        psMetadata *recipe = psMetadataLookupPtr(&mdok, config->recipes, RECIPE_NAME);
         assert(mdok && recipe);
 
-        if ((options->useVideoDark)&&(hasVideo)) {
-          if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
-            oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK");
-          } else {
-            dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK");
-          }
-        }
-        else {
-          if (psMetadataLookupBool(&mdok, recipe, "OLDDARK")) {
-            oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
-          } else {
-            dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
-          }
-        }
+        if ((options->useVideoDark) && (hasVideo))
+        {
+            if (psMetadataLookupBool(&mdok, recipe, "OLDDARK"))
+            {
+                oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.VIDEODARK");
+            }
+            else
+            {
+                dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.VIDEODARK");
+            }
+        }
+        else
+        {
+            if (psMetadataLookupBool(&mdok, recipe, "OLDDARK"))
+            {
+                oldDark = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.DARK");
+            }
+            else
+            {
+                dark = pmFPAfileThisCell(config->files, detview, "PPIMAGE.DARK");
+            }
+        }
     }
 
     // Bias and temperature-independent-dark subtraction are merged.
-    if (options->doBias) {
-        if (!pmBiasSubtract(input, bias, oldDark, view)) {
+    if (options->doBias)
+    {
+        if (!pmBiasSubtract(input, bias, oldDark, view))
+        {
             psError(PS_ERR_UNKNOWN, false, "Unable to subtract bias.");
             psFree(detview);
             return false;
         }
-        // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-    
+        // psLogMsg ("ppImage", 6, "apply bias: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
     // Weight on the basis of pixel value needs to be done after the overscan has been subtracted
-    if (options->doVarianceBuild) {
+    if (options->doVarianceBuild)
+    {
         // create the target mask and variance images
         psImage *noiseImage = NULL;
-        if (options->doNoiseMap) {
+        if (options->doNoiseMap)
+        {
             // XXX convert the noiseMap image to a binned image
             pmReadout *noiseMap = NULL;
             noiseMap = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.NOISEMAP");
-            noiseImage = psImageCopy (NULL, input->image, PS_TYPE_F32);
-            psImageInit (noiseImage, 0.0);
+            noiseImage = psImageCopy(NULL, input->image, PS_TYPE_F32);
+            psImageInit(noiseImage, 0.0);
 
             // XXX this works, but is not really quite right: the model shoud include the
             // offset information, we are not really getting exactly the right mapping from the
             // original file.
-            // CZW 2012-03-21: I do not believe this is true anymore.  In any case, this seems to be what
-            //                 ppImageReplaceBackground does to do sky subtraction.
+            // CZW 2012-03-21: I do not believe this is true anymore.  In any case, this seems to be what
+            //                 ppImageReplaceBackground does to do sky subtraction.
             psImageBinning *binning = psImageBinningAlloc();
             binning->nXruff = noiseMap->image->numCols;
@@ -220 +277 @@
             psImageBinningSetScale(binning, PS_IMAGE_BINNING_LEFT);
 
-            psImageUnbin (noiseImage, noiseMap->image, binning);
-            psFree (binning);
-        }
-
-        if (!pmReadoutGenerateVariance(input, noiseImage, true)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image.");
-            psFree(detview);
-            return false;
-        }
-        psFree (noiseImage);
-        // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-
-    if (options->doDark && dark) {
-        if (!pmDarkApply(input, dark, options->darkMask)) {
+            psImageUnbin(noiseImage, noiseMap->image, binning);
+            psFree(binning);
+        }
+
+        if (!pmReadoutGenerateVariance(input, noiseImage, true))
+        {
+            psError(PS_ERR_UNKNOWN, false, "Unable to generate a variance image.");
+            psFree(detview);
+            return false;
+        }
+        psFree(noiseImage);
+        // psLogMsg ("ppImage", 6, "generate variance: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
+    if (options->doDark && dark)
+    {
+        if (!pmDarkApply(input, dark, options->darkMask))
+        {
             psError(PS_ERR_UNKNOWN, false, "Unable to subtract dark.");
             psFree(detview);
             return false;
         }
-        // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-
-    if (options->doRemnance) {
+        // psLogMsg ("ppImage", 6, "apply dark: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
+    if (options->doRemnance)
+    {
         if (!pmRemnance(input, options->maskValue, options->lowMask,
-                        options->remnanceSize, options->remnanceThresh)) {
+                        options->remnanceSize, options->remnanceThresh))
+        {
             psError(PS_ERR_UNKNOWN, false, "Unable to mask remnance.");
             psFree(detview);
             return false;
         }
-        // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout"));
+        // psLogMsg ("ppImage", 6, "mask remnance: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
     // Shutter correction
-    if (options->doShutter) {
+    if (options->doShutter)
+    {
         pmReadout *shutter = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.SHUTTER");
-        if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config))) {
-            psFree(detview);
-            return false;
-        }
-        // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout"));
+        if (!pmShutterCorrectionApply(input, shutter, pmConfigMaskGet("FLAT", config)))
+        {
+            psFree(detview);
+            return false;
+        }
+        // psLogMsg ("ppImage", 6, "shutter correction: %f sec\n", psTimerMark ("detrend.readout"));
     }
 
     // Flat-field correction (no options used?)
-    if (options->doFlat) {
+    if (options->doFlat)
+    {
         pmReadout *flat = pmFPAfileThisReadout(config->files, detview, "PPIMAGE.FLAT");
-        if (!pmFlatField(input, flat, options->flatMask)) {
-            psFree(detview);
-            return false;
-        }
-        // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-
-/*     // Pattern noise correction */
-/*     if (options->doPattern) { */
-/*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
-/*                           options->patternThresh, options->patternMean, options->patternStdev, */
-/*                           options->maskValue, options->darkMask)) { */
-/*             psFree(detview); */
-/*             return false; */
-/*         } */
-/*     } */
+        if (!pmFlatField(input, flat, options->flatMask))
+        {
+            psFree(detview);
+            return false;
+        }
+        // psLogMsg ("ppImage", 6, "apply flat: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
+    /*     // Pattern noise correction */
+    /*     if (options->doPattern) { */
+    /*         if (!pmPatternRow(input, options->patternOrder, options->patternIter, options->patternRej, */
+    /*                           options->patternThresh, options->patternMean, options->patternStdev, */
+    /*                           options->maskValue, options->darkMask)) { */
+    /*             psFree(detview); */
+    /*             return false; */
+    /*         } */
+    /*     } */
 
     // Normalization by a single (known) constant
-    bool mdok;                          // Status of MD lookup
+    bool mdok; // Status of MD lookup
     float norm = psMetadataLookupF32(&mdok, config->arguments, "NORMALIZATION");
-    if (mdok && isfinite(norm) && norm != 1.0) {
+    if (mdok && isfinite(norm) && norm != 1.0)
+    {
         pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
-        psString comment = NULL;        // Comment to add
+        psString comment = NULL;              // Comment to add
         psStringAppend(&comment, "Normalization: %f", norm);
         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
@@ -293 +360 @@
 
         psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
-        // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-
-# if (1)
+        // psLogMsg ("ppImage", 6, "renormalize: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+
+#if (1)
     // Normalization by per-class values
     psMetadata *normlist = psMetadataLookupMetadata(&mdok, config->arguments, "NORMALIZATION.TABLE");
-    if (normlist) {
+    if (normlist)
+    {
         pmFPAfile *inputFile = psMetadataLookupPtr(&mdok, config->files, "PPIMAGE.INPUT");
 
         // get the menu of class IDs
         psMetadata *menu = psMetadataLookupMetadata(&mdok, inputFile->camera, "CLASSID");
-        if (!menu) {
+        if (!menu)
+        {
             psError(PS_ERR_IO, false, "Unable to find CLASSID metadata in camera configuration");
             psFree(detview);
@@ -311 +380 @@
         // get the rule for class_id for the desired class
         const char *rule = psMetadataLookupStr(&mdok, menu, options->normClass);
-        if (!rule) {
+        if (!rule)
+        {
             psError(PS_ERR_IO, false, "Unable to find NORM.CLASS value %s in CLASSID in camera configuration", options->normClass);
             psFree(detview);
@@ -318 +388 @@
         // get the class_id from the rule
         char *classID = pmFPAfileNameFromRule(rule, inputFile, view);
-        if (!classID) {
+        if (!classID)
+        {
             psError(PS_ERR_IO, false, "error converting CLASSID rule %s to name\n", rule);
             psFree(detview);
@@ -325 +396 @@
 
         // get normalization from the class_id
-        float norm = psMetadataLookupF32 (&mdok, normlist, classID);
-        if (!mdok) {
+        float norm = psMetadataLookupF32(&mdok, normlist, classID);
+        if (!mdok)
+        {
             psError(PS_ERR_IO, false, "failed to find class ID %s in normalization table\n", classID);
             psFree(detview);
@@ -333 +405 @@
 
         pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest
-        psString comment = NULL;        // Comment to add
+        psString comment = NULL;              // Comment to add
         psStringAppend(&comment, "Normalization: %f", norm);
         psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, "");
@@ -341 +413 @@
         psBinaryOp(input->image, input->image, "*", psScalarAlloc(norm, PS_TYPE_F32));
 
-        psFree (classID);
-
-        // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout"));
-    }
-# endif
-
-    if (options->doFringe) {
+        psFree(classID);
+
+        // psLogMsg ("ppImage", 6, "renormalize by class: %f sec\n", psTimerMark ("detrend.readout"));
+    }
+#endif
+
+    if (options->doFringe)
+    {
         pmCell *fringe = pmFPAfileThisCell(config->files, detview, "PPIMAGE.FRINGE");
-        if (!ppImageDetrendFringeMeasure(input, fringe, false, options)) {
-            psFree(detview);
-            return false;
-        }
-        // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout"));
+        if (!ppImageDetrendFringeMeasure(input, fringe, false, options))
+        {
+            psFree(detview);
+            return false;
+        }
+        // psLogMsg ("ppImage", 6, "measure fringe: %f sec\n", psTimerMark ("detrend.readout"));
     }
 

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;
 }
-

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

@@ -81 +81 @@
 // bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2);
 
-bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts);
+bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan);
 
 /// @}