Changeset 42653

Timestamp:

Apr 12, 2024, 2:49:40 PM (2 years ago)

Author:

eugene

Message:

added rudimentary 2D biassec mods

File:

• branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c (modified) (4 diffs)

branches/eam_branches/psModules.20240412/src/detrend/pmOverscan.c

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