r7939: Initial working version of fftw3 conversion -- tested only with pjrec

diff --git a/Makefile.in b/Makefile.in
index c6a9269f67f52ed01ef95af6817ddd1a8ec991da..be18d3ce013aac57d6ff48c8a3fb0622d95866e5 100644 (file)
--- a/Makefile.in
+++ b/Makefile.in
@@ -162,12 +162,12 @@ wxconfig = @wxconfig@
 wxlibs = @wxlibs@
 @INCLUDED_GETOPT_LONG_FALSE@EXTRA_DIRS1 = 
 @INCLUDED_GETOPT_LONG_TRUE@EXTRA_DIRS1 = getopt
+@HAVE_SGP_FALSE@EXTRA_DIRS2 = 
 
 @HAVE_SGP_TRUE@EXTRA_DIRS2 = libctgraphics
-@HAVE_SGP_FALSE@EXTRA_DIRS2 = 
+@HAVE_WXWINDOWS_FALSE@EXTRA_DIRS3 = 
 
 @HAVE_WXWINDOWS_TRUE@EXTRA_DIRS3 = src
-@HAVE_WXWINDOWS_FALSE@EXTRA_DIRS3 = 
 
 SUBDIRS = man libctsupport libctsim  html cgi-bin include $(EXTRA_DIRS1) $(EXTRA_DIRS2) $(EXTRA_DIRS3) tools helical
 

diff --git a/configure b/configure
index 8390af63a39c5c965828e4625844b25b6ea3983f..0a32152d8c3206bba8103ec096bf844199d03d39 100755 (executable)
--- a/configure
+++ b/configure
@@ -1600,7 +1600,7 @@ fi
 
 # Define the identity of the package.
  PACKAGE=ctsim
- VERSION=4.2.8
+ VERSION=4.3.0
 
 
 cat >>confdefs.h <<_ACEOF
@@ -6454,13 +6454,13 @@ _ACEOF
 
 fi
 
-echo "$as_me:$LINENO: checking for fftw_one in -lfftw" >&5
-echo $ECHO_N "checking for fftw_one in -lfftw... $ECHO_C" >&6
-if test "${ac_cv_lib_fftw_fftw_one+set}" = set; then
+echo "$as_me:$LINENO: checking for fftw_malloc in -lfftw3" >&5
+echo $ECHO_N "checking for fftw_malloc in -lfftw3... $ECHO_C" >&6
+if test "${ac_cv_lib_fftw3_fftw_malloc+set}" = set; then
   echo $ECHO_N "(cached) $ECHO_C" >&6
 else
   ac_check_lib_save_LIBS=$LIBS
-LIBS="-lfftw  $LIBS"
+LIBS="-lfftw3  $LIBS"
 cat >conftest.$ac_ext <<_ACEOF
 #line $LINENO "configure"
 /* confdefs.h.  */
@@ -6475,11 +6475,11 @@ extern "C"
 #endif
 /* We use char because int might match the return type of a gcc2
    builtin and then its argument prototype would still apply.  */
-char fftw_one ();
+char fftw_malloc ();
 int
 main ()
 {
-fftw_one ();
+fftw_malloc ();
   ;
   return 0;
 }
@@ -6496,19 +6496,19 @@ if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
   ac_status=$?
   echo "$as_me:$LINENO: \$? = $ac_status" >&5
   (exit $ac_status); }; }; then
-  ac_cv_lib_fftw_fftw_one=yes
+  ac_cv_lib_fftw3_fftw_malloc=yes
 else
   echo "$as_me: failed program was:" >&5
 sed 's/^/| /' conftest.$ac_ext >&5
 
-ac_cv_lib_fftw_fftw_one=no
+ac_cv_lib_fftw3_fftw_malloc=no
 fi
 rm -f conftest.$ac_objext conftest$ac_exeext conftest.$ac_ext
 LIBS=$ac_check_lib_save_LIBS
 fi
-echo "$as_me:$LINENO: result: $ac_cv_lib_fftw_fftw_one" >&5
-echo "${ECHO_T}$ac_cv_lib_fftw_fftw_one" >&6
-if test $ac_cv_lib_fftw_fftw_one = yes; then
+echo "$as_me:$LINENO: result: $ac_cv_lib_fftw3_fftw_malloc" >&5
+echo "${ECHO_T}$ac_cv_lib_fftw3_fftw_malloc" >&6
+if test $ac_cv_lib_fftw3_fftw_malloc = yes; then
   fftw=true;
 cat >>confdefs.h <<\_ACEOF
 #define HAVE_FFTW 1
@@ -9765,7 +9765,7 @@ if test "$zlib" = "true" ; then
   ctlibs_tools="$ctlibs_tools -lz"
 fi
 if test "$fftw" = "true" ; then
-  ctlibs_tools="$ctlibs_tools -lrfftw -lfftw"
+  ctlibs_tools="$ctlibs_tools -lfftw3"
 fi
 if test "$ctn" = "true"; then
   ctlibs_tools="$ctlibs_tools -lctn"

diff --git a/configure.ac b/configure.ac
index 77ebd09dc3c727e3b2f47409f5ed8bf0176cbc91..2ce392d076a09af503a529cc2e2accbd03bc9c69 100644 (file)
--- a/configure.ac
+++ b/configure.ac
@@ -6,7 +6,7 @@ dnl CDPATH=
 AC_INIT
 AC_CONFIG_SRCDIR([src/ctsim.cpp])
 AM_MAINTAINER_MODE
-AM_INIT_AUTOMAKE(ctsim,4.2.8)
+AM_INIT_AUTOMAKE(ctsim,4.3.0)
 AM_CONFIG_HEADER(config.h)
 
 dnl Checks for programs.
@@ -43,7 +43,7 @@ AC_CHECK_LIB(readline, main, [readline=true;
 wxwin=false
 AC_CHECK_LIB(wx_gtk-2.4, main, [wxwin=true; wx_gtk=true; AC_DEFINE(HAVE_WXWINDOWS,1,[wxwindows library])])
 AC_CHECK_LIB(wx_mac-2.4, main, [wxwin=true; wx_mac=true; AC_DEFINE(HAVE_WXWINDOWS,1,[wxwindows library])])
-AC_CHECK_LIB(fftw, fftw_one, [fftw=true; AC_DEFINE(HAVE_FFTW,1,[FFTW library])], [fftw=false])
+AC_CHECK_LIB(fftw3, fftw_malloc, [fftw=true; AC_DEFINE(HAVE_FFTW,1,[FFTW library])], [fftw=false])
 AC_CHECK_LIB(GL, main, [libgl=true], [libgl=false], [-L/usr/X11R6/lib -L/usr/X11R6/lib64 -lXt -lXext])
 AC_CHECK_LIB(pthread, main, [pthread=true], [pthread=false])
 
@@ -378,7 +378,7 @@ if test "$zlib" = "true" ; then
   ctlibs_tools="$ctlibs_tools -lz"
 fi
 if test "$fftw" = "true" ; then
-  ctlibs_tools="$ctlibs_tools -lrfftw -lfftw"
+  ctlibs_tools="$ctlibs_tools -lfftw3"
 fi
 if test "$ctn" = "true"; then
   ctlibs_tools="$ctlibs_tools -lctn"

diff --git a/debian/changelog b/debian/changelog
index 3a33e7c83d14eacb41c2f589bbf5258738e4f144..bb44dcbe8ebbb471ba23ade67c5de8f5636f3154 100644 (file)
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+ctsim (4.3.0-1) unstable; urgency=low
+
+  * Port from FFTW2 to incompatible FFTW3 library
+
+ -- Kevin M. Rosenberg <kmr@debian.org>  Sat,  4 Oct 2003 20:03:31 -0600
+
 ctsim (4.2.8-1) unstable; urgency=low
 
   * configure.ac: Support for compiling on AMD64

diff --git a/include/ct.h b/include/ct.h
index 0608d2f32e802938f178c65d024e38c3e6715bb1..9ee346bf018c1385518893e9f844a17d7443884b 100644 (file)
--- a/include/ct.h
+++ b/include/ct.h
@@ -153,8 +153,7 @@ extern "C" {
 
 
 #ifdef HAVE_FFTW
-#include <rfftw.h>
-#include <fftw.h>
+#include <fftw3.h>
 #define HAVE_FFT 1
 #endif
 

diff --git a/include/filter.h b/include/filter.h
index 48d13a0c03dd2a945e0aafd511f96fadb6ca8a7f..5974e9488cd0894bf1ab60b44daae02e8ed630e2 100644 (file)
--- a/include/filter.h
+++ b/include/filter.h
@@ -33,8 +33,7 @@
 #include "config.h"
 #endif
 #ifdef HAVE_FFTW
-#include <fftw.h>
-#include <rfftw.h>
+#include <fftw3.h>
 #endif
 
 

diff --git a/include/fourier.h b/include/fourier.h
index e1716f88fb9deef37de60be5ea9729a756eafce6..5215a3a25928ad645869be3c2eeadba9e9b2d875 100644 (file)
--- a/include/fourier.h
+++ b/include/fourier.h
@@ -26,6 +26,9 @@
 ******************************************************************************/
 
 #include <complex>
+#ifdef HAVE_FFTW
+#include <fftw3.h>
+#endif
 
 class ImageFile;
 
@@ -34,6 +37,11 @@ public:
     static void shuffleFourierToNaturalOrder (ImageFile& im);
     static void shuffleNaturalToFourierOrder (ImageFile& im);
 
+#ifdef HAVE_FFTW
+    static void shuffleFourierToNaturalOrder (fftw_complex* pc, const int n);
+    static void shuffleNaturalToFourierOrder (fftw_complex* pc, const int n);
+#endif
+    
 // Odd Number of Points
 //   Natural Frequency Order: -(n-1)/2...-1,0,1...(n-1)/2
 //   Fourier Frequency Order: 0, 1..(n-1)/2,-(n-1)/2...-1

diff --git a/include/procsignal.h b/include/procsignal.h
index 2c37de89ddc43ef17909606e599ffb50fca8dd27..14b8f6f85390139007efab727ca8e5ba73b40dfe 100644 (file)
--- a/include/procsignal.h
+++ b/include/procsignal.h
@@ -33,8 +33,7 @@
 #include "config.h"
 #endif
 #ifdef HAVE_FFTW
-#include <fftw.h>
-#include <rfftw.h>
+#include <fftw3.h>
 #endif
 
 #include <complex>
@@ -152,9 +151,9 @@ class ProcessSignal {
     static const int s_iFilterGenerationCount;
 
 #ifdef HAVE_FFTW
-    fftw_real* m_adRealFftInput, *m_adRealFftSignal;
-    rfftw_plan m_realPlanForward, m_realPlanBackward;
-    fftw_complex* m_adComplexFftInput, *m_adComplexFftSignal;
+    double *m_adRealFftInput, *m_adRealFftOutput, *m_adRealFftSignal, *m_adRealFftBackwardOutput;
+    fftw_plan m_realPlanForward, m_realPlanBackward;
+    fftw_complex *m_adComplexFftInput, *m_adComplexFftOutput, *m_adComplexFftSignal, *m_adComplexFftBackwardOutput;
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #endif
 

diff --git a/libctsim/fourier.cpp b/libctsim/fourier.cpp
index 83922d6c6af0530b89ab91100807ff952fea7f54..a3a62d984ef80966bfa847445be2597ecc78be9f 100644 (file)
--- a/libctsim/fourier.cpp
+++ b/libctsim/fourier.cpp
@@ -98,3 +98,84 @@ Fourier::shuffleNaturalToFourierOrder (ImageFile& im)
   }
   delete [] pRow;
 }
+
+#ifdef HAVE_FFTW
+void Fourier::shuffleNaturalToFourierOrder (fftw_complex* pVector, const int n) 
+{
+  fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
+  int i;
+  
+  if (isOdd(n)) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pTemp[0][0] = pVector[iHalfN][0];
+    pTemp[0][1] = pVector[iHalfN][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + 1][0] = pVector[i + 1 + iHalfN][0];
+      pTemp[i + 1][1] = pVector[i + 1 + iHalfN][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + iHalfN + 1][0] = pVector[i][0];
+      pTemp[i + iHalfN + 1][1] = pVector[i][1];
+    }
+  } else {     // Even
+    int iHalfN = n / 2;
+    pTemp[0][0] = pVector[iHalfN][0];
+    pTemp[0][1] = pVector[iHalfN][1];
+    for (i = 0; i < iHalfN - 1; i++) {
+      pTemp[i + 1][0] = pVector[i + iHalfN + 1][0];
+      pTemp[i + 1][1] = pVector[i + iHalfN + 1][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + iHalfN][0] = pVector[i][0];
+      pTemp[i + iHalfN][1] = pVector[i][1];
+    }
+  }
+  
+  for (i = 0; i < n; i++) {
+    pVector[i][0] = pTemp[i][0];
+    pVector[i][1] = pTemp[i][1];
+  }
+  fftw_free(pTemp);
+}
+
+void Fourier::shuffleFourierToNaturalOrder (fftw_complex* pVector, const int n)
+{
+  fftw_complex* pTemp = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * n));
+  int i;
+  if (isOdd(n)) { // Odd
+    int iHalfN = (n - 1) / 2;
+    
+    pTemp[iHalfN][0] = pVector[0][0];
+    pTemp[iHalfN][1] = pVector[0][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i + 1 + iHalfN][0] = pVector[i + 1][0];
+      pTemp[i + 1 + iHalfN][1] = pVector[i + 1][1];
+    }
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i][0] = pVector[i + iHalfN + 1][0];
+      pTemp[i][1] = pVector[i + iHalfN + 1][1];
+    }
+  } else {     // Even
+    int iHalfN = n / 2;
+    pTemp[iHalfN][0] = pVector[0][0];
+    pTemp[iHalfN][1] = pVector[0][1];
+    for (i = 0; i < iHalfN; i++) {
+      pTemp[i][0] = pVector[i + iHalfN][0];
+      pTemp[i][1] = pVector[i + iHalfN][1];
+    }
+    for (i = 0; i < iHalfN - 1; i++) {
+      pTemp[i + iHalfN + 1][0] = pVector[i+1][0];
+      pTemp[i + iHalfN + 1][1] = pVector[i+1][1];
+    }
+  }
+  
+  for (i = 0; i < n; i++) {
+    pVector[i][0] = pTemp[i][0];
+    pVector[i][1] = pTemp[i][1];
+  }
+  
+  fftw_free(pTemp);
+}
+#endif
+

diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp
index ef8b1e4bc7d262f4834aa549e2285b409a9a5568..ad13545d10bf918342bc12d28ebbaf3107e3592f 100644 (file)
--- a/libctsim/imagefile.cpp
+++ b/libctsim/imagefile.cpp
@@ -725,7 +725,7 @@ ImageFile::fft (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+  fftw_complex* in = static_cast<fftw_complex*> (fftw_malloc (sizeof(fftw_complex) * m_nx * m_ny));
   
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
@@ -734,17 +734,17 @@ ImageFile::fft (ImageFile& result) const
   unsigned int iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      in[iArray].re = vReal[ix][iy];
+      in[iArray][0] = vReal[ix][iy];
       if (isComplex())
-        in[iArray].im = vImag[ix][iy];
+        in[iArray][1] = vImag[ix][iy];
       else
-        in[iArray].im = 0;
+        in[iArray][1] = 0;
       iArray++;
     }
   }
   
-  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);  
-  fftwnd_one (plan, in, NULL);
+  fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);  
+  fftw_execute (plan);
   
   ImageFileArray vRealResult = result.getArray();
   ImageFileArray vImagResult = result.getImaginaryArray();
@@ -752,13 +752,13 @@ ImageFile::fft (ImageFile& result) const
   unsigned int iScale = m_nx * m_ny;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      vRealResult[ix][iy] = in[iArray].re / iScale;
-      vImagResult[ix][iy] = in[iArray].im / iScale;
+      vRealResult[ix][iy] = in[iArray][0] / iScale;
+      vImagResult[ix][iy] = in[iArray][1] / iScale;
       iArray++;
     }
   }
   delete in;
-  fftwnd_destroy_plan (plan);
+  fftw_destroy_plan (plan);
   
   Fourier::shuffleFourierToNaturalOrder (result);
   
@@ -796,32 +796,31 @@ ImageFile::ifft (ImageFile& result) const
   
   Fourier::shuffleNaturalToFourierOrder (result);
   
-  fftw_complex* in = new fftw_complex [m_nx * m_ny];
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx * m_ny));
   
   unsigned int iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      in[iArray].re = vRealResult[ix][iy];
-      in[iArray].im = vImagResult[ix][iy];
+      in[iArray][0] = vRealResult[ix][iy];
+      in[iArray][1] = vImagResult[ix][iy];
       iArray++;
     }
   }
   
-  fftwnd_plan plan = fftw2d_create_plan (m_nx, m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_plan plan = fftw_plan_dft_2d (m_nx, m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   
-  fftwnd_one (plan, in, NULL);
+  fftw_execute (plan);
   
   iArray = 0;
   for (ix = 0; ix < m_nx; ix++) {
     for (iy = 0; iy < m_ny; iy++) {
-      vRealResult[ix][iy] = in[iArray].re;
-      vImagResult[ix][iy] = in[iArray].im;
+      vRealResult[ix][iy] = in[iArray][0];
+      vImagResult[ix][iy] = in[iArray][1];
       iArray++;
     }
   }
-  fftwnd_destroy_plan (plan);
-  
-  delete in;
+  fftw_destroy_plan (plan);
+  fftw_free(in);
   
   return true;
 }
@@ -842,24 +841,24 @@ ImageFile::fftRows (ImageFile& result) const
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+  fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
 
-  fftw_complex* in = new fftw_complex [m_nx];
   std::complex<double>* pcRow = new std::complex<double> [m_nx];  
   for (unsigned int iy = 0; iy < m_ny; iy++) {
     unsigned int ix;
     for (ix = 0; ix < m_nx; ix++) {
-      in[ix].re = vReal[ix][iy];
+      in[ix][0] = vReal[ix][iy];
       if (isComplex())
-        in[ix].im = vImag[ix][iy];
+        in[ix][1] = vImag[ix][iy];
       else
-        in[ix].im = 0;
+        in[ix][1] = 0;
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (ix = 0; ix < m_nx; ix++)
-      pcRow[ix] = std::complex<double>(in[ix].re, in[ix].im);
+      pcRow[ix] = std::complex<double>(in[ix][0], in[ix][1]);
     
     Fourier::shuffleFourierToNaturalOrder (pcRow, m_nx);
     for (ix = 0; ix < m_nx; ix++) {
@@ -870,7 +869,7 @@ ImageFile::fftRows (ImageFile& result) const
   delete [] pcRow;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }     
@@ -888,12 +887,11 @@ ImageFile::ifftRows (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_nx];
-  
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
-  
-  fftw_plan plan = fftw_create_plan (m_nx, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+   
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_nx));
+  fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   std::complex<double>* pcRow = new std::complex<double> [m_nx];
   
   unsigned int ix, iy;
@@ -909,21 +907,21 @@ ImageFile::ifftRows (ImageFile& result) const
     Fourier::shuffleNaturalToFourierOrder (pcRow, m_nx);
     
     for (ix = 0; ix < m_nx; ix++) {
-      in[ix].re = pcRow[ix].real();
-      in[ix].im = pcRow[ix].imag();
+      in[ix][0] = pcRow[ix].real();
+      in[ix][1] = pcRow[ix].imag();
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (ix = 0; ix < m_nx; ix++) {
-      vReal[ix][iy] = in[ix].re;
-      vImag[ix][iy] = in[ix].im;
+      vReal[ix][iy] = in[ix][0];
+      vImag[ix][iy] = in[ix][1];
     }
   }
   delete [] pcRow;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }
@@ -944,24 +942,24 @@ ImageFile::fftCols (ImageFile& result) const
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_ny, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+  fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE);
 
   std::complex<double>* pcCol = new std::complex<double> [m_ny];  
-  fftw_complex* in = new fftw_complex [m_ny];
   for (unsigned int ix = 0; ix < m_nx; ix++) {
     unsigned int iy;
     for (iy = 0; iy < m_ny; iy++) {
-      in[iy].re = vReal[ix][iy];
+      in[iy][0] = vReal[ix][iy];
       if (isComplex())
-        in[iy].im = vImag[ix][iy];
+        in[iy][1] = vImag[ix][iy];
       else
-        in[iy].im = 0;
+        in[iy][1] = 0;
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (iy = 0; iy < m_ny; iy++)
-      pcCol[iy] = std::complex<double>(in[iy].re, in[iy].im);
+      pcCol[iy] = std::complex<double>(in[iy][0], in[iy][1]);
     
     Fourier::shuffleFourierToNaturalOrder (pcCol, m_ny);
     for (iy = 0; iy < m_ny; iy++) {
@@ -972,7 +970,7 @@ ImageFile::fftCols (ImageFile& result) const
   delete [] pcCol;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }
@@ -990,12 +988,11 @@ ImageFile::ifftCols (ImageFile& result) const
       return false;
   }
   
-  fftw_complex* in = new fftw_complex [m_ny];
-  
   ImageFileArrayConst vReal = getArray();
   ImageFileArrayConst vImag = getImaginaryArray();
   
-  fftw_plan plan = fftw_create_plan (m_ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* in = static_cast<fftw_complex*>(fftw_malloc(sizeof(fftw_complex) * m_ny));
+  fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE);
   std::complex<double>* pcCol = new std::complex<double> [m_ny];
   
   unsigned int ix, iy;
@@ -1011,21 +1008,21 @@ ImageFile::ifftCols (ImageFile& result) const
     Fourier::shuffleNaturalToFourierOrder (pcCol, m_ny);
     
     for (iy = 0; iy < m_ny; iy++) {
-      in[iy].re = pcCol[iy].real();
-      in[iy].im = pcCol[iy].imag();
+      in[iy][0] = pcCol[iy].real();
+      in[iy][1] = pcCol[iy].imag();
     }
     
-    fftw_one (plan, in, NULL);
+    fftw_execute (plan);
     
     for (iy = 0; iy < m_ny; iy++) {
-      vReal[ix][iy] = in[iy].re;
-      vImag[ix][iy] = in[iy].im;
+      vReal[ix][iy] = in[iy][0];
+      vImag[ix][iy] = in[iy][1];
     }
   }
   delete [] pcCol;
   
   fftw_destroy_plan (plan);
-  delete in;
+  fftw_free(in);
   
   return true;
 }

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index bdb9fa00626813ca37509601bd7eb7698d10c470..2206ae1838cf41ebd3f2e13eb34353fdb060670e 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -418,21 +418,26 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   }
   
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
-    m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    m_adRealFftInput = new fftw_real [ m_nFilterPoints ];
-    m_adRealFftSignal = new fftw_real [ m_nOutputPoints ];
+    m_adRealFftInput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+    m_adRealFftOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nFilterPoints));
+    m_realPlanForward = fftw_plan_r2r_1d (m_nFilterPoints, m_adRealFftInput, m_adRealFftOutput, FFTW_R2HC, FFTW_ESTIMATE);
+    m_adRealFftSignal = static_cast<double*>(fftw_malloc (sizeof(double) *  m_nOutputPoints));
+    m_adRealFftBackwardOutput = static_cast<double*>(fftw_malloc (sizeof(double) * m_nOutputPoints));
+    m_realPlanBackward = fftw_plan_r2r_1d (m_nOutputPoints, m_adRealFftSignal, m_adRealFftBackwardOutput, FFTW_HC2R, FFTW_ESTIMATE);
     for (i = 0; i < m_nFilterPoints; i++) 
       m_adRealFftInput[i] = 0;
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
-    m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD,  FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD,  FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    m_adComplexFftInput = new fftw_complex [ m_nFilterPoints ];
-    m_adComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
+    m_adComplexFftInput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+    m_adComplexFftOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nFilterPoints));
+    m_complexPlanForward = fftw_plan_dft_1d (m_nFilterPoints, m_adComplexFftInput, m_adComplexFftOutput, FFTW_FORWARD,  FFTW_ESTIMATE);
+    m_adComplexFftSignal = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+    m_adComplexFftBackwardOutput = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * m_nOutputPoints));
+    m_complexPlanBackward = fftw_plan_dft_1d (m_nOutputPoints, m_adComplexFftSignal, m_adComplexFftBackwardOutput, FFTW_BACKWARD,  FFTW_ESTIMATE);
+
     for (i = 0; i < m_nFilterPoints; i++) 
-      m_adComplexFftInput[i].re = m_adComplexFftInput[i].im = 0;
+      m_adComplexFftInput[i][0] = m_adComplexFftInput[i][1] = 0;
     for (i = 0; i < m_nOutputPoints; i++) 
-      m_adComplexFftSignal[i].re = m_adComplexFftSignal[i].im = 0;
+      m_adComplexFftSignal[i][0] = m_adComplexFftSignal[i][1] = 0;
   }
 #endif
   
@@ -448,14 +453,18 @@ ProcessSignal::~ProcessSignal (void)
   if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     fftw_destroy_plan(m_complexPlanForward);
     fftw_destroy_plan(m_complexPlanBackward);
-    delete [] m_adComplexFftInput;
-    delete [] m_adComplexFftSignal;
+    fftw_free (m_adComplexFftInput);
+    fftw_free (m_adComplexFftOutput);
+    fftw_free (m_adComplexFftSignal);
+    fftw_free (m_adComplexFftBackwardOutput);
   }
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
-    rfftw_destroy_plan(m_realPlanForward);
-    rfftw_destroy_plan(m_realPlanBackward);
-    delete [] m_adRealFftInput;
-    delete [] m_adRealFftSignal;
+    fftw_destroy_plan(m_realPlanForward);
+    fftw_destroy_plan(m_realPlanBackward);
+    fftw_free (m_adRealFftInput);
+    fftw_free (m_adRealFftOutput);
+    fftw_free (m_adRealFftSignal);
+    fftw_free (m_adRealFftBackwardOutput);
   }
 #endif
 }
@@ -595,35 +604,27 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     for (i = 0; i < m_nSignalPoints; i++)
       m_adRealFftInput[i] = input[i];
     
-    fftw_real* fftOutput = new fftw_real [ m_nFilterPoints ];
-    rfftw_one (m_realPlanForward, m_adRealFftInput, fftOutput);
+    fftw_execute (m_realPlanForward);
     for (i = 0; i < m_nFilterPoints; i++)
-      m_adRealFftSignal[i] = m_adFilter[i] * fftOutput[i];
-    delete [] fftOutput;
+      m_adRealFftSignal[i] = m_adFilter[i] * m_adRealFftOutput[i];
     for (i = m_nFilterPoints; i < m_nOutputPoints; i++)
             m_adRealFftSignal[i] = 0;
     
-    fftw_real* ifftOutput = new fftw_real [ m_nOutputPoints ];
-    rfftw_one (m_realPlanBackward, m_adRealFftSignal, ifftOutput);
+    fftw_execute (m_realPlanBackward);
     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++)
-      output[i] = ifftOutput[i];
-    delete [] ifftOutput;
+      output[i] = m_adRealFftBackwardOutput[i];
   } else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
     for (i = 0; i < m_nSignalPoints; i++)
-      m_adComplexFftInput[i].re = input[i];
+      m_adComplexFftInput[i][0] = input[i];
     
-    fftw_complex* fftOutput = new fftw_complex [ m_nFilterPoints ];
-    fftw_one (m_complexPlanForward, m_adComplexFftInput, fftOutput);
+    fftw_execute (m_complexPlanForward);
     for (i = 0; i < m_nFilterPoints; i++) {
-      m_adComplexFftSignal[i].re = m_adFilter[i] * fftOutput[i].re;
-      m_adComplexFftSignal[i].im = m_adFilter[i] * fftOutput[i].im;
+      m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
+      m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];
     }
-    delete [] fftOutput;
-    fftw_complex* ifftOutput = new fftw_complex [ m_nOutputPoints ];
-    fftw_one (m_complexPlanBackward, m_adComplexFftSignal, ifftOutput);
+    fftw_execute (m_complexPlanBackward);
     for (i = 0; i < m_nSignalPoints * m_iPreinterpolationFactor; i++) 
-      output[i] = ifftOutput[i].re;
-    delete [] ifftOutput;
+      output[i] = m_adComplexFftBackwardOutput[i][0];
   }
 #endif
   delete input;

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 871f62e6e52a4543df9262355872d3b22957f392..8674647146a8c4f64f160c4323ef7be974b304d9 100644 (file)
--- a/libctsim/projections.cpp
+++ b/libctsim/projections.cpp
@@ -1019,12 +1019,12 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
 
   int iInterpDet = static_cast<int>(static_cast<double>(sqrt(nx*nx+ny*ny)));
   int iNumInterpDetWithZeros = ProcessSignal::addZeropadFactor (iInterpDet, iZeropad);
-  double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.5);
+  double dProjScale = iInterpDet / (pProj->viewDiameter() * 0.05);
   double dZeropadRatio = static_cast<double>(iNumInterpDetWithZeros) / static_cast<double>(iInterpDet);
 
-  fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
+  fftw_complex* pcIn = static_cast<fftw_complex*> (fftw_malloc (sizeof(fftw_complex) * iNumInterpDetWithZeros));
+  fftw_plan plan = fftw_plan_dft_1d (iNumInterpDetWithZeros, pcIn, pcIn, FFTW_FORWARD, FFTW_ESTIMATE);
 
-  fftw_complex* pcIn = new fftw_complex [iNumInterpDetWithZeros];
   std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
   //double dInterpScale = (pProj->m_nDet-1) / static_cast<double>(iInterpDet-1);
   double dInterpScale = pProj->m_nDet / static_cast<double>(iInterpDet);
@@ -1040,28 +1040,30 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad
     LinearInterpolator<DetectorValue> projInterp (detval, pProj->m_nDet);
     for (int iDet = 0; iDet < iInterpDet; iDet++) {
       double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale;
-      pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dProjScale;
-      pcIn[iDet].im = 0;
+      pcIn[iDet][0] = projInterp.interpolate (dInterpPos) * dProjScale;
+      pcIn[iDet][1] = 0;
     }
 
     Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet);
     if (iZerosAdded > 0) {
-       for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--)
-        pcIn[iDet1+iZerosAdded] = pcIn[iDet1];
+      for (int iDet1 = iInterpDet -1; iDet1 >= iMidPoint; iDet1--) {
+        pcIn[iDet1+iZerosAdded][0] = pcIn[iDet1][0];
+        pcIn[iDet1+iZerosAdded][1] = pcIn[iDet1][1];
+      }
       for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++) 
-        pcIn[iDet2].re = pcIn[iDet2].im = 0;
+        pcIn[iDet2][0] = pcIn[iDet2][1] = 0;
     }
 
-    fftw_one (plan, pcIn, NULL);
+    fftw_execute (plan);
 
     ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
     for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
-      ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); 
+      ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD][0] * dFFTScale, pcIn[iD][1] * dFFTScale); 
     }
 
     Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], iNumInterpDetWithZeros);
   }
-  delete [] pcIn;
+  fftw_free(pcIn) ;
 
   fftw_destroy_plan (plan);  
   

diff --git a/src/views.cpp b/src/views.cpp
index 0bf911a79fe7d52582b4385260fc3fa5b66aac2e..1fe4bbfbbc2c71facb551f464b64bc7dec8a4b36 100644 (file)
--- a/src/views.cpp
+++ b/src/views.cpp
@@ -1560,19 +1560,19 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
   int ny = rIF.ny();
   
   if (v != NULL && yCursor < ny) {
-    fftw_complex* pcIn = new fftw_complex [nx];
+    fftw_complex* pcIn = static_cast<fftw_complex*>(fftw_malloc (sizeof(fftw_complex) * nx));
     
     int i;
     for (i = 0; i < nx; i++) {
-      pcIn[i].re = v[i][yCursor];
+      pcIn[i][0] = v[i][yCursor];
       if (rIF.isComplex())
-        pcIn[i].im = vImag[i][yCursor];
+        pcIn[i][1] = vImag[i][yCursor];
       else
-        pcIn[i].im = 0;
+        pcIn[i][1] = 0;
     }
     
-    fftw_plan plan = fftw_create_plan (nx, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    fftw_one (plan, pcIn, NULL);
+    fftw_plan plan = fftw_plan_dft_1d (nx, pcIn, pcIn, FFTW_FORWARD, FFTW_ESTIMATE);
+    fftw_execute (plan);
     fftw_destroy_plan (plan);
     
     double* pX = new double [nx];
@@ -1581,9 +1581,9 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
     double* pYMag = new double [nx];
     for (i = 0; i < nx; i++) {
       pX[i] = i;
-      pYReal[i] = pcIn[i].re / nx;
-      pYImag[i] = pcIn[i].im / nx;
-      pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
+      pYReal[i] = pcIn[i][0] / nx;
+      pYImag[i] = pcIn[i][1] / nx;
+      pYMag[i] = ::sqrt (pcIn[i][0] * pcIn[i][0] + pcIn[i][1] * pcIn[i][1]);
     }
     Fourier::shuffleFourierToNaturalOrder (pYReal, nx);
     Fourier::shuffleFourierToNaturalOrder (pYImag, nx);
@@ -1628,7 +1628,7 @@ ImageFileView::OnPlotFFTRow (wxCommandEvent& event)
     delete pYReal;
     delete pYImag;
     delete pYMag;
-    delete [] pcIn;
+    fftw_free(pcIn);
     
     if (theApp->getAskDeleteNewDocs())
       pPlotDoc->Modify (true);
@@ -1662,7 +1662,7 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
       pdTemp[i] = v[xCursor][i];
     Fourier::shuffleNaturalToFourierOrder (pdTemp, ny);
     for (i = 0; i < ny; i++) 
-      pcIn[i].re = pdTemp[i];
+      pcIn[i][0] = pdTemp[i];
     
     for (i = 0; i < ny; i++) {
       if (rIF.isComplex())
@@ -1672,10 +1672,10 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
     }
     Fourier::shuffleNaturalToFourierOrder (pdTemp, ny);
     for (i = 0; i < ny; i++)
-      pcIn[i].im = pdTemp[i];
+      pcIn[i][1] = pdTemp[i];
     
-    fftw_plan plan = fftw_create_plan (ny, FFTW_BACKWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM);
-    fftw_one (plan, pcIn, NULL);
+    fftw_plan plan = fftw_plan_dft_1d (ny, pcIn, pcIn, FFTW_BACKWARD, FFTW_ESTIMATE);
+    fftw_execute (plan);
     fftw_destroy_plan (plan);
     
     double* pX = new double [ny];
@@ -1684,9 +1684,9 @@ ImageFileView::OnPlotFFTCol (wxCommandEvent& event)
     double* pYMag = new double [ny];
     for (i = 0; i < ny; i++) {
       pX[i] = i;
-      pYReal[i] = pcIn[i].re / ny;
-      pYImag[i] = pcIn[i].im / ny;
-      pYMag[i] = ::sqrt (pcIn[i].re * pcIn[i].re + pcIn[i].im * pcIn[i].im);
+      pYReal[i] = pcIn[i][0] / ny;
+      pYImag[i] = pcIn[i][1] / ny;
+      pYMag[i] = ::sqrt (pcIn[i][0] * pcIn[i][0] + pcIn[i][1] * pcIn[i][1]);
     }
     
     PlotFileDocument* pPlotDoc = theApp->newPlotDoc();

diff --git a/tools/Makefile.in b/tools/Makefile.in
index 3c6a56823b0eb373910f8e3b56494818e98a4ab7..fd1f0d69fe2ea755acbc6f131ec3e9a44980bb36 100644 (file)
--- a/tools/Makefile.in
+++ b/tools/Makefile.in
@@ -163,9 +163,9 @@ wxlibs = @wxlibs@
 bin_PROGRAMS = @lamprograms@  ctsimtext 
 EXTRA_PROGRAMS = ctsimtext-lam
 INCLUDES = @my_includes@
+@HAVE_SGP_FALSE@SOURCE_DEPEND = ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
 
 @HAVE_SGP_TRUE@SOURCE_DEPEND = ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a
-@HAVE_SGP_FALSE@SOURCE_DEPEND = ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
 
 ctsimtext_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp phm2helix.cpp pjHinterp.cpp linogram.cpp
 ctsimtext_LDADD = @ctlibs@