/** @file  tst_psImageFFT.c
 *
 *  @brief Contains the tests for psFFT.[ch]
 *
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2004-12-15 03:57:03 $
 *
 *  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 (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); \
    } \
}

static psS32 testImageFFT(void);
static psS32 testImageRealImaginary(void);
static psS32 testImageComplex(void);
static psS32 testImageConjugate(void);
static psS32 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}
                          };

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

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

psS32 testImageFFT(void)
{
    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 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(PS_ERR_UNKNOWN, true,"FFT didn't produce complex values?");
        return 1;
    }
    if (img2->numCols != m || img2->numRows != n) {
        psError(PS_ERR_UNKNOWN, true,"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 (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 (n-1,0)
                if (! (col == 0 && (row == 1 || row == n-1))
                        && ! (row == 0 && (col==1 || col == n-1)) ) {
                    psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"FFT didn't produce complex values?");
        return 5;
    }

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

    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) {
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"FFT asked to make real result, but I got a %s type image?",typeStr);
        return 8;
    }

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

    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) {
                psError(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageFFT should return null for an invalid FFT direction.");
        return 11;
    }
    psFree(img4);

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

    return 0;
}

psS32 testImageRealImaginary(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* img3 = NULL;
    psImage* c64Img = NULL;
    psImage* c64Img2 = NULL;
    psImage* c64Img3 = NULL;
    psImage* ncImg = NULL;
    psImage* ncImg2 = NULL;
    psImage* ncImg3 = NULL;

    psU32 m = 128;
    psU32 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(PS_ERR_UNKNOWN, true,"psImageReal returned a NULL?");
        return 1;
    }
    if (img2->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN, true,"psImageReal returned a wrong type (%d)?",
                img2->type.type);
        return 2;
    }

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

    // 3. compare results to the real/imaginary components of input
    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) {
                psError(PS_ERR_UNKNOWN, true,"psImageReal didn't return the real portion at n=%d",
                        n);
                return 5;
            }
            if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN, true,"psImageImaginary didn't return the imag portion at n=%d",
                        n);
                return 6;
            }
        }
    }

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

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

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

    // 2. call psImageReal and psImageImaginary
    c64Img2 = psImageReal(c64Img2,c64Img);
    if (c64Img2 == NULL) {
        psError(PS_ERR_UNKNOWN, true,"psImageReal returned a NULL?");
        return 1;
    }
    if (c64Img2->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN, true,"psImageReal returned a wrong type (%d)?",
                img2->type.type);
        return 2;
    }

    c64Img3 = psImageImaginary(c64Img3,c64Img);
    if (c64Img3 == NULL) {
        psError(PS_ERR_UNKNOWN, true,"psImageImaginary returned a NULL?");
        return 3;
    }
    if (c64Img3->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN, true,"psImageImaginary returned a wrong type (%d)?",
                c64Img3->type.type);
        return 4;
    }

    // 3. compare results to the real/imaginary components of input
    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) {
                psError(PS_ERR_UNKNOWN, true,"psImageReal didn't return the real portion at n=%d",
                        n);
                return 5;
            }
            if (fabsf(img3Row[col] - col) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN, true,"psImageImaginary didn't return the imag portion at n=%d",
                        n);
                return 6;
            }
        }
    }

    GENIMAGE(ncImg,m,n,F32,row+col);
    ncImg2 = psImageReal(ncImg2,ncImg);
    ncImg3 = psImageImaginary(ncImg3,ncImg);
    if(ncImg2 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psImageReal returned NULL");
        return 40;
    }
    if(ncImg2->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psImageReal returned a wrong type");
        return 41;
    }
    if(ncImg3 == NULL) {
        psError(PS_ERR_UNKNOWN,true,"psImageImaginary returned NULL");
        return 42;
    }
    if(ncImg3->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psImageImaginary returned NULL");
        return 43;
    }
    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) {
                psError(PS_ERR_UNKNOWN,true,"psImageReal didn't return the real portion");
                return 45;
            }
            if(fabsf(ncImg3Row[col] - 0) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN,true,"psImageImaginary didn't return the imaginary portion");
                return 46;
            }
        }
    }
    psFree(ncImg);
    psFree(ncImg2);
    psFree(ncImg3);

    psFree(c64Img);
    psFree(c64Img2);
    psFree(c64Img3);

    // Perform psImageReal with null input
    if(psImageReal(NULL,NULL) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImageReal did not return null with null input");
        return 10;
    }

    // Perform psImageImaginary with null input
    if(psImageImaginary(NULL,NULL) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImageImaginary did not return null with null input");
        return 10;
    }

    return 0;
}

psS32 testImageComplex(void)
{
    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;

    /*
    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);
    GENIMAGE(c64Img,m,n,F64,row);
    GENIMAGE(c64Img2,m,n,F64,col);
    GENIMAGE(pImg,m,n,S16,row);
    GENIMAGE(pImg2,m,n,S16,col);

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

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

    // 4. call psImageReal and psImageImaginary on step 2 results (not needed, just use crealf/cimagf)
    // 5. compare step 4 results to input.
    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) {
                psError(PS_ERR_UNKNOWN, true,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
                        col,row,crealf(img3Row[col]),cimagf(img3Row[col]));
                return 2;
            }
            if (fabsf(crealf(c64Img3Row[col]) - row) > FLT_EPSILON ||
                    fabsf(cimagf(c64Img3Row[col]) - col) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN,true,"psImageComplex result is invalid");
                return 3;
            }
        }
    }

    // 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(PS_ERR_UNKNOWN, true,"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(PS_ERR_UNKNOWN, true,"psImageComplex returned a image though input sizes mismatched.");
        return 4;
    }

    // Perform psImageComplex with invalid type
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message.");
    pImg3 = psImageComplex(pImg3,pImg,pImg2);
    if ( pImg3 != NULL) {
        psError(PS_ERR_UNKNOWN,true,"psImageComplex did not return NULL");
        return 50;
    }

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

    // Perform psImageComplex with null input
    if(psImageComplex(NULL,NULL,NULL) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImageComplex did not return null with null input");
        return 20;
    }

    return 0;
}

psS32 testImageConjugate(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* c64Img = NULL;
    psImage* c64Img2 = NULL;
    psImage* pImg = NULL;
    psImage* pImg2 = NULL;

    psU32 m = 128;
    psU32 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);
    GENIMAGE(c64Img,m,n,C64,row + I*col);
    GENIMAGE(pImg,m,n,F32,row+col);

    // 2. call psImageConjugate
    img2 = psImageConjugate(img2,img);
    c64Img2 = psImageConjugate(c64Img2,c64Img);
    pImg2 = psImageConjugate(pImg2,pImg);

    // 3. verify result is psC32
    if (img2->type.type != PS_TYPE_C32) {
        psError(PS_ERR_UNKNOWN, true,"the psImageConjugate didn't return a C32 image.");
        return 1;
    }
    if (c64Img2->type.type != PS_TYPE_C64) {
        psError(PS_ERR_UNKNOWN,true,"psImageConjugate did not return a C64 image");
        return 10;
    }
    if (pImg2->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN,true,"psImageConjugate did not return a F32 image");
        return 11;
    }

    // 4. verify each value is conjugate of input (a+bi -> a-bi)
    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) {
                psError(PS_ERR_UNKNOWN, true,"psImageComplex result is invalid (%d,%d, %.2f+%.2fi)",
                        col,row,crealf(img2Row[col]),cimagf(img2Row[col]));
                return 2;
            }
            if (fabsf(crealf(c64Img2Row[col]) - row) > FLT_EPSILON ||
                    fabsf(cimagf(c64Img2Row[col]) + col) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN,true,"psImageComplex result is invalid");
                return 20;
            }
            if (fabsf(pImg2Row[col] - (row+col)) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN,true,"psImageComplex result is invalid");
                return 21;
            }
        }
    }

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

    // Perform psImageConjugate with null input
    if(psImageConjugate(NULL,NULL) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImageConjugate did not return null with null input");
        return 30;
    }

    return 0;
}

psS32 testImagePowerSpectrum(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psImage* c64Img = NULL;
    psImage* c64Img2 = NULL;
    psImage* pImg = NULL;

    psU32 m = 128;
    psU32 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);
    GENIMAGE(c64Img,m,n,C64,row+I*col);
    GENIMAGE(pImg,m,n,F32,row+col);

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

    // 3. verify result is psF32
    if (img2->type.type != PS_TYPE_F32) {
        psError(PS_ERR_UNKNOWN, true,"the type was not PS_TYPE_F32.");
        return 1;
    }
    if (c64Img2->type.type != PS_TYPE_F64) {
        psError(PS_ERR_UNKNOWN,true,"psImagePowerSpectrum did not return PS_TYPE_F64");
        return 10;
    }

    // 4. verify the values are the square of the absolute values of the original
    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) {
                psError(PS_ERR_UNKNOWN, true,"psImagePowerSpectrum result is invalid (%d,%d, %.2f vs %.2f)",
                        col,row,img2Row[col],power);
                return 2;
            }
            if (fabsf(c64Img2Row[col] - power64) > 2.0f*FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN,true,"psImagePowerSpectrum result is invalid");
                return 22;
            }
        }
    }

    psFree(img);
    psFree(img2);
    psFree(c64Img);
    psFree(c64Img2);

    // Perform psImagePowerSpectrum with invalid input
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate error message.");
    if(psImagePowerSpectrum(NULL,pImg) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImagePowerSpectrum did not return null with invalid type");
        return 41;
    }
    psFree(pImg);

    // Perform psImagePowerSpectrum with null input
    if(psImagePowerSpectrum(NULL,NULL) != NULL ) {
        psError(PS_ERR_UNKNOWN,true,"psImagePowerSpectrum did not return null with null input");
        return 40;
    }

    return 0;
}
