/** @file  tst_psImageFFT.c
 *
 *  @brief Contains the tests for psFFT.[ch]
 *
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.11 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-10-13 19:50:12 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 */

#include <math.h>
#include <float.h>

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

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

static int testImageFFT(void);
static int testImageRealImaginary(void);
static int testImageComplex(void);
static int testImageConjugate(void);
static int testImagePowerSpectrum(void);

testDescription tests[] = {
                              {
                                  testImageFFT,632,"psImageFFT",0,false
                              },
                              {
                                  testImageRealImaginary,633,"psImageRealImaginary",0,false
                              },
                              {
                                  testImageComplex,634,"psImageComplex",0,false
                              },
                              {
                                  testImageConjugate,635,"psImageConjugate",0,false
                              },
                              {
                                  testImagePowerSpectrum,636,"psImagePowerSpectrum",0,false
                              },
                              {
                                  NULL
                              }
                          };

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

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

int testImageFFT(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    unsigned int m = 128;
    unsigned int n = 64;
    psImage* img4 = NULL;
    psImage* img5 = NULL;

    /*
    1. assign a image to a radial sinisoid
    2. perform a forward transform
    3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    4. perform a 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)/16.0f*M_PI)+sinf((64.0f-col)/32.0f*M_PI));

    // 2. perform a forward transform
    img2 = psImageFFT(img2,img,PS_FFT_FORWARD);
    if (img2->type.type != PS_TYPE_C32) {
        psError(__func__,"FFT didn't produce complex values?");
        return 1;
    }
    if (img2->numCols != m || img2->numRows != n) {
        psError(__func__,"FFT didn't produce proper size result (%dx%d vs. expected %dx%d).",
                img2->numCols,img2->numRows,m,n);
        return 2;
    }

    // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
    for (unsigned int row=0;row<n;row++) {
        psC32* img2Row = img2->data.C32[row];
        for (unsigned int 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 (n-1,0)
                if (! (col == 0 && (row == 1 || row == n-1))
                        && ! (row == 0 && (col==1 || col == n-1)) ) {
                    psError(__func__,"Result invalid at %d,%d (%.2f)",col,row,mag);
                    return 3;
                }
            } else
                if ( (col == 0 && (row == 1 || row == n-1))
                        || (row == 0 && (col==1 || col == n-1)) ) {
                    psError(__func__,"Result invalid at %d,%d (%.2f)",col,row,mag);
                    return 4;
                }
        }
    }


    // 4. perform a reverse transform
    img3 = psImageFFT(img3,img2,PS_FFT_REVERSE);

    if (img3->type.type != PS_TYPE_C32) {
        psError(__func__,"FFT didn't produce complex values?");
        return 5;
    }

    if (img3->numCols != m || img3->numRows != n) {
        psError(__func__,"FFT didn't produce proper size result (%dx%d vs. expected %dx%d).",
                img3->numCols,img3->numRows,m,n);
        return 6;
    }

    for (unsigned int row=0;row<n;row++) {
        psC32* img3Row = img3->data.C32[row];
        psF32* imgRow = img->data.F32[row];
        for (unsigned int col=0;col<m;col++) {
            psF32 pixel = creal(img3Row[col])/m/n;
            if (fabsf(pixel-imgRow[col]) > 0.1) {
                psError(__func__,"Reverse FFT didn't give original image back (%d,%d %.2f vs %.2f)",
                        col,row,pixel,imgRow[col]);
                return 7;
            }
        }
    }

    // 4. perform a reverse transform to real result
    img3 = psImageFFT(img3,img2,PS_FFT_REVERSE|PS_FFT_REAL_RESULT);

    if (img3->type.type != PS_TYPE_F32) {
        char* typeStr;
        PS_TYPE_NAME(typeStr,img3->type.type)
        psError(__func__,"FFT asked to make real result, but I got a %s type image?",typeStr);
        return 8;
    }

    if (img3->numCols != m || img3->numRows != n) {
        psError(__func__,"FFT didn't produce proper size result (%dx%d vs. expected %dx%d).",
                img3->numCols,img3->numRows,m,n);
        return 9;
    }

    for (unsigned int row=0;row<n;row++) {
        psF32* img3Row = img3->data.F32[row];
        psF32* imgRow = img->data.F32[row];
        for (unsigned int col=0;col<m;col++) {
            psF32 pixel = img3Row[col]/m/n;
            if (fabsf(pixel-imgRow[col]) > 0.1) {
                psError(__func__,"Reverse FFT didn't give original image back (%d,%d %.2f vs %.2f)",
                        col,row,pixel,imgRow[col]);
                return 10;
            }
        }
    }

    // check if error occurs if FORWARD and REVERSE are both given.
    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error.");
    if (psImageFFT(NULL,img2,PS_FFT_REVERSE|PS_FFT_FORWARD) != NULL) {
        psError(__func__,"PS_FFT_REVERSE|PS_FFT_FORWARD option produced something?");
        return 11;
    }

    psLogMsg(__func__,PS_LOG_INFO,"Following should be an error.");
    if (psImageFFT(NULL,img2,PS_FFT_FORWARD|PS_FFT_REAL_RESULT) != NULL) {
        psError(__func__,"PS_FFT_FORWARD|PS_FFT_REAL_RESULT option produced something?");
        return 12;
    }

    /* Verify return null and program execution doesn't stop if input image is null */
    img4 = psImageFFT(NULL,NULL,PS_FFT_FORWARD);
    if (img4 != NULL ) {
        psError(__func__,"psImageFFT should return null for a null input image.");
        return 10;
    }

    /* Verify return null and program execution doesn't stop if input image is invalid direction */
    GENIMAGE(img4,8,8,S8,row+col);
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error for invalid direction.");
    img5 = psImageFFT(NULL,img4,PS_FFT_REAL_RESULT);
    if (img5 != NULL) {
        psError(__func__,"psImageFFT should return null for an invalid FFT direction.");
        return 11;
    }
    psFree(img4);

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

    return 0;
}

int testImageRealImaginary(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    unsigned int m = 128;
    unsigned int n = 64;

    /*
    1. create a C32 complex vector with distinctly different real and imaginary parts.
    2. call psImageReal and psImageImaginary
    3. compare results to the real/imaginary components of input
    */

    // 1. create a C32 complex vector with distinctly different real and imaginary parts.
    GENIMAGE(img,m,n,C32, row + I * col);

    // 2. call psImageReal and psImageImaginary
    img2 = psImageReal(img2,img);
    if (img2 == NULL) {
        psError(__func__,"psImageReal returned a NULL?");
        return 1;
    }
    if (img2->type.type != PS_TYPE_F32) {
        psError(__func__,"psImageReal returned a wrong type (%d)?",
                img2->type.type);
        return 2;
    }

    img3 = psImageImaginary(img3,img);
    if (img3 == NULL) {
        psError(__func__,"psImageImaginary returned a NULL?");
        return 3;
    }
    if (img3->type.type != PS_TYPE_F32) {
        psError(__func__,"psImageImaginary returned a wrong type (%d)?",
                img3->type.type);
        return 4;
    }

    // 3. compare results to the real/imaginary components of input
    for (unsigned int row=0;row<n;row++) {
        psF32* img2Row = img2->data.F32[row];
        psF32* img3Row = img3->data.F32[row];
        for (unsigned int col=0;col<m;col++) {
            if (fabsf(img2Row[col] - row) > FLT_EPSILON) {
                psError(__func__,"psImageReal didn't return the real portion at n=%d",
                        n);
                return 5;
            }
            if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
                psError(__func__,"psImageImaginary didn't return the imag portion at n=%d",
                        n);
                return 6;
            }
        }
    }

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

    return 0;
}

int testImageComplex(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    unsigned int m = 128;
    unsigned int n = 64;

    /*
    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.
    */

    // 1. create two unique psF32 vectors of the same size
    GENIMAGE(img,m,n,F32,row);
    GENIMAGE(img2,m,n,F32,col);

    // 2. call psImageComplex
    img3 = psImageComplex(img3,img,img2);

    // 3. verify that the result is a psC32
    if (img3->type.type != PS_TYPE_C32) {
        psError(__func__,"Image Type from psImageComplex is not complex? (%d)",
                img3->type.type);
        return 1;
    }

    // 4. call psImageReal and psImageImaginary on step 2 results (not needed, just use crealf/cimagf)
    // 5. compare step 4 results to input.
    for (unsigned int row=0;row<n;row++) {
        psC32* img3Row = img3->data.C32[row];
        for (unsigned int col=0;col<m;col++) {
            if (fabsf(crealf(img3Row[col]) - row) > FLT_EPSILON ||
                    fabsf(cimagf(img3Row[col]) - col) > FLT_EPSILON) {
                psError(__func__,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
                        col,row,crealf(img3Row[col]),cimagf(img3Row[col]));
                return 2;
            }
        }
    }

    // 6. create a psF32 and a psF64 image of the same size
    img2 = psImageRecycle(img2,m,n,PS_TYPE_F64);

    // 7. call psImageComplex
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error (type mismatch).");
    img3 = psImageComplex(img3,img,img2);
    // 8. verify that an appropriate error occurred. (this partially has to be done via inspection)
    if (img3 != NULL) {
        psError(__func__,"psImageComplex returned a image though input types mismatched.");
        return 3;
    }

    // 9. create two psf32 vectors of different sizes
    img2 = psImageRecycle(img2,m/2,n,PS_TYPE_F32);

    // 10. call psImageComplex
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error (size mismatch).");
    img3 = psImageComplex(img3,img,img2);

    // 11. verify thet an appropriate error occurred.
    if (img3 != NULL) {
        psError(__func__,"psImageComplex returned a image though input sizes mismatched.");
        return 4;
    }

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

    return 0;
}

int testImageConjugate(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    unsigned int m = 128;
    unsigned int n = 64;

    /*
    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)
    */

    // 1. create a psC32 with unique real and imaginary values.
    GENIMAGE(img,m,n,C32, row + I * col);

    // 2. call psImageConjugate
    img2 = psImageConjugate(img2,img);

    // 3. verify result is psC32
    if (img2->type.type != PS_TYPE_C32) {
        psError(__func__,"the psImageConjugate didn't return a C32 image.");
        return 1;
    }

    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    for (unsigned int row=0;row<n;row++) {
        psC32* img2Row = img2->data.C32[row];
        for (unsigned int col=0;col<m;col++) {
            if (fabsf(crealf(img2Row[col]) - row) > FLT_EPSILON ||
                    fabsf(cimagf(img2Row[col]) + col) > FLT_EPSILON) {
                psError(__func__,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
                        col,row,crealf(img2Row[col]),cimagf(img2Row[col]));
                return 2;
            }
        }
    }

    psFree(img);
    psFree(img2);

    return 0;
}

int testImagePowerSpectrum(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    unsigned int m = 128;
    unsigned int n = 64;

    /*
    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
    */

    // 1. create a psC32 vector with unique real and imaginary components
    GENIMAGE(img,m,n,C32, row + I * col);

    // 2. call psImagePowerSpectrum
    img2 = psImagePowerSpectrum(img2,img);

    // 3. verify result is psF32
    if (img2->type.type != PS_TYPE_F32) {
        psError(__func__,"the type was not PS_TYPE_F32.");
        return 1;
    }

    // 4. verify the values are the square of the absolute values of the original
    for (unsigned int row=0;row<n;row++) {
        psC32* imgRow = img->data.C32[row];
        psF32* img2Row = img2->data.F32[row];
        for (unsigned int col=0;col<m;col++) {
            psF32 power = cabs(imgRow[col]);
            power *= power/n/n/m/m;

            if (fabsf(img2Row[col] - power) > 2.0f*FLT_EPSILON) {
                psError(__func__,"psImagePowerSpectrum result is invalid (%d,%d, %.2f vs %.2f)",
                        col,row,img2Row[col],power);
                return 2;
            }
        }
    }

    psFree(img);
    psFree(img2);

    return 0;
}
