/** @file tst_pmSubtractBias.c
 *
 *  @brief Contains the tests for pmSubtractBias.c:
 *
 * test00: This code will subtract full bias frames from the input image.
 * test01: Multiple overscan regions, calculate a scalar statistic for
 *  each, then subtract from the input image.
 * test02: Calculate a column overscan vector and subtract it from each
 *  column in the input image.
 * test03: Calculate a row overscan vector and subtract it from each
 *  row in the input image.
 *
 *  @author GLG, MHPCC
 *
 *  @version $Revision: 1.15 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-01-10 22:13:39 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

#include "psTest.h"
#include "pslib.h"
#include "pmSubtractBias.h"
static int test00(void);
static int test01(void);
static int test02(void);
static int test03(void);
static int test04(void);
static int testX(void);
testDescription tests[] = {
                              {test00, 000, "pmSubtractBias", 0, false},
                              {test01, 000, "pmSubtractBias", 0, false},
                              {test02, 000, "pmSubtractBias", 0, false},
                              {test03, 000, "pmSubtractBias", 0, false},
                              {test04, 000, "pmSubtractBias", 0, false},
                              {testX,  000, "pmSubtractBias", 0, false},
                              {NULL}
                          };


int main(int argc, char* argv[])
{
    return !runTestSuite(stderr, "Test Point Driver", tests, argc, argv);
}

#define NUM_ROWS 8
#define NUM_COLS 8
/******************************************************************************
doSubtractBiasFullFrame(): a sample psReadout as well as a bias image are
created and the bias image is subtracted from the psReadout.
 *****************************************************************************/
int doSubtractBiasFullFrame(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    //    psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImage1);
    //    psReadout *myBias = psReadoutAlloc(numCols, numRows, tmpImage2);
    psReadout *myReadout = psReadoutAlloc();
    psReadout *myBias = psReadoutAlloc();
    myReadout->image = tmpImage1;
    myBias->image = tmpImage2;


    printPositiveTestHeader(stdout, "pmSubtractBias", "doSubtractBiasFullFrame");
    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            myBias->image->data.F32[i][j] = 1.0;
        }
    }

    myReadout = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, NULL,
                               0, PM_FIT_NONE, myBias);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 1.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(myBias);
    printFooter(stdout, "pmSubtractBias", "doSubtractBiasFullFrame", true);
    return(testStatus);
}


int test00( void )
{
    int testStatus = 0;

    testStatus |= doSubtractBiasFullFrame(1, 1);
    testStatus |= doSubtractBiasFullFrame(NUM_COLS, 1);
    testStatus |= doSubtractBiasFullFrame(1, NUM_ROWS);
    testStatus |= doSubtractBiasFullFrame(NUM_COLS, NUM_ROWS);
    return(testStatus);
}

/******************************************************************************
doSubtractFullOverscans(): a sample psReadout as well as several overscan
images of the same size are created.  The overscan images are then subtracted
from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscans(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    //    psReadout *myReadout = psReadoutAlloc(numCols, numRows, tmpImage1);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;


    printPositiveTestHeader(stdout, "pmSubtractBias", "doSubtractFullOverscans");
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ALL, stat,
                               0, PM_FIT_NONE, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);

    printFooter(stdout, "pmSubtractBias", "doSubtractFullOverscans", true);
    return(testStatus);
}

int test01( void )
{
    int testStatus = 0;

    testStatus |= doSubtractFullOverscans(1, 1);
    testStatus |= doSubtractFullOverscans(1, NUM_ROWS);
    testStatus |= doSubtractFullOverscans(NUM_COLS, 1);
    testStatus |= doSubtractFullOverscans(NUM_COLS, NUM_ROWS);

    return(testStatus);
}

/******************************************************************************
doSubtractFullOverscans(): a sample psReadout as well as several overscan
images of the same size are created.  The overscan images are collected
pixel-by-pixel then subtracted column-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanColumns(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");

    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_COLUMNS, stat,
                               0, PM_FIT_NONE, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);

    printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    return(testStatus);
}
/******************************************************************************
doSubtractFullOverscans(): a sample psReadout as well as several overscan
images of the same size are created.  The overscan images are collected
pixel-by-pixel then subtracted column-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanColumnsPoly(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");

    myReadout = pmSubtractBias(myReadout, myPoly, list, PM_OVERSCAN_COLUMNS, stat,
                               0, PM_FIT_POLYNOMIAL, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);
    psFree(myPoly);

    printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    return(testStatus);
}
/******************************************************************************
doSubtractFullOverscansSmall(): a sample psReadout as well as several overscan
images of smaller size are created.  The overscan images are collected
pixel-by-pixel then subtracted column-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanColumnsSmall(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psImage *tmpImage2 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols/2, numRows/2, PS_TYPE_F32);
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }
    for (i=0;i<numRows/2;i++) {
        for (j=0;j<numCols/2;j++) {
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");

    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_COLUMNS, stat,
                               0, PM_FIT_POLYNOMIAL, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);

    printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    return(testStatus);
}

int test02( void )
{
    int testStatus = 0;

    testStatus |= doSubtractFullOverscanColumns(1, 1);
    testStatus |= doSubtractFullOverscanColumns(1, NUM_ROWS);
    testStatus |= doSubtractFullOverscanColumns(NUM_COLS, 1);
    testStatus |= doSubtractFullOverscanColumns(NUM_COLS, NUM_ROWS);
    /* These tests do not make sense until the SDRS is clarified.
        testStatus |= doSubtractFullOverscanColumnsSmall(1, 1);
        testStatus |= doSubtractFullOverscanColumnsSmall(1, NUM_ROWS);
        testStatus |= doSubtractFullOverscanColumnsSmall(NUM_COLS, 1);
        testStatus |= doSubtractFullOverscanColumnsSmall(NUM_COLS, NUM_ROWS);
    */

    return(testStatus);
}

/******************************************************************************
doSubtractFullOverscans(): a sample psReadout as well as several overscan
images of the same size are created.  The overscan images are collected
pixel-by-pixel then subtracted row-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanRows(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Row Overscans");

    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ROWS, stat,
                               0, PM_FIT_NONE, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);

    printFooter(stdout, "pmSubtractBias", "Row Overscans", true);
    return(testStatus);
}

/******************************************************************************
doSubtractFullOverscansSmall(): a sample psReadout as well as several overscan
images of smaller size are created.  The overscan images are collected
pixel-by-pixel then subtracted row-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanRowsSmall(int numCols, int numRows)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;
    psS32 OSnumRows = numRows-1;
    if (OSnumRows == 0) {
        OSnumRows = 1;
    }
    psS32 OSnumCols = numCols-1;
    if (OSnumCols == 0) {
        OSnumCols = 1;
    }
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(OSnumCols, OSnumRows, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psList *list;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }
    for (i=0;i<OSnumRows;i++) {
        for (j=0;j<OSnumCols;j++) {
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Row Overscans");

    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ROWS, stat,
                               0, PM_FIT_SPLINE, NULL);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);

    printFooter(stdout, "pmSubtractBias", "Row Overscans", true);
    return(testStatus);
}


int test03( void )
{
    int testStatus = 0;

    testStatus |= doSubtractFullOverscanRows(1, 1);
    testStatus |= doSubtractFullOverscanRows(1, NUM_ROWS);
    testStatus |= doSubtractFullOverscanRows(NUM_COLS, 1);
    testStatus |= doSubtractFullOverscanRows(NUM_COLS, NUM_ROWS);
    testStatus |= doSubtractFullOverscanRowsSmall(1, 1);
    testStatus |= doSubtractFullOverscanRowsSmall(1, NUM_ROWS);
    testStatus |= doSubtractFullOverscanRowsSmall(NUM_COLS, 1);
    testStatus |= doSubtractFullOverscanRowsSmall(NUM_COLS, NUM_ROWS);

    return(testStatus);
}



int doSubtractOverscansTestInputCases(int numCols, int numRows)
{
    int i;
    int j;
    int testStatus = 0;
    psImage *tmpImage1 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psImage *tmpImage2Short = psImageAlloc(numCols-1, numRows-1, PS_TYPE_F32);
    psImage *tmpImage3Short = psImageAlloc(numCols-1, numRows-1, PS_TYPE_F32);
    psImage *tmpImage4Short = psImageAlloc(numCols-1, numRows-1, PS_TYPE_F32);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    psReadout *rc = NULL;
    psList *list;
    psList *listShort;
    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    psImage *tmpImage5 = psImageAlloc(numCols, numRows, PS_TYPE_F32);
    psReadout *myBias = psReadoutAlloc();
    myBias->image = tmpImage5;
    printPositiveTestHeader(stdout, "pmSubtractBias", "Testing input parameter error conditions");

    psImage *tmpImage5ShortRows = psImageAlloc(numCols, numRows-1, PS_TYPE_F32);
    psReadout *myBiasShortRows = psReadoutAlloc();
    myBiasShortRows->image = tmpImage5ShortRows;
    psImage *tmpImage5ShortCols = psImageAlloc(numCols-1, numRows, PS_TYPE_F32);
    psReadout *myBiasShortCols = psReadoutAlloc();
    myBiasShortCols->image = tmpImage5ShortCols;

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
            myBias->image->data.F32[i][j] = 1.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    for (i=0;i<numRows-1;i++) {
        for (j=0;j<numCols-1;j++) {
            tmpImage2Short->data.F32[i][j] = 3.0;
            tmpImage3Short->data.F32[i][j] = 4.0;
            tmpImage4Short->data.F32[i][j] = 5.0;
        }
    }
    listShort = psListAlloc(tmpImage2Short);
    psListAdd(listShort, PS_LIST_HEAD, tmpImage3Short);
    psListAdd(listShort, PS_LIST_HEAD, tmpImage4Short);
    for (i=0;i<numRows-1;i++) {
        for (j=0;j<numCols;j++) {
            myBiasShortRows->image->data.F32[i][j] = 1.0;
        }
    }
    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols-1;j++) {
            myBiasShortCols->image->data.F32[i][j] = 1.0;
        }
    }


    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_ALL.  Should generate error.\n");
    rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_ALL, stat, 0, PM_FIT_NONE, NULL);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
    }

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_ROWS.  Should generate error.\n");
    rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_ROWS, stat, 0, PM_FIT_NONE, NULL);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_COLUMNS.  Should generate error.\n");
    rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_COLUMNS, stat, 0, PM_FIT_NONE, NULL);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with non-NULL overscan list and PM_OVERSCAN_NONE.  Should generate warning.\n");
    rc = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_NONE, stat, 0, PM_FIT_NONE, myBias);

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            psF32 expect = ((float) (i + j)) - 1.0;
            psF32 actual = rc->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }

            // Restore myReadout for next test.
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }

    /* XXX: This does not seem to be a requirement.
        printf("------------------------------------------------------------------\n");
        printf("Calling pmSubtractBias() with NULL overscan list and PM_OVERSCAN_NONE.  Should generate warning.\n");
        rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, stat,
                            0, PM_FIT_NONE, myBias);
     
        for (i=0;i<numRows;i++) {
            for (j=0;j<numCols;j++) {
                psF32 expect = ((float) (i + j)) - 1.0;
                psF32 actual = rc->image->data.F32[i][j];
                if (FLT_EPSILON < fabs(expect - actual)) {
                    printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                    testStatus = 1;
                }
     
                // Restore myReadout for next test.
                myReadout->image->data.F32[i][j] = (float) (i + j);
            }
        }
    */

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with PM_OVERSCAN_NONE and PM_FIT_POLYNOMIAL.  Should generate Warning.\n");
    rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, stat, 0, PM_FIT_POLYNOMIAL, myBias);
    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            psF32 expect = ((float) (i + j)) - 1.0;
            psF32 actual = rc->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }

            // Restore myReadout for next test.
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }


    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with PM_OVERSCAN_ALL and PM_FIT_SPLINE.  Should generate Warning.\n");
    rc = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, stat, 0, PM_FIT_SPLINE, myBias);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            psF32 expect = ((float) (i + j)) - 1.0;
            psF32 actual = rc->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }

            // Restore myReadout for next test.
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }


    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with multiple stats->options.  Should generate Warning.\n");
    stat->options|= PS_STAT_SAMPLE_MEDIAN;
    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ALL, stat,
                               0, PM_FIT_NONE, NULL);
    for (i=0;i<numRows;i++) {
        for (j=0;j<numCols;j++) {
            psF32 expect = ((float) (i + j)) - 12.0;
            psF32 actual = rc->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }
    stat->options = PS_STAT_SAMPLE_MEAN;

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() undersize overscans (PM_OVERSCAN_ROWS).  Should generate Warning.\n");
    rc = pmSubtractBias(myReadout, NULL, listShort, PM_OVERSCAN_ROWS, stat,
                        0, PM_FIT_NONE, NULL);
    /*
        for (i=0;i<numRows;i++) {
            for (j=0;j<numCols;j++) {
                psF32 expect = ((float) (i + j)) - 12.0;
                psF32 actual = rc->image->data.F32[i][j];
                if (FLT_EPSILON < fabs(expect - actual)) {
                    printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                    testStatus = 1;
                }
            }
        }
    */

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() undersize overscans (PM_OVERSCAN_COLUMNS).  Should generate Warning.\n");
    rc = pmSubtractBias(myReadout, NULL, listShort, PM_OVERSCAN_COLUMNS, stat,
                        0, PM_FIT_NONE, NULL);
    /*
        for (i=0;i<numRows;i++) {
            for (j=0;j<numCols;j++) {
                psF32 expect = ((float) (i + j)) - 12.0;
                psF32 actual = rc->image->data.F32[i][j];
                if (FLT_EPSILON < fabs(expect - actual)) {
                    printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                    testStatus = 1;
                }
            }
        }
    */

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() undersize bias image (short rows).  Should generate Error.\n");
    myReadout = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, NULL,
                               0, PM_FIT_NONE, myBiasShortRows);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }


    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() undersize bias image (short columns).  Should generate Error.\n");
    myReadout = pmSubtractBias(myReadout, NULL, NULL, PM_OVERSCAN_NONE, NULL,
                               0, PM_FIT_NONE, myBiasShortCols);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with bogus PM_FIT.  Should generate Error.\n");
    myReadout = pmSubtractBias(myReadout, NULL, list, PM_OVERSCAN_ROWS, stat,
                               0, 54321, NULL);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    printf("------------------------------------------------------------------\n");
    printf("Calling pmSubtractBias() with bogus overScanAxis.  Should generate Error.\n");
    myReadout = pmSubtractBias(myReadout, NULL, list, 54321, stat,
                               0, PM_FIT_NONE, NULL);
    if (rc != myReadout) {
        printf("TEST ERROR: pmSubtractBias() did not return input psReadout.\n");
        testStatus = false;
        psFree(rc);
    }

    /*
        for (i=0;i<numRows;i++) {
            for (j=0;j<numCols;j++) {
                psF32 expect = ((float) (i + j)) - 12.0;
                psF32 actual = rc->image->data.F32[i][j];
                if (FLT_EPSILON < fabs(expect - actual)) {
                    printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                    testStatus = 1;
                }
            }
        }
    */

    printf("------------------------------------------------------------------\n");
    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(tmpImage2Short);
    psFree(tmpImage3Short);
    psFree(tmpImage4Short);
    psFree(myBias);
    psFree(myBiasShortRows);
    psFree(myBiasShortCols);
    psFree(stat);
    psFree(list);
    psFree(listShort);

    printFooter(stdout, "pmSubtractBias", "Testing input parameter error conditions", true);
    return(testStatus);
}

int test04( void )
{
    int testStatus = 0;

    testStatus |= doSubtractOverscansTestInputCases(NUM_COLS, NUM_ROWS);
    return(testStatus);
}

/******************************************************************************
doSubtractFullOverscansSmall(): a sample psReadout as well as several overscan
images of smaller size are created.  The overscan images are collected
pixel-by-pixel then subtracted column-wise from the psReadout.
 *****************************************************************************/
int doSubtractFullOverscanColumnsGeneric(int imageNumCols,
        int imageNumRows,
        int overscanNumCols,
        int overscanNumRows,
        pmOverscanAxis overscanaxis,
        pmFit fit,
        psS32 nBin)
{
    int i;
    int j;
    float actual;
    float expect;
    int testStatus = 0;

    psImage *tmpImage1 = psImageAlloc(imageNumCols, imageNumRows, PS_TYPE_F32);
    //    psReadout *myReadout = psReadoutAlloc(imageNumCols, imageNumRows, tmpImage1);
    psReadout *myReadout = psReadoutAlloc();
    myReadout->image = tmpImage1;
    for (i=0;i<imageNumRows;i++) {
        for (j=0;j<imageNumCols;j++) {
            myReadout->image->data.F32[i][j] = (float) (i + j);
        }
    }

    psStats *stat = psStatsAlloc(PS_STAT_SAMPLE_MEAN);
    psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD);
    psSpline1D *mySpline = NULL;


    /*
        if (overscanNumRows <= 0) {
            overscanNumRows = 1;
        }
        if (overscanNumCols <= 0) {
            overscanNumCols = 1;
        }
    */
    psImage *tmpImage2 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    psImage *tmpImage3 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    psImage *tmpImage4 = psImageAlloc(overscanNumCols, overscanNumRows, PS_TYPE_F32);
    psList *list;
    for (i=0;i<overscanNumRows;i++) {
        for (j=0;j<overscanNumCols;j++) {
            tmpImage2->data.F32[i][j] = 3.0;
            tmpImage3->data.F32[i][j] = 4.0;
            tmpImage4->data.F32[i][j] = 5.0;
        }
    }
    list = psListAlloc(tmpImage2);
    psListAdd(list, PS_LIST_HEAD, tmpImage3);
    psListAdd(list, PS_LIST_HEAD, tmpImage4);

    printPositiveTestHeader(stdout, "pmSubtractBias", "Column Overscans");

    if (fit == PM_FIT_NONE) {
        myReadout = pmSubtractBias(myReadout, NULL, list, overscanaxis, stat,
                                   nBin, fit, NULL);
    } else if (fit == PM_FIT_POLYNOMIAL) {
        myReadout = pmSubtractBias(myReadout, myPoly, list, overscanaxis, stat,
                                   nBin, fit, NULL);
    } else if (fit == PM_FIT_SPLINE) {
        myReadout = pmSubtractBias(myReadout, mySpline, list, overscanaxis, stat,
                                   nBin, fit, NULL);
    }

    for (i=0;i<imageNumRows;i++) {
        for (j=0;j<imageNumCols;j++) {
            expect = ((float) (i + j)) - 12.0;
            actual = myReadout->image->data.F32[i][j];
            if (FLT_EPSILON < fabs(expect - actual)) {
                printf("TEST ERROR: image[%d][%d] is %f, should be %f\n", i, j, actual, expect);
                testStatus = 1;
            }
        }
    }


    psFree(myReadout);
    psFree(tmpImage2);
    psFree(tmpImage3);
    psFree(tmpImage4);
    psFree(stat);
    psFree(list);
    psFree(myPoly);
    psFree(mySpline);

    printFooter(stdout, "pmSubtractBias", "Column Overscans", true);
    return(testStatus);
}

int testX( void )
{
    int testStatus = 0;

    // imageNumCols, imageNumRows, overscanNumCols, overscanNumRows,
    // overscanaxis, fit, nBin

    //
    // Overscan images are same size, no fit, bin factor is 1.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_NONE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_NONE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_NONE, 1);

    //
    // Overscan images are same size, no fit, bin factor is 2.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_NONE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_NONE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_NONE, 2);
    //
    // Overscan images are too small, spline fit, bin factor is 1.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_ALL,
                  PM_FIT_SPLINE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_SPLINE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_SPLINE, 1);

    //
    // Overscan images are too small, spline fit, bin factor is 2.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_ALL,
                  PM_FIT_SPLINE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_SPLINE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS-1, NUM_ROWS-1,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_SPLINE, 2);

    //
    // Overscan images are same size, spline fit, bin factor is 1.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_SPLINE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_SPLINE, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_SPLINE, 1);
    //
    // Overscan images are same size, spline fit, bin factor is 2.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_SPLINE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_SPLINE, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_SPLINE, 2);

    //
    // Overscan images are same size, polynomial fit, bin factor is 1.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_POLYNOMIAL, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_POLYNOMIAL, 1);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_POLYNOMIAL, 1);
    //
    // Overscan images are same size, polynomial fit, bin factor is 2.
    //
    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ALL,
                  PM_FIT_POLYNOMIAL, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_COLUMNS,
                  PM_FIT_POLYNOMIAL, 2);

    testStatus |= doSubtractFullOverscanColumnsGeneric(NUM_COLS, NUM_ROWS,
                  NUM_COLS, NUM_ROWS,
                  PM_OVERSCAN_ROWS,
                  PM_FIT_POLYNOMIAL, 2);

    return(testStatus);
}

