Changeset 1983 for trunk/psLib/src/image/psImageFFT.c

Timestamp:

Oct 6, 2004, 11:31:30 AM (22 years ago)

Author:

desonia

Message:

work on FFT and convolution. Fourier-domain convolution is still
broken.

File:

• trunk/psLib/src/image/psImageFFT.c (modified) (6 diffs)

trunk/psLib/src/image/psImageFFT.c

@@ -5 +5 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-23 18:31:49 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-10-06 21:31:30 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21 +21 @@
 #include "psLogMsg.h"
 #include "psImageExtraction.h"
+#include "psImageIO.h"
 
 #include "psImageErrors.h"
@@ -41 +42 @@
     }
 
+    if ( ((direction & PS_FFT_FORWARD) != 0) ) {
+        if ((direction & PS_FFT_REVERSE) != 0) {
+            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
+                       PS_ERR_BAD_PARAMETER_VALUE, true,
+                       PS_ERRORTEXT_psImageFFT_FORWARD_REVERSE);
+            psFree(out);
+            return NULL;
+        }
+        if ((direction & PS_FFT_REAL_RESULT) != 0) {
+            psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
+                       PS_ERR_BAD_PARAMETER_VALUE, true,
+                       PS_ERRORTEXT_psImageFFT_REAL_FORWARD_NOTSUPPORTED);
+            psFree(out);
+            return NULL;
+        }
+    } else if ((direction & PS_FFT_REVERSE) == 0) {
+        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
+                   PS_ERR_BAD_PARAMETER_VALUE, true,
+                   PS_ERRORTEXT_psImageFFT_NO_DIRECTION_OPTION);
+        psFree(out);
+        return NULL;
+    }
+
     type = in->type.type;
-
-    if ((type != PS_TYPE_F32) && (type != PS_TYPE_C32)) {
-        char* typeStr;
-        PS_TYPE_NAME(typeStr,type);
-        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
-                   PS_ERR_BAD_PARAMETER_TYPE, true,
-                   PS_ERRORTEXT_psImageFFT_IMAGE_TYPE_UNSUPPORTED,
-                   typeStr);
-        psFree(out);
-        return NULL;
-    }
-
-    if (type != PS_TYPE_C32 && direction == PS_FFT_REVERSE) {
-        char* typeStr;
-        PS_TYPE_NAME(typeStr,type);
-        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
-                   PS_ERR_BAD_PARAMETER_TYPE, true,
-                   PS_ERRORTEXT_psImageFFT_REVERSE_NOT_COMPLEX,
-                   typeStr);
-        psFree(out);
-        return NULL;
-
-    }
-
-    if (type != PS_TYPE_F32 && direction == PS_FFT_FORWARD) {
-        char* typeStr;
-        PS_TYPE_NAME(typeStr,type);
-        psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
-                   PS_ERR_BAD_PARAMETER_TYPE, true,
-                   PS_ERRORTEXT_psImageFFT_FORWARD_NOT_REAL,
-                   typeStr);
-        psFree(out);
-        return NULL;
-    }
 
     /* make sure the system-level wisdom information is imported. */
@@ -86 +76 @@
     numCols = in->numCols;
 
-    // n.b. FFTW can perform a in-place transform at the same or faster than out-of-place.
+    // n.b. FFTW can perform a in-place transform at the same rate or faster than out-of-place.
+    int sign = ((direction & PS_FFT_FORWARD) != 0) ? FFTW_FORWARD : FFTW_BACKWARD;
     out = psImageCopy(out, in, PS_TYPE_C32);
     plan = fftwf_plan_dft_2d(numCols, numRows,
-                             (fftwf_complex *) out->data.C32[0],
-                             (fftwf_complex *) out->data.C32[0], direction, PS_FFTW_PLAN_RIGOR);
-
-    /* check if a plan exists now -- if not, it is a real problem */
+                             out->data.C32[0],
+                             out->data.C32[0],
+                             sign,
+                             PS_FFTW_PLAN_RIGOR);
+
+    /* check if a plan exists now -- if not, it is a real problem at this point */
     if (plan == NULL) {
         psErrorMsg(PS_ERRORNAME_DOMAIN "psImageFFT",
@@ -105 +98 @@
 
     fftwf_destroy_plan(plan);
+
+    if (direction & PS_FFT_REAL_RESULT) {
+        // n.b., we do this instead of using fftwf_plan_dft_c2r because that
+        // plan requires a half-image, which would require a image reordering
+        // that is not as simple as performing a normal complex transform and
+        // then taking the real part of the result.  If performance here
+        // becomes an issue, the use of fftwf_plan_dft_r2c should be considered
+        // as well as fftwf_plan_dft_c2r.
+        psImage* realOut = psImageCopy(NULL,out,PS_TYPE_F32);
+        psFree(out);
+        out = realOut;
+    }
 
     return out;
@@ -424 +429 @@
                 real = creal(inRow[col]);
                 imag = cimag(inRow[col]);
-                outRow[col] = real * real + imag * imag / numElementsSquared;
+                outRow[col] = (real * real + imag * imag) / numElementsSquared;
             }
         }