Index: /trunk/psLib/test/image/tst_psImageFFT.c
===================================================================
--- /trunk/psLib/test/image/tst_psImageFFT.c	(revision 1013)
+++ /trunk/psLib/test/image/tst_psImageFFT.c	(revision 1014)
@@ -6,6 +6,6 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-06-12 02:17:42 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-06-12 19:46:04 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -28,19 +28,15 @@
 
 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}
                           };
@@ -73,7 +69,5 @@
 
     // 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");
+    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
@@ -84,19 +78,18 @@
     }
 
-    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;
+    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);
+                }
+            } 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);
             }
         }
@@ -109,14 +102,18 @@
         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
+
+    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 gime original image back (%d,%d %.2f vs %.2f)",
+                        col,row,pixel,imgRow[col]);
+                return 5;
+            }
+        }
+    }
+
     psImageFree(img);
     psImageFree(img2);
@@ -125,5 +122,5 @@
     return 0;
 }
-#if 0
+
 int testImageRealImaginary(void)
 {
@@ -131,4 +128,6 @@
     psImage* img2 = NULL;
     psImage* img3 = NULL;
+    unsigned int m = 128;
+    unsigned int n = 64;
 
     /*
@@ -139,9 +138,5 @@
 
     // 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);
-    }
+    GENIMAGE(img,m,n,C32, row + I * col);
 
     // 2. call psImageReal and psImageImaginary
@@ -169,16 +164,18 @@
 
     // 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;
+    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;
+            }
         }
     }
@@ -196,4 +193,6 @@
     psImage* img2 = NULL;
     psImage* img3 = NULL;
+    unsigned int m = 128;
+    unsigned int n = 64;
 
     /*
@@ -214,12 +213,6 @@
 
     // 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);
-    }
+    GENIMAGE(img,m,n,F32,row);
+    GENIMAGE(img2,m,n,F32,col);
 
     // 2. call psImageComplex
@@ -228,5 +221,5 @@
     // 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)",
+        psError(__func__,"Image Type from psImageComplex is not complex? (%d)",
                 img3->type.type);
         return 1;
@@ -235,16 +228,18 @@
     // 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);
+    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
@@ -253,17 +248,18 @@
     // 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.");
+        psError(__func__,"psImageComplex returned a image though input types mismatched.");
         return 3;
     }
 
     // 9. create two psf32 vectors of different sizes
-    img2 = psImageRecycle(img2,PS_TYPE_F32,200);
+    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. (actually, it isn't an error...)
-    if (img3->n != 100) {
-        psError(__func__,"psImageComplex returned a vector though input sizes mismatched.");
+    // 11. verify thet an appropriate error occurred.
+    if (img3 != NULL) {
+        psError(__func__,"psImageComplex returned a image though input sizes mismatched.");
         return 4;
     }
@@ -280,4 +276,6 @@
     psImage* img = NULL;
     psImage* img2 = NULL;
+    unsigned int m = 128;
+    unsigned int n = 64;
 
     /*
@@ -289,9 +287,5 @@
 
     // 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);
-    }
+    GENIMAGE(img,m,n,C32, row + I * col);
 
     // 2. call psImageConjugate
@@ -300,16 +294,19 @@
     // 3. verify result is psC32
     if (img2->type.type != PS_TYPE_C32) {
-        psError(__func__,"the psImageConjugate didn't return a C32 vector");
+        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 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;
-        };
+    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;
+            }
+        }
     }
 
@@ -324,5 +321,6 @@
     psImage* img = NULL;
     psImage* img2 = NULL;
-    psF32 val;
+    unsigned int m = 128;
+    unsigned int n = 64;
 
     /*
@@ -334,9 +332,5 @@
 
     // 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);
-    }
+    GENIMAGE(img,m,n,C32, row + I * col);
 
     // 2. call psImagePowerSpectrum
@@ -350,39 +344,22 @@
 
     // 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
+    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;
+            }
+        }
+    }
+
+    psImageFree(img);
+    psImageFree(img2);
+
+    return 0;
+}