/** @file  tst_psImageManip.c
 *
 *  @brief Contains the tests for psImageManip.[ch]
 *
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-04-07 20:27:42 $
 *
 *  Copyright 2004-2005 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_strict.h"
#include "psType.h"

static psS32 testImageClip(void);
static psS32 testImageClipNAN(void);
static psS32 testImageClipComplexRegion(void);
static psS32 testImageOverlay(void);
static psS32 testImageRebin(void);
static psS32 testImageRoll(void);
static psS32 testImageRotate(void);
static psS32 testImageShift(void);
static psS32 testImageShiftCase(psS32 cols, psS32 rows, float colShift,float rowShift);
static psS32 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}
                          };

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

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


psS32 testImageClip(void)
{
    psImage* img = NULL;
    psU32 c = 128;
    psU32 r = 256;
    psF64 min;
    psF64 max;
    psS32 numClipped = 0;
    psS32 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 (psU32 row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (psU32 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 (psU32 row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (psU32 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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Expected %d clips, but got %d", \
                numClipped,retVal); \
        return 2; \
    } \
    psFree(img);

    #define testImageClipByComplexType(datatype) \
    img = psImageAlloc(c,r,PS_TYPE_##datatype); \
    for (psU32 row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (psU32 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 (psU32 row=0;row<r;row++) { \
        ps##datatype* imgRow = img->data.datatype[row]; \
        for (psU32 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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Expected zero return for clips when max < min.");
        return 4;
    }

    return 0;
}

psS32 testImageClipNAN(void)
{
    psImage* img = NULL;
    psU32 c = 128;
    psU32 r = 256;
    psS32 numClipped = 0;
    psS32 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(PS_ERR_UNKNOWN, true,"Pixel value is not as expected (%f vs %f) at %d,%d", \
                        imgRow[col],value,col,row); \
                return 1; \
            } \
        } \
    } \
    if (retVal != numClipped) { \
        psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Expected zero return for clip of invalid image type.");
        return 4;
    }
    psFree(img);

    return 0;
}

psS32 testImageClipComplexRegion(void)
{
    psImage* img = NULL;
    psU32 c = 1024;
    psU32 r = 2048;
    psS32 numClipped = 0;
    psS32 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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Expected zero return for clip of invalid image type.");
        return 84;
    }
    psFree(img);

    return 0;
}

psS32 testImageOverlay(void)
{

    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    psImage* img4 = NULL;
    psU32 c = 128;
    psU32 r = 256;
    psS32 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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned zero with %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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageOverlaySection did not return "
                "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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageOverlaySection did not return "
                "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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned non-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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned non-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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned non-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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned nonzero eventhough "
                " overlay image type is different than input image.");
        return 10;
    }

    /*
    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(PS_ERR_UNKNOWN, true,"psImageOverlaySection returned zero when "
                "checking divide-by-zero.");
        return 12;
    }

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

    return 0;
}

static psS32 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 (psS32 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 (psS32 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 (psS32 row = 0; row<4; row++) { \
        psF32* meanTruthRow = meanTruth->data.F32[row]; \
        psF32* meanTruthWMaskRow = meanTruthWMask->data.F32[row]; \
        for (psS32 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(PS_ERR_UNKNOWN, true,"psImageRebin returned a NULL pointer!?"); \
        return 1; \
    } \
    if (out->numRows != 4 || out->numCols != 4) { \
        psError(PS_ERR_UNKNOWN, true,"psImageRebin didn't produce the proper size image " \
                "(%d x %d).", \
                out->numCols, out->numRows); \
        return 2; \
    } \
    for (psS32 row = 0; row<4; row++) { \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        psF32* truthRow = meanTruth->data.F32[row]; \
        for (psS32 col = 0; col<4; col++) { \
            if (fabsf((float)outRow[col]-(float)truthRow[col]) > FLT_EPSILON) { \
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row = 0; row<4; row++) { \
        ps##DATATYPE* outRow = out3->data.DATATYPE[row]; \
        psF32* truthRow = meanTruthWMask->data.F32[row]; \
        for ( psS32 col = 0; col<4; col++) { \
            if(abs((psS32)outRow[col]-(psS32)truthRow[col]) > FLT_EPSILON) { \
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageRebin didn't recycle a psImage properly!?"); \
        return 7; \
    } \
    if (out == NULL) { \
        psError(PS_ERR_UNKNOWN, true,"psImageRebin returned a NULL pointer!?"); \
        return 4; \
    } \
    if (out->numRows != 6 || out->numCols != 6) { \
        psError(PS_ERR_UNKNOWN, true,"psImageRebin didn't produce the proper size image " \
                "(%d x %d).", \
                out->numCols, out->numRows); \
        return 5; \
    } \
    for (psS32 row = 0; row<6; row++) { \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        psF32* truthRow = maxTruth->data.F32[row]; \
        for (psS32 col = 0; col<6; col++) { \
            if (fabsf((float)outRow[col]-(float)truthRow[col]) > FLT_EPSILON) { \
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageRebin returned an image though the stats "
                "options was PS_STAT_SAMPLE_MEAN+PS_STAT_MAX!?");
        return 13;
    }

    psFree(in);

    return 0;
}

static psS32 testImageRoll(void)
{

    psImage* in;
    psImage* out;
    psImage* out2;
    psS32 rows = 64;
    psS32 cols = 64;
    psS32 rows1 = 8;
    psS32 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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            inRow[col] = (psF32)row+(psF32)col/1000.0f;
        }
    }

    out = psImageRoll(NULL,in,0,0);
    for (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[(col+cols/4) % cols] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageRoll didn't recycle my out psImage!?");
        return 1;
    }

    out = psImageRoll(out,in,0,rows/4);
    for (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[(col+cols/4) % cols] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[row];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[(col+(cols-cols/4)) % cols] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows-rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[col] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows;row++) {
        psF32* inRow = in->data.F32[(row+rows-rows/4)%rows];
        psF32* outRow = out->data.F32[row];
        for (psS32 col=0;col<cols;col++) {
            if (inRow[(col+cols-cols/4) % cols] != outRow[col]) {
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<rows1;row++) { \
        ps##DATATYPE* inRow = in->data.DATATYPE[row]; \
        for (psS32 col=0;col<cols1;col++) { \
            inRow[col] = (ps##DATATYPE)row+(ps##DATATYPE)col; \
        } \
    } \
    \
    out = psImageRoll(NULL,in,rows1/4,cols1/4); \
    for (psS32 row=0;row<rows1;row++) { \
        ps##DATATYPE* inRow = in->data.DATATYPE[(row+rows1/4)%rows1]; \
        ps##DATATYPE* outRow = out->data.DATATYPE[row]; \
        for (psS32 col=0;col<cols1;col++) { \
            if (inRow[(col+cols1/4)%cols1] != outRow[col]) { \
                psError(PS_ERR_UNKNOWN, true,"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;
}

psS32 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* fBiOut = NULL;
    psImage* sBiOut = NULL;
    psImage* fTruth = NULL;
    psImage* sTruth = NULL;
    psImage* fBiTruth = NULL;
    psImage* sBiTruth = NULL;
    psS32 rows = 64;
    psS32 cols = 64;
    psImage* fImg = psImageAlloc(cols,rows,PS_TYPE_F32);
    psImage* sImg = psImageAlloc(cols,rows,PS_TYPE_S16);

    for(psS32 row=0;row<rows;row++) {
        psF32* fRow = fImg->data.F32[row];
        psS16* sRow = sImg->data.S16[row];
        for (psS32 col=0;col<cols;col++) {
            fRow[col] = (psF32)(row)+(psF32)(col)/100.0f;
            sRow[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");
    ;
    psS32 index = 0;
    psBool fail = false;
    psF32 radianRot;
    for (psS32 rot=-180;rot<=180;rot+=45) {
        psImage* oldOut = fOut;
        psImage* oldBiOut = fBiOut;

        if (rot == 90) {
            radianRot = PS_PI_2;
        } else if (rot == -90) {
            radianRot = PS_PI+PS_PI_2;
        } else if (rot == 180 || rot == -180) {
            radianRot = PS_PI;
        } else {
            radianRot = ((float)rot)*PS_PI/180.0;
        }

        fOut = psImageRotate(fOut,fImg,radianRot,-1.0,PS_INTERPOLATE_FLAT);
        fBiOut = psImageRotate(fBiOut,fImg,radianRot,-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(PS_ERR_UNKNOWN, true,"the output recycle functionality failed");
            return 2;
        }
        if (oldBiOut != NULL && oldBiOut != fBiOut) {
            psError(PS_ERR_UNKNOWN, true,"the output recycle functionality failed");
            return 4;
        }
        sOut = psImageRotate(sOut,sImg,radianRot,-1.0,PS_INTERPOLATE_FLAT);
        sBiOut = psImageRotate(sBiOut,sImg,radianRot,-1.0,PS_INTERPOLATE_BILINEAR);
        if (! psImageWriteSection(fOut,0,0,0,NULL,index,"temp/fOut.fits") ) {
            psError(PS_ERR_UNKNOWN, true,"Can not write to temp/fOut.fits, so why continue!?");
            return 20;
        }
        if (! psImageWriteSection(sOut,0,0,0,NULL,index,"temp/sOut.fits") ) {
            psError(PS_ERR_UNKNOWN, true,"Can not write to temp/sOut.fits, so why continue!?");
            return 21;
        }
        if (! psImageWriteSection(fBiOut,0,0,0,NULL,index,"temp/fBiOut.fits") ) {
            psError(PS_ERR_UNKNOWN, true,"Can not write to temp/fBiOut.fits, so why continue!?");
            return 40;
        }
        if (! psImageWriteSection(sBiOut,0,0,0,NULL,index,"temp/sBiOut.fits") ) {
            psError(PS_ERR_UNKNOWN, true,"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,"temp/fOut.fits");
        sTruth = psImageReadSection(sTruth,0,0,0,0,0,NULL,index,"temp/sOut.fits");
        fBiTruth = psImageReadSection(fBiTruth,0,0,0,0,0,NULL,index,"temp/fBiOut.fits");
        sBiTruth = psImageReadSection(sBiTruth,0,0,0,0,0,NULL,index,"temp/sBiOut.fits");
        if (fTruth == NULL) {
            psError(PS_ERR_UNKNOWN, true,"verified psF32 image failed to be read (%d deg. rotation)",
                    rot);
            fail = true;
        } else {
            if (fTruth->numRows != fOut->numRows || fTruth->numCols != fOut->numCols) {
                psError(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<fTruth->numRows;row++) {
                    psF32* truthRow = fTruth->data.F32[row];
                    psF32* outRow = fOut->data.F32[row];
                    for (psS32 col=0;col<fTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(PS_ERR_UNKNOWN, true,"Float Image mismatch (%f vs %f) at %d,%d.",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

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


        if (fBiTruth == NULL) {
            psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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 (psS32 row=0;row<fBiTruth->numRows;row++) {
                    psF32* truthRow = fBiTruth->data.F32[row];
                    psF32* outRow = fBiOut->data.F32[row];
                    for (psS32 col=0;col<fBiTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(PS_ERR_UNKNOWN, true,"Float Image mismatch (%f vs %f) at %d,%d. BILINEAR",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

        if (sBiTruth == NULL) {
            psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Rotated psS16 image size did not match truth "
                        "image for %d deg rotation. BILINEAR",rot);
                fail = true;
            } else {
                for (psS32 row=0;row<sBiTruth->numRows;row++) {
                    psS16* truthRow = sBiTruth->data.S16[row];
                    psS16* outRow = sBiOut->data.S16[row];
                    for (psS32 col=0;col<sBiTruth->numCols;col++) {
                        if (fabsf(truthRow[col]-outRow[col]) > 1) {
                            psError(PS_ERR_UNKNOWN, true,"Short Image mismatch (%d vs %d) "
                                    "at %d,%d. BILINEAR",
                                    outRow[col], truthRow[col],col,row);
                            fail = true;
                        }
                    }
                }
            }
        }

        index++;
    }

    if (fail) {
        psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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);

    return 0;
}

static psS32 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.

    */

    psS32 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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageShift didn't return NULL given an invalid interpolation mode.");
        return 12;
    }

    psFree(fImg);

    return 0;
}

static psS32 testImageShiftCase(psS32 cols,
                                psS32 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(psS32 row=0;row<rows;row++) {
        psF32* fRow = fImg->data.F32[row];
        psS16* sRow = sImg->data.S16[row];
        for (psS32 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(psS32 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 (psS32 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(PS_ERR_UNKNOWN, true,"Float image not shifted correctly at %d,%d (%g vs %g)",
                        col,row,fRow[col],fValue);
                return 1;
            }
            if (sRow[col] != sValue) {
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"Float image not shifted correctly at %d,%d (%g vs %g)",
                        col,row,fBiRow[col],fBiValue);
                return 1;
            }
            if (sBiRow[col] != sBiValue) {
                psError(PS_ERR_UNKNOWN, true,"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 psS32 testImageResample(void)
{

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

    psImage* image = psImageAlloc(cols,rows,PS_TYPE_F32);
    for(psS32 row=0;row<rows;row++) {
        psF32* imageRow = image->data.F32[row];
        for (psS32 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(psS32 row=0;row<result->numRows;row++) {
        for (psS32 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,-1);

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

    // Verify that that an invalid image type is handled properly
    psErrorClear();
    psImage* invImage = psImageAlloc(cols,rows,PS_TYPE_BOOL);
    memset(invImage->rawDataBuffer,0,cols*rows*sizeof(psBool)); // make sure the image is of all NULLs
    result = psImageResample(result,invImage,2,PS_INTERPOLATE_FLAT);
    if (result != NULL) {
        psLogMsg(__func__,PS_LOG_ERROR,"return was not NULL, as expected.");
        return 20;
    }
    err = psErrorLast();
    if (err->code != PS_ERR_BAD_PARAMETER_TYPE) {
        psLogMsg(__func__,PS_LOG_ERROR,"error message was not appropriate type.");
        return 21;
    }
    psFree(err);

    psFree(image);
    psFree(result);
    psFree(invImage);

    return 0;
}

