/** @file  pmSubtractBias.c
 *
 *  This file will contain a module which will subtract the detector bias
 *  in place from an input image.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2006-02-06 21:03:25 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */
/*****************************************************************************/
/* INCLUDE FILES                                                             */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include <string.h>
#include "pslib.h"
#if HAVE_CONFIG_H
#include <config.h>
#endif
#include "pmSubtractBias.h"

/*****************************************************************************/
/* DEFINE STATEMENTS                                                         */
/*****************************************************************************/
// XXX: put these in psConstants.h
void PS_POLY1D_PRINT(
    psPolynomial1D *poly)
{
    printf("-------------- PS_POLY1D_PRINT() --------------\n");
    printf("poly->nX is %d\n", poly->nX);
    for (psS32 i = 0 ; i < (1 + poly->nX) ; i++) {
        printf("poly->coeff[%d] is %f\n", i, poly->coeff[i]);
    }
}

void PS_PRINT_SPLINE(psSpline1D *mySpline)
{
    printf("-------------- PS_PRINT_SPLINE() --------------\n");
    printf("mySpline->n is %d\n", mySpline->n);
    for (psS32 i = 0 ; i < mySpline->n ; i++) {
        PS_POLY1D_PRINT(mySpline->spline[i]);
    }
    PS_VECTOR_PRINT_F32(mySpline->knots);
}

#define PS_IMAGE_PRINT_F32_HIDEF(NAME) \
printf("======== printing %s ========\n", #NAME); \
for (int i = 0 ; i < (NAME)->numRows ; i++) { \
    for (int j = 0 ; j < (NAME)->numCols ; j++) { \
        printf("%.5f ", (NAME)->data.F32[i][j]); \
    } \
    printf("\n"); \
}\

/*****************************************************************************/
/* TYPE DEFINITIONS                                                          */
/*****************************************************************************/

/*****************************************************************************/
/* GLOBAL VARIABLES                                                          */
/*****************************************************************************/
psS32 currentId = 0;                // XXX: remove
psS32 memLeaks = 0;                 // XXX: remove
//PRINT_MEMLEAKS(8); XXX
/*****************************************************************************/
/* FILE STATIC VARIABLES                                                     */
/*****************************************************************************/

/*****************************************************************************/
/* FUNCTION IMPLEMENTATION - LOCAL                                           */
/*****************************************************************************/

/******************************************************************************
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.
*****************************************************************************/
static pmReadout *SubtractFrame(
    pmReadout *in,
    const pmReadout *bias)
{
    // XXX: When did the ->row0 and ->col0 offsets get coded?
    for (psS32 i=0;i<in->image->numRows;i++) {
        for (psS32 j=0;j<in->image->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];
            }
        }
    }

    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;i<in->image->numRows;i++) {
            for (psS32 j=0;j<in->image->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 (psS32 i=0;i<in->image->numRows;i++) {
            for (psS32 j=0;j<in->image->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);
}

/******************************************************************************
ImageSubtractScalar(): subtract a scalar from the input image.
 
XXX: Is there a psLib function for this?
 *****************************************************************************/
static psImage *ImageSubtractScalar(
    psImage *image,
    psF32 scalar)
{
    for (psS32 i=0;i<image->numRows;i++) {
        for (psS32 j=0;j<image->numCols;j++) {
            image->data.F32[i][j]-= scalar;
        }
    }
    return(image);
}

/******************************************************************************
GenNewStatOptions(): this routine will take as input the options member of the
stat data structure, determine if multiple options have been specified, issue
a warning message if so, and return the highest priority option (according to
the order of the if-statements in this code).  The higher priority options are
listed lower in the code.
 *****************************************************************************/
static psStatsOptions GenNewStatOptions(const psStats *stat)
{
    psS32 numOptions = 0;
    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) {
            opt = PS_STAT_ROBUST_MEDIAN;
        }
        numOptions++;
    }

    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) {
            opt = PS_STAT_CLIPPED_MEAN;
        }
        numOptions++;
    }

    if (stat->options & PS_STAT_SAMPLE_MEDIAN) {
        if (numOptions == 0) {
            opt = PS_STAT_SAMPLE_MEDIAN;
        }
        numOptions++;
    }

    if (stat->options & PS_STAT_SAMPLE_MEAN) {
        numOptions++;
        opt = PS_STAT_SAMPLE_MEAN;
    }


    if (numOptions == 0) {
        psError(PS_ERR_UNKNOWN,true, "No allowable statistics options have been specified.\n");
    }
    if (numOptions != 1) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmSubtractBias.c: GenNewStatOptions(): Too many statistics options have been specified\n");
    }
    return(opt);
}

/******************************************************************************
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);
    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);
}

/******************************************************************************
ScaleOverscanVector(): this routine takes as input an arbitrary vector,
creates a new vector of length n, and fills the new vector with the
interpolated values of the old vector.  The type of interpolation is:
    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.
 *****************************************************************************/
static psVector *ScaleOverscanVector(
    psVector *overscanVector,
    psS32 n,
    void *fitSpec,
    pmFit fit)
{
    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
            "---- ScaleOverscanVector() begin (%d -> %d) ----\n", overscanVector->n, n);

    if (NULL == overscanVector) {
        return(overscanVector);
    }

    // Allocate the new vector.
    psVector *newVec = psVectorAlloc(n, PS_TYPE_F32);

    //
    // 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));
    }
    psF32 x;

    if (fit == PM_FIT_POLYNOMIAL) {
        // Fit a polynomial to the old overscan vector.
        psPolynomial1D *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");
            return(NULL);
        }

        // For each element of the new vector, convert the x-ordinate to that
        // 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;i<n;i++) {
            x = ((psF32) i) * ((psF32) overscanVector->n) / ((psF32) n);
            newVec->data.F32[i] = psPolynomial1DEval(myPoly, x);
        }
    } else if (fit == PM_FIT_SPLINE) {
        //
        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
        // splines endpoints will be from 0.0 to overscanVector->n-1.  Must scale
        // properly when doing the spline eval.
        //
        psSpline1D *mySpline = psVectorFitSpline1D(NULL, overscanVector);
        if (mySpline == NULL) {
            psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to the psVector.\n");
            return(NULL);
        }

        // For each element of the new vector, convert the x-ordinate to that
        // of the old vector, use the fitted spline to determine the
        // interpolated value at that point, and set the new vector.
        for (psS32 i=0;i<n;i++) {
            // Scale to [0 : overscanVector->n - 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);
    } else {
        psError(PS_ERR_UNKNOWN, true, "unknown fit type.  Returning NULL.\n");
        psFree(newVec);
        return(NULL);
    }

    psTrace(".psModule.pmSubtracBias.ScaleOverscanVector", 4,
            "---- ScaleOverscanVector() exit ----\n");
    return(newVec);
}

/******************************************************************************
 *****************************************************************************/
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);
}

/******************************************************************************
psListLength(list): determine the length of a psList.
 
XXX: Put this elsewhere.
 *****************************************************************************/
static psS32 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 = psListLength(bias);
    psVector *statsAll = psVectorAlloc(numOverscanImages, PS_TYPE_F32);
    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);

    //
    // 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;i<overscanImage->numRows;i++) {
        for (psS32 j=0;j<overscanImage->numCols;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);

    //
    // 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;i<overscanImage->numCols;i++) {
        for (psS32 j=0;j<overscanImage->numRows;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);

    //
    // 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 = psListLength(bias);
    psVector **overscanVectors = (psVector **) psAlloc(numOverscanImages * sizeof(psVector *));
    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);
    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);

        for (psS32 i=0;i<numBins;i++) {
            for(psS32 j=0;j<nBinOrig;j++) {
                if (overscanVector->n > ((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;i<numBins;i++) {
            overscanVector->data.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);
    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);
    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;i<inImg->numCols;i++) {
            for (psS32 j=0;j<inImg->numRows;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;i<inImg->numRows;i++) {
            for (psS32 j=0;j<inImg->numCols;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)
{
    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)) {
                return(in);
            }
            psFree(myStats);
        } 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);
            }

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


