From d42d3d062dd1aca92b5a2552a1f474aab0bee610 Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Sun, 5 Oct 2003 05:47:22 +0000 Subject: [PATCH] r7939: Initial working version of fftw3 conversion -- tested only with pjrec --- Makefile.in | 4 +- configure | 26 +++++----- configure.ac | 6 +-- debian/changelog | 6 +++ include/ct.h | 3 +- include/filter.h | 3 +- include/fourier.h | 8 +++ include/procsignal.h | 9 ++-- libctsim/fourier.cpp | 81 ++++++++++++++++++++++++++++++ libctsim/imagefile.cpp | 105 +++++++++++++++++++-------------------- libctsim/procsignal.cpp | 69 ++++++++++++------------- libctsim/projections.cpp | 24 +++++---- src/views.cpp | 34 ++++++------- tools/Makefile.in | 2 +- 14 files changed, 236 insertions(+), 144 deletions(-) diff --git a/Makefile.in b/Makefile.in index c6a9269..be18d3c 100644 --- 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 8390af6..0a32152 100755 --- 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 77ebd09..2ce392d 100644 --- 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 3a33e7c..bb44dcb 100644 --- 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 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 0608d2f..9ee346b 100644 --- a/include/ct.h +++ b/include/ct.h @@ -153,8 +153,7 @@ extern "C" { #ifdef HAVE_FFTW -#include -#include +#include #define HAVE_FFT 1 #endif diff --git a/include/filter.h b/include/filter.h index 48d13a0..5974e94 100644 --- a/include/filter.h +++ b/include/filter.h @@ -33,8 +33,7 @@ #include "config.h" #endif #ifdef HAVE_FFTW -#include -#include +#include #endif diff --git a/include/fourier.h b/include/fourier.h index e1716f8..5215a3a 100644 --- a/include/fourier.h +++ b/include/fourier.h @@ -26,6 +26,9 @@ ******************************************************************************/ #include +#ifdef HAVE_FFTW +#include +#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 2c37de8..14b8f6f 100644 --- a/include/procsignal.h +++ b/include/procsignal.h @@ -33,8 +33,7 @@ #include "config.h" #endif #ifdef HAVE_FFTW -#include -#include +#include #endif #include @@ -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 83922d6..a3a62d9 100644 --- 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_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_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 ef8b1e4..ad13545 100644 --- 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_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_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_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* pcRow = new std::complex [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(in[ix].re, in[ix].im); + pcRow[ix] = std::complex(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_malloc(sizeof(fftw_complex) * m_nx)); + fftw_plan plan = fftw_plan_dft_1d (m_nx, in, in, FFTW_BACKWARD, FFTW_ESTIMATE); std::complex* pcRow = new std::complex [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_malloc(sizeof(fftw_complex) * m_ny)); + fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_FORWARD, FFTW_ESTIMATE); std::complex* pcCol = new std::complex [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(in[iy].re, in[iy].im); + pcCol[iy] = std::complex(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_malloc(sizeof(fftw_complex) * m_ny)); + fftw_plan plan = fftw_plan_dft_1d (m_ny, in, in, FFTW_BACKWARD, FFTW_ESTIMATE); std::complex* pcCol = new std::complex [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 bdb9fa0..2206ae1 100644 --- 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(fftw_malloc (sizeof(double) * m_nFilterPoints)); + m_adRealFftOutput = static_cast(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(fftw_malloc (sizeof(double) * m_nOutputPoints)); + m_adRealFftBackwardOutput = static_cast(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_malloc (sizeof(fftw_complex) * m_nFilterPoints)); + m_adComplexFftOutput = static_cast(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_malloc (sizeof(fftw_complex) * m_nOutputPoints)); + m_adComplexFftBackwardOutput = static_cast(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 871f62e..8674647 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -1019,12 +1019,12 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad int iInterpDet = static_cast(static_cast(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(iNumInterpDetWithZeros) / static_cast(iInterpDet); - fftw_plan plan = fftw_create_plan (iNumInterpDetWithZeros, FFTW_FORWARD, FFTW_IN_PLACE | FFTW_ESTIMATE | FFTW_USE_WISDOM); + fftw_complex* pcIn = static_cast (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** ppcDetValue = new std::complex* [pProj->m_nView]; //double dInterpScale = (pProj->m_nDet-1) / static_cast(iInterpDet-1); double dInterpScale = pProj->m_nDet / static_cast(iInterpDet); @@ -1040,28 +1040,30 @@ Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad LinearInterpolator 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 [iNumInterpDetWithZeros]; for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) { - ppcDetValue[iView][iD] = std::complex (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale); + ppcDetValue[iView][iD] = std::complex (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 0bf911a..1fe4bbf 100644 --- 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_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 3c6a568..fd1f0d6 100644 --- 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@ -- 2.34.1