Index: /trunk/psLib/test/dataManip/Makefile
===================================================================
--- /trunk/psLib/test/dataManip/Makefile	(revision 1010)
+++ /trunk/psLib/test/dataManip/Makefile	(revision 1011)
@@ -3,6 +3,6 @@
 ##  Makefile:   test/sysUtils
 ##
-##  $Revision: 1.13 $  $Name: not supported by cvs2svn $
-##  $Date: 2004-06-11 20:36:04 $
+##  $Revision: 1.14 $  $Name: not supported by cvs2svn $
+##  $Date: 2004-06-12 02:17:42 $
 ##
 ##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -44,5 +44,7 @@
  tst_psImageStats01 \
  tst_psImageStats02 \
- tst_psImageIO
+ tst_psImageIO \
+ tst_psVectorFFT \
+ tst_psImageFFT
 
 DEPENDENCIES = $(addprefix builddir/,$(addsuffix .d,$(TARGET)))
Index: /trunk/psLib/test/image/tst_psImageFFT.c
===================================================================
--- /trunk/psLib/test/image/tst_psImageFFT.c	(revision 1011)
+++ /trunk/psLib/test/image/tst_psImageFFT.c	(revision 1011)
@@ -0,0 +1,388 @@
+/** @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: 2004-06-12 02:17:42 $
+ *
+ *  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);
+#if 0
+static int testImageRealImaginary(void);
+static int testImageComplex(void);
+static int testImageConjugate(void);
+static int testImagePowerSpectrum(void);
+#endif
+
+testDescription tests[] = {
+                              {testImageFFT,"600-testImageFFT",0},
+                              #if 0
+                              {testImageRealImaginary,"601-testImageRealImaginary",0},
+                              {testImageComplex,"602-testImageComplex",0},
+                              {testImageConjugate,"603-testImageConjugate",0},
+                              {testImagePowerSpectrum,"604-testImagePowerSpectrum",0},
+                              #endif
+                              {NULL}
+                          };
+
+int main()
+{
+    psLogSetLevel(PS_LOG_INFO);
+
+    if (! runTestSuite(stderr,"psFFT",tests)) {
+        psAbort(__FILE__,"One or more tests failed");
+    }
+    return 0;
+}
+
+int testImageFFT(void)
+{
+    psImage* img = NULL;
+    psImage* img2 = NULL;
+    psImage* img3 = NULL;
+    unsigned int m = 128;
+    unsigned int n = 64;
+
+    /*
+    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(sqrt((n/2-row)*(n/2-row))/n*2*M_PI));
+    GENIMAGE(img,m,n,F32, 1);
+    psImageWriteSection(img,0,0,0,NULL,0,"test.fits");
+
+    // 2. perform a forward transform
+    img2 = psImageFFT(NULL,img,PS_FFT_FORWARD);
+    if (img2->type.type != PS_TYPE_C32) {
+        psError(__func__,"FFT didn't produce complex values?");
+        return 1;
+    }
+
+    img3 = psImageReal(img3,img2);
+    psImageWriteSection(img2,0,0,0,NULL,1,"test.fits");
+
+    #if 0
+    // 3. verify that the only significant component cooresponds to the freqency of the input in step 1.
+    for (unsigned int n = 0; n<100; n++) {
+        if (n==1 || n==99) {
+            if (fabsf(cabsf(img2->data.C32[n]) - 50.0f) > 0.1f) {
+                psError(__func__,"FFT didn't work for vector (n=%d)",n);
+                return 2;
+            }
+        } else {
+            if (fabsf(cabsf(img2->data.C32[n])) > 0.1f) {
+                psError(__func__,"FFT didn't work for vector (n=%d)",n);
+                return 3;
+            }
+        }
+    }
+
+    // 4. perform a reverse transform
+    img3 = psImageFFT(NULL,img2,PS_FFT_REVERSE);
+    if (img3->type.type != PS_TYPE_C32) {
+        psError(__func__,"FFT didn't produce complex values?");
+        return 4;
+    }
+    for (unsigned int n = 0; n<100; n++) {
+        psF32 val =sinf((psF32)n / 50.0f * M_PI);
+        psF32 vecVal = crealf(img3->data.C32[n])/100;
+        if (fabsf(vecVal - val) > 0.1f) {
+            psError(__func__,"Reverse FFT didn't give me the original vector back (n=%d) (%.2f vs %.2f)",
+                    n,vecVal,val);
+            return 5;
+        }
+    }
+    #endif
+    psImageFree(img);
+    psImageFree(img2);
+    psImageFree(img3);
+
+    return 0;
+}
+#if 0
+int testImageRealImaginary(void)
+{
+    psImage* img = NULL;
+    psImage* img2 = NULL;
+    psImage* img3 = NULL;
+
+    /*
+    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.
+    img=psImageAlloc(100,PS_TYPE_C32);
+    img->n = img->nalloc;
+    for (unsigned int n = 0; n<100; n++) {
+        img->data.C32[n] = n + I * (n*2);
+    }
+
+    // 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 n = 0; n<100; n++) {
+        psF32 r = n;
+        psF32 i = (n*2);
+        if (fabsf(img2->data.F32[n] -r) > FLT_EPSILON) {
+            psError(__func__,"psImageReal didn't return the real portion at n=%d",
+                    n);
+            return 5;
+        }
+        if (fabsf(img3->data.F32[n] -i) > FLT_EPSILON) {
+            psError(__func__,"psImageImaginary didn't return the real portion at n=%d",
+                    n);
+            return 6;
+        }
+    }
+
+    psImageFree(img);
+    psImageFree(img2);
+    psImageFree(img3);
+
+    return 0;
+}
+
+int testImageComplex(void)
+{
+    psImage* img = NULL;
+    psImage* img2 = NULL;
+    psImage* img3 = NULL;
+
+    /*
+    1. create two unique psF32 vectors of the same size
+    2. call psImageComplex
+    3. verify that the result is a psC32
+    4. call psImageReal and psImageImaginary 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
+    img=psImageAlloc(100,PS_TYPE_F32);
+    img2=psImageAlloc(100,PS_TYPE_F32);
+    img->n = img->nalloc;
+    img2->n = img2->nalloc;
+    for (unsigned int n = 0; n<100; n++) {
+        img->data.F32[n] = n;
+        img2->data.F32[n] = (n*2);
+    }
+
+    // 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__,"Vector 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 n = 0; n<100; n++) {
+        if (fabsf(crealf(img3->data.C32[n]) - n) > FLT_EPSILON ||
+                fabsf(cimagf(img3->data.C32[n]) - (n*2)) > FLT_EPSILON) {
+            psError(__func__,"psImageComplex result is invalid (n=%d, %.2f+%.2fi)",
+                    n,crealf(img3->data.C32[n]),cimagf(img3->data.C32[n]));
+            return 2;
+        };
+    }
+
+
+    // 6. create a psF32 and a psF64 vector of the same size
+    img2 = psImageRecycle(img2,PS_TYPE_F64, 100);
+
+    // 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 vector though input types mismatched.");
+        return 3;
+    }
+
+    // 9. create two psf32 vectors of different sizes
+    img2 = psImageRecycle(img2,PS_TYPE_F32,200);
+
+    // 10. call psImageComplex
+    img3 = psImageComplex(img3,img,img2);
+
+    // 11. verify thet an appropriate error occurred. (actually, it isn't an error...)
+    if (img3->n != 100) {
+        psError(__func__,"psImageComplex returned a vector though input sizes mismatched.");
+        return 4;
+    }
+
+    psImageFree(img);
+    psImageFree(img2);
+    psImageFree(img3);
+
+    return 0;
+}
+
+int testImageConjugate(void)
+{
+    psImage* img = NULL;
+    psImage* img2 = NULL;
+
+    /*
+    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.
+    img=psImageAlloc(100,PS_TYPE_C32);
+    img->n = img->nalloc;
+    for (unsigned int n = 0; n<100; n++) {
+        img->data.C32[n] = n + I * (n*2);
+    }
+
+    // 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 vector");
+        return 1;
+    }
+
+    // 4. verify each value is conjugate of input (a+bi -> a-bi)
+    for (unsigned int n = 0; n<100; n++) {
+        if (fabsf(crealf(img->data.C32[n]) - crealf(img2->data.C32[n])) > FLT_EPSILON ||
+                fabsf(cimagf(img->data.C32[n]) + cimagf(img2->data.C32[n])) > FLT_EPSILON) {
+            psError(__func__,"psImageComplex result is invalid (n=%d, %.2f+%.2fi)",
+                    n,crealf(img2->data.C32[n]),cimagf(img2->data.C32[n]));
+            return 2;
+        };
+    }
+
+    psImageFree(img);
+    psImageFree(img2);
+
+    return 0;
+}
+
+int testImagePowerSpectrum(void)
+{
+    psImage* img = NULL;
+    psImage* img2 = NULL;
+    psF32 val;
+
+    /*
+    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
+    img=psImageAlloc(100,PS_TYPE_C32);
+    img->n = img->nalloc;
+    for (unsigned int n = 0; n<100; n++) {
+        img->data.C32[n] = n + I * sinf(((psF32)n) / 50.f * M_PI);
+    }
+
+    // 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
+    //   (ADD specifies something else)
+    //   P_0 = |C_0|^2/N^2
+    //   P_j = (|C_j|^2+|C_N-j|^2)/N^2
+    //   P_N/2 = |C_N/2|^2/N^2
+    //  where j = 1,2,...,(N/2-1)
+
+    val = cabsf(img->data.C32[0])*cabsf(img->data.C32[0])/100/100;
+    if (fabsf(img2->data.F32[0] - val) > FLT_EPSILON) {
+        psError(__func__,"psImagePowerSpectrum result is invalid (n=0, %.2f %.2f)",
+                img2->data.F32[0],val);
+        return 2;
+    };
+
+    for (unsigned int n = 1; n<50; n++) {
+        val = ( cabsf(img->data.C32[n])*cabsf(img->data.C32[n])+
+                cabsf(img->data.C32[100-n])*cabsf(img->data.C32[100-n]) ) /100/100;
+
+        if (fabsf(val - img2->data.F32[n]) > FLT_EPSILON) {
+            psError(__func__,"psImagePowerSpectrum result is invalid (n=%d, %.2f %.2f)",
+                    n,img2->data.F32[n],val);
+            return 2;
+        };
+    }
+
+    val = cabsf(img->data.C32[50])*cabsf(img->data.C32[50])/100/100;
+    if (fabsf(img2->data.F32[50] - val) > FLT_EPSILON) {
+        psError(__func__,"psImagePowerSpectrum result is invalid (n=50, %.2f %.2f)",
+                img2->data.F32[0],val);
+        return 2;
+    };
+
+    psImageFree(img);
+    psImageFree(img2);
+
+    return 0;
+}
+#endif