Index: trunk/psModules/src/imsubtract/pmSubtractBias.c =================================================================== --- trunk/psModules/src/imsubtract/pmSubtractBias.c (revision 6511) +++ trunk/psModules/src/imsubtract/pmSubtractBias.c (revision 6873) @@ -1,2 +1,7 @@ +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// XXX WARNING: I have completely replaced this file with an OLD VERSION (that works) instead of the +// one that was being worked on. +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + /** @file pmSubtractBias.c * @@ -6,28 +11,28 @@ * @author GLG, MHPCC * - * @version $Revision: 1.11 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-04 01:01:33 $ + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-17 18:10:08 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii * */ -/*****************************************************************************/ -/* INCLUDE FILES */ -/*****************************************************************************/ -#include -#include -#include -#include "pslib.h" + #if HAVE_CONFIG_H #include #endif + +#include #include "pmSubtractBias.h" -/*****************************************************************************/ -/* DEFINE STATEMENTS */ -/*****************************************************************************/ +#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2 +#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3 + + +#define MAX(a,b) ((a) > (b) ? (a) : (b)) +#define MIN(a,b) ((a) < (b) ? (a) : (b)) + + // XXX: put these in psConstants.h -void PS_POLY1D_PRINT( - psPolynomial1D *poly) +void PS_POLY1D_PRINT(psPolynomial1D *poly) { printf("-------------- PS_POLY1D_PRINT() --------------\n"); @@ -57,105 +62,98 @@ }\ -/*****************************************************************************/ -/* TYPE DEFINITIONS */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* GLOBAL VARIABLES */ -/*****************************************************************************/ -psS32 currentId = 0; // XXX: remove -psS32 memLeaks = 0; // XXX: remove -//PRINT_MEMLEAKS(8); XXX -/*****************************************************************************/ -/* FILE STATIC VARIABLES */ -/*****************************************************************************/ - -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - LOCAL */ -/*****************************************************************************/ + +void overscanOptionsFree(pmOverscanOptions *options) +{ + psFree(options->stat); + psFree(options->poly); + psFree(options->spline); +} + +pmOverscanOptions *pmOverscanOptionsAlloc(bool single, pmFit fitType, unsigned int order, psStats *stat) +{ + pmOverscanOptions *opts = psAlloc(sizeof(pmOverscanOptions)); + psMemSetDeallocator(opts, (psFreeFunc)overscanOptionsFree); + + // Inputs + opts->single = single; + opts->fitType = fitType; + opts->order = order; + opts->stat = psMemIncrRefCounter(stat); + + // Outputs + opts->poly = NULL; + opts->spline = NULL; + + return opts; +} + /****************************************************************************** -psSubtractFrame(): this routine will take as input the pmReadout for the input -image and a pmReadout for the bias image. The bias image is subtracted in -place from the input image. We assume that sizes and types are checked -elsewhere. - -XXX: Verify that the image and readout offsets are being used the right way. - -XXX: Ensure that it does the correct thing with image size. +psSubtractFrame(): this routine will take as input a readout for the input +image and a readout for the bias image. The bias image is subtracted in +place from the input image. *****************************************************************************/ -static pmReadout *SubtractFrame( - pmReadout *in, - const pmReadout *bias) -{ - // XXX: When did the ->row0 and ->col0 offsets get coded? - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - bias->image->data.F32[i+in->row0-bias->row0][j+in->col0-bias->col0]; - - if ((in->mask != NULL) && (bias->mask != NULL)) { - (in->mask->data.U8[i][j])|= - bias->mask->data.U8[i+in->row0-bias->row0][j+in->col0-bias->col0]; +static bool SubtractFrame(pmReadout *in,// Input readout + const pmReadout *sub, // Readout to be subtracted from input + float scale // Scale to apply before subtracting + ) +{ + assert(in); + assert(sub); + + psImage *inImage = in->image; // The input image + psImage *inMask = in->mask; // The input mask + psImage *subImage = sub->image; // The image to be subtracted + psImage *subMask = sub->mask; // The mask for the subtraction image + + // Offsets of the cells + int x0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.X0"); + int y0in = psMetadataLookupS32(NULL, in->parent->concepts, "CELL.Y0"); + int x0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.X0"); + int y0sub = psMetadataLookupS32(NULL, sub->parent->concepts, "CELL.Y0"); + + if ((inImage->numCols + x0in - x0sub) > subImage->numCols) { + psError(PS_ERR_UNKNOWN, true, "Image does not have enough columns for subtraction.\n"); + return false; + } + if ((inImage->numRows + y0in - y0sub) > subImage->numRows) { + psError(PS_ERR_UNKNOWN, true, "Image does not have enough rows for subtraction.\n"); + return false; + } + + if (scale == 1.0) { + for (int i = 0; i < inImage->numRows; i++) { + for (int j = 0; j < inImage->numCols; j++) { + inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub]; + if (inMask && subMask) { + inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; + } } } - } - - return(in); -} - - -/****************************************************************************** -psSubtractDarkFrame(): this routine will take as input the pmReadout for the -input image and a pmReadout for the dark image. The dark image is scaled and -subtracted in place from the input image. - -XXX: Verify that the image and readout offsets are being used the right way. - -XXX: Ensure that it does the correct thing with image size. -*****************************************************************************/ -static pmReadout *SubtractDarkFrame( - pmReadout *in, - const pmReadout *dark, - psF32 scale) -{ - // XXX: When did the ->row0 and ->col0 offsets get coded? - if (fabs(scale) > FLT_EPSILON) { - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - (scale * dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]); - - if ((in->mask != NULL) && (dark->mask != NULL)) { - (in->mask->data.U8[i][j])|= - dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; + } else { + for (int i = 0; i < inImage->numRows; i++) { + for (int j = 0; j < inImage->numCols; j++) { + inImage->data.F32[i][j] -= subImage->data.F32[i+y0in-y0sub][j+x0in-x0sub] * scale; + if (inMask && subMask) { + inMask->data.U8[i][j] |= subMask->data.U8[i+y0in-y0sub][j+x0in-x0sub]; } } } - } else { - for (psS32 i=0;iimage->numRows;i++) { - for (psS32 j=0;jimage->numCols;j++) { - in->image->data.F32[i][j]-= - dark->image->data.F32[i+in->row0-dark->row0][j+in->col0-dark->col0]; - - if ((in->mask != NULL) && (dark->mask != NULL)) { - (in->mask->data.U8[i][j])|= - dark->mask->data.U8[i+in->row0-dark->row0][j+in->col0-dark->col0]; - } - } - } - } - - return(in); -} - + } + + return true; +} + + +#if 0 /****************************************************************************** ImageSubtractScalar(): subtract a scalar from the input image. -XXX: Is there a psLib function for this? +XXX: Use a psLib function for this. + +XXX: This should *****************************************************************************/ -static psImage *ImageSubtractScalar( - psImage *image, - psF32 scalar) +static psImage *ImageSubtractScalar(psImage *image, + psF32 scalar) { for (psS32 i=0;inumRows;i++) { @@ -166,4 +164,5 @@ return(image); } +#endif /****************************************************************************** @@ -179,12 +178,4 @@ psStatsOptions opt = 0; - /* - if (stat->options & PS_STAT_ROBUST_MODE) { - if (numOptions == 0) { - opt = PS_STAT_ROBUST_MODE; - } - numOptions++; - } - */ if (stat->options & PS_STAT_ROBUST_MEDIAN) { if (numOptions == 0) { @@ -194,11 +185,4 @@ } - if (stat->options & PS_STAT_FITTED_MEAN) { - if (numOptions == 0) { - opt = PS_STAT_FITTED_MEAN; - } - numOptions++; - } - if (stat->options & PS_STAT_CLIPPED_MEAN) { if (numOptions == 0) { @@ -222,5 +206,5 @@ if (numOptions == 0) { - psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n"); + psError(PS_ERR_UNKNOWN,true, "No statistics options have been specified.\n"); } if (numOptions != 1) { @@ -231,97 +215,7 @@ } -/****************************************************************************** -Polynomial1DCopy(): This private function copies the members of the existing -psPolynomial1D "in" into the existing psPolynomial1D "out". The previous -members of the existing psPolynomial1D "out" are psFree'ed. - *****************************************************************************/ -static psBool Polynomial1DCopy( - psPolynomial1D *out, - psPolynomial1D *in) -{ - psFree(out->coeff); - psFree(out->coeffErr); - psFree(out->mask); - - out->type = in->type; - out->nX = in->nX; - - out->coeff = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->coeff[i] = in->coeff[i]; - } - - out->coeffErr = (psF64 *) psAlloc((in->nX + 1) * sizeof(psF64)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->coeffErr[i] = in->coeffErr[i]; - } - - out->mask = (psMaskType *) psAlloc((in->nX + 1) * sizeof(psMaskType)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->nX + 1) ; i++) { - out->mask[i] = in->mask[i]; - } - - return(true); -} - -/****************************************************************************** -Polynomial1DDup(): This private function duplicates and then returns the input -psPolynomial1D "in". - *****************************************************************************/ -static psPolynomial1D *Polynomial1DDup( - psPolynomial1D *in) -{ - psPolynomial1D *out = psPolynomial1DAlloc(in->type, in->nX); - Polynomial1DCopy(out, in); - return(out); -} - - -/****************************************************************************** -SplineCopy(): This private function copies the members of the existing -psSpline in into the existing psSpline out. - *****************************************************************************/ -static psBool SplineCopy( - psSpline1D *out, - psSpline1D *in) -{ - PS_ASSERT_PTR_NON_NULL(out, false); - PS_ASSERT_PTR_NON_NULL(in, false); - - for (psS32 i = 0 ; i < out->n ; i++) { - psFree(out->spline[i]); - } - psFree(out->spline); - psFree(out->knots); - psFree(out->p_psDeriv2); - - out->n = in->n; - out->spline = (psPolynomial1D **) psAlloc(in->n * sizeof(psPolynomial1D *)); - for (psS32 i = 0 ; i < in->n ; i++) { - out->spline[i] = Polynomial1DDup(in->spline[i]); - } - - // XXX: use psVectorCopy if they get it working. - out->knots = psVectorAlloc(in->knots->n, in->knots->type.type); - out->knots->n = out->knots->nalloc; - for (psS32 i = 0 ; i < in->knots->n ; i++) { - out->knots->data.F32[i] = in->knots->data.F32[i]; - } - /* - out->knots = psVectorCopy(out->knots, in->knots, in->knots->type.type); - */ - - out->p_psDeriv2 = (psF32 *) psAlloc((in->n + 1) * sizeof(psF32)); - // XXX: use memcpy - for (psS32 i = 0 ; i < (in->n + 1) ; i++) { - out->p_psDeriv2[i] = in->p_psDeriv2[i]; - } - - return(true); -} - + + +#if 0 /****************************************************************************** ScaleOverscanVector(): this routine takes as input an arbitrary vector, @@ -330,14 +224,16 @@ PM_FIT_POLYNOMIAL: fit a polynomial to the entire input vector data. PM_FIT_SPLINE: fit splines to the input vector data. -The resulting spline or polynomial is set in the fitSpec argument. +XXX: Doesn't it make more sense to do polynomial interpolation on a few +elements of the input vector, rather than fit a polynomial to the entire +vector? *****************************************************************************/ -static psVector *ScaleOverscanVector( - psVector *overscanVector, - psS32 n, - void *fitSpec, - pmFit fit) +static psVector *ScaleOverscanVector(psVector *overscanVector, + psS32 n, + void *fitSpec, + pmFit fit) { psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4, "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n); + // PS_VECTOR_PRINT_F32(overscanVector); if (NULL == overscanVector) { @@ -347,21 +243,25 @@ // Allocate the new vector. psVector *newVec = psVectorAlloc(n, PS_TYPE_F32); - newVec->n = newVec->nalloc; + // // If the new vector is the same size as the old, simply copy the data. // if (n == overscanVector->n) { - return(psVectorCopy(newVec, overscanVector, PS_TYPE_F32)); - } + for (psS32 i = 0 ; i < n ; i++) { + newVec->data.F32[i] = overscanVector->data.F32[i]; + } + return(newVec); + } + psPolynomial1D *myPoly; + psSpline1D *mySpline; psF32 x; - + psS32 i; if (fit == PM_FIT_POLYNOMIAL) { // Fit a polynomial to the old overscan vector. - psPolynomial1D *myPoly = (psPolynomial1D *) fitSpec; + myPoly = (psPolynomial1D *) fitSpec; PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY1D(myPoly, NULL); myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the psVector.\n"); + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(1): Could not fit a polynomial to the psVector.\n"); return(NULL); } @@ -370,9 +270,17 @@ // of the old vector, use the fitted polynomial to determine the // interpolated value at that point, and set the new vector. - for (psS32 i=0;in) / ((psF32) n); newVec->data.F32[i] = psPolynomial1DEval(myPoly, x); } } else if (fit == PM_FIT_SPLINE) { + psS32 mustFreeSpline = 0; + // Fit a spline to the old overscan vector. + mySpline = (psSpline1D *) fitSpec; + // XXX: Does it make any sense to have a psSpline argument? + if (mySpline == NULL) { + mustFreeSpline = 1; + } + // // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL, @@ -380,26 +288,26 @@ // properly when doing the spline eval. // - psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector); + // mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL); + mySpline = psVectorFitSpline1D(NULL, overscanVector); if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n"); + psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector()(2): Could not fit a spline to the psVector.\n"); return(NULL); } + // PS_PRINT_SPLINE(mySpline); // For each element of the new vector, convert the x-ordinate to that - // of the old vector, use the fitted spline to determine the + // of the old vector, use the fitted polynomial to determine the // interpolated value at that point, and set the new vector. - for (psS32 i=0;in - 1] x = ((psF32) i) * ((psF32) (overscanVector->n-1)) / ((psF32) n); newVec->data.F32[i] = psSpline1DEval(mySpline, x); } - - psSpline1D *ptrSpline = (psSpline1D *) fitSpec; - if (ptrSpline != NULL) { - // Copy the resulting spline fit into ptrSpline. - PS_ASSERT_SPLINE(ptrSpline, NULL); - SplineCopy(ptrSpline, mySpline); - } - psFree(mySpline); + if (mustFreeSpline ==1) { + psFree(mySpline); + } + // PS_VECTOR_PRINT_F32(newVec); + + } else { psError(PS_ERR_UNKNOWN, true, "unknown fit type. Returning NULL.\n"); @@ -413,882 +321,257 @@ } +#endif + +// Produce an overscan vector from an array of pixels +static psVector *overscanVector(pmOverscanOptions *overscanOpts, // Overscan options + const psArray *pixels, // Array of vectors containing the pixel values + psStats *myStats // Statistic to use in reducing the overscan + ) +{ + // 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_U8); // 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.U8[i] = 0; + ordinate->data.F32[i] = 2.0*(float)i/(float)pixels->n - 1.0; // Scale to [-1,1] + psVectorStats(myStats, values, NULL, NULL, 0); + double reducedVal = NAN; // Result of statistics + if (! p_psGetStatValue(myStats, &reducedVal)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result " + "of statistics on row %d.\n", i); + return NULL; + } + reduced->data.F32[i] = reducedVal; + } 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.U8[i] = 1; + } + } + + // Fit the overscan, if required + switch (overscanOpts->fitType) { + case PM_FIT_NONE: + // No fitting --- that's easy. + 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); + psFree(reduced); + reduced = 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); + psFree(reduced); + reduced = 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 requires an + // input spline + overscanOpts->spline = psVectorFitSpline1D(reduced, ordinate); + psFree(reduced); + reduced = psSpline1DEvalVector(overscanOpts->spline, ordinate); + break; + default: + psError(PS_ERR_UNKNOWN, true, "Unknown value for the fitting type: %d\n", overscanOpts->fitType); + return NULL; + break; + } + + psFree(ordinate); + psFree(mask); + + return reduced; +} + + + /****************************************************************************** +XXX: The SDRS does not specify type support. F32 is implemented here. *****************************************************************************/ -static psS32 GetOverscanSize( - psImage *inImg, - pmOverscanAxis overScanAxis) -{ - if (overScanAxis == PM_OVERSCAN_ROWS) { - return(inImg->numCols); - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - return(inImg->numRows); - } else if (overScanAxis == PM_OVERSCAN_ALL) { - return(1); - } - return(0); -} - -/****************************************************************************** -GetOverscanAxis(in) this private routine determines the appropiate overscan -axis from the parent cell metadata. - -XXX: Verify the READDIR corresponds with my overscan axis. - *****************************************************************************/ -static pmOverscanAxis GetOverscanAxis(pmReadout *in) -{ - psBool rc; - if ((in->parent != NULL) && (in->parent->concepts)) { - psS32 dir = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.READDIR"); - if (rc == true) { - if (dir == 1) { - return(PM_OVERSCAN_ROWS); - } else if (dir == 2) { - return(PM_OVERSCAN_COLUMNS); - } else if (dir == 3) { - return(PM_OVERSCAN_ALL); - } - } - } - - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): could not determine CELL.READDIR from in->parent metadata. Setting overscan axis to PM_OVERSCAN_NONE.\n"); - return(PM_OVERSCAN_NONE); -} - -/****************************************************************************** -my_psListLength(list): determine the length of a psList. - -XXX: Put this elsewhere. - -XXX: psList.h now has a version of this function. Use that instead. - *****************************************************************************/ - -static psS32 my_psListLength( - psList *list) -{ - psS32 length = 0; - psListElem *tmpElem = (psListElem *) list->head; - while (NULL != tmpElem) { - tmpElem = tmpElem->next; - length++; - } - return(length); -} - -/****************************************************************************** -Note: this isn't needed anymore as of psModule SDRS 12-09. - *****************************************************************************/ -static psBool OverscanReducePixel( - psImage *in, - psList *bias, - psStats *myStats) -{ - PS_ASSERT_PTR_NON_NULL(in, NULL); - PS_ASSERT_PTR_NON_NULL(bias, NULL); - PS_ASSERT_PTR_NON_NULL(bias->head, NULL); - PS_ASSERT_PTR_NON_NULL(myStats, NULL); - - // Allocate a psVector with one element per overscan image. - psS32 numOverscanImages = my_psListLength(bias); - psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32); - statsAll->n = statsAll->nalloc; - psListElem *tmpOverscan = (psListElem *) bias->head; - psS32 i = 0; - psF64 statValue; - // - // We loop through each overscan image, calculating the specified - // statistic on that image. - // - while (NULL != tmpOverscan) { - psImage *myOverscanImage = (psImage *) tmpOverscan->data; - - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, NULL); - myStats = psImageStats(myStats, myOverscanImage, NULL, (psMaskType)0xffffffff); - if (myStats == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - statsAll->data.F32[i] = statValue; - i++; - tmpOverscan = tmpOverscan->next; - } - - // - // We reduce the individual stats for each overscan image to - // a single psF32. - // - myStats = psVectorStats(myStats, statsAll, NULL, NULL, 0); - if (myStats == NULL) { - psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - psFree(statsAll); - return(false); - } - - // - // Subtract the result and return. - // - ImageSubtractScalar(in, statValue); - psFree(statsAll); - return(in); -} - -/****************************************************************************** -ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces -a single psImage to a column by combining all pixels from each row into a -single pixel via requested statistic in myStats. - *****************************************************************************/ -static psVector *ReduceOverscanImageToCol( - psImage *overscanImage, - psStats *myStats) -{ - psF64 statValue; - psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); - psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); - tmpRow->n = tmpRow->nalloc; - tmpCol->n = tmpCol->nalloc; - - // - // For each row, we store all pixels in that row into a temporary psVector, - // then we run psVectorStats() on that vector. - // - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - tmpRow->data.F32[j] = overscanImage->data.F32[i][j]; - } - - psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); - return(NULL); - } - - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); - return(NULL); - } - - tmpCol->data.F32[i] = (psF32) statValue; - } - psFree(tmpRow); - - return(tmpCol); -} - -/****************************************************************************** -ReduceOverscanImageToCol(overscanImage, myStats): This private routine reduces -a single psImage to a row by combining all pixels from each column into a -single pixel via requested statistic in myStats. - *****************************************************************************/ -static psVector *ReduceOverscanImageToRow( - psImage *overscanImage, - psStats *myStats) -{ - psF64 statValue; - psVector *tmpRow = psVectorAlloc(overscanImage->numCols, PS_TYPE_F32); - psVector *tmpCol = psVectorAlloc(overscanImage->numRows, PS_TYPE_F32); - tmpRow->n = tmpRow->nalloc; - tmpCol->n = tmpCol->nalloc; - - // - // For each column, we store all pixels in that column into a temporary psVector, - // then we run psVectorStats() on that vector. - // - for (psS32 i=0;inumCols;i++) { - for (psS32 j=0;jnumRows;j++) { - tmpCol->data.F32[j] = overscanImage->data.F32[j][i]; - } - - psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, true, "psVectorStats() could not perform requested statistical operation. Returning in image.\n"); - return(NULL); - } - - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, true, "p_psGetStatValue() could not determine result from requested statistical operation. Returning in image.\n"); - return(NULL); - } - - tmpRow->data.F32[i] = (psF32) statValue; - } - psFree(tmpCol); - - return(tmpRow); -} - -/****************************************************************************** -OverscanReduce(vecSize, bias, myStats): This private routine takes a psList of -overscan images (in bias) and reduces them to a single psVector via the -specified psStats struct. The vector is then scaled to the length or the -row/column in inImg. - *****************************************************************************/ -static psVector* OverscanReduce( - psImage *inImg, - pmOverscanAxis overScanAxis, - psList *bias, - void *fitSpec, - pmFit fit, - psStats *myStats) -{ - if ((overScanAxis != PM_OVERSCAN_ROWS) && (overScanAxis != PM_OVERSCAN_COLUMNS)) { - psError(PS_ERR_UNKNOWN, true, "overScanAxis must be PM_OVERSCAN_ROWS or PM_OVERSCAN_COLUMNS\n"); - return(NULL); - } - PS_ASSERT_PTR_NON_NULL(inImg, NULL); - PS_ASSERT_PTR_NON_NULL(bias, NULL); - PS_ASSERT_PTR_NON_NULL(bias->head, NULL); - PS_ASSERT_PTR_NON_NULL(myStats, NULL); - // - // Allocate a psVector for the output of this routine. - // - psS32 vecSize = GetOverscanSize(inImg, overScanAxis); - psVector *overscanVector = psVectorAlloc(vecSize, PS_TYPE_F32); - overscanVector->n = overscanVector->nalloc; - - // - // Allocate an array of psVectors with one psVector per element of the - // final oversan column vector. These psVectors will be used with - // psStats to reduce the multiple elements from each overscan column - // vector to a single final column vector. - // - psS32 numOverscanImages = my_psListLength(bias); - psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *)); - // (*overscanVectors)->n = (*overscanVectors)->nalloc; - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - overscanVectors[i] = NULL; - } - - // - // We iterate through the list of overscan images. For each image, - // we reduce it to a single column or row. Save the overscan vector - // in overscanVectors[]. - // - psListElem *tmpOverscan = (psListElem *) bias->head; - psS32 overscanID = 0; - while (tmpOverscan != NULL) { - psImage *tmpOverscanImage = (psImage *) tmpOverscan->data; - if (overScanAxis == PM_OVERSCAN_ROWS) { - overscanVectors[overscanID] = ReduceOverscanImageToRow(tmpOverscanImage, myStats); - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - overscanVectors[overscanID] = ReduceOverscanImageToCol(tmpOverscanImage, myStats); - } - - tmpOverscan = tmpOverscan->next; - overscanID++; - } - - // - // For each overscan vector, if necessary, we scale that column or - // row to vecSize. Note: we should have already ensured that the - // fit is poly or spline. - // - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psVector *tmpOverscanVector = overscanVectors[i]; - - if (tmpOverscanVector->n != vecSize) { - overscanVectors[i] = ScaleOverscanVector(tmpOverscanVector, vecSize, fitSpec, fit); - if (overscanVectors[i] == NULL) { - psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpOverscanVector); - return(NULL); - } - psFree(tmpOverscanVector); - } - } - - // - // We collect all elements in the overscan vectors for the various - // overscan images into a single psVector (tmpVec). Then we call - // psStats on that vector to determine the final values for the - // overscan vector. - // - psVector *tmpVec = psVectorAlloc(numOverscanImages, PS_TYPE_F32); - tmpVec->n = tmpVec->nalloc; - psF64 statValue; - for (psS32 i = 0 ; i < vecSize ; i++) { - // Collect the i-th elements from each overscan vector into a single vector. - for (psS32 j = 0 ; j < numOverscanImages ; j++) { - tmpVec->data.F32[j] = overscanVectors[j]->data.F32[i]; - } - - if (NULL == psVectorStats(myStats, tmpVec, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpVec); - return(NULL); - } - if (false == p_psGetStatValue(myStats, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - psFree(tmpVec); - return(NULL); - } - - overscanVector->data.F32[i] = (psF32) statValue; - } - - // - // We're done. Free the intermediate overscan vectors. - // - psFree(tmpVec); - for (psS32 i = 0 ; i < numOverscanImages ; i++) { - psFree(overscanVectors[i]); - } - psFree(overscanVectors); - - // - // Return the computed overscanVector - // - return(overscanVector); -} - -/****************************************************************************** -RebinOverscanVector(overscanVector, nBinOrig, myStats): this private routine -takes groups of nBinOrig elements in the input vector, combines them into a -single pixel via myStats and psVectorStats(), and then outputs a vector of -those pixels. - *****************************************************************************/ -static psS32 RebinOverscanVector( - psVector *overscanVector, - psS32 nBinOrig, - psStats *myStats) -{ - psF64 statValue; - psS32 nBin; - if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) { - psS32 numBins = 1+((overscanVector->n)/nBinOrig); - psVector *myBin = psVectorAlloc(numBins, PS_TYPE_F32); - psVector *binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32); - myBin->n = myBin->nalloc; - binVec->n = binVec->nalloc; - - for (psS32 i=0;in > ((i*nBinOrig)+j)) { - binVec->data.F32[j] = overscanVector->data.F32[(i*nBinOrig)+j]; - } else { - // XXX: we get here if nBinOrig does not evenly divide - // the overscanVector vector. This is the last bin. Should - // we change the binVec->n to acknowledge that? - binVec->n = j; - } - } - psStats *rc = psVectorStats(myStats, binVec, NULL, NULL, 0); - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation. Returning in image.\n"); - return(-1); - } - if (false == p_psGetStatValue(rc, &statValue)) { - psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning in image.\n"); - return(-1); - } - myBin->data.F32[i] = statValue; - } - - // Change the effective size of overscanVector. - overscanVector->n = numBins; - for (psS32 i=0;idata.F32[i] = myBin->data.F32[i]; - } - psFree(binVec); - psFree(myBin); - nBin = nBinOrig; - } else { - nBin = 1; - } - - return(nBin); -} - -/****************************************************************************** -FitOverscanVectorAndUnbin(inImg, overscanVector, overScanAxis, fitSpec, fit, -nBin): this private routine fits a psPolynomial or psSpline to the overscan -vector. It then creates a new vector, with a size determined by the input -image, evaluates the psPolynomial or psSpline at each element in that vector, -then returns that vector. - *****************************************************************************/ -static psVector *FitOverscanVectorAndUnbin( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis, - void *fitSpec, - pmFit fit, - psS32 nBin) -{ - psPolynomial1D* myPoly = NULL; - psSpline1D *mySpline = NULL; - // - // Fit a polynomial or spline to the overscan vector. - // - if (fit == PM_FIT_POLYNOMIAL) { - myPoly = (psPolynomial1D *) fitSpec; - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_POLY1D(myPoly, NULL); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, overscanVector, NULL, NULL); - if (myPoly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector. Returning NULL.\n"); - return(NULL); - } - } else if (fit == PM_FIT_SPLINE) { - mySpline = psVectorFitSpline1D(NULL, overscanVector); - if (mySpline == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector. Returning NULL.\n"); - return(NULL); - } - if (fitSpec != NULL) { - // Copy the resulting spline fit into fitSpec. - psSpline1D *ptrSpline = (psSpline1D *) fitSpec; - PS_ASSERT_SPLINE(ptrSpline, NULL); - SplineCopy(ptrSpline, mySpline); - } - } - - // - // Evaluate the poly/spline at each pixel in the overscan row/column. - // - psS32 vecSize = GetOverscanSize(inImg, overScanAxis); - psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); - newVec->n = newVec->nalloc; - if ((nBin > 1) && (nBin < overscanVector->n)) { - for (psS32 i = 0 ; i < vecSize ; i++) { - if (fit == PM_FIT_POLYNOMIAL) { - newVec->data.F32[i] = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin)); - } else if (fit == PM_FIT_SPLINE) { - newVec->data.F32[i] = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin)); - } - } - } else { - for (psS32 i = 0 ; i < vecSize ; i++) { - if (fit == PM_FIT_POLYNOMIAL) { - newVec->data.F32[i] = psPolynomial1DEval(myPoly, (psF32) i); - } else if (fit == PM_FIT_SPLINE) { - newVec->data.F32[i] = psSpline1DEval(mySpline, (psF32) i); - } - } - } - - psFree(mySpline); - psFree(overscanVector); - return(newVec); -} - - - -/****************************************************************************** -UnbinOverscanVector(inImg, overscanVector, overScanAxis, nBin): this private -routine takes a psVector overscanVector that was previously binned by a factor -of nBin, and then expands it to its original size, duplicated elements nBin -times for each element in the input vector overscanVector. - *****************************************************************************/ -static psVector *UnbinOverscanVector( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis, - psS32 nBin) -{ - psS32 vecSize = 0; - - if (overScanAxis == PM_OVERSCAN_ROWS) { - vecSize = inImg->numCols; - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - vecSize = inImg->numRows; - } - - psVector *newVec = psVectorAlloc(vecSize, PS_TYPE_F32); - newVec->n = newVec->nalloc; - for (psS32 i = 0 ; i < vecSize ; i++) { - newVec->data.F32[i] = overscanVector->data.F32[i/nBin]; - } - - psFree(overscanVector); - return(newVec); -} - - -/****************************************************************************** -SubtractVectorFromImage(inImg, overscanVector, overScanAxis): this private -routine subtracts the overscanVector column-wise or row-wise from inImg. - *****************************************************************************/ -static psImage *SubtractVectorFromImage( - psImage *inImg, - psVector *overscanVector, - pmOverscanAxis overScanAxis) -{ - // - // Subtract overscan vector row-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_ROWS) { - for (psS32 i=0;inumCols;i++) { - for (psS32 j=0;jnumRows;j++) { - inImg->data.F32[j][i]-= overscanVector->data.F32[i]; - } - } - } - - // - // Subtract overscan vector column-wise from the image. - // - if (overScanAxis == PM_OVERSCAN_COLUMNS) { - for (psS32 i=0;inumRows;i++) { - for (psS32 j=0;jnumCols;j++) { - inImg->data.F32[i][j]-= overscanVector->data.F32[i]; - } - } - } - - return(inImg); -} - - - -typedef enum { - PM_ERROR_NO_SUBTRACTION, - PM_WARNING_NO_SUBTRACTION, - PM_ERROR_NO_BIAS_SUBTRACT, - PM_WARNING_NO_BIAS_SUBTRACT, - PM_ERROR_NO_DARK_SUBTRACT, - PM_WARNING_NO_DARK_SUBTRACT, - PM_OKAY -} pmSubtractBiasAssertStatus; -/****************************************************************************** -AssertCodeOverscan(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for overscan subtraction. - *****************************************************************************/ -pmSubtractBiasAssertStatus AssertCodeOverscan( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - - PS_ASSERT_READOUT_NON_NULL(in, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_READOUT_NON_EMPTY(in, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - - if (overscan == true) { - pmOverscanAxis overScanAxis = GetOverscanAxis(in); - PS_ASSERT_PTR_NON_NULL(stat, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_PTR_NON_NULL(in->bias, PM_ERROR_NO_SUBTRACTION); - PS_ASSERT_PTR_NON_NULL(in->bias->head, PM_ERROR_NO_SUBTRACTION); - // - // Check the type, size of each bias image. - // - psListElem *tmpOverscan = (psListElem *) in->bias->head; - psS32 numOverscans = 0; - while (NULL != tmpOverscan) { - numOverscans++; - psImage *myOverscanImage = (psImage *) tmpOverscan->data; - PS_ASSERT_IMAGE_TYPE(myOverscanImage, PS_TYPE_F32, PM_ERROR_NO_SUBTRACTION); - // XXX: Get this right with the rows and columns. - if (overScanAxis == PM_OVERSCAN_ROWS) { - if (myOverscanImage->numRows != in->image->numRows) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image (# %d) has %d rows, input image has %d rows\n", - numOverscans, myOverscanImage->numCols, in->image->numRows); - if (fit == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - } else if (overScanAxis == PM_OVERSCAN_COLUMNS) { - if (myOverscanImage->numCols != in->image->numCols) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan image (# %d) has %d columns, input image has %d columns\n", - numOverscans, myOverscanImage->numCols, in->image->numCols); - if (fit == PM_FIT_NONE) { - psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vectors. Set fit to PM_FIT_POLYNOMIAL or PM_FIT_SPLINE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - } else if (overScanAxis != PM_OVERSCAN_ALL) { - psError(PS_ERR_UNKNOWN, true, "Must specify and overscan axis.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - tmpOverscan = tmpOverscan->next; - } - } else { - if (fit != PM_FIT_NONE) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): overscan is FALSE and fit is not PM_FIT_NONE.\n"); - return(PM_WARNING_NO_SUBTRACTION); - } - } - - // XXX: I do not like the following spec since it's useless to specify - // a psSpline as the fitSpec. - if (0) { - if ((fitSpec == NULL) && - ((fit != PM_FIT_NONE) || (overscan == true))) { - psError(PS_ERR_UNKNOWN, true, "fitSpec is NULL and fit is not PM_FIT_NONE or overscan is TRUE.\n"); - return(PM_ERROR_NO_SUBTRACTION); - } - } - - return(PM_OKAY); -} - -/****************************************************************************** -AssertCodeBias(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for bias subtraction. - *****************************************************************************/ -static pmSubtractBiasAssertStatus AssertCodeBias( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - if ((in->image->numRows + in->row0 - bias->row0) > bias->image->numRows) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough rows. Returning in image\n"); - return(PM_ERROR_NO_BIAS_SUBTRACT); - } - if ((in->image->numCols + in->col0 - bias->col0) > bias->image->numCols) { - psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns. Returning in image\n"); - return(PM_ERROR_NO_BIAS_SUBTRACT); - } - - if (bias != NULL) { - PS_ASSERT_READOUT_NON_EMPTY(bias, PM_ERROR_NO_BIAS_SUBTRACT); - PS_ASSERT_READOUT_TYPE(bias, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); - } - return(PM_OKAY); -} - -/****************************************************************************** -AssertCodeDark(....) this private routine verifies that the various input -parameters to pmSubtractBias() are correct for dark subtraction. - *****************************************************************************/ -pmSubtractBiasAssertStatus AssertCodeDark( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBinOrig, - const pmReadout *bias, - const pmReadout *dark) -{ - if ((in->image->numRows + in->row0 - dark->row0) > dark->image->numRows) { - psError(PS_ERR_UNKNOWN, true, "dark image does not have enough rows. Returning in image\n"); - return(PM_ERROR_NO_DARK_SUBTRACT); - } - if ((in->image->numCols + in->col0 - dark->col0) > dark->image->numCols) { - psError(PS_ERR_UNKNOWN, true, "dark image does not have enough columns. Returning in image\n"); - return(PM_ERROR_NO_DARK_SUBTRACT); - } - - if (dark != NULL) { - PS_ASSERT_READOUT_NON_EMPTY(dark, PM_ERROR_NO_DARK_SUBTRACT); - PS_ASSERT_READOUT_TYPE(dark, PS_TYPE_F32, PM_ERROR_NO_DARK_SUBTRACT); - } - return(PM_OKAY); -} - -/****************************************************************************** -p_psDetermineTrimmedImage(): global routine: determines the region of the -input pmReadout which will be operated on by the various detrend modules. It -does a metadata fetch on "CELL.TRIMSEC" for the parent cell of the pmReadout. - -Use it this way: - PS_WARN_PTR_NON_NULL(in->parent); - if (in->parent != NULL) { - PS_WARN_PTR_NON_NULL(in->parent->concepts); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - -XXX: Create a pmUtils.c file and put this routine there. - *****************************************************************************/ -psImage *p_psDetermineTrimmedImage(pmReadout *in) -{ - if ((in->parent == NULL) || (in->parent->concepts == NULL)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata (NULL).\n"); - return(in->image); - } - - psBool rc = false; - psImage *trimmedImg = NULL; - psRegion *trimRegion = psMetadataLookupPtr(&rc, in->parent->concepts, - "CELL.TRIMSEC"); - if (rc == false) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: could not determine CELL.TRIMSEC from parent cell Metadata.\n"); - trimmedImg = in->image; - } else { - trimmedImg = psImageSubset(in->image, *trimRegion); - } - - return(trimmedImg); -} - - -/****************************************************************************** -pmSubtractBias(....): see SDRS for complete specification. - -XXX: Code and assert type support: U16, S32, F32. -XXX: Add trace messages. - *****************************************************************************/ -pmReadout *pmSubtractBias( - pmReadout *in, - void *fitSpec, - pmFit fit, - bool overscan, - psStats *stat, - int nBin, - const pmReadout *bias, - const pmReadout *dark) +pmReadout *pmSubtractBias(pmReadout *in, pmOverscanOptions *overscanOpts, + const pmReadout *bias, const pmReadout *dark) { psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, "---- pmSubtractBias() begin ----\n"); - // - // Check input parameters, generate warnings and errors. - // - if (PM_OKAY != AssertCodeOverscan(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - return(in); - } - // - // Determine trimmed image from metadata. - // - psImage *trimmedImg = p_psDetermineTrimmedImage(in); - - // - // Subtract overscan frames if necessary. - // - if (overscan == true) { - pmOverscanAxis overScanAxis = GetOverscanAxis(in); - // - // Create a psStats data structure and determine the highest - // priority stats option. - // - psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); - if (stat != NULL) { - myStats->options = GenNewStatOptions(stat); - } - - // - // Reduce overscan images to a single pixel, then subtract. - // This code is no longer required as of SDRS 12-09. - // - if (overScanAxis == PM_OVERSCAN_ALL) { - if (false == OverscanReducePixel(trimmedImg, in->bias, myStats)) { + PS_ASSERT_READOUT_NON_NULL(in, NULL); + PS_ASSERT_READOUT_NON_EMPTY(in, NULL); + PS_ASSERT_READOUT_TYPE(in, PS_TYPE_F32, NULL); + + psImage *image = in->image; // The input image + + // Overscan processing + if (overscanOpts) { + // 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(in); + } + + psList *overscans = in->bias; // List of the overscan images + + psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN); // A new psStats, to avoid clobbering original + myStats->options = GenNewStatOptions(overscanOpts->stat); + + // Reduce all overscan pixels to a single value + if (overscanOpts->single) { + psVector *pixels = psVectorAlloc(0, PS_TYPE_F32); + pixels->n = 0; + 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); + // XXX Reimplement with memcpy + for (int i = 0; i < overscan->numRows; i++) { + for (int j = 0; j < overscan->numCols; j++) { + pixels->data.F32[index++] = overscan->data.F32[i][j]; + } + } + + } + psFree(iter); + + (void)psVectorStats(myStats, pixels, NULL, NULL, 0); + double reduced = NAN; // Result of statistics + if (! p_psGetStatValue(myStats, &reduced)) { + psError(PS_ERR_UNKNOWN, false, "p_psGetStatValue(): could not determine result from requested statistical operation. Returning input image.\n"); return(in); } - psFree(myStats); + (void)psBinaryOp(image, image, "-", psScalarAlloc((float)reduced, PS_TYPE_F32)); } else { - // - // Reduce the overscan images to a single overscan vector. - // - psVector *overscanVector = OverscanReduce(in->image, overScanAxis, - in->bias, fitSpec, - fit, myStats); - if (overscanVector == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not reduce overscan images to a single overscan vector. Returning in image\n"); - psFree(myStats); - return(in); + + // We do the regular overscan subtraction + + bool readRows = psMetadataLookupBool(NULL, in->parent->concepts, "CELL.READDIR");// Read direction + + if (readRows) { + // The read direction is rows + psArray *pixels = psArrayAlloc(image->numRows); // Array of vectors containing pixels + for (int i = 0; i < pixels->n; i++) { + psVector *values = psVectorAlloc(0, PS_TYPE_F32); + values->n = 0; + 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))) { + int diff = image->row0 - overscan->row0; // Offset between the two regions + for (int i = MAX(0,diff); i < MIN(image->numRows, overscan->numRows + diff); i++) { + // i is row on overscan + // XXX Reimplement with memcpy + psVector *values = pixels->data[i]; + int index = values->n; // Index in the vector + values = psVectorRealloc(values, values->n + overscan->numCols); + for (int j = 0; j < overscan->numCols; j++) { + values->data.F32[index++] = overscan->data.F32[i][j]; + } + values->n += overscan->numCols; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = overscanVector(overscanOpts, pixels, myStats); + psFree(pixels); + if (! reduced) { + return in; + } + + // 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]; + } + } + psFree(reduced); + + } else { + // 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); + values->n = 0; + 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))) { + int diff = image->col0 - overscan->col0; // Offset between the two regions + for (int i = MAX(0,diff); i < MIN(image->numCols, overscan->numCols + diff); i++) { + // i is column on overscan + // XXX Reimplement with memcpy + 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[i][j]; + } + values->n += overscan->numRows; + pixels->data[i] = values; // Update the pointer in case it's moved + } + } + psFree(iter); + + // Reduce the overscans + psVector *reduced = overscanVector(overscanOpts, pixels, myStats); + psFree(pixels); + if (! reduced) { + return in; + } + + // Subtract column by column + for (int i = 0; i < image->numCols; i++) { + for (int j = 0; j < image->numRows; j++) { + image->data.F32[j][i] -= reduced->data.F32[i]; + } + } + psFree(reduced); } - - // - // Rebin the overscan vector if necessary. - // - psS32 newBin = RebinOverscanVector(overscanVector, nBin, myStats); - if (newBin < 0) { - psError(PS_ERR_UNKNOWN, false, "Could rebin the overscan vector. Returning in image\n"); - psFree(myStats); - return(in); - } - - // - // If necessary, fit a psPolynomial or psSpline to the overscan vector. - // Then, unbin the overscan vector to appropriate length for the in image. - // - if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) { - overscanVector = FitOverscanVectorAndUnbin(trimmedImg, overscanVector, overScanAxis, fitSpec, fit, newBin); - if (overscanVector == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit the polynomial or spline to the overscan vector. Returning in image\n"); - psFree(myStats); - return(in); - } - } else { - overscanVector = UnbinOverscanVector(trimmedImg, overscanVector, overScanAxis, newBin); - } - - // - // Subtract the overscan vector from the input image. - // - SubtractVectorFromImage(trimmedImg, overscanVector, overScanAxis); - psFree(myStats); - psFree(overscanVector); - } - } - - // - // Perform bias subtraction if necessary. - // - if (bias != NULL) { - if (PM_OKAY == AssertCodeBias(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - SubtractFrame(in, bias); - } - } - - // - // Perform dark subtraction if necessary. - // - if (dark != NULL) { - if (PM_OKAY == AssertCodeDark(in, fitSpec, fit, overscan, stat, nBin, bias, dark)) { - psBool rc; - psF32 scale = 0.0; - if (in->parent != NULL) { - scale = psMetadataLookupS32(&rc, in->parent->concepts, "CELL.DARKTIME"); - if (rc == false) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: pmSubtractBias.(): could not determine CELL.FARKTIME from in->parent metadata.\n"); - } - } - SubtractDarkFrame(in, dark, scale); - } - } - - // - // All done. - // - psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4, - "---- pmSubtractBias() exit ----\n"); - return(in); -} - - + } + psFree(myStats); + } // End of overscan subtraction + + // Bias frame subtraction + if (bias) { + SubtractFrame(in, bias, 1.0); + } + + if (dark) { + // Get the scaling + float inTime = psMetadataLookupF32(NULL, in->parent->concepts, "CELL.DARKTIME"); + float darkTime = psMetadataLookupF32(NULL, dark->parent->concepts, "CELL.DARKTIME"); + SubtractFrame(in, dark, inTime/darkTime); + } + + return in; +} + +