/** @file  pmSubtractBias.c
 *
 *  This file will contain a module which will subtract the detector bias
 *  in place from an input image.
 *
 *  @author George Gusciora, MHPCC
 *
 *  @version $Revision: 1.30 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-01-17 19:58:18 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */
#if HAVE_CONFIG_H
#include <config.h>
#endif

#include "pmSubtractBias.h"

#define PM_SUBTRACT_BIAS_POLYNOMIAL_ORDER 2
#define PM_SUBTRACT_BIAS_SPLINE_ORDER 3

/******************************************************************************
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.
 *****************************************************************************/
psReadout *psSubtractFrame(psReadout *in,
                           const psReadout *bias)
{
    psS32 i;
    psS32 j;

    if (bias == NULL) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmSubtractBias.c: psSubtractFrame(): bias frame is NULL.  Returning original image.\n");
        return(in);
    }

    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(in);
    }
    if (((in->image)->numCols + in->row0 - bias->row0) > (bias->image)->numCols) {
        psError(PS_ERR_UNKNOWN,true, "bias image does not have enough columns.  Returning in image\n");
        return(in);
    }

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

/******************************************************************************
p_psImageSubtractScalar(): subtract a scalar from the input image.
 *****************************************************************************/
psImage *p_psImageSubtractScalar(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.
 *****************************************************************************/
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_ROBUST_MEAN) {
        if (numOptions == 0) {
            opt = PS_STAT_ROBUST_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 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);
}

/******************************************************************************
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.
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?
 *****************************************************************************/
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) {
        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.
        myPoly = (psPolynomial1D *) fitSpec;
        PS_POLY_CHECK_NULL(myPoly, NULL);
        myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL);
        if (myPoly == NULL) {
            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): 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 (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) {
        psS32 mustFreeSpline = 0;
        // Fit a spline to the old overscan vector.
        mySpline = (psSpline1D *) fitSpec;
        if (mySpline == NULL) {
            mustFreeSpline = 1;
        }

        //
        // NOTE: Since the X arg in the psVectorFitSpline1D() function is NULL,
        // splines enpoints will be from 0.0 to overscanVector->n-1.  Must scale
        // properly when doing the spline eval.
        //
        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
        if (mySpline == NULL) {
            psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): 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 polynomial to determine the
        // interpolated value at that point, and set the new vector.
        for (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);
        }
        if (mustFreeSpline ==1) {
            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);
}

/******************************************************************************
XXX: The SDRS does not specify type support.  F32 is implemented here.
 *****************************************************************************/
psReadout *pmSubtractBias(psReadout *in,
                          void *fitSpec,
                          const psList *overscans,
                          pmOverscanAxis overScanAxis,
                          psStats *stat,
                          psS32 nBinOrig,
                          pmFit fit,
                          const psReadout *bias)
{
    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
            "---- pmSubtractBias() begin ----\n");
    PS_READOUT_CHECK_NULL(in, NULL);
    PS_READOUT_CHECK_EMPTY(in, NULL);
    PS_READOUT_CHECK_TYPE(in, PS_TYPE_F32, NULL);

    //
    // If the overscans != NULL, then check the type of each image.
    //
    if (overscans != NULL) {
        psListElem *tmpOverscan = (psListElem *) overscans->head;
        while (NULL != tmpOverscan) {
            psImage *myOverscanImage = (psImage *) tmpOverscan->data;
            PS_IMAGE_CHECK_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
            tmpOverscan = tmpOverscan->next;
        }
    }

    if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
        psError(PS_ERR_UNKNOWN,true, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE).  Returning in image\n");
        return(in);
    }

    // Check for an unallowable pmFit.
    if ((fit != PM_OVERSCAN_NONE) &&
            (fit != PM_OVERSCAN_ROWS) &&
            (fit != PM_OVERSCAN_COLUMNS) &&
            (fit != PM_OVERSCAN_ALL)) {
        psError(PS_ERR_UNKNOWN, true, "fit is unallowable (%d).  Returning in image.\n", fit);
        return(in);
    }
    // Check for an unallowable pmOverscanAxis.
    if ((overScanAxis != PM_OVERSCAN_NONE) &&
            (overScanAxis != PM_OVERSCAN_ROWS) &&
            (overScanAxis != PM_OVERSCAN_COLUMNS) &&
            (overScanAxis != PM_OVERSCAN_ALL)) {
        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).  Returning in image.\n", overScanAxis);
        return(in);
    }
    psS32 i;
    psS32 j;
    psS32 numBins = 0;
    static psVector *overscanVector = NULL;
    psVector *tmpRow = NULL;
    psVector *tmpCol = NULL;
    psVector *myBin = NULL;
    psVector *binVec = NULL;
    psListElem *tmpOverscan = NULL;
    double statValue;
    psImage *myOverscanImage = NULL;
    psPolynomial1D *myPoly = NULL;
    psSpline1D *mySpline = NULL;
    psS32 nBin;

    //
    //  Create a static stats data structure and determine the highest
    //  priority stats option.
    //
    static psStats *myStats = NULL;
    if (myStats == NULL) {
        myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
        p_psMemSetPersistent(myStats, true);
    }
    if (stat != NULL) {
        myStats->options = GenNewStatOptions(stat);
    }




    if (overScanAxis == PM_OVERSCAN_NONE) {
        if (fit != PM_FIT_NONE) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE, and fit does not equal NONE.  Proceeding to full fram subtraction.\n");
        }

        if (overscans != NULL) {
            psLogMsg(__func__, PS_LOG_WARN,
                     "WARNING: pmSubtractBias.(): overScanAxis equals NONE and overscans does not equal NULL.  Proceeding to full fram subtraction.\n");
        }
        return(psSubtractFrame(in, bias));
    }

    if ((overScanAxis == PM_OVERSCAN_ALL) && (fit != PM_FIT_NONE)) {
        psLogMsg(__func__, PS_LOG_WARN,
                 "WARNING: pmSubtractBias.(): overScanAxis equals ALL, and fit does not equal NONE.  Proceeding with the rest of the module.\n");
    }


    //
    // We subtract each overscan region from the image data.
    // If we get here we know that overscans != NULL.
    //

    if (overScanAxis == PM_OVERSCAN_ALL) {
        tmpOverscan = (psListElem *) overscans->head;
        while (NULL != tmpOverscan) {
            myOverscanImage = (psImage *) tmpOverscan->data;

            PS_IMAGE_CHECK_TYPE(myOverscanImage, PS_TYPE_F32, NULL);
            psStats *rc = psImageStats(myStats, myOverscanImage, NULL, 0xffffffff);
            if (rc == NULL) {
                psError(PS_ERR_UNKNOWN, false, "psImageStats(): could not perform requested statistical operation.  Returning in image.\n");
                return(in);
            }
            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");
                return(in);
            }
            p_psImageSubtractScalar(in->image, statValue);

            tmpOverscan = tmpOverscan->next;
        }
        return(in);
    }

    // This check is redundant with above code.
    if (!((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS))) {
        psError(PS_ERR_UNKNOWN, true, "overScanAxis is unallowable (%d).\nReturning in image.\n", overScanAxis);
        return(in);
    }

    tmpOverscan = (psListElem *) overscans->head;
    while (NULL != tmpOverscan) {
        myOverscanImage = (psImage *) tmpOverscan->data;

        if (overScanAxis == PM_OVERSCAN_ROWS) {
            if (myOverscanImage->numCols != (in->image)->numCols) {
                psLogMsg(__func__, PS_LOG_WARN,
                         "WARNING: pmSubtractBias.(): overscan image has %d columns, input image has %d columns\n",
                         myOverscanImage->numCols, in->image->numCols);
            }

            // We create a row vector and subtract this vector from image.
            // XXX: Is there a better way to extract a psVector from a psImage without
            // having to copy every element in that vector?
            overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
            for (i=0;i<overscanVector->n;i++) {
                overscanVector->data.F32[i] = 0.0;
            }
            tmpRow = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);

            // For each column of the input image, loop through every row,
            // collect the pixel in that row, then performed the specified
            // statistical op on those pixels.  Store this in overscanVector.
            for (i=0;i<myOverscanImage->numCols;i++) {
                for (j=0;j<myOverscanImage->numRows;j++) {
                    tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
                }
                psStats *rc = psVectorStats(myStats, tmpRow, NULL, NULL, 0);
                if (rc == NULL) {
                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
                    return(in);
                }
                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(in);
                }
                overscanVector->data.F32[i] = statValue;
            }
            psFree(tmpRow);

            // Scale the overscan vector to the size of the input image.
            if (overscanVector->n != in->image->numCols) {
                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
                    psVector *newVec = ScaleOverscanVector(overscanVector,
                                                           in->image->numCols,
                                                           fitSpec, fit);
                    if (newVec == NULL) {
                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
                        return(in);
                    }
                    psFree(overscanVector);
                    overscanVector = newVec;
                } else {
                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
                    psFree(overscanVector);
                    return(in);
                }
            }
        }

        if (overScanAxis == PM_OVERSCAN_COLUMNS) {
            if (myOverscanImage->numRows != (in->image)->numRows) {
                psLogMsg(__func__, PS_LOG_WARN,
                         "WARNING: pmSubtractBias.(): overscan image has %d rows, input image has %d rows\n",
                         myOverscanImage->numRows, in->image->numRows);
            }

            // We create a column vector and subtract this vector from image.
            overscanVector = psVectorAlloc(myOverscanImage->numRows, PS_TYPE_F32);
            for (i=0;i<overscanVector->n;i++) {
                overscanVector->data.F32[i] = 0.0;
            }
            tmpCol = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);

            // For each row of the input image, loop through every column,
            // collect the pixel in that row, then performed the specified
            // statistical op on those pixels.  Store this in overscanVector.
            for (i=0;i<myOverscanImage->numRows;i++) {
                for (j=0;j<myOverscanImage->numCols;j++) {
                    tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
                }
                psStats *rc = psVectorStats(myStats, tmpCol, NULL, NULL, 0);
                if (rc == NULL) {
                    psError(PS_ERR_UNKNOWN, false, "psVectorStats(): could not perform requested statistical operation.  Returning in image.\n");
                    return(in);
                }
                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(in);
                }
                overscanVector->data.F32[i] = statValue;
            }
            psFree(tmpCol);

            // Scale the overscan vector to the size of the input image.
            if (overscanVector->n != in->image->numRows) {
                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
                    psVector *newVec = ScaleOverscanVector(overscanVector,
                                                           in->image->numRows,
                                                           fitSpec, fit);
                    if (newVec == NULL) {
                        psError(PS_ERR_UNKNOWN, false, "ScaleOverscanVector(): could not scale the overscan vector.  Returning in image.\n");
                        return(in);
                    }
                    psFree(overscanVector);
                    overscanVector = newVec;
                } else {
                    psError(PS_ERR_UNKNOWN, true, "Don't know how to scale the overscan vector.  Set fit to PM_FIT_SPLINE or PM_FIT_POLYNOMIAL.  Returning in image.\n");
                    psFree(overscanVector);
                    return(in);
                }
            }
        }

        //
        // Re-bin the overscan vector (change its length).
        //
        // Only if nBinOrig > 1.
        if ((nBinOrig > 1) && (nBinOrig < overscanVector->n)) {
            numBins = 1+((overscanVector->n)/nBinOrig);
            myBin = psVectorAlloc(numBins, PS_TYPE_F32);
            binVec = psVectorAlloc(nBinOrig, PS_TYPE_F32);

            for (i=0;i<numBins;i++) {
                for(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(in);
                }
                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(in);
                }
                myBin->data.F32[i] = statValue;
            }

            // Change the effective size of overscanVector.
            overscanVector->n = numBins;
            for (i=0;i<numBins;i++) {
                overscanVector->data.F32[i] = myBin->data.F32[i];
            }
            psFree(binVec);
            psFree(myBin);
            nBin = nBinOrig;
        } else {
            nBin = 1;
        }

        // At this point the number of data points in overscanVector should be
        // equal to the number of rows/columns (whatever is appropriate) in the
        // image divided by numBins.
        //


        //
        // This doesn't seem right.  The only way to do a spline fit is if,
        // by SDRS requirements, fitSpec is not-NULL>  But in order for it
        // to be non-NULL, someone must have called psSpline1DAlloc() with
        // the min, max, and number of splines.
        //
        if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
            //
            // Fit a polynomial or spline to the overscan vector.
            //
            if (fit == PM_FIT_POLYNOMIAL) {
                myPoly = (psPolynomial1D *) fitSpec;
                myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL);
                if (myPoly == NULL) {
                    psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to overscan vector.  Returning in image.\n");
                    psFree(overscanVector);
                    return(in);
                }
            } else if (fit == PM_FIT_SPLINE) {
                mySpline = (psSpline1D *) fitSpec;
                mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, NULL);
                if (mySpline == NULL) {
                    psError(PS_ERR_UNKNOWN, false, "Could not fit a spline to overscan vector.  Returning in image.\n");
                    psFree(overscanVector);
                    return(in);
                }
            }

            //
            // Subtract fitted overscan vector row-wise from the image.
            //
            if (overScanAxis == PM_OVERSCAN_ROWS) {
                for (i=0;i<(in->image)->numCols;i++) {
                    psF32 tmpF32 = 0.0;
                    if (fit == PM_FIT_POLYNOMIAL) {
                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
                    } else if (fit == PM_FIT_SPLINE) {
                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
                    }
                    for (j=0;j<(in->image)->numRows;j++) {
                        (in->image)->data.F32[j][i]-= tmpF32;
                    }
                }
            }

            //
            // Subtract fitted overscan vector column-wise from the image.
            //
            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
                for (i=0;i<(in->image)->numRows;i++) {
                    psF32 tmpF32 = 0.0;
                    if (fit == PM_FIT_POLYNOMIAL) {
                        tmpF32 = psPolynomial1DEval(myPoly, ((psF32) i) / ((psF32) nBin));
                    } else if (fit == PM_FIT_SPLINE) {
                        tmpF32 = psSpline1DEval(mySpline, ((psF32) i) / ((psF32) nBin));
                    }

                    for (j=0;j<(in->image)->numCols;j++) {
                        (in->image)->data.F32[i][j]-= tmpF32;
                    }
                }
            }
        } else {
            //
            // If we get here, then no polynomials were fit to the overscan
            // vector.  We simply subtract it, taking into account binning,
            // from the image.
            //

            //
            // Subtract overscan vector row-wise from the image.
            //
            if (overScanAxis == PM_OVERSCAN_ROWS) {
                for (i=0;i<(in->image)->numCols;i++) {
                    for (j=0;j<(in->image)->numRows;j++) {
                        (in->image)->data.F32[j][i]-= overscanVector->data.F32[i/nBin];
                    }
                }
            }

            //
            // Subtract overscan vector column-wise from the image.
            //
            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
                for (i=0;i<(in->image)->numRows;i++) {
                    for (j=0;j<(in->image)->numCols;j++) {
                        (in->image)->data.F32[i][j]-= overscanVector->data.F32[i/nBin];
                    }
                }
            }
        }

        psFree(overscanVector);

        tmpOverscan = tmpOverscan->next;
    }

    psTrace(".psModule.pmSubtracBias.pmSubtractBias", 4,
            "---- pmSubtractBias() exit ----\n");

    if (bias != NULL) {
        return(psSubtractFrame(in, bias));
    }
    return(in);
}


