/** @file  tst_psImageManip.c
 *
 *  @brief Contains the tests for psImageManip.[ch]
 *
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.29 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-14 01:22:59 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */


#include <math.h>
#include <float.h>
#include <string.h>
#include <stdlib.h>
#include <string.h>                    // for memset
#include <sys/stat.h>
#include <sys/types.h>

#include "psTest.h"
#include "pslib.h"
#include "psType.h"

static int testImageClip(void);
static int testImageClipNAN(void);
static int testImageClipComplexRegion(void);
static int testImageOverlay(void);
static int testImageRebin(void);
static int testImageRoll(void);
static int testImageRotate(void);
static int testImageShift(void);
static int testImageShiftCase(int cols, int rows, float colShift,float rowShift);
static int testImageResample(void);

testDescription tests[] = {
                              {testImageClip,571,"psImageClip",0,false},
                              {testImageClipNAN,572,"psImageClipNAN",0,false},
                              {testImageClipComplexRegion,673,"psImageClipComplexRegion",0,false},
                              {testImageOverlay,573,"psImageOverlay",0,false},
                              {testImageRebin,559,"psImageRebin",0,false},
                              {testImageRoll,562,"psImageRoll",0,false},
                              {testImageRotate,560,"psImageRotate",0,false},
                              {testImageShift,561,"psImageShift",0,false},
                              {testImageResample,743,"psImageResample",0,false},
                              {NULL}
                          };

int main(int argc, char* argv[])
{
    psLogSetLevel(PS_LOG_INFO);

    return ! runTestSuite(stderr,"psImage",tests,argc,argv);
}


int testImageClip(void)
{
    psImage* img = NULL;
    unsigned int c = 128;
    unsigned int r = 256;
    psF64 min;
    psF64 max;
    int numClipped = 0;
    int retVal;

    psLogMsg(__func__,PS_LOG_INFO,
             "psImageClip shall limit the minimum and maximum data value within a psImage structure");

    /*

        psImageClip shall limit the minimum and maximum data value within a
        psImage structure to a specified min and max value.

        Verify the returned integer is equal to the number of pixels clipped,
        if the input psImage structure contains known values and input parameters
        min and max have know values.

        Verify the psImage structure specified by the input parameter input is
        modified to contain the expected values, if the input psImage structure
        contains known values, min and max are specified and vmin and vmax
        parameters are known.

        Verify the retuned integer is zero, psImage structure input is unmodified
        and program executions doesn't stop, if input parameter psImage structure
        pointer is null.

        Verify the retuned integer is zero, psImage structure input is unmodified
        and program executions doesn't stop, if input parameter min is larger than max.
    */

    // create image
    #define testImageClipByType(datatype) \
    img = psImageAlloc(c,r,PS_TYPE_##datatype); \
    for (unsigned int row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned int col=0;col<c;col++) { \
            imgRow[col] = (ps##datatype)(row+col); \
        } \
    } \
    min = (psF64)r/2.0; \
    max = (psF64)r; \
    \
    retVal = psImageClip(img,min,(double)PS_MIN_##datatype,max,(double)PS_MAX_##datatype); \
    \
    numClipped = 0; \
    for (unsigned int row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned int col=0;col<c;col++) { \
            ps##datatype value = (ps##datatype)(row+col); \
            if (value < min) { \
                numClipped++; \
                value = PS_MIN_##datatype; \
            } else if (value > max) { \
                numClipped++; \
                value = PS_MAX_##datatype; \
            } \
            if (fabsf(imgRow[col]-value) > FLT_EPSILON) { \
                psError(__func__,"Pixel value is not as expected (%g vs %g) at %u,%u", \
                        (psF64)imgRow[col],(psF64)value,col,row); \
                return 1; \
            } \
        } \
    } \
    if (retVal != numClipped) { \
        psError(__func__,"Expected %d clips, but got %d", \
                numClipped,retVal); \
        return 2; \
    } \
    psFree(img);

    #define testImageClipByComplexType(datatype) \
    img = psImageAlloc(c,r,PS_TYPE_##datatype); \
    for (unsigned int row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned int col=0;col<c;col++) { \
            imgRow[col] = (ps##datatype)(row+I*col); \
        } \
    } \
    min = (float)r/2.0f; \
    max = (float)r; \
    \
    retVal = psImageClip(img,min,-1.0f,max,-2.0f); \
    \
    numClipped = 0; \
    for (unsigned int row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned int col=0;col<c;col++) { \
            ps##datatype value = row+I*col; \
            if (cabs(value) < min) { \
                numClipped++; \
                value = -1.0f; \
            } else if (cabs(value) > max) { \
                numClipped++; \
                value = -2.0f; \
            } \
            if (fabsf(creal(imgRow[col])-creal(value)) > FLT_EPSILON || \
                    fabsf(cimag(imgRow[col])-cimag(value)) > FLT_EPSILON) { \
                psError(__func__,"Pixel value is not as expected (%.2f+%.2fi vs %.2f+%.2fi) at %u,%u", \
                        creal(imgRow[col]),cimag(imgRow[col]),creal(value),cimag(value),col,row); \
                return 1; \
            } \
        } \
    } \
    if (retVal != numClipped) { \
        psError(__func__,"Expected %d clips, but got %d", \
                numClipped,retVal); \
        return 2; \
    } \
    psFree(img);

    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of F64 imagery");
    testImageClipByType(F64);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of F32 imagery");
    testImageClipByType(F32);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of S32 imagery");
    testImageClipByType(S32);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of S16 imagery");
    testImageClipByType(S16);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of S8 imagery");
    testImageClipByType(S8);
    //    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of U32 imagery");
    //    testImageClipByType(U32);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of U16 imagery");
    testImageClipByType(U16);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of U8 imagery");
    testImageClipByType(U8);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of C32 imagery");
    testImageClipByComplexType(C32);
    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping of C64 imagery");
    testImageClipByComplexType(C64);

    // Verify the retuned integer is zero, psImage structure input is unmodified
    // and program executions doesn't stop, if input parameter psImage structure
    // pointer is null.
    retVal = psImageClip(NULL,min,-1.0f,max,-2.0f);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clips of a NULL image.");
        return 3;
    }

    // Verify the retuned integer is zero, psImage structure input is unmodified
    // and program executions doesn't stop, if input parameter min is larger than max.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error (max<min)");
    retVal = psImageClip(img,max,-1.0f,min,-2.0f);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clips when max < min.");
        return 4;
    }

    // Verify program execution doesn't stop if the input image type is something
    // other than U8, U16, S8, S16, F32, F64, C32, C64.
    img = psImageAlloc(c,r,PS_TYPE_U32);
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error (invalid type)");
    retVal = psImageClip(img,min,-1.0f,max,-2.0f);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clip of invalid image type.");
        return 5;
    }
    psFree(img);

    return 0;
}

int testImageClipNAN(void)
{
    psImage* img = NULL;
    unsigned int c = 128;
    unsigned int r = 256;
    int numClipped = 0;
    int retVal;

    psLogMsg(__func__,PS_LOG_INFO,
             "psImageClipNaN shall modified pixel values of NaN with a specified value");

    /*
        psImageClipNaN shall modify a psImage structure with pixel values set
        to NaN to a value specified as an input parameter.

        Verify the returned integer is equal to the number of pixels modified
        and the psImage is modified at locations where NaN pixels where
        located to the value specified in the input parameter value.

        Verify the returned integer is zero and program execution doesn't stop,
        if the input parameter psImage structure pointer is null.
    */

    // create image
    #define testImageClipNaNByType(datatype) \
    img = psImageAlloc(c,r,PS_TYPE_##datatype); \
    for (unsigned row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned col=0;col<c;col++) { \
            if (row == col) { \
                imgRow[col] = NAN; \
            } else if (row+1 == col) { \
                imgRow[col] = INFINITY; \
            } else { \
                imgRow[col] = (ps##datatype)(row+col); \
            } \
        } \
    } \
    \
    retVal = psImageClipNaN(img,-1.0f); \
    \
    numClipped = 0; \
    for (unsigned row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned col=0;col<c;col++) { \
            ps##datatype value = (ps##datatype)(row+col); \
            if ( (row == col) || (row+1 == col) ) { \
                numClipped++; \
                value = -1.0; \
            } \
            if (fabsf(imgRow[col]-value) > FLT_EPSILON) { \
                psError(__func__,"Pixel value is not as expected (%f vs %f) at %d,%d", \
                        imgRow[col],value,col,row); \
                return 1; \
            } \
        } \
    } \
    if (retVal != numClipped) { \
        psError(__func__,"Expected %d clips, but got %d", \
                numClipped,retVal); \
        return 2; \
    } \
    psFree(img);

    testImageClipNaNByType(F32);
    testImageClipNaNByType(F64);
    testImageClipNaNByType(C32);
    testImageClipNaNByType(C64);

    // Verify the retuned integer is zero, psImage structure input is unmodified
    // and program executions doesn't stop, if input parameter psImage structure
    // pointer is null.
    retVal = psImageClipNaN(NULL,-1.0f);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clips of a NULL image.");
        return 3;
    }

    // Verify program execution doesn't stop if the input image type is something
    // other than F32, F64, C32, C64.
    img = psImageAlloc(c,r,PS_TYPE_S32);
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error (invalid type)");
    retVal = psImageClipNaN(img,2.0f);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clip of invalid image type.");
        return 4;
    }
    psFree(img);

    return 0;
}

int testImageClipComplexRegion(void)
{
    psImage* img = NULL;
    unsigned int c = 1024;
    unsigned int r = 2048;
    int numClipped = 0;
    int retVal;

    psLogMsg(__func__,PS_LOG_INFO,
             "psImageClipNaN shall modified pixel values of NaN with a specified value");

    /*
    1. Create a complex image with a wide range of complex values

    2. call psImageClipComplexRegion with min and max where there is at least 
       2 pixels in the image above) that is:
        a) real(p) < real(min) && complex(p) < complex(min),
        b) real(p) < real(min) && complex(min) < complex(p) < complex(max)
        c) real(min) < real(p) < real(max) && complex(p) < complex(min)
        d) real(min) < real(p) < real(max) && complex(min) < complex(p) < complex(max)
        e) real(p) > real(max) && complex(min) < complex(p) < complex(max)
        f) real(pmin) < real(p) < real(max) && complex(p) > complex(max)
        g) real(p) > real(max) && complex(p) > complex(max)
        h) real(p) < real(min) && complex(p) > complex(max)
        i) real(p) > real(max) && complex(p) < complex(min)

    3. verify that All pixels in case (a), (b), and (c) have the value vmin

    4. verify that all pixels in case (d) are unchanged from input

    5. verify that all pixels in case (e), (f), (g), (h), and (i) have the 
       value vmax

    */

    #define testImageClipComplexByType(datatype,MIN,MAX) /* datatype must be complex */ \
    /* create image */ \
    img = psImageAlloc(c,r,PS_TYPE_##datatype); \
    for (unsigned row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned col=0;col<c;col++) { \
            imgRow[col] = row+I*col; \
        } \
    } \
    \
    retVal = psImageClipComplexRegion(img,MIN,-1.0-1.0*I,MAX,-2.0-2.0*I); \
    \
    numClipped = 0; \
    for (unsigned row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (unsigned col=0;col<c;col++) { \
            ps##datatype value = (ps##datatype)(row+I*col); \
            if ( (row > creal(MAX)) || (col > cimag(MAX)) ) { \
                numClipped++; \
                value = -2.0-2.0*I; \
            } else if ((row < creal(MIN)) || (col < cimag(MIN)) ) { \
                numClipped++; \
                value = -1.0-1.0*I; \
            } \
            if (cabs(imgRow[col]-value) > FLT_EPSILON) { \
                psError(__func__,"Pixel value is not as expected (%g%+gi vs %g%+gi) at %d,%d", \
                        creal(imgRow[col]),cimag(imgRow[col]),creal(value),cimag(value),col,row); \
                return 1; \
            } \
        } \
    } \
    if (retVal != numClipped) { \
        psError(__func__,"Expected %d clips, but got %d", \
                numClipped,retVal); \
        return 2; \
    } \
    psFree(img);

    complex double min = ((double)r)/5.0+I*((double)c)/4.0;
    complex double max = ((double)r)/3.0+I*((double)c)/2.0;

    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping at %g%+gi to %g%+gi for psC32",
             creal(min),cimag(min),creal(max),cimag(max));

    testImageClipComplexByType(C32,min,max);

    psLogMsg(__func__,PS_LOG_INFO,"Testing clipping at %g%+gi to %g%+gi for psC64",
             creal(min),cimag(min),creal(max),cimag(max));
    testImageClipComplexByType(C64,min,max);

    //  6. Call psImageClipComplexRegion with NULL input parameter; should error
    //     but not stop execution.

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(NULL,0,0,0,0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clips of a NULL image.");
        return 3;
    }


    img = psImageAlloc(c,r,PS_TYPE_C32);
    for (unsigned row=0;row<r;row++) {
        psC32* imgRow = img->data.C32[row];
        for (unsigned col=0;col<c;col++) {
            imgRow[col] = row+I*col;
        }
    }

    //  7. Call psImageClipComplexRegion with min > max; should error and return 0,
    //     but not stop execution

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,10.0+I*1.0,-1.0,5.0+5.0*I,-2.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for creal(min)>creal(max).");
        return 3;
    }
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,1.0+I*10.0,-1.0,5.0+5.0*I,-2.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for cimag(min)>cimag(max).");
        return 3;
    }
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,10.0+I*10.0,-1.0,5.0+5.0*I,-2.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for min>max.");
        return 3;
    }

    //  8. Call psImageClipComplexRegion with the follow vmin/vmax values; each
    //     should error and return 0, but not stop execution
    //      a) vmin < datatype region's minimum
    //      b) vmax < datatype region's minimum
    //      c) vmin > datatype region's maximum
    //      d) vmax > datatype region's maximum

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      -2.0*(double)FLT_MAX,
                                      5.0+5.0*I,
                                      0.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmin not in datatype range.");
        return 80;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      2.0*(double)FLT_MAX,
                                      5.0+5.0*I,
                                      0.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmin not in datatype range.");
        return 81;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      -2.0*(double)FLT_MAX*I,
                                      5.0+5.0*I,
                                      0.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmin not in datatype range.");
        return 82;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      2.0*(double)FLT_MAX*I,
                                      5.0+5.0*I,
                                      0.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmin not in datatype range.");
        return 83;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      0.0,
                                      5.0+5.0*I,
                                      -2.0*(double)FLT_MAX);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmax not in datatype range.");
        return 84;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      0.0,
                                      5.0+5.0*I,
                                      2.0*(double)FLT_MAX);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmax not in datatype range.");
        return 85;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      0.0,
                                      5.0+5.0*I,
                                      -2.0*(double)FLT_MAX*I);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmax not in datatype range.");
        return 87;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error:");
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      0.0,
                                      5.0+5.0*I,
                                      2.0*(double)FLT_MAX*I);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for vmax not in datatype range.");
        return 88;
    }


    // now check if vmin > vmax is OK
    for (unsigned row=0;row<r;row++) {
        psC32* imgRow = img->data.C32[row];
        for (unsigned col=0;col<c;col++) {
            imgRow[col] = row+I*col;
        }
    }
    retVal = psImageClipComplexRegion(img,
                                      1.0+I*1.0,
                                      10.0,
                                      5.0+5.0*I,
                                      0.0);
    if (retVal == 0) {
        psError(__func__,"Didn't expect zero return for vmin > vmax.");
        return 83;
    }


    psFree(img);
    img = NULL;

    //  9. Call psImageClipComplexRegion with the max value out of datatype's
    //     range; should clip as expected (see step 1-5). Repeat with min value
    //     out of datatype's range.

    testImageClipComplexByType(C32,-(double)FLT_MAX*2.0-I*(double)FLT_MAX*2.0,10.0+I*10.0);
    testImageClipComplexByType(C32,10.0+I*10.0,(double)FLT_MAX*2.0+I*(double)FLT_MAX*2.0);

    // Verify program execution doesn't stop if the input image type is something
    // other than C32, C64.
    img = psImageAlloc(c,r,PS_TYPE_S32);
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error (invalid type)");
    retVal = psImageClipComplexRegion(img,2.0,10.0,5.0,0.0);
    if (retVal != 0) {
        psError(__func__,"Expected zero return for clip of invalid image type.");
        return 84;
    }
    psFree(img);

    return 0;
}

int testImageOverlay(void)
{

    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    psImage* img4 = NULL;
    unsigned int c = 128;
    unsigned int r = 256;
    int retVal;

    /*
    psImageSectionOverlay shall modified pixel values in a psImage structure to
    be equal to the value of the originaldata and an overlay image with a
    specified operation. Valid operations include =, +, -, *, /.

    Verify the returned integer is zero
    and the input parameter psImage structure is modified at the specified
    location and range with the given overlay image and the specified
    function. Cases should include all the valid operations. Comparison of
    expected values should include a delta to allow for testing on
    different platforms.

    */



    #define testOverlayTypeOP(DATATYPE,OP,OPSTRING) \
    img = psImageAlloc(c,r,PS_TYPE_##DATATYPE); \
    for (unsigned row=0;row<r;row++) { \
        ps##DATATYPE* imgRow = img->data.DATATYPE[row]; \
        for (unsigned col=0;col<c;col++) { \
            imgRow[col] = 6.0; \
        } \
    } \
    img2 = psImageAlloc(c/2,r/2,PS_TYPE_##DATATYPE); \
    for (unsigned row=0;row<r/2;row++) { \
        ps##DATATYPE* img2Row = img2->data.DATATYPE[row]; \
        for (unsigned col=0;col<c/2;col++) { \
            img2Row[col] = 2.0; \
        } \
    } \
    retVal = psImageOverlaySection(img,img2,c/4,r/4,OPSTRING); \
    if (retVal != 0) { \
        psError(__func__,"psImageOverlaySection returned non-zero %s op",OPSTRING); \
        return 1; \
    } \
    for (unsigned row=0;row<r;row++) { \
        ps##DATATYPE* imgRow = img->data.DATATYPE[row]; \
        ps##DATATYPE* img2Row = img2->data.DATATYPE[row]; \
        for (unsigned col=0;col<c;col++) { \
            ps##DATATYPE val = 6.0; \
            if ( ! (row < r/4 || row >= r/2+r/4 || col < c/4 || col >= c/2+c/4)) { \
                val OP 2.0; \
            } \
            if (fabsf(imgRow[col] - val) > FLT_EPSILON) { \
                psError(__func__,"Value incorrect at %d,%d (%.2f vs %.2f for %s)", \
                        col,row,imgRow[col],val,OPSTRING); \
                return 2; \
            } \
            if (row < r/2 && col < c/2 && fabsf(img2Row[col] - 2.0) > FLT_EPSILON) { \
                psError(__func__,"Overlay modified at %d,%d (%.2f for %s)", \
                        col,row,img2Row[col],OPSTRING); \
                return 2; \
            } \
        } \
    } \
    psFree(img); \
    psFree(img2);

    #define testOverlayType(DATATYPE) \
    testOverlayTypeOP(DATATYPE,+=,"+"); \
    testOverlayTypeOP(DATATYPE,-=,"-"); \
    testOverlayTypeOP(DATATYPE,*=,"*");\
    testOverlayTypeOP(DATATYPE,/=,"/");\
    testOverlayTypeOP(DATATYPE,=,"=");

    testOverlayType(C64);
    testOverlayType(C32);
    testOverlayType(F64);
    testOverlayType(F32);
    testOverlayType(S16);
    testOverlayType(S8);
    testOverlayType(U16);
    testOverlayType(U8);

    /*
    Verify the returned integer is equal to non-zero and the input psImage structure
    is unmodified, if the overlay specified is not within the data range of the
    input psImage structure.
    */

    img = psImageAlloc(c,r,PS_TYPE_F32);
    for (unsigned row=0;row<r;row++) {
        psF32* imgRow = img->data.F32[row];
        for (unsigned col=0;col<c;col++) {
            imgRow[col] = 6.0f;
        }
    }
    img2 = psImageAlloc(c,r,PS_TYPE_F32);
    for (unsigned row=0;row<r;row++) {
        psF32* img2Row = img2->data.F32[row];
        for (unsigned col=0;col<c;col++) {
            img2Row[col] = 2.0f;
        }
    }
    img3 = psImageAlloc(c,r,PS_TYPE_S64);
    for (unsigned row=0;row<r;row++) {
        psS64* img3Row = img3->data.S64[row];
        for (unsigned col=0;col<c;col++) {
            img3Row[col] = 6.0f;
        }
    }
    img4 = psImageAlloc(c,r,PS_TYPE_S64);
    for (unsigned row=0;row<r;row++) {
        psS64* img4Row = img4->data.S64[row];
        for (unsigned col=0;col<c;col++) {
            img4Row[col] = 2.0f;
        }
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should error as overlay isn't "
             "within image boundaries");
    retVal = psImageOverlaySection(img,img2,c/4,r/4,"+");
    if (retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero even though "
                "overlay too big");
        return 3;
    }
    for (unsigned row=0;row<r;row++) {
        psF32* imgRow = img->data.F32[row];
        for (unsigned col=0;col<c;col++) {
            if (imgRow[col] != 6.0f) {
                psError(__func__,"Input image modified when overlay size too big");
                return 4;
            }
        }
    }

    /*
    Verify the returned integer is equal to non-zero, the input psImage
    structure is unmodified and program execution doesn't stop, if the
    overlay specified is null.
    */

    psLogMsg(__func__,PS_LOG_INFO,"Following should error as overlay is NULL");
    retVal = psImageOverlaySection(img,NULL,c/4,r/4,"+");
    if (retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero even though "
                "overlay too big");
        return 5;
    }
    for (unsigned row=0;row<r;row++) {
        psF32* imgRow = img->data.F32[row];
        for (unsigned col=0;col<c;col++) {
            if (imgRow[col] != 6.0f) {
                psError(__func__,"Input image modified when overlay NULL");
                return 6;
            }
        }
    }

    /*
    Verify the returned integer is equal to non-zero and program execution
    doesn't stop, if the input parameter image is null.
    */

    psLogMsg(__func__,PS_LOG_INFO,"Following should error as image input is NULL");
    retVal = psImageOverlaySection(NULL,img2,c/4,r/4,"+");
    if (retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero even though "
                "overlay too big");
        return 7;
    }

    /*
    Verify the return integer is equal to non-zero and program execution
    doesn't stop, if the specified operator is not =,+,-,*,/
    */

    psLogMsg(__func__,PS_LOG_INFO,"Following should error as operator is invalid");
    retVal = psImageOverlaySection(img,img2,0,0,"$");
    if (retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero even though "
                "overlay operator is invalid");
        return 8;
    }

    /*
    Verify the return integer is equal to non-zero and program execution
    doesn't stop, if the specified operator is NULL
    */

    psLogMsg(__func__,PS_LOG_INFO,"Following should error as operator is invalid");
    retVal = psImageOverlaySection(img,img2,0,0,NULL);
    if(retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero even though "
                "overlay operator is NULL");
        return 9;
    }

    /*
    Verify the return integer is equal to non-zero and program execution 
    doesn't stop, if overlay image is a different type than the input image
    */
    psLogMsg(__func__,PS_LOG_INFO,"Following should error as overlay is "
             "a different type");
    retVal = psImageOverlaySection(img,img3,0,0,"+");
    if(retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero eventhough "
                " overlay image type is different than input image.");
        return 10;
    }

    /*
    Verify the return integer is equal to non-zero and program execution
    doesn't stop, if invalid input and overlay images
    */
    psLogMsg(__func__,PS_LOG_INFO,"Following should error as overlay and "
             "image are not valid for the psImageOverlaySection function.");
    retVal = psImageOverlaySection(img3,img4,0,0,"+");
    if(retVal == 0) {
        psError(__func__,"psImageOverlaySection returned zero eventhough "
                "overlay and image type are invalid type.");
        return 11;
    }

    /*
    Verify program execution doen't stop, if the overly image contains
    zero values with division operation is specified.
    */
    for (unsigned row=0;row<r;row++) {
        psF32* img2Row = img2->data.F32[row];
        for (unsigned col=0;col<c;col++) {
            img2Row[col] = 0.0f;
        }
    }
    retVal = psImageOverlaySection(img,img2,0,0,"/");
    if (retVal != 0) {
        psError(__func__,"psImageOverlaySection returned non-zero when "
                "checking divide-by-zero.");
        return 12;
    }

    psFree(img);
    psFree(img2);
    psFree(img3);
    psFree(img4);

    return 0;
}

static int testImageRebin(void)
{

    /*
    This function shall generate a rescaled version of a psImage structure 
    derived from a specified statistics method.
    */

    psImage* in = NULL;
    psImage* out = NULL;
    psImage* out2 = NULL;
    psImage* out3 = NULL;
    psImage* mask = NULL;
    psImage* meanTruth = NULL;
    psImage* meanTruthWMask = NULL;
    psImage* maxTruth = NULL;
    psStats stats;

    /*
    Verify the returned psImage structure contains expected values, if the 
    input parameter input contains known data, the input scale is a known 
    value with a known statistical method specified in stats. Cases should 
    include at least two different scales and statistical methods. Comparison 
    of expected values should include a delta to allow testing on different 
    platforms.
    */

    #define testRebinType(DATATYPE)  \
    in = psImageAlloc(16,16,PS_TYPE_##DATATYPE); \
    mask = psImageAlloc(16,16,PS_TYPE_U8); \
    meanTruth = psImageAlloc(4,4,PS_TYPE_F32); \
    meanTruthWMask = psImageAlloc(4,4,PS_TYPE_F32); \
    maxTruth = psImageAlloc(6,6,PS_TYPE_F32); \
    memset(meanTruth->data.F32[0],0,sizeof(psF32)*4*4); \
    memset(meanTruthWMask->data.F32[0],0,sizeof(psF32)*4*4); \
    memset(maxTruth->data.F32[0],0,sizeof(psF32)*6*6); \
    for (int row = 0; row<16; row++) { \
        ps##DATATYPE* inRow = in->data.DATATYPE[row]; \
        psF32* meanTruthRow = meanTruth->data.F32[row/4]; \
        psF32* meanTruthWMaskRow = meanTruthWMask->data.F32[row/4]; \
        psF32* maxTruthRow = maxTruth->data.F32[row/3]; \
        psU8* maskRow = mask->data.U8[row]; \
        for (int col = 0; col<16; col++) { \
            if(col != 15) { \
                maskRow[col] = 0; \
            } else { \
                maskRow[col] = 1; \
            } \
            inRow[col] = row + col; \
            meanTruthRow[col/4] += row + col; \
            if (maxTruthRow[col/3] < row + col) { \
                maxTruthRow[col/3] = row+col; \
            } \
            if(maskRow[col] == 0 ) { \
                meanTruthWMaskRow[col/4] += row + col; \
            } \
        } \
    } \
    for (int row = 0; row<4; row++) { \
        psF32* meanTruthRow = meanTruth->data.F32[row]; \
        psF32* meanTruthWMaskRow = meanTruthWMask->data.F32[row]; \
        for (int col = 0; col<4; col++) { \
            meanTruthRow[col] /= 16; \
            if ( col == 3 ) { \
                meanTruthWMaskRow[col] /= 12; \
            } else { \
                meanTruthWMaskRow[col] /= 16; \
            } \
        } \
    } \
    stats.options = PS_STAT_SAMPLE_MEAN; \
    out = psImageRebin(NULL,in,NULL,0,4,&stats); \
    if (out == NULL) { \
        psError(__func__,"psImageRebin returned a NULL pointer!?"); \
        return 1; \
    } \
    if (out->numRows != 4 || out->numCols != 4) { \
        psError(__func__,"psImageRebin didn't produce the proper size image " \
                "(%d x %d).", \
                out->numCols, out->numRows); \
        return 2; \
    } \
    for (int row = 0; row<4; row++) { \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        psF32* truthRow = meanTruth->data.F32[row]; \
        for (int col = 0; col<4; col++) { \
            if (fabsf((float)outRow[col]-(float)truthRow[col]) > FLT_EPSILON) { \
                psError(__func__,"psImageRebin didn't produce the proper mean " \
                        "result at (%d,%d) [%f vs %f].", \
                        col,row,outRow[col],truthRow[col]); \
                return 3; \
            } \
        } \
    } \
    stats.options = PS_STAT_SAMPLE_MEAN; \
    out3 = psImageRebin(NULL,in,mask,1,4,&stats); \
    for (int row = 0; row<4; row++) { \
        ps##DATATYPE* outRow = out3->data.DATATYPE[row]; \
        psF32* truthRow = meanTruthWMask->data.F32[row]; \
        for ( int col = 0; col<4; col++) { \
            if(abs((int)outRow[col]-(int)truthRow[col]) > FLT_EPSILON) { \
                psError(__func__,"psImageRebin with mask didn't produce the proper mean " \
                        "result at (%d,%d) [%f vs %f].", \
                        col,row,outRow[col],truthRow[col]); \
                return 3; \
            } \
        } \
    } \
    stats.options = PS_STAT_MAX; \
    out2 = psImageRebin(out,in,NULL,0,3,&stats); \
    if (out != out2) { \
        psError(__func__,"psImageRebin didn't recycle a psImage properly!?"); \
        return 7; \
    } \
    if (out == NULL) { \
        psError(__func__,"psImageRebin returned a NULL pointer!?"); \
        return 4; \
    } \
    if (out->numRows != 6 || out->numCols != 6) { \
        psError(__func__,"psImageRebin didn't produce the proper size image " \
                "(%d x %d).", \
                out->numCols, out->numRows); \
        return 5; \
    } \
    for (int row = 0; row<6; row++) { \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        psF32* truthRow = maxTruth->data.F32[row]; \
        for (int col = 0; col<6; col++) { \
            if (fabsf((float)outRow[col]-(float)truthRow[col]) > FLT_EPSILON) { \
                psError(__func__,"psImageRebin didn't produce the proper " \
                        "max result at (%d,%d) [%f vs %f].", \
                        col,row,outRow[col],truthRow[col]); \
                return 6; \
            } \
        } \
    } \
    psFree(in); \
    psFree(out); \
    psFree(out3); \
    psFree(mask); \
    psFree(meanTruth); \
    psFree(meanTruthWMask); \
    psFree(maxTruth);

    testRebinType(F32);
    testRebinType(F64);
    testRebinType(U16);
    testRebinType(S8);

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input image type is not supported.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for unsupported type.");
    in = psImageAlloc(16,16,PS_TYPE_U8);
    mask = psImageAlloc(16,16,PS_TYPE_F32);
    stats.options = PS_STAT_SAMPLE_MEAN;
    out = psImageRebin(NULL,in,NULL,0,4,&stats);
    if(out != NULL) {
        psError(__func__,"psImageRebin return an image eventhough the "
                "type is not handled.");
        return 14;
    }
    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the mask type is not U8
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for invallid mask type.");
    out = psImageRebin(NULL,in,mask,1,4,&stats);
    if(out != NULL) {
        psError(__func__,"psImageRebin return an image eventhough the "
                "mask is not the correct type.");
        return 17;
    }
    psFree(mask);
    psFree(in);

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input parameter input is null.

    out2 = psImageRebin(NULL,NULL,NULL,0,1,&stats);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the input was "
                "NULL!?");
        return 8;
    }

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input parameter scale is less than or equal to zero.
    in = psImageAlloc(16, 16, PS_TYPE_F32);
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for scale < 0.");
    out2 = psImageRebin(NULL,in,NULL,0,0,&stats);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the scale was "
                "zero!?");
        return 9;
    }

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input parameter stats is null.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for stats null.");
    out2 = psImageRebin(NULL,in,NULL,0,1,NULL);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the stats was "
                "NULL!?");
        return 10;
    }

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input parameter psStats structure member options
    // is zero or any value which doesn't correspond to a valid statistical
    // method.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for stats options 0.");
    stats.options = 0;
    out2 = psImageRebin(NULL,in,NULL,0,1,&stats);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the stats "
                "options was zero!?");
        return 11;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for stat options use range.");
    stats.options = PS_STAT_USE_RANGE;
    out2 = psImageRebin(NULL,in,NULL,0,1,&stats);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the stats "
                "options was PS_STAT_USE_RANGE!?");
        return 12;
    }

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the input parameter psStats structure member options
    // specifies more than one valid statistical method.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for stats with multiple options.");
    stats.options = PS_STAT_SAMPLE_MEAN + PS_STAT_MAX;
    out2 = psImageRebin(NULL,in,NULL,0,1,&stats);

    if (out2 != NULL) {
        psError(__func__,"psImageRebin returned an image though the stats "
                "options was PS_STAT_SAMPLE_MEAN+PS_STAT_MAX!?");
        return 13;
    }

    psFree(in);

    return 0;
}

static int testImageRoll(void)
{

    psImage* in;
    psImage* out;
    psImage* out2;
    int rows = 64;
    int cols = 64;
    int rows1 = 8;
    int cols1 = 8;

    /*
     The function psImageRoll shall generate a new psImage structure by 
     rolling the input image the correponding number of pixels in the vertical
     and/or horizontal direction. The image output image shall be the same size
     as the input image. Values which roll off the image are wrapped to the 
     other side.

     Verify the returned psImage structure contains expected values, if the 
     input image contains known values and the roll performed is known. 
     Cases should include no roll, vertical roll, horizontal roll and 
     combination vertical/horizontal rolls. Positive and negative rolls 
     should be performed.
    */

    in = psImageAlloc(cols,rows,PS_TYPE_F32);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        for (int col=0;col<cols;col++) {
            inRow[col] = (psF32)row+(psF32)col/1000.0f;
        }
    }

    out = psImageRoll(NULL,in,0,0);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=0, dy=0.",
                        col,row,inRow[col],outRow[col]);
                return 3;
            }
        }
    }

    out2 = psImageRoll(out,in,cols/4,0);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[(col+cols/4) % cols] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=cols/4, dy=0.",
                        col,row,inRow[(col+cols/4) % cols],outRow[col]);
                return 4;
            }
        }
    }

    // Verify the returned psImage structure pointer is equal to the input
    // parameter out if provided.
    if (out2 != out) {
        psError(__func__,"psImageRoll didn't recycle my out psImage!?");
        return 1;
    }

    out = psImageRoll(out,in,0,rows/4);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=0, dy=rows/4.",
                        col,row,inRow[col],outRow[col]);
                return 5;
            }
        }
    }

    out = psImageRoll(out,in,cols/4,rows/4);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[(col+cols/4) % cols] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=cols/4, dy=rows/4.",
                        col,row,inRow[(col+cols/4) % cols],outRow[col]);
                return 6;
            }
        }
    }

    out = psImageRoll(out,in,-cols/4,0);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[(col+(cols-cols/4)) % cols] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=-cols/4, dy=0.",
                        col,row,inRow[(col+(cols-cols/4)) % cols],outRow[col]);
                return 7;
            }
        }
    }

    out = psImageRoll(out,in,0,-rows/4);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows-rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=0, dy=-rows/4.",
                        col,row,inRow[col],outRow[col]);
                return 8;
            }
        }
    }

    out = psImageRoll(out,in,-cols/4,-rows/4);
    for (int row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows-rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (int col=0;col<cols;col++) {
            if (inRow[(col+cols-cols/4) % cols] != outRow[col]) {
                psError(__func__,"psImageRoll didn't produce expected result "
                        "at %d,%d (%f vs %f) for dx=cols/4, dy=rows/4.",
                        col,row,inRow[(col+cols-cols/4) % cols],outRow[col]);
                return 9;
            }
        }
    }


    // Verify the returned psImage structure pointer is null and program
    // execution doesn't stop, if input parameter input is null.
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error.");
    out2 = psImageRoll(NULL,NULL,0,0);
    if (out2 != NULL) {
        psError(__func__,"psImageRoll did not return NULL though input image was NULL!?");
        return 2;
    }

    psFree(in);
    psFree(out);

    #define testRollType(DATATYPE) \
    in = psImageAlloc(rows1,cols1,PS_TYPE_##DATATYPE); \
    \
    for (int row=0;row<rows1;row++) { \
        ps##DATATYPE* inRow = in->data.DATATYPE[row]; \
        for (int col=0;col<cols1;col++) { \
            inRow[col] = (ps##DATATYPE)row+(ps##DATATYPE)col; \
        } \
    } \
    \
    out = psImageRoll(NULL,in,rows1/4,cols1/4); \
    for (int row=0;row<rows1;row++) { \
        ps##DATATYPE* inRow = in->data.DATATYPE[(row+rows1/4)%rows1]; \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        for (int col=0;col<cols1;col++) { \
            if (inRow[(col+cols1/4)%cols1] != outRow[col]) { \
                psError(__func__,"psImageRoll didn't produce expected result " \
                        "at %d,%d (%f vs %f) for dx=0, dy=0.", \
                        col,row,(float)inRow[col],(float)outRow[col]); \
                return 3; \
            } \
        } \
    } \
    psFree(in); \
    psFree(out);

    testRollType(U8);
    testRollType(U16);
    testRollType(S8);
    testRollType(S16);
    testRollType(F64);
    testRollType(C32);
    testRollType(C64);

    return 0;
}

int testImageRotate(void)
{
    /*

    This function shall calculate a new psImage structure based upon the 
    rotation of a given psImage structure. The center of rotation shall be the 
    center pixel of the input image.

    The following steps of the testpoint are done manually via inspection of 
    temp/fOut.fits & temp/sOut.fits.

        * Verify the returned psImage structure contains expected values, if 
          the input parameter psImage contains known values. Cases should 
          include rotations of 0, 45, 90, 135, 180, 225, 270, 315, 360 and at leat one 
          other arbitrary angle. Cases of the input image should include image
          with a center pixel and an image without a center pixel.
        * Verify the returned psImage structure contains pixels set to exposed value, if
          the rotation and input psImage to not correspond to the output image.

    */

    psImage* fOut = NULL;
    psImage* sOut = NULL;
    psImage* s32Out = NULL;
    psImage* fBiOut = NULL;
    psImage* sBiOut = NULL;
    psImage* fTruth = NULL;
    psImage* sTruth = NULL;
    psImage* fBiTruth = NULL;
    psImage* sBiTruth = NULL;
    int rows = 64;
    int cols = 64;
    psImage* fImg = psImageAlloc(cols,rows,PS_TYPE_F32);
    psImage* sImg = psImageAlloc(cols,rows,PS_TYPE_S16);
    psImage* s32Img = psImageAlloc(cols,rows,PS_TYPE_S32);

    for(int row=0;row<rows;row++) {
        psF32* fRow = fImg->data.F32[row];
        psS16* sRow = sImg->data.S16[row];
        psS32* s32Row = s32Img->data.S32[row];
        for (int col=0;col<cols;col++) {
            fRow[col] = (psF32)(row)+(psF32)(col)/100.0f;
            sRow[col] = row-2*col;
            s32Row[col] = row-2*col;
        }
    }

    // since interpolation is involved, etc., the simplist way to verify things
    // is to verify the results manually and bless it for automated comparison
    // thereafter


    // write results of various rotates to a file and verify with truth images
    mkdir("temp",0777);
    remove
        ("temp/fOut.fits")
        ;
    remove
        ("temp/sOut.fits")
        ;
    remove
        ("temp/fBiOut.fits")
        ;
    remove
        ("temp/sBiOut.fits");
    ;
    int index = 0;
    bool fail = false;
    for (int rot=-180;rot<=180;rot+=45) {
        psImage* oldOut = fOut;
        psImage* oldBiOut = fBiOut;
        fOut = psImageRotate(fOut,fImg,rot,-1.0,PS_INTERPOLATE_FLAT);
        fBiOut = psImageRotate(fBiOut,fImg,rot,-1.0,PS_INTERPOLATE_BILINEAR);
        // Verify the returned psImage structure is equal to the input
        // parameter out if provided.
        if (oldOut != NULL && oldOut != fOut) {
            psError(__func__,"the output recycle functionality failed");
            return 2;
        }
        if (oldBiOut != NULL && oldBiOut != fBiOut) {
            psError(__func__,"the output recycle functionality failed");
            return 4;
        }
        sOut = psImageRotate(sOut,sImg,rot,-1.0,PS_INTERPOLATE_FLAT);
        sBiOut = psImageRotate(sBiOut,sImg,rot,-1.0,PS_INTERPOLATE_BILINEAR);
        if (! psImageWriteSection(fOut,0,0,0,NULL,index,"temp/fOut.fits") ) {
            psError(__func__,"Can not write to temp/fOut.fits, so why continue!?");
            return 20;
        }
        if (! psImageWriteSection(sOut,0,0,0,NULL,index,"temp/sOut.fits") ) {
            psError(__func__,"Can not write to temp/sOut.fits, so why continue!?");
            return 21;
        }
        if (! psImageWriteSection(fBiOut,0,0,0,NULL,index,"temp/fBiOut.fits") ) {
            psError(__func__,"Can not write to temp/fBiOut.fits, so why continue!?");
            return 40;
        }
        if (! psImageWriteSection(sBiOut,0,0,0,NULL,index,"temp/sBiOut.fits") ) {
            psError(__func__,"Can not write to temp/sBiOut.fits, so why continue!?");
            return 41;
        }

        // now, let's compare this with the verified file
        fTruth = psImageReadSection(fTruth,0,0,0,0,0,NULL,index,"verified/fOut.fits");
        sTruth = psImageReadSection(sTruth,0,0,0,0,0,NULL,index,"verified/sOut.fits");
        fBiTruth = psImageReadSection(fBiTruth,0,0,0,0,0,NULL,index,"verified/fBiOut.fits");
        sBiTruth = psImageReadSection(sBiTruth,0,0,0,0,0,NULL,index,"verified/sBiOut.fits");
        if (fTruth == NULL) {
            psError(__func__,"verified psF32 image failed to be read (%d deg. rotation)",
                    rot);
            fail = true;
        } else {
            if (fTruth->numRows != fOut->numRows || fTruth->numCols != fOut->numCols) {
                psError(__func__,"Rotated float image size did not match truth "
                        "image for %d deg rotation (%dx%d vs %dx%d).",
                        rot,fOut->numCols,fOut->numRows,fTruth->numCols,fTruth->numRows);
                fail = true;
            } else {
                for (int row=0;row<fTruth->numRows;row++) {
                    psF32* truthRow = fTruth->data.F32[row];
                    psF32* outRow = fOut->data.F32[row];
                    for (int col=0;col<fTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(__func__,"Float Image mismatch (%f vs %f) at %d,%d.",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

        if (sTruth == NULL) {
            psError(__func__,"verified psS16 image failed to be read "
                    "(%d deg. rotation)",rot);
            fail = true;
        } else {
            if (sTruth->numRows != sOut->numRows ||
                    sTruth->numCols != sOut->numCols) {
                psError(__func__,"Rotated short image size did not match truth "
                        "image for %d deg rotation.",rot);
                fail = true;
            } else {
                for (int row=0;row<sTruth->numRows;row++) {
                    psS16* truthRow = sTruth->data.S16[row];
                    psS16* outRow = sOut->data.S16[row];
                    for (int col=0;col<sTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(__func__,"Short Image mismatch (%d vs %d) "
                                    "at %d,%d.",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }


        if (fBiTruth == NULL) {
            psError(__func__,"verified psF32 Bi image failed to be read (%d deg. rotation)",
                    rot);
            fail = true;
        } else {
            if (fBiTruth->numRows != fBiOut->numRows || fBiTruth->numCols != fBiOut->numCols) {
                psError(__func__,"Rotated float image size did not match truth "
                        "image for %d deg rotation (%dx%d vs %dx%d). BILINEAR",
                        rot,fBiOut->numCols,fBiOut->numRows,fBiTruth->numCols,fBiTruth->numRows);
                fail = true;
            } else {
                for (int row=0;row<fBiTruth->numRows;row++) {
                    psF32* truthRow = fBiTruth->data.F32[row];
                    psF32* outRow = fBiOut->data.F32[row];
                    for (int col=0;col<fBiTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(__func__,"Float Image mismatch (%f vs %f) at %d,%d. BILINEAR",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

        if (sBiTruth == NULL) {
            psError(__func__,"verified psS16 image failed to be read "
                    "(%d deg. rotation) BILINEAR",rot);
            fail = true;
        } else {
            if (sBiTruth->numRows != sBiOut->numRows ||
                    sBiTruth->numCols != sBiOut->numCols) {
                psError(__func__,"Rotated short image size did not match truth "
                        "image for %d deg rotation. BILINEAR",rot);
                fail = true;
            } else {
                for (int row=0;row<sBiTruth->numRows;row++) {
                    psS16* truthRow = sBiTruth->data.S16[row];
                    psS16* outRow = sBiOut->data.S16[row];
                    for (int col=0;col<sBiTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(__func__,"Short Image mismatch (%d vs %d) "
                                    "at %d,%d. BILINEAR",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

        // Verify psImage structure returned is NULL and program execution doesn't stop
        // if the input image type is not supported. Do not check for 0 degree rotation
        // since an imgae copy will be used.
        if (rot != 0 ) {
            psLogMsg(__func__,PS_LOG_INFO,"Following should be an error invalid type for rotate.");
            s32Out = psImageRotate(s32Out,s32Img,rot,-1.0,PS_INTERPOLATE_FLAT);
            if( s32Out != NULL ) {
                psError(__func__,"psImage returned is not NULL for invalid image type.");
                return 5;
            }
        }

        index++;
    }

    if (fail) {
        psError(__func__,"One or more images didn't match truth or truth did "
                "not exist.");
        return 10;
    }


    // Verify the returned psImage structure pointer is null and program
    // execution doesn't stop, if the input parameter input is null.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error");
    fOut = psImageRotate(fOut,NULL,0,0,PS_INTERPOLATE_FLAT);
    if (fOut != NULL) {
        psError(__func__,"NULL wasn't returned though the input image was NULL.");
        return 3;
    }

    // Verify the returned psImage structure pointer is null and program
    // execution doesn't stop, if the specified interpolation mode is invalid
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for invalid "
             "interpolation type.");
    fOut = psImageRotate(fOut, fImg, 33, 0, -1);
    if (fOut != NULL) {
        psError(__func__,"NULL wasn't returned though the interpolation mode "
                "is invalid.");
        return 4;
    }

    psFree(sOut);
    psFree(fImg);
    psFree(sImg);
    psFree(fTruth);
    psFree(sTruth);
    psFree(sBiOut);
    psFree(fBiTruth);
    psFree(sBiTruth);
    psFree(fBiOut);
    psFree(s32Img);
    psFree(s32Out);

    return 0;
}

static int testImageShift(void)
{
    /* psImageShift:

       This functions shall generate a new psImage structure by shifting the 
       input psImage structure a specified number of pixels in the horizontal 
       and/or vertical directions.

       Verify the returned psImage structure contains expected values, if the 
       input psImage structure contains known values and a know shift in the 
       vertical and/or horizontal directions. Cases should include no shift, 
       vertical only(up,down), horizontal only(right, left), and combination 
       shift. Cases should include fractional shifts. Comparison of expected 
       values should include a delta to allow for testing on different 
       platforms.

       Verify the returned psImage structure contains values for pixels not in 
       the original image set to the input parameter exposed.

    */

    int retVal=0;

    // integer shift
    retVal |= testImageShiftCase(64,128,0.0f,0.0f);
    retVal |= testImageShiftCase(64,128,0.0f,16.0f);
    retVal |= testImageShiftCase(64,128,0.0f,-16.0f);
    retVal |= testImageShiftCase(64,128,32.0f,0.0f);
    retVal |= testImageShiftCase(64,128,-32.0f,0.0f);
    retVal |= testImageShiftCase(64,128,32.0f,16.0f);
    retVal |= testImageShiftCase(64,128,32.0f,-16.0f);
    retVal |= testImageShiftCase(64,128,-32.0f,16.0f);
    retVal |= testImageShiftCase(64,128,-32.0f,-16.0f);

    if (retVal != 0) {
        return retVal;
    }

    // fractional shift
    retVal |= testImageShiftCase(64,128,0.0f,16.4f);
    retVal |= testImageShiftCase(64,128,0.0f,-16.4f);
    retVal |= testImageShiftCase(64,128,32.7f,0.0f);
    retVal |= testImageShiftCase(64,128,-32.7f,0.0f);
    retVal |= testImageShiftCase(64,128,32.6f,16.2f);
    retVal |= testImageShiftCase(64,128,32.6f,-16.2f);
    retVal |= testImageShiftCase(64,128,-32.6f,16.2f);
    retVal |= testImageShiftCase(64,128,-32.6f,-16.2f);

    if (retVal != 0) {
        return retVal;
    }

    /*
       Verify the returned psImage structure pointer is equal to the input 
       parameter out if provided.
    */
    psImage* fImg = psImageAlloc(32,32,PS_TYPE_F32);
    psImage* fRecycle = psImageAlloc(32,32,PS_TYPE_F32);
    psImage* fOut = psImageShift(fRecycle, fImg, 8,8, NAN, PS_INTERPOLATE_FLAT);

    if (fRecycle != fOut) {
        psError(__func__,"psImageShift didn't recycle my image?");
        return 10;
    }

    /*
       Verify the returned psImage structure pointer is null and program 
       execution doesn't stop, if the input psImage structure pointer is null.
    */
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error...");
    fOut = psImageShift(fOut,NULL,8,8,NAN,PS_INTERPOLATE_FLAT);
    if (fOut != NULL) {
        psError(__func__,"psImageShift didn't return NULL given a NULL input image.");
        return 11;
    }

    // Verify the returned psImage structure is null and program execution
    // doesn't stop, if the specified interpolation mode is invalid.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error for invalid interpolation mode.");
    fOut = psImageShift(fOut,fImg,8,8,NAN,-1);
    if (fOut != NULL ) {
        psError(__func__,"psImageShift didn't return NULL given an invalid interpolation mode.");
        return 12;
    }

    psFree(fImg);

    return 0;
}

static int testImageShiftCase(int cols,
                              int rows,
                              float colShift,
                              float rowShift)
{
    psImage* fOut = NULL;
    psImage* sOut = NULL;

    psLogMsg(__func__,PS_LOG_INFO,"Testing psImageShift with a %dx%d image for "
             "a shift of %g,%g.",cols,rows,colShift,rowShift);

    psImage* fImg = psImageAlloc(cols,rows,PS_TYPE_F32);
    psImage* sImg = psImageAlloc(cols,rows,PS_TYPE_S16);
    psImage* fBiOut = psImageAlloc(cols,rows,PS_TYPE_F32);
    psImage* sBiOut = psImageAlloc(cols,rows,PS_TYPE_S16);

    for(int row=0;row<rows;row++) {
        psF32* fRow = fImg->data.F32[row];
        psS16* sRow = sImg->data.S16[row];
        for (int col=0;col<cols;col++) {
            fRow[col] = (psF32)(row)+(psF32)(col)/100.0f;
            sRow[col] = row-2*col;
        }
    }

    fOut = psImageShift(fOut, fImg, colShift, rowShift, NAN, PS_INTERPOLATE_FLAT);
    sOut = psImageShift(sOut, sImg, colShift, rowShift, -1, PS_INTERPOLATE_FLAT);
    fBiOut = psImageShift(fBiOut, fImg, colShift, rowShift, NAN, PS_INTERPOLATE_BILINEAR);
    sBiOut = psImageShift(sBiOut, sImg, colShift, rowShift, -1, PS_INTERPOLATE_BILINEAR);

    for(int row=0;row<rows;row++) {
        psF32* fRow = fOut->data.F32[row];
        psS16* sRow = sOut->data.S16[row];
        psF32* fBiRow = fBiOut->data.F32[row];
        psS16* sBiRow = sBiOut->data.S16[row];

        for (int col=0;col<cols;col++) {
            psF32 fValue = psImagePixelInterpolate(fImg,col+colShift,
                                                   row+rowShift,NULL,0,NAN,PS_INTERPOLATE_FLAT);
            psS16 sValue = (psS16)psImagePixelInterpolate(sImg,col+colShift,
                           row+rowShift,NULL,0,-1,PS_INTERPOLATE_FLAT);

            psF32 fBiValue = psImagePixelInterpolate(fImg,col+colShift,
                             row+rowShift,NULL,0,NAN,PS_INTERPOLATE_BILINEAR);
            psS16 sBiValue = (psS16)psImagePixelInterpolate(sImg,col+colShift,
                             row+rowShift,NULL,0,-1,PS_INTERPOLATE_BILINEAR);

            if (fabsf(fRow[col] - fValue) > FLT_EPSILON) {
                psError(__func__,"Float image not shifted correctly at %d,%d (%g vs %g)",
                        col,row,fRow[col],fValue);
                return 1;
            }
            if (sRow[col] != sValue) {
                psError(__func__,"Short image not shifted correctly at %d,%d (%d vs %d)",
                        col,row,sRow[col],sValue);
                return 2;
            }
            if (fabsf(fBiRow[col] - fBiValue) > FLT_EPSILON) {
                psError(__func__,"Float image not shifted correctly at %d,%d (%g vs %g)",
                        col,row,fBiRow[col],fBiValue);
                return 1;
            }
            if (sBiRow[col] != sBiValue) {
                psError(__func__,"Short image not shifted correctly at %d,%d (%d vs %d)",
                        col,row,sBiRow[col],sBiValue);
                return 2;
            }
        }
    }

    psFree(fImg);
    psFree(sImg);
    psFree(fOut);
    psFree(sOut);
    psFree(fBiOut);
    psFree(sBiOut);

    return 0;
}

static int testImageResample(void)
{

    int rows = 60;
    int cols = 80;
    psImage* result = NULL;
    int scale = 4;
    psErr* err;

    psImage* image = psImageAlloc(cols,rows,PS_TYPE_F32);
    for(int row=0;row<rows;row++) {
        psF32* imageRow = image->data.F32[row];
        for (int col=0;col<cols;col++) {
            imageRow[col] = row+2*col;
        }
    }
    result = psImageCopy(NULL,image,PS_TYPE_F64);
    psImage* orig = result;
    result = psImageResample(result,image,scale,PS_INTERPOLATE_FLAT);

    if (result == NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "NULL return unexpected");
        return 1;
    }

    if (result != orig) {
        psLogMsg(__func__,PS_LOG_ERROR,"failure to recycle image.");
        return 2;
    }

    if (result->type.type != PS_TYPE_F32) {
        psLogMsg(__func__,PS_LOG_ERROR,"unexpected type");
        return 3;
    }


    if (result->numCols != image->numCols*scale ||
            result->numRows != image->numRows*scale) {
        psLogMsg(__func__,PS_LOG_ERROR,"The size of the result is %dx%d, but %dx%d was expected.",
                 result->numCols,result->numRows,
                 image->numCols*scale, image->numRows*scale);
        return 4;
    }

    psF32 truthValue;
    for(int row=0;row<result->numRows;row++) {
        for (int col=0;col<result->numCols;col++) {
            truthValue = psImagePixelInterpolate(image,
                                                 (float)col/(float)scale,(float)row/(float)scale,
                                                 NULL,0,-1,PS_INTERPOLATE_FLAT);
            if (fabs(truthValue - result->data.F32[row][col]) > FLT_EPSILON) {
                psLogMsg(__func__,PS_LOG_ERROR,"value bad at (%d,%d).  Got %g, expected %g.",
                         col,row,result->data.F32[row][col], truthValue);
                return 5;
            }
        }
    }

    // verify that image=null is handled properly.
    psErrorClear();
    result = psImageResample(result,NULL,scale,PS_INTERPOLATE_FLAT);

    if (result != NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "return was not NULL, as expected.");
        return 6;
    }
    err = psErrorLast();
    if (err->code != PS_ERR_BAD_PARAMETER_NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "error message was not appropriate type.");
        return 7;
    }
    psFree(err);

    // verify that scale < 1 is handled properly
    psErrorClear();
    result = psImageResample(result,image,0,PS_INTERPOLATE_FLAT);

    if (result != NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "return was not NULL, as expected.");
        return 8;
    }
    err = psErrorLast();
    if (err->code != PS_ERR_BAD_PARAMETER_VALUE) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "error message was not appropriate type.");
        return 9;
    }
    psFree(err);

    // verify that invalid interpolation mode is handled properly
    psErrorClear();
    result = psImageResample(result,image,2,0);

    if (result != NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "return was not NULL, as expected.");
        return 10;
    }
    err = psErrorLast();
    if (err->code != PS_ERR_BAD_PARAMETER_VALUE) {
        psLogMsg(__func__,PS_LOG_ERROR,
                 "error message was not appropriate type.");
        return 11;
    }
    psFree(err);

    psFree(image);
    psFree(result);

    return 0;
}

