/** @file  tst_psImageFFT.c
 *
 *  @brief Contains the tests for psFFT.[ch]
 *
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2007-01-30 21:32:56 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */
#include <stdio.h>
#include <string.h>
#include <pslib.h>
#include "tap.h"
#include "pstap.h"

#define GENIMAGE(img,c,r,TYP, valueFcn) \
img = psImageAlloc(c,r,PS_TYPE_##TYP); \
for (psU32 row=0;row<r;row++) { \
    ps##TYP* imgRow = img->data.TYP[row]; \
    for (psU32 col=0;col<c;col++) { \
        imgRow[col] = (ps##TYP)(valueFcn); \
    } \
}

psS32 main(psS32 argc, char* argv[])
{
    psLogSetFormat("HLNM");
    psLogSetLevel(PS_LOG_INFO);
    plan_tests(68);

    // psImageFFT(void)
    {
        psMemId id = psMemGetId();
        psImage* img = NULL;
        psImage* img2 = NULL;
        psImage* img3 = NULL;
        psU32 m = 128;
        psU32 n = 64;
        psImage* img4 = NULL;
        psImage* img5 = NULL;

        // 1. assign a image to a radial sinisoid
        // 2. perform forward transform
        // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
        // 4. perform reverse transform
        // 5. compare to original (should be equal to within a reasonable error)

        // 1. assign a image to a radial sinisoid
        GENIMAGE(img,m,n,F32, sinf((32.0f-row)/32.0f*M_PI)+sinf((64.0f-col)/64.0f*M_PI));

        // 2. perform forward transform
        img2 = psImageFFT(img2,img,PS_FFT_FORWARD);
        ok(img2 != NULL, "psImageFFT() returned non-NULL");
        skip_start(img2 == NULL, 1, "Skipping tests because psImageFFT() returned NULL");
        ok(img2->type.type == PS_TYPE_C32, "FFT produced complex values");
        ok(img2->numCols == m && img2->numRows == n, "FFT produced proper size result (%dx%d vs. expected %dx%d).",
           img2->numCols,img2->numRows,m,n);

        // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psC32* img2Row = img2->data.C32[row];
            for (psU32 col=0;col<m;col++) {
                psF32 mag = cabsf(img2Row[col])/m/n;
                if (mag > 0.1f) {
                    // must be (0,1) or (0,n-1) or (1,0) or (m-1,0)
                    if (! (col == 0 && (row == 1 || row == n-1))
                            && ! (row == 0 && (col==1 || col == m-1)) ) {
                        diag("Result incorrect at %d,%d (%.2f)",col,row,mag);
                        errorFlag = true;
                    }
                } else {
                    if ( (col == 0 && (row == 1 || row == n-1))
                            || (row == 0 && (col==1 || col == m-1)) ) {
                        diag("Result incorrect at %d,%d (%.2f)",col,row,mag);
                        errorFlag = true;
                    }
                }
            }
        }
        ok(!errorFlag, "FFT produced correct data values");



        // 4. perform reverse transform
        img3 = psImageFFT(img3,img2,PS_FFT_REVERSE);
        ok(img3 != NULL, "psImageFFT() returned non-NULL");
        skip_start(img3 == NULL, 1, "Skipping tests because psImageFFT() returned NULL");
        ok(img3->type.type == PS_TYPE_C32, "FFT produced complex values");

        ok(img3->numCols == m && img3->numRows == n, "FFT didt produce proper size result (%dx%d vs. expected %dx%d).",
           img3->numCols,img3->numRows,m,n);
















        errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psC32* img3Row = img3->data.C32[row];
            psF32* imgRow = img->data.F32[row];
            for (psU32 col=0;col<m;col++) {
                psF32 pixel = creal(img3Row[col])/m/n;
                if (fabsf(pixel-imgRow[col]) > 0.1) {
                    diag("Reverse FFT gave original image back (%d,%d %.2f vs %.2f)",
                         col,row,pixel,imgRow[col]);
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "FFT produced correct data values");
        ok(!errorFlag, "FFT produced correct data values");



        // 4. perform reverse transform to real result
        img3 = psImageFFT(img3,img2,PS_FFT_REVERSE|PS_FFT_REAL_RESULT);
        ok(img3 != NULL, "psImageFFT() returned non-NULL");
        skip_start(img3 == NULL, 1, "Skipping tests because psImageFFT() returned NULL");
        ok(img3->type.type == PS_TYPE_F32, "FFT asked to make real result");
        ok(img3->numCols == m && img3->numRows == n, "FFT produced proper size result (%dx%d vs. expected %dx%d).",
           img3->numCols,img3->numRows,m,n);

        errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psF32* img3Row = img3->data.F32[row];
            psF32* imgRow = img->data.F32[row];
            for (psU32 col=0;col<m;col++) {
                psF32 pixel = img3Row[col]/m/n;
                if (fabsf(pixel-imgRow[col]) > 0.1) {
                    diag("Reverse FFT gave original image back (%d,%d %.2f vs %.2f)",
                         col,row,pixel,imgRow[col]);
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "FFT produced correct data values");
        ok(!errorFlag, "FFT produced correct data values");



        // check if error occurs if FORWARD and REVERSE are both given.
        // Following should be an error
        // XXX: werify error
        ok(psImageFFT(NULL,img2,PS_FFT_REVERSE|PS_FFT_FORWARD) == NULL,
           "PS_FFT_REVERSE|PS_FFT_FORWARD returned NULL");

        // Following should be an error
        // XXX: werify error
        ok(psImageFFT(NULL,img2,PS_FFT_FORWARD|PS_FFT_REAL_RESULT) == NULL,
           "PS_FFT_FORWARD|PS_FFT_REAL_RESULT returned NULL");

        // Verify return null and program execution doesn't stop if input image is null
        img4 = psImageFFT(NULL,NULL,PS_FFT_FORWARD);
        ok(img4 == NULL, "psImageFFT should return null for a null input image");

        // Verify return null and program execution doesn't stop if input image is incorrect direction
        // Following should generate error for incorrect direction
        // XXX: werify error
        GENIMAGE(img4,8,8,S8,row+col);
        img5 = psImageFFT(NULL,img4,PS_FFT_REAL_RESULT);
        ok(img5 == NULL, "psImageFFT should return null for an incorrect FFT direction");

        psFree(img4);
        psFree(img);
        psFree(img2);
        psFree(img3);
        skip_end();
        skip_end();
        skip_end();
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // 1. create a C64 complex vector with distinctly different real and imaginary parts.
    // 2. call psImageReal and psImageImaginary
    // 3. compare results to the real/imaginary components of input
    {
        psMemId id = psMemGetId();
        psImage* c64Img = NULL;
        psImage* c64Img2 = NULL;
        psImage* c64Img3 = NULL;
        psU32 m = 128;
        psU32 n = 64;

        // XXX: What is I?
        GENIMAGE(c64Img,m,n,C64, row + I * col);
        c64Img2 = psImageReal(c64Img2,c64Img);
        ok(c64Img2 != NULL, "psImageReal() returned non-NULL");
        ok(c64Img2->type.type == PS_TYPE_F64, "psImageReal() returned the correct type");
        c64Img3 = psImageImaginary(c64Img3,c64Img);
        ok(c64Img3 != NULL, "psImageImaginary() returned non-NULL");
        ok(c64Img3->type.type == PS_TYPE_F64, "psImageImaginary() returned the correct type");

        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psF64* img2Row = c64Img2->data.F64[row];
            psF64* img3Row = c64Img3->data.F64[row];
            for (psU32 col=0;col<m;col++) {
                if (fabsf(img2Row[col] - row) > FLT_EPSILON) {
                    diag("psImageReal didn't return the real portion at n=%d", n);
                    errorFlag = true;
                }
                if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
                    diag("psImageImaginary didn't return the imag portion at n=%d", n);
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psImageReal() and psImageImaginary() returned the correct data");
        psFree(c64Img);
        psFree(c64Img2);
        psFree(c64Img3);
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }



    {
        psMemId id = psMemGetId();
        psImage* ncImg = NULL;
        psImage* ncImg2 = NULL;
        psImage* ncImg3 = NULL;
        psU32 m = 128;
        psU32 n = 64;

        GENIMAGE(ncImg,m,n,F32,row+col);
        ncImg2 = psImageReal(ncImg2,ncImg);
        ncImg3 = psImageImaginary(ncImg3,ncImg);
        ok(ncImg2 != NULL, "psImageReal() returned non-NULL");
        ok(ncImg2->type.type == PS_TYPE_F32, "psImageReal() returned the correct type");
        ok(ncImg3 != NULL, "psImageImaginary() returned non-NULL");
        ok(ncImg3->type.type == PS_TYPE_F32, "psImageImaginary() returned the correct type");

        bool errorFlag = false;
        for(psU32 row=0; row<n; row++)
        {
            psF32* ncImg2Row = ncImg2->data.F32[row];
            psF32* ncImg3Row = ncImg3->data.F32[row];
            for(psU32 col=0; col<m; col++) {
                if(fabsf(ncImg2Row[col] - (row+col)) > FLT_EPSILON) {
                    diag("psImageReal didn't return the real portion");
                    errorFlag = true;
                }
                if(fabsf(ncImg3Row[col] - 0) > FLT_EPSILON) {
                    diag("psImageImaginary didn't return the imaginary portion");
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psImageReal() and psImageImaginary() returned the correct data");
        psFree(ncImg);
        psFree(ncImg2);
        psFree(ncImg3);
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }



    // Perform psImageReal with null input
    {
        psMemId id = psMemGetId();
        ok(psImageReal(NULL,NULL) == NULL, "psImageReal() returned NULL with NULL inputs");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // Perform psImageImaginary with null input
    {
        psMemId id = psMemGetId();
        ok(psImageImaginary(NULL,NULL) == NULL, "psImageImaginary() returned NULL with NULL inputs");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testImageComplex()
    // 1. create two unique psF32 vectors of the same size
    // 2. call psImageComplex
    // 3. verify that the result is a psC32
    // 4. use crealf and cimagf on step 2 results
    // 5. compare step 4 results to input.
    // 6. create a psF32 and a psF64 vector of the same size
    // 7. call psImageComplex
    // 8. verify that an appropriate error occurred.
    // 9. create two psf32 vectors of different sizes
    // 10. call psImageComplex
    // 11. verify thet an appropriate error occurred.
    {
        psMemId id = psMemGetId();
        psImage* img = NULL;
        psImage* img2 = NULL;
        psImage* img3 = NULL;
        psImage* c64Img = NULL;
        psImage* c64Img2 = NULL;
        psImage* c64Img3 = NULL;
        psImage* pImg = NULL;
        psImage* pImg2 = NULL;
        psImage* pImg3 = NULL;
        psU32 m = 128;
        psU32 n = 64;
        GENIMAGE(img,m,n,F32,row);
        GENIMAGE(img2,m,n,F32,col);
        GENIMAGE(c64Img,m,n,F64,row);
        GENIMAGE(c64Img2,m,n,F64,col);
        GENIMAGE(pImg,m,n,S16,row);
        GENIMAGE(pImg2,m,n,S16,col);
        img3 = psImageComplex(img3,img,img2);
        ok(img3 != NULL, "psImageComplex() returned non-NULL");
        ok(img3->type.type == PS_TYPE_C32, "psImageComplex() returned the correct type");
        c64Img3 = psImageComplex(c64Img3,c64Img,c64Img2);
        ok(c64Img3 != NULL, "psImageComplex() returned non-NULL");
        ok(c64Img3->type.type == PS_TYPE_C64, "psImageComplex() returned the correct type");
        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psC32* img3Row = img3->data.C32[row];
            psC64* c64Img3Row = c64Img3->data.C64[row];
            for (psU32 col=0;col<m;col++) {
                if (fabsf(crealf(img3Row[col]) - row) > FLT_EPSILON ||
                        fabsf(cimagf(img3Row[col]) - col) > FLT_EPSILON) {
                    diag("psImageComplex result is incorrect (%d,%d, %.2f+%.2fi)",
                         col,row,crealf(img3Row[col]),cimagf(img3Row[col]));
                    errorFlag = true;
                }
                if (fabsf(crealf(c64Img3Row[col]) - row) > FLT_EPSILON ||
                        fabsf(cimagf(c64Img3Row[col]) - col) > FLT_EPSILON) {
                    diag("psImageComplex result is incorrect");
                    errorFlag = true;

                }
            }
        }
        ok(!errorFlag, "psImageComplex() returned correct results");

        img2 = psImageRecycle(img2,m,n,PS_TYPE_F64);

        // Following should be an error (type mismatch)
        // Verify that an appropriate error occurred
        img3 = psImageComplex(img3,img,img2);
        ok(img3 == NULL, "psImageComplex() returned NULL when input types mismatched.");

        // Following should be an error (size mismatch)
        img2 = psImageRecycle(img2,m/2,n,PS_TYPE_F32);
        img3 = psImageComplex(img3,img,img2);
        ok(img3 == NULL, "psImageComplex() returned a NULL when input sizes mismatched");

        // Perform psImageComplex with incorrect type
        // Following should generate an error message.");
        pImg3 = psImageComplex(pImg3,pImg,pImg2);
        ok(pImg3 == NULL, "psImageComplex returned NULL with incorrect type input");

        psFree(img);
        psFree(img2);
        psFree(img3);
        psFree(c64Img);
        psFree(c64Img2);
        psFree(c64Img3);
        psFree(pImg);
        psFree(pImg2);

        // Perform psImageComplex with null input
        ok(psImageComplex(NULL,NULL,NULL) == NULL, "psImageComplex() returned NULL with null input");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // 1. create a psC32 with unique real and imaginary values.
    // 2. call psImageConjugate
    // 3. verify result is psC32
    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    // testImageConjugate()
    {
        psMemId id = psMemGetId();
        psImage* img = NULL;
        psImage* img2 = NULL;
        psImage* c64Img = NULL;
        psImage* c64Img2 = NULL;
        psImage* pImg = NULL;
        psImage* pImg2 = NULL;
        psU32 m = 128;
        psU32 n = 64;
        GENIMAGE(img,m,n,C32, row + I * col);
        GENIMAGE(c64Img,m,n,C64,row + I*col);
        GENIMAGE(pImg,m,n,F32,row+col);


        img2 = psImageConjugate(img2,img);
        ok(img2 != NULL, "psImageConjugate() returned non-NULL");
        c64Img2 = psImageConjugate(c64Img2,c64Img);
        ok(c64Img2 != NULL, "psImageConjugate() returned non-NULL");
        pImg2 = psImageConjugate(pImg2,pImg);
        ok(pImg2 != NULL, "psImageConjugate() returned non-NULL");

        ok(img2->type.type == PS_TYPE_C32, "psImageConjugate() returned the correct type");
        ok(c64Img2->type.type == PS_TYPE_C64, "psImageConjugate() returned the correct type");
        ok(pImg2->type.type == PS_TYPE_F32, "psImageConjugate() returned the correct type");

        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psC32* img2Row = img2->data.C32[row];
            psC64* c64Img2Row = c64Img2->data.C64[row];
            psF32* pImg2Row = pImg2->data.F32[row];
            for (psU32 col=0;col<m;col++) {
                if (fabsf(crealf(img2Row[col]) - row) > FLT_EPSILON ||
                        fabsf(cimagf(img2Row[col]) + col) > FLT_EPSILON) {
                    diag("psImageComplex result is incorrect (%d,%d, %.2f+%.2fi)",
                         col,row,crealf(img2Row[col]),cimagf(img2Row[col]));
                    errorFlag = true;
                }
                if (fabsf(crealf(c64Img2Row[col]) - row) > FLT_EPSILON ||
                        fabsf(cimagf(c64Img2Row[col]) + col) > FLT_EPSILON) {
                    diag("psImageComplex result is incorrect");
                    errorFlag = true;
                }
                if (fabsf(pImg2Row[col] - (row+col)) > FLT_EPSILON) {
                    diag("psImageComplex result is incorrect");
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psImageConjugate() generated the correct data values");
        psFree(img);
        psFree(img2);
        psFree(c64Img);
        psFree(c64Img2);
        psFree(pImg);
        psFree(pImg2);

        // Perform psImageConjugate with null input
        ok(psImageConjugate(NULL,NULL) == NULL, "psImageConjugate() returned NULL with NULL input");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }

    // testImagePowerSpectrum()
    // 1. create a psC32 vector with unique real and imaginary components
    // 2. call psImagePowerSpectrum
    // 3. verify result is psF32
    // 4. verify the values are the square of the absolute values of the original
    {
        psMemId id = psMemGetId();
        psImage* img = NULL;
        psImage* img2 = NULL;
        psImage* c64Img = NULL;
        psImage* c64Img2 = NULL;
        psImage* pImg = NULL;
        psU32 m = 128;
        psU32 n = 64;

        GENIMAGE(img,m,n,C32, row + I * col);
        GENIMAGE(c64Img,m,n,C64,row+I*col);
        GENIMAGE(pImg,m,n,F32,row+col);

        img2 = psImagePowerSpectrum(img2,img);
        c64Img2 = psImagePowerSpectrum(c64Img2, c64Img);
        ok(img2 != NULL, "psImagePowerSpectrum() returned non-NULL");
        ok(c64Img2 != NULL, "psImagePowerSpectrum() returned non-NULL");
        ok(img2->type.type == PS_TYPE_F32, "psImagePowerSpectrum() returned the correct type");
        ok(c64Img2->type.type == PS_TYPE_F64, "psImagePowerSpectrum() returned the correct type");

        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psC32* imgRow = img->data.C32[row];
            psF32* img2Row = img2->data.F32[row];
            psC64* c64ImgRow = c64Img->data.C64[row];
            psF64* c64Img2Row = c64Img2->data.F64[row];
            for (psU32 col=0;col<m;col++) {
                psF32 power = cabs(imgRow[col]);
                psF64 power64 = cabs(c64ImgRow[col]);
                power *= power/n/n/m/m;
                power64 *= power64/n/n/m/m;

                if (fabsf(img2Row[col] - power) > 2.0f*FLT_EPSILON) {
                    diag("psImagePowerSpectrum result is incorrect (%d,%d, %.2f vs %.2f)",
                         col,row,img2Row[col],power);
                    errorFlag = true;
                }
                if (fabsf(c64Img2Row[col] - power64) > 2.0f*FLT_EPSILON) {
                    diag("psImagePowerSpectrum result is incorrect");
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psImagePowerSPectrum() generated the correct data values");
        psFree(img);
        psFree(img2);
        psFree(c64Img);
        psFree(c64Img2);

        // Perform psImagePowerSpectrum with incorrect input
        // Following should generate error message
        // XXX: Verify error
        ok(psImagePowerSpectrum(NULL,pImg) == NULL, "psImagePowerSpectrum() returned NULL with incorrect type");
        psFree(pImg);

        // Perform psImagePowerSpectrum with NULL input
        ok(psImagePowerSpectrum(NULL,NULL) == NULL, "psImagePowerSpectrum() returned NULL with NULL input");
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


    // testImageRealImaginary()
    // 1. create a C32 complex image with distinctly different real and imaginary parts
    // 2. call psImageReal and psImageImaginary
    // 3. compare results to the real/imaginary components of input
    // XXX: If we put this block above the previous 2 blocks, then we get memory errors
    {
        psMemId id = psMemGetId();
        psImage* img = NULL;
        psImage* img3 = NULL;
        psU32 m = 128;
        psU32 n = 64;
        GENIMAGE(img,m,n,C32, row + I * col);
        psImage* img2 = psImageReal(img2,img);
        ok(img2 != NULL, "psImageReal returned non-NULL");
        skip_start(img2 == NULL, 4, "Skipping tests because psImageReal() returned NULL");
        ok(img2->type.type == PS_TYPE_F32, "psImageReal returned the correct type");
        img3 = psImageImaginary(img3,img);
        ok(img3 != NULL, "psImageImaginary() returned non-NULL");
        skip_start(img3 == NULL, 2, "Skipping tests because psImageImaginary() returned NULL");
        ok(img3->type.type == PS_TYPE_F32, "psImageImaginary() returned the correct type");

        bool errorFlag = false;
        for (psU32 row=0;row<n;row++)
        {
            psF32* img2Row = img2->data.F32[row];
            psF32* img3Row = img3->data.F32[row];
            for (psU32 col=0;col<m;col++) {
                if (fabsf(img2Row[col] - row) > FLT_EPSILON) {
                    diag("psImageReal didn't return the real portion at n=%d", n);
                    errorFlag = true;
                }
                if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
                    diag("psImageImaginary didn't return the imag portion at n=%d", n);
                    errorFlag = true;
                }
            }
        }
        ok(!errorFlag, "psImageReal(), psImageImaginary() returned the correct data");
        skip_end();
        skip_end();
        psFree(img);
        psFree(img2);
        psFree(img3);
        ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks");
    }


}
