Index: trunk/psModules/src/detrend/pmOverscan.c =================================================================== --- trunk/psModules/src/detrend/pmOverscan.c (revision 42888) +++ trunk/psModules/src/detrend/pmOverscan.c (revision 42889) @@ -17,537 +17,807 @@ #include "pmOverscan.h" -#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends +#define SMOOTH_NSIGMA 4.0 // Number of Gaussian sigma the smoothing kernel extends + +static void pmOverscanOptionsFree(pmOverscanOptions *options); +static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options); +bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced); + +bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan) +{ + + assert(input); + + if (overscanOpts == NULL) + return true; // no overscan subtraction requested + + pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest + psImage *image = input->image; + + // check for 'soft bias' (simple, fixed offset to be subtracted) + if (overscanOpts->constant) + { + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", overscanOpts->value); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + // NOTE psBinaryOp frees arg2 if it is a scalar + (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32)); + + return true; + } + + // we are performing a statitical analysis of the overscan region + + // Check for an unallowable pmFit. + if (overscanOpts->primary) + { + if (overscanOpts->primary->fitType != PM_FIT_NONE && + overscanOpts->primary->fitType != PM_FIT_POLY_ORD && + overscanOpts->primary->fitType != PM_FIT_POLY_CHEBY && + overscanOpts->primary->fitType != PM_FIT_SPLINE) + { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", + overscanOpts->primary->fitType); + return false; + } + } + if (overscanOpts->secondary) + { + if (overscanOpts->secondary->fitType != PM_FIT_NONE && + overscanOpts->secondary->fitType != PM_FIT_POLY_ORD && + overscanOpts->secondary->fitType != PM_FIT_POLY_CHEBY && + overscanOpts->secondary->fitType != PM_FIT_SPLINE) + { + psError(PS_ERR_UNKNOWN, true, "Invalid fit type (2D) (%d). Returning original image.\n", + overscanOpts->secondary->fitType); + return false; + } + } + + psList *overscans = input->bias; // List of the overscan images + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) + { + + // extract overscan pixels to a single vector + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + int index = pixels->n; // Index + pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); + pixels->n += overscan->numRows * overscan->numCols; + for (int i = 0; i < overscan->numRows; i++) + { + memcpy(&pixels->data.F32[index], overscan->data.F32[i], + overscan->numCols * sizeof(psF32)); + index += overscan->numCols; + } + } + psFree(iter); + + // statistic to be calculated + psStatsOptions statistic = psStatsSingleOption(overscanOpts->primary->stat->options); // Statistic to use + if (!statistic) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", + overscanOpts->primary->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + if (!psVectorStats(stats, pixels, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + return false; + } + psFree(pixels); + double reduced = psStatsGetValue(stats, statistic); // Result of statistics + + psString comment = NULL; // Comment to add + psStringAppend(&comment, "Overscan value: %f", reduced); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + + // write metadata header value + // XXX EAM : this could / should write the stdev of the overscan region + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", reduced); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); + + psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32); + psBinaryOp(image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use + + // subtract the measured value from each overscan region as well + iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use + } + psFree(iter); + psFree(reducedScalar); + + // EAM 2022.03.29 : if the calculated overscan value is below the threshold, + // declare the readout dead and mask + + if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid)) + { + fprintf(stderr, "bad overscan (1) %f, masking readout\n", reduced); + psImage *mask = input->mask; + for (int y = 0; y < mask->numRows; y++) + { + for (int x = 0; x < mask->numCols; x++) + { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; + } + } + } + + psFree(stats); + return true; + } + + bool mdok = false; + + // We are performing a row-by-row overscan subtraction + int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction + + if ((cellreaddir != 1) && (cellreaddir != 2)) + { + psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n"); + return false; + } + + float chi2 = NAN; // chi^2 from fit + + // adjust operation depending on the read direction : need to re-org pixels for columns + if (!doTwoDOverscan && (cellreaddir == 1)) + { + // The read direction is rows + psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) + { + pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + // XXX HG: investigate whether the y-overscan is being used or the x-overscan + while ((overscan = psListGetAndIncrement(iter))) + { + + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->row0 - image->row0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) + { + int j = i - diff; + // i is row on image + // j is row on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numCols); + values->n += overscan->numCols; + // XXX double-check the range of values->n here + memcpy(&values->data.F32[index], overscan->data.F32[j], + overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + index += overscan->numCols; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels, false); + psFree(pixels); + if (!reduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Subtract row by row + for (int i = 0; i < image->numRows; i++) + { + for (int j = 0; j < image->numCols; j++) + { + image->data.F32[i][j] -= reduced->data.F32[i]; + } + } + + // subtract from the overscan regions + { + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. + int diff = overscan->row0 - image->row0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) + { + int j = i - diff; + // i is row on image + // j is row on overscan + for (int k = 0; k < overscan->numCols; k++) + { + overscan->data.F32[j][k] -= reduced->data.F32[j]; + } + } + } + psFree(iter); + } + psFree(reduced); + } + + if (!doTwoDOverscan && (cellreaddir == 2)) + { + // The read direction is columns + psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) + { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + pixels->data[i] = values; + } + + // Pull the pixels out into the vectors + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. pixels->data represents + // the image row pixels. + int diff = overscan->col0 - image->col0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) + { + int iFixed = i - diff; + // i is column on image + // iFixed is column on overscan + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numRows); + for (int j = 0; j < overscan->numRows; j++) + { + values->data.F32[index++] = overscan->data.F32[j][iFixed]; + } + values->n += overscan->numRows; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = pmOverscanVector(&chi2, overscanOpts->primary, pixels, false); + psFree(pixels); + if (!reduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); + return false; + } + + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, reduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Subtract column by column + for (int j = 0; j < image->numRows; j++) + { + for (int i = 0; i < image->numCols; i++) + { + image->data.F32[j][i] -= reduced->data.F32[i]; + } + } + + // subtract from the overscan regions + { + psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator + psImage *overscan = NULL; // Overscan image from iterator + while ((overscan = psListGetAndIncrement(iter))) + { + // the overscan and image might not be aligned. + int diff = overscan->col0 - image->col0; // Offset between the two regions + for (int i = PS_MAX(0, diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) + { + int j = i - diff; + // i is col on image + // j is col on overscan + for (int k = 0; i < overscan->numRows; k++) + { + overscan->data.F32[k][j] -= reduced->data.F32[j]; + } + } + } + psFree(iter); + } + + psFree(reduced); + } + + // 2D bias subtraction with x-dir readout direction: the + // bias is constructed by combining a 1D pattern in the + // readout direction from the top overscan region and a second + // 1D pattern in the cross direction from the overscan + if (doTwoDOverscan && (cellreaddir == 1)) + { + psAssert(overscanOpts->secondary, "2D overscan subtraction requires OVERSCAN.2D parameters"); + + //parse 2D overscan BISASEC region information + psRegion yRegion = overscanOpts->secondary->biassecslow; + psRegion xRegion = overscanOpts->secondary->biassecfast; + + // Cut out relevant region image. Need to use the parent uncut image + psImage *parent = (psImage *) image->parent; + + // the serial (fast readout) direction is columns (x-dir) + psImage *yscan = psImageSubset(parent, yRegion); // overscan region spanning all rows + psImage *xscan = psImageSubset(parent, xRegion); // overscan region spanning all columns + + // Extract the y-dir overscan vector. The overscan and image might not be aligned: + // diff represents the offset between the rows in the image data and the overscan. + // pixels->data represents the image row pixels. For example, the image region may be + // inset in the y-direction but the overscan could cover the entire y-range + + // The read direction is rows + psArray *yscanPixels = psArrayAlloc(yscan->numRows); // Array of vectors containing pixels + for (int i = 0; i < yscanPixels->n; i++) + { + yscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // XXX this code allows multiple yscans to be appended, but this does not + // match the concept of how they are assigned above: biassec[0] = yscan + // int yDiff = yscan->row0 - image->row0; // Offset between the two regions + for (int i = 0; i < yscanPixels->n; i++) + { + psVector *values = yscanPixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + yscan->numCols); + values->n += yscan->numCols; + // XXX double-check the range of values->n here + memcpy(&values->data.F32[index], yscan->data.F32[i], + yscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); + yscanPixels->data[i] = values; // Update the pointer in case it's moved + } + + // Extract the x-dir overscan vector. The overscan and image might not be aligned: + // diff represents the offset between the rows in the image data and the overscan. + // pixels->data represents the image row pixels. For example, the image region may be + // inset in the x-direction but the overscan could cover the entire y-range + + // Extract the top region as a vector of the columns + psArray *xscanPixels = psArrayAlloc(xscan->numCols); // Array of vectors containing pixels + for (int i = 0; i < xscanPixels->n; i++) + { + xscanPixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); + } + + // int xDiff = xscan->col0 - image->col0; // Offset between the two regions + for (int ix = 0; ix < xscanPixels->n; ix++) + { + psVector *values = xscanPixels->data[ix]; + values = psVectorRealloc(values, xscan->numRows); + values->n = xscan->numRows; + for (int iy = 0; iy < xscan->numRows; iy++) + { + values->data.F32[iy] = xscan->data.F32[iy][ix]; + } + xscanPixels->data[ix] = values; // Update the pointer in case it's moved + } + + // Reduce the overscans + // XXX need to save 2 different chi-square values + psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels, true); + psFree(yscanPixels); + if (!yReduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate y-dir overscan vector.\n"); + return false; + } + if (!pmOverscanUpdateHeaderVector(input, hdu, overscanOpts, yReduced)) + { + psError(PS_ERR_UNKNOWN, false, "failure to update header"); + return false; + } + + // Reduce the overscans + // XXX need to save 2 different chi-square values + psVector *xReduced = pmOverscanVector(&chi2, overscanOpts->secondary, xscanPixels, true); + psFree(xscanPixels); + if (!xReduced) + { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate x-dir overscan vector.\n"); + return false; + } + + // subtract the 2D bias from the image + if (yscan->col0 >= xscan->col0 && yscan->col0 + yscan->numCols <= xscan->col0 + xscan->numCols) + { + // define how to normalize the xReduced vector + // slice the part overlapping withe yoverscan from xReduced + long startIndex = yscan->col0 - xscan->col0; + long normSize = xReduced->n - startIndex; // Calculate the size of the new vector + // Allocate the new vector + psVector *normVector = psVectorAlloc(normSize, PS_TYPE_F32); + // Copy the data from xReduced starting from startIndex to the end + for (long i = startIndex; i < xReduced->n; i++) + { + normVector->data.F32[i - startIndex] = xReduced->data.F32[i]; + } + // Now normVector holds the sliced data; then, compute its robust mean + psStatsOptions statistic = PS_STAT_ROBUST_MEDIAN; + psStats *stats = psStatsAlloc(statistic); + if (!psVectorStats(stats, normVector, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure robust median as the normalization of xReduced"); + psFree(stats); + psFree(normVector); + return NULL; + } + float xReducedNormalized = psStatsGetValue(stats, statistic); + psFree(stats); + psFree(normVector); + // XXX apply the 2D bias correction here + int yDiff = yscan->row0 - image->row0; // y offset between the science and the yoverscan region + int xDiff = xscan->col0 - image->col0; // x offset between the science and the xoverscan region + for (int i = 0; i < image->numRows; i++) + { + for (int j = 0; j < image->numCols; j++) + { + int iy = i + yDiff; + int jx = j + xDiff; + image->data.F32[i][j] -= yReduced->data.F32[iy] - xReducedNormalized + xReduced->data.F32[jx]; + } + } + } + else + { + psError(PS_ERR_UNKNOWN, true, "x dimension of yscan is not fully contained by xscan\n"); + return false; + } + + // subtract the y-dir vector from the y-dir overscan regions (why?) + { + // the overscan and image might not be aligned. + // int diff = yscan->row0 - image->row0; // Offset between the two regions + for (int i = 0; i < yscan->numRows; i++) + { + for (int j = 0; j < yscan->numCols; j++) + { + yscan->data.F32[i][j] -= yReduced->data.F32[i]; + } + } + } + + // subtract the x-dir vector from the x-dir overscan regions (why?) + { + // the overscan and image might not be aligned. + // int diff = xscan->col0 - image->col0; // Offset between the two regions + for (int i = 0; i < xscan->numCols; i++) + { + // int j = i - diff; + // i is column on image (aligned with xReduced) + // j is column on xscan + for (int j = 0; j < xscan->numRows; j++) + { + xscan->data.F32[j][i] -= xReduced->data.F32[i]; + } + } + } + psFree(xReduced); + psFree(yReduced); + } + // pmOverscanUpdateHeader (hdu, overscanOpts, chi2); + return true; + +} // End of overscan subtraction static void pmOverscanOptionsFree(pmOverscanOptions *options) { - psFree(options->stat); - psFree(options->poly); - psFree(options->spline); + psFree(options->primary); + psFree(options->secondary); } -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat, - int boxcar, float gauss) +static void pmOverscanStatOptionsFree(pmOverscanStatOptions *options) { - pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); - psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree); - - // Inputs - opts->single = single; - opts->constant = false; - opts->fitType = fitType; - opts->order = order; - opts->stat = psMemIncrRefCounter(stat); - - opts->minValid = 0.0; // default value if not defined - opts->maxValid = (float) 0x10000; // default value if not defined - opts->maskVal = 0x0001; // default value if not defined - - // Smoothing - opts->boxcar = boxcar; - opts->gauss = gauss; - - // Outputs - opts->poly = NULL; - opts->spline = NULL; - - return opts; + psFree(options->stat); + psFree(options->poly); + psFree(options->spline); } +// Globally defined +pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void) +{ + pmOverscanStatOptions *opts = psAlloc(sizeof(pmOverscanStatOptions)); + psMemSetDeallocator(opts, (psFreeFunc)pmOverscanStatOptionsFree); + + // Inputs + opts->fitType = PM_FIT_NONE; + opts->order = 0; + opts->stat = NULL; + + // Smoothing + opts->boxcar = 0; + opts->gauss = 0.0; + + // Outputs + opts->poly = NULL; + opts->spline = NULL; + + return opts; +} + +// Globally defined +pmOverscanOptions *pmOverscanOptionsAlloc(void) +{ + pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); + psMemSetDeallocator(opts, (psFreeFunc)pmOverscanOptionsFree); + + // Inputs + opts->single = false; + opts->constant = false; + opts->TwoD = false; + + opts->value = 0.0; + + opts->minValid = 0.0; // default value if not defined + opts->maxValid = (float)0x10000; // default value if not defined + opts->maskVal = 0x0001; // default value if not defined + + // stat options + opts->primary = NULL; + opts->secondary = NULL; + + return opts; +} + // Produce an overscan vector from an array of pixels -psVector *pmOverscanVector(float *chi2, // chi^2 from fit - pmOverscanOptions *overscanOpts, // Overscan options - const psArray *pixels, // Array of vectors containing the pixel values - psStats *myStats // Statistic to use in reducing the overscan - ) +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanStatOptions *overscanOpts, // Overscan statistic options + const psArray *pixels, // Array of vectors containing the pixel values + bool robust) // Use robust statistics for boxcar smoothing { - assert(overscanOpts); - assert(pixels); - assert(myStats); - - psStatsOptions statistic = psStatsSingleOption(myStats->options); // Statistic to use - assert(statistic != 0); - - // Reduce the overscans - psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row - psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate - psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting - - for (int i = 0; i < pixels->n; i++) { - psVector *values = pixels->data[i]; // Vector with overscan values - if (values->n > 0) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; - ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] - if (!psVectorStats(myStats, values, NULL, NULL, 0)) { + assert(overscanOpts); + assert(pixels); + + // statisctic to be calculated + psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use + if (!statistic) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", overscanOpts->stat); + return false; + } + psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original + + // Reduce the overscans + psVector *reduced = psVectorAlloc(pixels->n, PS_TYPE_F32); // Overscan for each row + psVector *ordinate = psVectorAlloc(pixels->n, PS_TYPE_F32); // Ordinate + psVector *mask = psVectorAlloc(pixels->n, PS_TYPE_VECTOR_MASK); // Mask for fitting + + for (int i = 0; i < pixels->n; i++) + { + psVector *values = pixels->data[i]; // Vector with overscan values + if (values->n > 0) + { + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; + ordinate->data.F32[i] = 2.0 * (float)i / (float)pixels->n - 1.0; // Scale to [-1,1] + if (!psVectorStats(stats, values, NULL, NULL, 0)) + { + psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); + goto escape; + } + reduced->data.F32[i] = psStatsGetValue(stats, statistic); + } + else + { + if (overscanOpts->fitType == PM_FIT_NONE) + { + psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the image, and no fit is requested.\n"); + goto escape; + } + else + { + // We'll fit this one out + mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; + } + } + } + // Smooth the reduced vector + if (overscanOpts->boxcar > 0) + { + psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar, robust); // Smoothed vector + psFree(reduced); + reduced = smoothed; + } + if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) + { + if (overscanOpts->boxcar > 0) + { + psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); + } + psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); + psFree(reduced); + reduced = smoothed; + } + + // Fit the overscan, if required + psVector *fitted = NULL; // Fitted overscan values + switch (overscanOpts->fitType) + { + case PM_FIT_NONE: + // No fitting --- that's easy. + fitted = psMemIncrRefCounter(reduced); + break; + case PM_FIT_POLY_ORD: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) + { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (!overscanOpts->poly) + { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_POLY_CHEBY: + if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || + overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) + { + psFree(overscanOpts->poly); + overscanOpts->poly = NULL; + } + if (!overscanOpts->poly) + { + overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); + } + psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); + fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); + break; + case PM_FIT_SPLINE: + + // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes + // a knot for every input point. it needs an argument like 'number of knots' for the + // output spline. EAM: still true 2023.01.22 + + // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); + // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); + psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n"); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); + goto escape; + } + + if (chi2) + { + *chi2 = 0.0; // chi^2 (sort of) + for (int i = 0; i < reduced->n; i++) + { + *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); + } + } + + psFree(reduced); + psFree(ordinate); + psFree(mask); + psFree(stats); + return fitted; + +escape: + psFree(reduced); + psFree(ordinate); + psFree(mask); + psFree(stats); + return NULL; +} + +// XXX EAM 2024.04.13 : this function seems poorly conceived. it is replacing the stats from +// the reduced overscan vector with the 0-order element from the polynomial fit. Not clear is +// adding any useful information +bool pmOverscanUpdateHeader(pmHDU *hdu, pmOverscanStatOptions *overscanOpts, float chi2) +{ + + psString comment = NULL; // Comment to add + + switch (overscanOpts->fitType) + { + case PM_FIT_POLY_ORD: + case PM_FIT_POLY_CHEBY: + { + psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); + psPolynomial1D *poly = overscanOpts->poly; // The polynomial + for (int i = 0; i < poly->nX; i++) + { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + comment = NULL; + + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", poly->coeff[0]); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", poly->coeffErr[0]); + break; + } + case PM_FIT_SPLINE: + { + /* + psSpline1D *spline = overscanOpts->spline; // The spline + for (int i = 0; i < spline->n; i++) { + psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); + psPolynomial1D *poly = spline->spline[i]; // i-th polynomial + for (int j = 0; j < poly->nX; j++) { + psStringAppend(&comment, "%.1f ", poly->coeff[i]); + } + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + comment, ""); + psFree(comment); + comment = NULL; + } + */ + // write metadata header value + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, + "Overscan value", NAN); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, + "Overscan stdev", NAN); + break; + } + case PM_FIT_NONE: + break; + default: + psAbort("Should never get here!!!\n"); + } + return true; +} + +// generate stats of overscan vector for header +// reduced: 1D vector with overscan stats +bool pmOverscanUpdateHeaderVector(pmReadout *input, pmHDU *hdu, pmOverscanOptions *overscanOpts, psVector *reduced) +{ + + psString comment = NULL; // Comment to add + psStats *vectorStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + if (!psVectorStats(vectorStats, reduced, NULL, NULL, 0)) + { psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); return false; - } - reduced->data.F32[i] = psStatsGetValue(myStats, statistic); - } else if (overscanOpts->fitType == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "The overscan is not supplied for all points on the " - "image, and no fit is requested.\n"); - return NULL; - } else { - // We'll fit this one out - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 1; - } - } - - // Smooth the reduced vector - if (overscanOpts->boxcar > 0) { - psVector *smoothed = psVectorBoxcar(NULL, reduced, overscanOpts->boxcar); // Smoothed vector - psFree(reduced); - reduced = smoothed; - } - if (isfinite(overscanOpts->gauss) && overscanOpts->gauss > 0) { - if (overscanOpts->boxcar > 0) { - psWarning("Gaussian smoothing the boxcar smoothed overscan --- you asked for it."); - } - psVector *smoothed = psVectorSmooth(NULL, reduced, overscanOpts->gauss, SMOOTH_NSIGMA); - psFree(reduced); - reduced = smoothed; - } - - // Fit the overscan, if required - psVector *fitted = NULL; // Fitted overscan values - switch (overscanOpts->fitType) { - case PM_FIT_NONE: - // No fitting --- that's easy. - fitted = psMemIncrRefCounter(reduced); - break; - case PM_FIT_POLY_ORD: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_ORD)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_POLY_CHEBY: - if (overscanOpts->poly && (overscanOpts->poly->nX != overscanOpts->order || - overscanOpts->poly->type != PS_POLYNOMIAL_CHEB)) { - psFree(overscanOpts->poly); - overscanOpts->poly = NULL; - } - if (! overscanOpts->poly) { - overscanOpts->poly = psPolynomial1DAlloc(PS_POLYNOMIAL_CHEB, overscanOpts->order); - } - psVectorFitPolynomial1D(overscanOpts->poly, mask, 1, reduced, NULL, ordinate); - fitted = psPolynomial1DEvalVector(overscanOpts->poly, ordinate); - break; - case PM_FIT_SPLINE: - - // XXX I don't think psSpline1D is up to scratch yet --- it has no mask, and it assumes - // a knot for every input point. it needs an argument like 'number of knots' for the - // output spline. EAM: still true 2023.01.22 - - // overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); - // fitted = psSpline1DEvalVector(overscanOpts->spline, ordinate); - psError(PS_ERR_UNKNOWN, true, "Spline overscan fitting is broken\n"); - break; - default: - psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); - return NULL; - break; - } - - if (chi2) { - *chi2 = 0.0; // chi^2 (sort of) - for (int i = 0; i < reduced->n; i++) { - *chi2 += PS_SQR(fitted->data.F32[i] - reduced->data.F32[i]); - } - } - - psFree(reduced); - psFree(ordinate); - psFree(mask); - - return fitted; -} - -bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2) { - - psString comment = NULL; // Comment to add - - switch (overscanOpts->fitType) { - case PM_FIT_POLY_ORD: - case PM_FIT_POLY_CHEBY: { - psStringAppend(&comment, "Overscan fit (chi2: %.2f): ", chi2); - psPolynomial1D *poly = overscanOpts->poly; // The polynomial - for (int i = 0; i < poly->nX; i++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - comment = NULL; - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, - "Overscan value", poly->coeff[0]); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, - "Overscan stdev", poly->coeffErr[0]); - break; - } - case PM_FIT_SPLINE: { - /* - psSpline1D *spline = overscanOpts->spline; // The spline - for (int i = 0; i < spline->n; i++) { - psStringAppend(&comment, "Overscan fit (chi2: %.2f) %d:", chi2, i); - psPolynomial1D *poly = spline->spline[i]; // i-th polynomial - for (int j = 0; j < poly->nX; j++) { - psStringAppend(&comment, "%.1f ", poly->coeff[i]); - } - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - comment, ""); - psFree(comment); - comment = NULL; - } - */ - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, - "Overscan value", NAN); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, - "Overscan stdev", NAN); - break; - } - case PM_FIT_NONE: - break; - default: - psAbort("Should never get here!!!\n"); - } - return true; -} - -bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts) { - - assert (input); - - if (overscanOpts == NULL) return true; // no overscan subtraction requested - - pmHDU *hdu = pmHDUFromReadout(input); // HDU of interest - psImage *image = input->image; - - // check for 'soft bias' (simple, fixed offset to be subtracted) - if (overscanOpts->constant) { - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", - overscanOpts->value); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); - - // NOTE psBinaryOp frees arg2 if it is a scalar - (void)psBinaryOp(input->image, input->image, "-", psScalarAlloc((float)overscanOpts->value, PS_TYPE_F32)); - - return true; - } - - // we are performing a statitical analysis of the overscan region - - // Check for an unallowable pmFit. - if (overscanOpts->fitType != PM_FIT_NONE && overscanOpts->fitType != PM_FIT_POLY_ORD && - overscanOpts->fitType != PM_FIT_POLY_CHEBY && overscanOpts->fitType != PM_FIT_SPLINE) { - psError(PS_ERR_UNKNOWN, true, "Invalid fit type (%d). Returning original image.\n", - overscanOpts->fitType); - return false; - } - - psList *overscans = input->bias; // List of the overscan images - - psStatsOptions statistic = psStatsSingleOption(overscanOpts->stat->options); // Statistic to use - if (statistic == 0) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Multiple or no statistics options set: %p\n", - overscanOpts->stat); - return false; - } - psStats *stats = psStatsAlloc(statistic); // A new psStats, to avoid clobbering original - - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Subtracting overscan (stat %x; type %x; order %d)", - statistic, overscanOpts->fitType, overscanOpts->order); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, - comment, ""); - psTrace ("psModules.detrend", 4, "%s\n", comment); - psFree(comment); - - // Reduce all overscan pixels to a single value - if (overscanOpts->single) { - psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - int index = pixels->n; // Index - pixels = psVectorRealloc(pixels, pixels->n + overscan->numRows * overscan->numCols); - pixels->n += overscan->numRows * overscan->numCols; - for (int i = 0; i < overscan->numRows; i++) { - memcpy(&pixels->data.F32[index], overscan->data.F32[i], - overscan->numCols * sizeof(psF32)); - index += overscan->numCols; - } - } - psFree(iter); - - if (!psVectorStats(stats, pixels, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - psFree(pixels); - double reduced = psStatsGetValue(stats, statistic); // Result of statistics - - psString comment = NULL; // Comment to add - psStringAppend(&comment, "Overscan value: %f", reduced); + } + psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); psFree(comment); // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan value", - reduced); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", NAN); - - psScalar *reducedScalar = psScalarAlloc(reduced, PS_TYPE_F32); - psBinaryOp (image, image, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use - - // subtract the measured value from each overscan region as well - iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - psBinaryOp(overscan, overscan, "-", psMemIncrRefCounter(reducedScalar)); // NOTE: psBinaryOp frees arg2 if it a scalar, so we need to bump to re-use - } - psFree(iter); - psFree(reducedScalar); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean); + psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev); // EAM 2022.03.29 : if the calculated overscan value is below the threshold, // declare the readout dead and mask - if ((reduced < overscanOpts->minValid) || (reduced > overscanOpts->maxValid)) { - fprintf (stderr, "bad overscan (1) %f, masking readout\n", reduced); - psImage *mask = input->mask; - for (int y = 0; y < mask->numRows; y++) { - for (int x = 0; x < mask->numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; - } - } - } - - psFree(stats); + if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) + { + fprintf(stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean); + psImage *mask = input->mask; + for (int y = 0; y < mask->numRows; y++) + { + for (int x = 0; x < mask->numCols; x++) + { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; + } + } + } + + psFree(vectorStats); return true; - } - - bool mdok = false; - - // We are performing a row-by-row overscan subtraction - int cellreaddir = psMetadataLookupS32(&mdok, input->parent->concepts, "CELL.READDIR"); // Read direction - if ((cellreaddir != 1) && (cellreaddir != 2)) { - psError(PS_ERR_UNKNOWN, true, "CELL.READDIR must be 1 (rows) or 2 (cols)\n"); - return false; - } - - float chi2 = NAN; // chi^2 from fit - - // adjust operation depending on the read direction : need to re-org pixels for columns - if (cellreaddir == 1) { - // The read direction is rows - psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - pixels->data[i] = psVectorAlloc(0, PS_TYPE_F32); - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image - // j is row on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numCols); - values->n += overscan->numCols; - memcpy(&values->data.F32[index], overscan->data.F32[j], - overscan->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); - index += overscan->numCols; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - psFree(stats); - return false; - } - - // generate stats of overscan vector for header - { - psString comment = NULL; // Comment to add - psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev); - - // EAM 2022.03.29 : if the calculated overscan value is below the threshold, - // declare the readout dead and mask - - if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) { - fprintf (stderr, "bad overscan (2) %f, masking readout\n", vectorStats->sampleMean); - psImage *mask = input->mask; - for (int y = 0; y < mask->numRows; y++) { - for (int x = 0; x < mask->numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; - } - } - } - - psFree (vectorStats); - } - - // Subtract row by row - for (int i = 0; i < image->numRows; i++) { - for (int j = 0; j < image->numCols; j++) { - image->data.F32[i][j] -= reduced->data.F32[i]; - } - } - - // subtract from the overscan regions - { - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. - int diff = overscan->row0 - image->row0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numRows, overscan->numRows + diff); i++) { - int j = i - diff; - // i is row on image - // j is row on overscan - for (int k = 0; k < overscan->numCols; k++) { - overscan->data.F32[j][k] -= reduced->data.F32[j]; - } - } - } - psFree(iter); - } - psFree(reduced); - } - if (cellreaddir == 2) { - // The read direction is columns - psArray *pixels = psArrayAlloc(image->numCols); // Array of vectors containing pixels - for (int i = 0; i < pixels->n; i++) { - psVector *values = psVectorAlloc(0, PS_TYPE_F32); - pixels->data[i] = values; - } - - // Pull the pixels out into the vectors - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. pixels->data represents - // the image row pixels. - int diff = overscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { - int iFixed = i - diff; - // i is column on image - // iFixed is column on overscan - psVector *values = pixels->data[i]; - int index = values->n; // Index in the vector - values = psVectorRealloc(values, values->n + overscan->numRows); - for (int j = 0; j < overscan->numRows; j++) { - values->data.F32[index++] = overscan->data.F32[j][iFixed]; - } - values->n += overscan->numRows; - pixels->data[i] = values; // Update the pointer in case it's moved - } - } - psFree(iter); - - // Reduce the overscans - psVector *reduced = pmOverscanVector(&chi2, overscanOpts, pixels, stats); - psFree(pixels); - if (! reduced) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to generate overscan vector.\n"); - psFree(stats); - return false; - } - - // generate stats of overscan vector for header - { - psString comment = NULL; // Comment to add - psStats *vectorStats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - if (!psVectorStats (vectorStats, reduced, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - psStringAppend(&comment, "Mean Overscan value: %f", vectorStats->sampleMean); - psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); - psFree(comment); - - // write metadata header value - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_VAL", PS_META_REPLACE, "Overscan mean", vectorStats->sampleMean); - psMetadataAddF32(hdu->header, PS_LIST_TAIL, "OVER_SIG", PS_META_REPLACE, "Overscan stdev", vectorStats->sampleStdev); - - // EAM 2022.03.29 : if the calculated overscan value is below the threshold, - // declare the readout dead and mask - - if ((vectorStats->sampleMean < overscanOpts->minValid) || (vectorStats->sampleMean > overscanOpts->maxValid)) { - fprintf (stderr, "bad overscan (3) %f, masking readout\n", vectorStats->sampleMean); - psImage *mask = input->mask; - for (int y = 0; y < mask->numRows; y++) { - for (int x = 0; x < mask->numCols; x++) { - mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= overscanOpts->maskVal; - } - } - } - psFree (vectorStats); - } - - // Subtract column by column - for (int j = 0; j < image->numRows; j++) { - for (int i = 0; i < image->numCols; i++) { - image->data.F32[j][i] -= reduced->data.F32[i]; - } - } - - // subtract from the overscan regions - { - psListIterator *iter = psListIteratorAlloc(overscans, PS_LIST_HEAD, false); // Iterator - psImage *overscan = NULL; // Overscan image from iterator - while ((overscan = psListGetAndIncrement(iter))) { - // the overscan and image might not be aligned. - int diff = overscan->col0 - image->col0; // Offset between the two regions - for (int i = PS_MAX(0,diff); i < PS_MIN(image->numCols, overscan->numCols + diff); i++) { - int j = i - diff; - // i is col on image - // j is col on overscan - for (int k = 0; i < overscan->numRows; k++) { - overscan->data.F32[k][j] -= reduced->data.F32[j]; - } - } - } - psFree(iter); - } - - psFree(reduced); - } - - pmOverscanUpdateHeader (hdu, overscanOpts, chi2); - psFree(stats); - - return true; - -} // End of overscan subtraction - +} Index: trunk/psModules/src/detrend/pmOverscan.h =================================================================== --- trunk/psModules/src/detrend/pmOverscan.h (revision 42888) +++ trunk/psModules/src/detrend/pmOverscan.h (revision 42889) @@ -18,9 +18,10 @@ /// Type of fit to perform -typedef enum { - PM_FIT_NONE, ///< No fit - PM_FIT_POLY_ORD, ///< Fit ordinary polynomial - PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial - PM_FIT_SPLINE ///< Fit cubic splines +typedef enum +{ + PM_FIT_NONE, ///< No fit + PM_FIT_POLY_ORD, ///< Fit ordinary polynomial + PM_FIT_POLY_CHEBY, ///< Fit Chebyshev polynomial + PM_FIT_SPLINE ///< Fit cubic splines } pmFit; @@ -35,42 +36,56 @@ { // Inputs - bool single; ///< Reduce all overscan regions to a single value? - bool constant; ///< use a supplied constant value (do not measure region) - float value; ///< supplied value if needed (per above) - pmFit fitType; ///< Type of fit to overscan - unsigned int order; ///< Order of polynomial, or number of spline pieces - psStats *stat; ///< Statistic to use when reducing the minor direction - int boxcar; ///< Boxcar smoothing radius - float gauss; ///< Gaussian smoothing sigma - float minValid; ///< if overscan is too low, readout is dead : mask - float maxValid; ///< if overscan is too high, readout is dead : mask - psImageMaskType maskVal; ///< Mask value to give dead readouts + pmFit fitType; ///< Type of fit to overscan + unsigned int order; ///< Order of polynomial, or number of spline pieces + psStats *stat; ///< Statistic to use when reducing the minor direction + int boxcar; ///< Boxcar smoothing radius + float gauss; ///< Gaussian smoothing sigma + psRegion biassecslow; ///< The 2D slow bias section + psRegion biassecfast; ///< The 2D fast bias section - // Outputs - psPolynomial1D *poly; ///< Result of polynomial fit - psSpline1D *spline; ///< Result of spline fit -} -pmOverscanOptions; + + // These are used to carry the fitted functions, but are only used to generate + // an updated reduced vector + psPolynomial1D *poly; ///< Result of polynomial fit + psSpline1D *spline; ///< Result of spline fit +} pmOverscanStatOptions; + +/// Options for overscan subtraction +/// +/// The overscan subtraction may be performed by reducing all overscan regions to a single value (e.g., if +/// there is no structure); or the overscan may be fit perpendicular to the read direction (usually the +/// columns) with a particular functional form; or a single value may be subtracted for each read/scan without +/// fitting (if the structure defies characterisation). In any case, statistics are required to reduce +/// multiple values to a single value (either for the scan, or for the entire overscan regions). +typedef struct +{ + // Inputs + bool single; ///< Reduce all overscan regions to a single value? + bool constant; ///< use a supplied constant value (do not measure region) + bool TwoD; ///< use a supplied constant value (do not measure region) + float value; ///< supplied value if needed (per above) + float minValid; ///< if overscan is too low, readout is dead : mask + float maxValid; ///< if overscan is too high, readout is dead : mask + psImageMaskType maskVal; ///< Mask value to give dead readouts + + pmOverscanStatOptions *primary; + pmOverscanStatOptions *secondary; +} pmOverscanOptions; /// Allocator for overscan options -pmOverscanOptions *pmOverscanOptionsAlloc(bool single, ///< Reduce all overscan regions to a single value? - pmFit fitType, ///< Type of fit to overscan - unsigned int order, ///< Order of polynomial, or number of splines - psStats *stat, ///< Statistic to use - int boxcar, ///< Boxcar smoothing radius - float gauss ///< Gaussian smoothing sigma - ); +pmOverscanOptions *pmOverscanOptionsAlloc(void); -psVector *pmOverscanVector(float *chi2, // chi^2 from fit - pmOverscanOptions *overscanOpts, // Overscan options - const psArray *pixels, // Array of vectors containing the pixel values - psStats *myStats // Statistic to use in reducing the overscan - ); +pmOverscanStatOptions *pmOverscanStatOptionsAlloc(void); -bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2); +psVector *pmOverscanVector(float *chi2, // chi^2 from fit + pmOverscanStatOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + bool robust // use robust statistics for boxcar smoothing +); -bool pmOverscanSubtract (pmReadout *input, pmOverscanOptions *overscanOpts); +// bool pmOverscanUpdateHeader (pmHDU *hdu, pmOverscanOptions *overscanOpts, float chi2); + +bool pmOverscanSubtract(pmReadout *input, pmOverscanOptions *overscanOpts, bool doTwoDOverscan); /// @} #endif - Index: trunk/psModules/src/detrend/pmPattern.c =================================================================== --- trunk/psModules/src/detrend/pmPattern.c (revision 42888) +++ trunk/psModules/src/detrend/pmPattern.c (revision 42889) @@ -239,4 +239,5 @@ # define READNOISE 10 /* arbitrary number */ float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + if(isnan(sigma)) sigma=20.; float delta = PS_MIN (thresh * sigma, 2*nominalAmplitude); float lower = stats->robustMedian - delta; // Lower bound for data @@ -246,7 +247,13 @@ // if the noise from the background is too large float significance = nominalAmplitude / sigma; + significance = abs(significance); // XXX EAM : arbitrary number if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + return true; + } + //Add some safeties in case of bad fits + if (isnan(significance) || (background < 0.0)) { psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); return true; @@ -554,4 +561,5 @@ # define READNOISE 10 /* arbitrary number */ float sigma = sqrt(stats->robustMedian + PS_SQR(READNOISE)); + if(isnan(sigma)) sigma=20.; // XXX EAM : somewhat arbitrary number float delta = PS_MIN (thresh * sigma, 5*nominalAmplitude); float lower = stats->robustMedian - delta; // Lower bound for data @@ -561,5 +569,6 @@ // if the noise from the background is too large float significance = nominalAmplitude / sigma; - + significance = abs(significance); + // XXX EAM : arbitrary number if (isfinite(nominalAmplitude) && (significance < 1.0/6.0)) { @@ -569,4 +578,12 @@ return true; } + //Add some safeties in case of bad fits + if (isnan(significance) || (background < 0.0)) { + psLogMsg("ppImage", PS_LOG_INFO, "Skipping row pattern correction due to bad fit for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); + psFree(stats); + psFree(rng); + return true; + } + psLogMsg("ppImage", PS_LOG_INFO, "Performing row pattern correction for %s, %s, stats: %f - %f - %f : %f %f %f\n", chipName, cellName, lower, background, upper, sigma, nominalAmplitude, significance); # endif