/** @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.16 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-15 02:20:39 $
 *
 *  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)
{
    int i;
    int j;

    if (bias == NULL) {
        printf("WARNING: bias frame is NULL\n");
        return(in);
    }

    if (((in->image)->numRows + in->row0 - bias->row0) > (bias->image)->numRows) {
        psError(__func__, "bias image does not have enough rows\n");
    }
    if (((in->image)->numCols + in->row0 - bias->row0) > (bias->image)->numCols) {
        psError(__func__, "bias image does not have enough columns\n");
    }

    // XXX: Verify that this is the correct way to use offsets.
    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.
 
// XXX: Is there already a function which does this?
 *****************************************************************************/
psImage *p_psImageSubtractScalar(psImage *image,
                                 float scalar)
{
    int i;
    int j;

    for (i=0;i<image->numRows;i++) {
        for (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).
 *****************************************************************************/
psStatsOptions GenNewStatOptions(const psStats *stat)
{
    int numOptions = 0;
    psStatsOptions opt = 0;

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

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

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

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

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

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

    if (numOptions == 0) {
        psError(__func__, "No statistics options have been specified\n");
    }
    if (numOptions != 1) {
        printf("WARNING: 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,
                              int n,
                              void *fitSpec,
                              pmFit fit)
{
    printf("ScaleOverscanVector(%d)\n", n);
    if (NULL == overscanVector) {
        return(overscanVector);
    }

    psVector *newVec = psVectorAlloc(n, PS_TYPE_F32);
    psPolynomial1D *myPoly;
    psSpline1D *mySpline;
    float x;
    int i;

    if (fit == PM_FIT_POLYNOMIAL) {
        // Fit a polynomial to the old overscan vector.
        myPoly = (psPolynomial1D *) fitSpec;
        myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, 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 = ((float) i) * ((float) overscanVector->n) / ((float) n);
            newVec->data.F32[i] = psPolynomial1DEval(x, myPoly);
        }
    } else if (fit == PM_FIT_SPLINE) {
        // Fit a spline to the old overscan vector.
        mySpline = (psSpline1D *) fitSpec;
        mySpline = psVectorFitSpline1D(mySpline, NULL, overscanVector, 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 = ((float) i) * ((float) overscanVector->n) / ((float) n);
            newVec->data.F32[i] = psSpline1DEval(mySpline, x);
        }
    } else {
        psError(__func__, "unknown fit type\n");
    }

    return(newVec);
}

psReadout *pmSubtractBias(psReadout *in,
                          void *fitSpec,
                          const psList *overscans,
                          pmOverscanAxis overScanAxis,
                          psStats *stat,
                          int nBinOrig,
                          pmFit fit,
                          const psReadout *bias)
{
    int i;
    int j;
    int 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;
    psSpline1D *mySpline;
    int nBin;

    PS_CHECK_NULL_PTR_RETURN_NULL(in);

    //
    //  Determine if multiple options are specified in stat.
    //
    if (stat != NULL) {
        stat->options = GenNewStatOptions(stat);
    }

    if ((overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)) {
        psError(__func__, "(overscans == NULL) && (overScanAxis != PM_OVERSCAN_NONE)\n");
        return(in);
    }

    if (((overScanAxis == PM_OVERSCAN_NONE) ||
            (overScanAxis == PM_OVERSCAN_ALL)) &&
            (fit == PM_FIT_NONE)) {
        printf("WARNING: overScanAxis equals NONE or ALL, and fit equals NONE\n");
        //XXX: Then what?
    }

    if (overScanAxis == PM_OVERSCAN_NONE) {
        if (overscans != NULL) {
            printf("WARNING: overScanAxis equals NONE and overscans does not equal NULL\n");
        }
        return(psSubtractFrame(in, bias));
    }

    //
    // We subtract each overscan region from the image data.
    //
    if (overScanAxis == PM_OVERSCAN_ALL) {
        tmpOverscan = (psListElem *) overscans->head;
        while (NULL != tmpOverscan) {
            myOverscanImage = (psImage *) tmpOverscan->data;
            psImageStats(stat, myOverscanImage, NULL, 0xffffffff);
            p_psGetStatValue(stat, &statValue);
            p_psImageSubtractScalar(in->image, statValue);

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


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

        // XXX: Is there a better way to extract a psVector from a psImage without
        // having to copy every element in that vector?
        // XXX: If we get here, do we know that overScanAxis == PM_OVERSCAN_ROWS
        // or PM_OVERSCAN_COL?

        if ((overScanAxis == PM_OVERSCAN_ROWS) || (overScanAxis == PM_OVERSCAN_COLUMNS)) {
            if (overScanAxis == PM_OVERSCAN_ROWS) {
                // We create a row vector and subtract this vector from image.
                overscanVector = psVectorAlloc(myOverscanImage->numCols, PS_TYPE_F32);
                for (i=0;i<overscanVector->n;i++) {
                    overscanVector->data.F32[i] = 0.0;
                }
                tmpRow = psVectorAlloc((in->image)->numRows, PS_TYPE_F32);

                for (i=0;i<myOverscanImage->numCols;i++) {
                    for (j=0;j<myOverscanImage->numRows;j++) {
                        tmpRow->data.F32[j] = myOverscanImage->data.F32[j][i];
                    }
                    stat = psVectorStats(stat, tmpRow, NULL, 0);
                    p_psGetStatValue(stat, &statValue);
                    overscanVector->data.F32[i] = statValue;
                }
                psFree(tmpRow);
                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
                    psVector *newVec = ScaleOverscanVector(overscanVector, (in->image)->numCols, fitSpec, fit);
                    psFree(overscanVector);
                    overscanVector = newVec;
                }
            }

            if (overScanAxis == PM_OVERSCAN_COLUMNS) {
                // 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((in->image)->numCols, PS_TYPE_F32);

                for (i=0;i<myOverscanImage->numRows;i++) {
                    for (j=0;j<myOverscanImage->numCols;j++) {
                        tmpCol->data.F32[j] = myOverscanImage->data.F32[i][j];
                    }
                    stat = psVectorStats(stat, tmpCol, NULL, 0);
                    p_psGetStatValue(stat, &statValue);
                    overscanVector->data.F32[i] = statValue;
                }

                psFree(tmpCol);
                if ((fit == PM_FIT_POLYNOMIAL) || (fit == PM_FIT_SPLINE)) {
                    psVector *newVec = ScaleOverscanVector(overscanVector, (in->image)->numRows, fitSpec, fit);
                    psFree(overscanVector);
                    overscanVector = newVec;
                }
            }

            //
            // Re-bin the overscan vector (change its length).
            //
            if ((nBinOrig > 0) && (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;
                        }
                    }
                    stat = psVectorStats(stat, binVec, NULL, 0);
                    p_psGetStatValue(stat, &statValue);
                    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;
            }

            //
            // Fit a polynomial to the overscan vector.
            //
            if (!((fitSpec == NULL) || (fit == PM_FIT_NONE))) {
                if (fit == PM_FIT_POLYNOMIAL) {
                    myPoly = (psPolynomial1D *) fitSpec;
                    myPoly = psVectorFitPolynomial1D(myPoly, NULL, overscanVector, NULL);
                    for (i=0;i<numBins;i++) {
                        overscanVector->data.F32[i] = psPolynomial1DEval((float) i, myPoly);
                    }

                } else if (fit == PM_FIT_SPLINE) {
                    mySpline = (psSpline1D *) fitSpec;
                    // XXX: The spline functions do not currently work if the number of
                    // data points and the number of splines are not equal.
                }
            }

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


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


