From a8ba12a8c971de1d8cb3ef1c3a7d2d9fcf45affa Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Tue, 4 Jul 2000 18:33:35 +0000 Subject: [PATCH] r129: *** empty log message *** --- ChangeLog | 4 +- acconfig.h | 2 + cgi-bin/ctsim.cgi.in | 3 +- cgi-bin/ctsim.conf.in | 2 + config.h.in | 2 + configure | 260 +++++++++++++++++++++--------------- configure.in | 5 + html/simulate.html.in | 5 + include/ct.h | 6 +- include/filter.h | 12 +- libctgraphics/ezplot.cpp | 6 +- libctgraphics/ezsupport.cpp | 8 +- libctsim/filter.cpp | 228 ++++++++++++++++++++++++------- libctsim/projections.cpp | 25 ++-- src/pjrec.cpp | 11 +- 15 files changed, 395 insertions(+), 184 deletions(-) diff --git a/ChangeLog b/ChangeLog index 023c1c7..17dfc1f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,9 +1,9 @@ -1.9.9 - 6/29/00 +1.9.9 - 7/04/00 Fixed const issue with ImageFileArray Fixed Array2dFile::labelsCopy() Added copy constructor and assignment for Array2dFileLabel class Added Timer to if-2.cpp and ifinfo.cpp - Added beginning of frequency (FFT) based filter to SignalFilter + Added beginning of frequency-based (DFT & FFT) filter to SignalFilter 1.9.8 - 6/27/2000 Rewrote Array2dFile class to be non-templated diff --git a/acconfig.h b/acconfig.h index 69da21f..8d68fc0 100644 --- a/acconfig.h +++ b/acconfig.h @@ -62,3 +62,5 @@ #undef HAVE_X11 #undef HAVE_WXWINDOWS + +#undef HAVE_FFTW diff --git a/cgi-bin/ctsim.cgi.in b/cgi-bin/ctsim.cgi.in index 06d16e8..54a38a6 100755 --- a/cgi-bin/ctsim.cgi.in +++ b/cgi-bin/ctsim.cgi.in @@ -47,6 +47,7 @@ $error .= "Projection RotAngle must be between 0.1 and 2
" if ($PJ_RotAngle < #my $IR_Ny = FilterToNumber($in{'IR_Ny'}); my $IR_Nx = $Phantom_Nx; my $IR_Ny = $Phantom_Ny; +my $IR_FilterMethod = FilterMetaChars($in{'IR_FilterMethod'}); my $IR_Filter = FilterMetaChars($in{'IR_Filter'}); my $IR_Filter_Param = FilterToNumber($in{'IR_Filter_Param'}); my $IR_Interp = FilterMetaChars($in{'IR_Interp'}); @@ -107,7 +108,7 @@ $phm2if_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2if-lam" if $MPI; my $gp_cmd = "$phm2if_ver $phantom_fname $Phantom_Nx $Phantom_Ny --phantom $Phantom_Name --nsample $Phantom_NSample"; my $pj_cmd = "$phm2pj_ver $pj_fname $PJ_NDet $PJ_NRot --phantom $Phantom_Name --nray $PJ_NRay --rotangle $PJ_RotAngle"; my $pj_if_cmd = "$::bindir/pj2if $pj_fname $pj_if_fname"; -my $pjrec_cmd = "$pjrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj"; +my $pjrec_cmd = "$pjrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj --filter-method $IR_FilterMethod"; my $diff_cmd = "$diff_ver $phantom_fname $ir_fname $diff_fname --comp"; my $compare_cmd = "$ifinfo_ver $phantom_fname $ir_fname"; diff --git a/cgi-bin/ctsim.conf.in b/cgi-bin/ctsim.conf.in index e6fd9b4..2b09f38 100755 --- a/cgi-bin/ctsim.conf.in +++ b/cgi-bin/ctsim.conf.in @@ -6,3 +6,5 @@ $::jobdir = "@webdatadir@"; $::datadir = "@webdatadir@"; $::url_datadir = "@webdataurl@"; $::mpi_enable = "@mpienable@"; + +1; diff --git a/config.h.in b/config.h.in index 377a0fd..1a9d08d 100644 --- a/config.h.in +++ b/config.h.in @@ -54,6 +54,8 @@ #undef HAVE_WXWINDOWS +#undef HAVE_FFTW + /* The number of bytes in a double. */ #undef SIZEOF_DOUBLE diff --git a/configure b/configure index ccbce8f..1c3e0c3 100755 --- a/configure +++ b/configure @@ -1989,10 +1989,50 @@ else dmalloc=false fi +echo $ac_n "checking for main in -lfftw""... $ac_c" 1>&6 +echo "configure:1994: checking for main in -lfftw" >&5 +ac_lib_var=`echo fftw'_'main | sed 'y%./+-%__p_%'` +if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then + echo $ac_n "(cached) $ac_c" 1>&6 +else + ac_save_LIBS="$LIBS" +LIBS="-lfftw $LIBS" +cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then + rm -rf conftest* + eval "ac_cv_lib_$ac_lib_var=yes" +else + echo "configure: failed program was:" >&5 + cat conftest.$ac_ext >&5 + rm -rf conftest* + eval "ac_cv_lib_$ac_lib_var=no" +fi +rm -f conftest* +LIBS="$ac_save_LIBS" + +fi +if eval "test \"`echo '$ac_cv_lib_'$ac_lib_var`\" = yes"; then + echo "$ac_t""yes" 1>&6 + fftw=true; cat >> confdefs.h <<\EOF +#define HAVE_FFTW 1 +EOF + +else + echo "$ac_t""no" 1>&6 +fftw=false +fi + if test "$zlib" = "true" ; then echo $ac_n "checking for main in -lpng""... $ac_c" 1>&6 -echo "configure:1996: checking for main in -lpng" >&5 +echo "configure:2036: checking for main in -lpng" >&5 ac_lib_var=`echo png'_'main | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -2000,14 +2040,14 @@ else ac_save_LIBS="$LIBS" LIBS="-lpng $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2051: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -2034,7 +2074,7 @@ fi fi echo $ac_n "checking how to run the C preprocessor""... $ac_c" 1>&6 -echo "configure:2038: checking how to run the C preprocessor" >&5 +echo "configure:2078: checking how to run the C preprocessor" >&5 # On Suns, sometimes $CPP names a directory. if test -n "$CPP" && test -d "$CPP"; then CPP= @@ -2049,13 +2089,13 @@ else # On the NeXT, cc -E runs the code through the compiler's parser, # not just through cpp. cat > conftest.$ac_ext < Syntax Error EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2059: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2099: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then : @@ -2066,13 +2106,13 @@ else rm -rf conftest* CPP="${CC-cc} -E -traditional-cpp" cat > conftest.$ac_ext < Syntax Error EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2076: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2116: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then : @@ -2083,13 +2123,13 @@ else rm -rf conftest* CPP="${CC-cc} -nologo -E" cat > conftest.$ac_ext < Syntax Error EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2093: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2133: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then : @@ -2114,12 +2154,12 @@ fi echo "$ac_t""$CPP" 1>&6 echo $ac_n "checking for ANSI C header files""... $ac_c" 1>&6 -echo "configure:2118: checking for ANSI C header files" >&5 +echo "configure:2158: checking for ANSI C header files" >&5 if eval "test \"`echo '$''{'ac_cv_header_stdc'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < #include @@ -2127,7 +2167,7 @@ else #include EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2131: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2171: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then rm -rf conftest* @@ -2144,7 +2184,7 @@ rm -f conftest* if test $ac_cv_header_stdc = yes; then # SunOS 4.x string.h does not declare mem*, contrary to ANSI. cat > conftest.$ac_ext < EOF @@ -2162,7 +2202,7 @@ fi if test $ac_cv_header_stdc = yes; then # ISC 2.0.2 stdlib.h does not declare free, contrary to ANSI. cat > conftest.$ac_ext < EOF @@ -2183,7 +2223,7 @@ if test "$cross_compiling" = yes; then : else cat > conftest.$ac_ext < #define ISLOWER(c) ('a' <= (c) && (c) <= 'z') @@ -2194,7 +2234,7 @@ if (XOR (islower (i), ISLOWER (i)) || toupper (i) != TOUPPER (i)) exit(2); exit (0); } EOF -if { (eval echo configure:2198: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext} && (./conftest; exit) 2>/dev/null +if { (eval echo configure:2238: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext} && (./conftest; exit) 2>/dev/null then : else @@ -2221,17 +2261,17 @@ for ac_hdr in fcntl.h unistd.h getopt.h sys/fcntl.h setjmp.h stdarg.h stddef.h s do ac_safe=`echo "$ac_hdr" | sed 'y%./+-%__p_%'` echo $ac_n "checking for $ac_hdr""... $ac_c" 1>&6 -echo "configure:2225: checking for $ac_hdr" >&5 +echo "configure:2265: checking for $ac_hdr" >&5 if eval "test \"`echo '$''{'ac_cv_header_$ac_safe'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2235: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2275: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then rm -rf conftest* @@ -2259,12 +2299,12 @@ done echo $ac_n "checking for working const""... $ac_c" 1>&6 -echo "configure:2263: checking for working const" >&5 +echo "configure:2303: checking for working const" >&5 if eval "test \"`echo '$''{'ac_cv_c_const'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_compile) 2>&5; }; then +if { (eval echo configure:2357: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then rm -rf conftest* ac_cv_c_const=yes else @@ -2334,12 +2374,12 @@ EOF fi echo $ac_n "checking for off_t""... $ac_c" 1>&6 -echo "configure:2338: checking for off_t" >&5 +echo "configure:2378: checking for off_t" >&5 if eval "test \"`echo '$''{'ac_cv_type_off_t'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < #if STDC_HEADERS @@ -2367,12 +2407,12 @@ EOF fi echo $ac_n "checking for size_t""... $ac_c" 1>&6 -echo "configure:2371: checking for size_t" >&5 +echo "configure:2411: checking for size_t" >&5 if eval "test \"`echo '$''{'ac_cv_type_size_t'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < #if STDC_HEADERS @@ -2400,12 +2440,12 @@ EOF fi echo $ac_n "checking whether struct tm is in sys/time.h or time.h""... $ac_c" 1>&6 -echo "configure:2404: checking whether struct tm is in sys/time.h or time.h" >&5 +echo "configure:2444: checking whether struct tm is in sys/time.h or time.h" >&5 if eval "test \"`echo '$''{'ac_cv_struct_tm'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < #include @@ -2413,7 +2453,7 @@ int main() { struct tm *tp; tp->tm_sec; ; return 0; } EOF -if { (eval echo configure:2417: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then +if { (eval echo configure:2457: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then rm -rf conftest* ac_cv_struct_tm=time.h else @@ -2435,12 +2475,12 @@ fi echo $ac_n "checking for vprintf""... $ac_c" 1>&6 -echo "configure:2439: checking for vprintf" >&5 +echo "configure:2479: checking for vprintf" >&5 if eval "test \"`echo '$''{'ac_cv_func_vprintf'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2507: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_vprintf=yes" else @@ -2487,12 +2527,12 @@ fi if test "$ac_cv_func_vprintf" != yes; then echo $ac_n "checking for _doprnt""... $ac_c" 1>&6 -echo "configure:2491: checking for _doprnt" >&5 +echo "configure:2531: checking for _doprnt" >&5 if eval "test \"`echo '$''{'ac_cv_func__doprnt'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2559: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func__doprnt=yes" else @@ -2542,12 +2582,12 @@ fi for ac_func in strtod strtol snprintf htonl do echo $ac_n "checking for $ac_func""... $ac_c" 1>&6 -echo "configure:2546: checking for $ac_func" >&5 +echo "configure:2586: checking for $ac_func" >&5 if eval "test \"`echo '$''{'ac_cv_func_$ac_func'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2614: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_$ac_func=yes" else @@ -2595,12 +2635,12 @@ fi done echo $ac_n "checking for basename""... $ac_c" 1>&6 -echo "configure:2599: checking for basename" >&5 +echo "configure:2639: checking for basename" >&5 if eval "test \"`echo '$''{'ac_cv_func_basename'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2667: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_basename=yes" else @@ -2643,12 +2683,12 @@ else fi echo $ac_n "checking for setjmp""... $ac_c" 1>&6 -echo "configure:2647: checking for setjmp" >&5 +echo "configure:2687: checking for setjmp" >&5 if eval "test \"`echo '$''{'ac_cv_func_setjmp'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2715: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_setjmp=yes" else @@ -2694,12 +2734,12 @@ if test "${OSTYPE}" = "cygwin" ; then getopt_long=false else echo $ac_n "checking for getopt_long""... $ac_c" 1>&6 -echo "configure:2698: checking for getopt_long" >&5 +echo "configure:2738: checking for getopt_long" >&5 if eval "test \"`echo '$''{'ac_cv_func_getopt_long'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:2766: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_getopt_long=yes" else @@ -2769,7 +2809,7 @@ fi fi echo $ac_n "checking whether to enable verbose warnings""... $ac_c" 1>&6 -echo "configure:2773: checking whether to enable verbose warnings" >&5 +echo "configure:2813: checking whether to enable verbose warnings" >&5 # Check whether --enable-verbose-warnings or --disable-verbose-warnings was given. if test "${enable_verbose_warnings+set}" = set; then enableval="$enable_verbose_warnings" @@ -2819,7 +2859,7 @@ fi echo $ac_n "checking debug""... $ac_c" 1>&6 -echo "configure:2823: checking debug" >&5 +echo "configure:2863: checking debug" >&5 # Check whether --enable-debug or --disable-debug was given. if test "${enable_debug+set}" = set; then enableval="$enable_debug" @@ -2865,7 +2905,7 @@ if test "$withval" != "no" ; then fi echo $ac_n "checking for LAM MPI installation""... $ac_c" 1>&6 -echo "configure:2869: checking for LAM MPI installation" >&5 +echo "configure:2909: checking for LAM MPI installation" >&5 for testlamdir in "." $trylamdir /usr/local /usr/local/lam /usr /usr/lam /opt /opt/lam ; do if test -x "$testlamdir/bin/hcc" ; then LDFLAGS="$LDFLAGS -L$testlamdir/lib" @@ -2883,7 +2923,7 @@ else fi echo $ac_n "checking for web access""... $ac_c" 1>&6 -echo "configure:2887: checking for web access" >&5 +echo "configure:2927: checking for web access" >&5 # Check whether --with-cgibin-dir or --without-cgibin-dir was given. if test "${with_cgibin_dir+set}" = set; then withval="$with_cgibin_dir" @@ -2973,7 +3013,7 @@ fi # Uses ac_ vars as temps to allow command line to override cache and checks. # --without-x overrides everything else, but does not touch the cache. echo $ac_n "checking for X""... $ac_c" 1>&6 -echo "configure:2977: checking for X" >&5 +echo "configure:3017: checking for X" >&5 # Check whether --with-x or --without-x was given. if test "${with_x+set}" = set; then @@ -3035,12 +3075,12 @@ if test "$ac_x_includes" = NO; then # First, try using that file with no special directory specified. cat > conftest.$ac_ext < EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:3044: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:3084: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then rm -rf conftest* @@ -3109,14 +3149,14 @@ if test "$ac_x_libraries" = NO; then ac_save_LIBS="$LIBS" LIBS="-l$x_direct_test_library $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3160: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* LIBS="$ac_save_LIBS" # We can link X programs with no special library path. @@ -3222,17 +3262,17 @@ else case "`(uname -sr) 2>/dev/null`" in "SunOS 5"*) echo $ac_n "checking whether -R must be followed by a space""... $ac_c" 1>&6 -echo "configure:3226: checking whether -R must be followed by a space" >&5 +echo "configure:3266: checking whether -R must be followed by a space" >&5 ac_xsave_LIBS="$LIBS"; LIBS="$LIBS -R$x_libraries" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3276: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_nospace=yes else @@ -3248,14 +3288,14 @@ rm -f conftest* else LIBS="$ac_xsave_LIBS -R $x_libraries" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3299: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_space=yes else @@ -3287,7 +3327,7 @@ rm -f conftest* # libraries were built with DECnet support. And karl@cs.umb.edu says # the Alpha needs dnet_stub (dnet does not exist). echo $ac_n "checking for dnet_ntoa in -ldnet""... $ac_c" 1>&6 -echo "configure:3291: checking for dnet_ntoa in -ldnet" >&5 +echo "configure:3331: checking for dnet_ntoa in -ldnet" >&5 ac_lib_var=`echo dnet'_'dnet_ntoa | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3295,7 +3335,7 @@ else ac_save_LIBS="$LIBS" LIBS="-ldnet $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3350: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3328,7 +3368,7 @@ fi if test $ac_cv_lib_dnet_dnet_ntoa = no; then echo $ac_n "checking for dnet_ntoa in -ldnet_stub""... $ac_c" 1>&6 -echo "configure:3332: checking for dnet_ntoa in -ldnet_stub" >&5 +echo "configure:3372: checking for dnet_ntoa in -ldnet_stub" >&5 ac_lib_var=`echo dnet_stub'_'dnet_ntoa | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3336,7 +3376,7 @@ else ac_save_LIBS="$LIBS" LIBS="-ldnet_stub $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3391: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3376,12 +3416,12 @@ fi # The nsl library prevents programs from opening the X display # on Irix 5.2, according to dickey@clark.net. echo $ac_n "checking for gethostbyname""... $ac_c" 1>&6 -echo "configure:3380: checking for gethostbyname" >&5 +echo "configure:3420: checking for gethostbyname" >&5 if eval "test \"`echo '$''{'ac_cv_func_gethostbyname'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3448: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_gethostbyname=yes" else @@ -3425,7 +3465,7 @@ fi if test $ac_cv_func_gethostbyname = no; then echo $ac_n "checking for gethostbyname in -lnsl""... $ac_c" 1>&6 -echo "configure:3429: checking for gethostbyname in -lnsl" >&5 +echo "configure:3469: checking for gethostbyname in -lnsl" >&5 ac_lib_var=`echo nsl'_'gethostbyname | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3433,7 +3473,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lnsl $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3488: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3474,12 +3514,12 @@ fi # -lsocket must be given before -lnsl if both are needed. # We assume that if connect needs -lnsl, so does gethostbyname. echo $ac_n "checking for connect""... $ac_c" 1>&6 -echo "configure:3478: checking for connect" >&5 +echo "configure:3518: checking for connect" >&5 if eval "test \"`echo '$''{'ac_cv_func_connect'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3546: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_connect=yes" else @@ -3523,7 +3563,7 @@ fi if test $ac_cv_func_connect = no; then echo $ac_n "checking for connect in -lsocket""... $ac_c" 1>&6 -echo "configure:3527: checking for connect in -lsocket" >&5 +echo "configure:3567: checking for connect in -lsocket" >&5 ac_lib_var=`echo socket'_'connect | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3531,7 +3571,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lsocket $X_EXTRA_LIBS $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3586: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3566,12 +3606,12 @@ fi # gomez@mi.uni-erlangen.de says -lposix is necessary on A/UX. echo $ac_n "checking for remove""... $ac_c" 1>&6 -echo "configure:3570: checking for remove" >&5 +echo "configure:3610: checking for remove" >&5 if eval "test \"`echo '$''{'ac_cv_func_remove'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3638: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_remove=yes" else @@ -3615,7 +3655,7 @@ fi if test $ac_cv_func_remove = no; then echo $ac_n "checking for remove in -lposix""... $ac_c" 1>&6 -echo "configure:3619: checking for remove in -lposix" >&5 +echo "configure:3659: checking for remove in -lposix" >&5 ac_lib_var=`echo posix'_'remove | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3623,7 +3663,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lposix $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3678: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3658,12 +3698,12 @@ fi # BSDI BSD/OS 2.1 needs -lipc for XOpenDisplay. echo $ac_n "checking for shmat""... $ac_c" 1>&6 -echo "configure:3662: checking for shmat" >&5 +echo "configure:3702: checking for shmat" >&5 if eval "test \"`echo '$''{'ac_cv_func_shmat'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3730: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_func_shmat=yes" else @@ -3707,7 +3747,7 @@ fi if test $ac_cv_func_shmat = no; then echo $ac_n "checking for shmat in -lipc""... $ac_c" 1>&6 -echo "configure:3711: checking for shmat in -lipc" >&5 +echo "configure:3751: checking for shmat in -lipc" >&5 ac_lib_var=`echo ipc'_'shmat | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3715,7 +3755,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lipc $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3770: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3759,7 +3799,7 @@ fi # libraries we check for below, so use a different variable. # --interran@uluru.Stanford.EDU, kb@cs.umb.edu. echo $ac_n "checking for IceConnectionNumber in -lICE""... $ac_c" 1>&6 -echo "configure:3763: checking for IceConnectionNumber in -lICE" >&5 +echo "configure:3803: checking for IceConnectionNumber in -lICE" >&5 ac_lib_var=`echo ICE'_'IceConnectionNumber | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3767,7 +3807,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lICE $X_EXTRA_LIBS $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3822: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3815,7 +3855,7 @@ X_BASIC_LIBS="-lXext -lX11" our_saved_LDFLAGS="$LDFLAGS" LDFLAGS="$X_LIBS $LDFLAGS" echo $ac_n "checking for XtToolkitThreadInitialize in -lXt""... $ac_c" 1>&6 -echo "configure:3819: checking for XtToolkitThreadInitialize in -lXt" >&5 +echo "configure:3859: checking for XtToolkitThreadInitialize in -lXt" >&5 ac_lib_var=`echo Xt'_'XtToolkitThreadInitialize | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3823,7 +3863,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lXt $X_PRE_LIBS $X_BASIC_LIBS $X_EXTRA_LIBS $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3878: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3862,7 +3902,7 @@ LDFLAGS="$our_saved_LDFLAGS" our_saved_LDFLAGS="$LDFLAGS" LDFLAGS="$X_LIBS $LDFLAGS" echo $ac_n "checking for XdbeQueryExtension in -lXext""... $ac_c" 1>&6 -echo "configure:3866: checking for XdbeQueryExtension in -lXext" >&5 +echo "configure:3906: checking for XdbeQueryExtension in -lXext" >&5 ac_lib_var=`echo Xext'_'XdbeQueryExtension | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3870,7 +3910,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lXext -lX11 "$X_EXTRA_LIBS" $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3925: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3905,7 +3945,7 @@ else fi echo $ac_n "checking for XmbufQueryExtension in -lXext""... $ac_c" 1>&6 -echo "configure:3909: checking for XmbufQueryExtension in -lXext" >&5 +echo "configure:3949: checking for XmbufQueryExtension in -lXext" >&5 ac_lib_var=`echo Xext'_'XmbufQueryExtension | sed 'y%./+-%__p_%'` if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 @@ -3913,7 +3953,7 @@ else ac_save_LIBS="$LIBS" LIBS="-lXext -lX11 "$X_EXTRA_LIBS" $LIBS" cat > conftest.$ac_ext <&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then +if { (eval echo configure:3968: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* eval "ac_cv_lib_$ac_lib_var=yes" else @@ -3954,17 +3994,17 @@ for ac_hdr in X11/extensions/Xdbe.h X11/extensions/multibuf.h do ac_safe=`echo "$ac_hdr" | sed 'y%./+-%__p_%'` echo $ac_n "checking for $ac_hdr""... $ac_c" 1>&6 -echo "configure:3958: checking for $ac_hdr" >&5 +echo "configure:3998: checking for $ac_hdr" >&5 if eval "test \"`echo '$''{'ac_cv_header_$ac_safe'+set}'`\" = set"; then echo $ac_n "(cached) $ac_c" 1>&6 else cat > conftest.$ac_ext < EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:3968: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:4008: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } ac_err=`grep -v '^ *+' conftest.out | grep -v "^conftest.${ac_ext}\$"` if test -z "$ac_err"; then rm -rf conftest* @@ -4006,7 +4046,7 @@ my_includes="$my_includes -I../include -I.." echo $ac_n "checking interactive graphics""... $ac_c" 1>&6 -echo "configure:4010: checking interactive graphics" >&5 +echo "configure:4050: checking interactive graphics" >&5 if test "$no_x" != "yes" ; then cat >> confdefs.h <<\EOF #define HAVE_X11 1 @@ -4045,6 +4085,10 @@ if test "$wx_gtk" = "true" ; then ctlibs_tools="$ctlibs_tools -lwx_gtk" fi +if test "$fftw" = "true" ; then + ctlibs_tools="$ctlibs_tools -lrfftw -lfftw" +fi + LDFLAGS="$LDFLAGS -L../libctsupport -L../libctsim" ctlibs="$ctlibs_base -lctsim $ctlibs_graphics -lctsupport $ctlibs_tools" diff --git a/configure.in b/configure.in index 5f43602..e96baa4 100644 --- a/configure.in +++ b/configure.in @@ -80,6 +80,7 @@ AC_CHECK_LIB(g2, main, [g2=true], [g2=false]) AC_CHECK_LIB(wx_gtk, main, [wxwin=true; wx_gtk=true; AC_DEFINE(HAVE_WXWINDOWS)], [wxwin=false]) AC_CHECK_LIB(wx_msw, main, [wxwin=true; wx_msw=true; AC_DEFINE(HAVE_WXWINDOWS)], [wxwin=false]) AC_CHECK_LIB(dmallocxx, main, [dmalloc=true], [dmalloc=false]) +AC_CHECK_LIB(fftw, main, [fftw=true; AC_DEFINE(HAVE_FFTW)], [fftw=false]) if test "$zlib" = "true" ; then AC_CHECK_LIB(png, main, [png=true ; AC_DEFINE(HAVE_PNG)], [png=false]) @@ -309,6 +310,10 @@ if test "$wx_gtk" = "true" ; then ctlibs_tools="$ctlibs_tools -lwx_gtk" fi +if test "$fftw" = "true" ; then + ctlibs_tools="$ctlibs_tools -lrfftw -lfftw" +fi + dnl Setting projet libraries and includes LDFLAGS="$LDFLAGS -L../libctsupport -L../libctsim" ctlibs="$ctlibs_base -lctsim $ctlibs_graphics -lctsupport $ctlibs_tools" diff --git a/html/simulate.html.in b/html/simulate.html.in index 018a838..632b4ef 100644 --- a/html/simulate.html.in +++ b/html/simulate.html.in @@ -34,6 +34,11 @@ Number of Rays
(samples) per detector: as a multiple of PI:

Image Reconstruction

+Filter Method:
+Convolution
+Fourier
+FFT (Fast Fourier
+

Interpolation:
Linear
Nearest Neighbor
diff --git a/include/ct.h b/include/ct.h index 54dd3cf..5cbb361 100644 --- a/include/ct.h +++ b/include/ct.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ct.h,v 1.24 2000/07/02 18:21:39 kevin Exp $ +** $Id: ct.h,v 1.25 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -109,6 +109,10 @@ extern "C" { #include /* Standard ints on Linux */ #endif +#ifdef HAVE_FFTW +#include +#endif + #include #include #include diff --git a/include/filter.h b/include/filter.h index 7a50a7d..48f8210 100644 --- a/include/filter.h +++ b/include/filter.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.h,v 1.7 2000/07/03 11:02:06 kevin Exp $ +** $Id: filter.h,v 1.8 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -106,8 +106,12 @@ class SignalFilter { static void finiteFourierTransform (const float input[], complex output[], const int n, const int direction); + static void finiteFourierTransform (const complex input[], complex output[], const int n, const int direction); + void finiteFourierTransform (const float input[], complex output[], const int direction) const; + void finiteFourierTransform (const complex input[], complex output[], const int direction) const; + bool fail(void) const {return m_fail;} const string& failMessage(void) const {return m_failMessage;} @@ -129,6 +133,8 @@ class SignalFilter { static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param); + static void dotProduct (const double v1[], const complex v2[], complex output[], const int n); + private: double m_bw; int m_nFilterPoints; @@ -140,6 +146,10 @@ class SignalFilter { double* m_vecFilter; double* m_vecFourierCosTable; double* m_vecFourierSinTable; + complex* m_complexVecFilter; +#ifdef HAVE_FFTW + fftw_plan m_planForward, m_planBackward; +#endif bool m_fail; string m_failMessage; diff --git a/libctgraphics/ezplot.cpp b/libctgraphics/ezplot.cpp index 81af0dd..666f759 100644 --- a/libctgraphics/ezplot.cpp +++ b/libctgraphics/ezplot.cpp @@ -6,7 +6,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezplot.cpp,v 1.3 2000/06/20 17:54:51 kevin Exp $ +** $Id: ezplot.cpp,v 1.4 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -99,8 +99,8 @@ static int y_fldwid, y_frac; /* in fraction of number, used for printf() */ static double xtl_wid, ytl_wid; /* length of ticks labels in NDC */ static double tl_height; /* height of tick labels in NDC */ -static char x_numfmt[10]; /* format to print x tick labels */ -static char y_numfmt[10]; /* format to print y tick labels */ +static char x_numfmt[20]; /* format to print x tick labels */ +static char y_numfmt[20]; /* format to print y tick labels */ /*---------------------------------------------------------------------- diff --git a/libctgraphics/ezsupport.cpp b/libctgraphics/ezsupport.cpp index 8945787..de495e9 100644 --- a/libctgraphics/ezsupport.cpp +++ b/libctgraphics/ezsupport.cpp @@ -2,7 +2,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ezsupport.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $ +** $Id: ezsupport.cpp,v 1.3 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -151,7 +151,7 @@ make_numfmt (char *fmtstr, int *fldwid, int *nfrac, double minval, double maxval wid = 2 + frac + expon; if (minval < 0. || maxval < 0.) ++wid; - snprintf (fmtstr, sizeof(fmtstr), "%%%d.%dle", wid, frac); + sprintf (fmtstr, "%s%d%s%d%s", "%", wid, ".", frac, "g"); } else { /* use fixed format */ wid = static_cast(trunc(logt)) + 1; if (wid < 1) @@ -168,9 +168,9 @@ make_numfmt (char *fmtstr, int *fldwid, int *nfrac, double minval, double maxval frac = *nfrac; wid += 1 + frac; - snprintf (fmtstr, sizeof(fmtstr), "%%%d.%dlf", wid, frac); + sprintf (fmtstr, "%s%d%s%d%s", "%", wid, ".", frac, "f"); } - + *fldwid = wid; *nfrac = frac; } diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 3215d6a..27d1829 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: filter.cpp,v 1.7 2000/07/03 11:02:06 kevin Exp $ +** $Id: filter.cpp,v 1.8 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -37,14 +37,14 @@ * int filt_type Type of filter wanted * double bw Bandwidth of filter * double filterMin, filterMax Filter limits - * int n Number of points in signal + * int nSignalPoints Number of points in signal * double param General input parameter to filters * int domain FREQUENCY or SPATIAL domain wanted * int numint Number if intervals for calculating discrete inverse fourier xform * for spatial domain filters. For ANALYTIC solutions, use numint = 0 */ -SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int n, double param, const char* domainName, int numIntegral = 0) +SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int numIntegral = 0) { m_vecFilter = NULL; m_vecFourierCosTable = NULL; @@ -70,12 +70,12 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName m_failMessage += domainName; return; } - init (m_idFilter, m_idFilterMethod, bw, signalIncrement, n, param, m_idDomain, numIntegral); + init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, numIntegral); } -SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domainID, int numIntegral = 0) +SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int numIntegral = 0) { - init (filterID, filterMethodID, bw, signalIncrement, n, param, domainID, numIntegral); + init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, numIntegral); } SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0) @@ -105,7 +105,7 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub } void -SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domainID, int numint) +SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int numint) { m_bw = bw; m_idFilter = filterID; @@ -119,53 +119,96 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID m_nameDomain = convertDomainIDToName (m_idDomain); m_nameFilterMethod = convertFilterMethodIDToName (m_idFilterMethod); m_fail = false; - m_nSignalPoints = n; - m_nFilterPoints = 2 * m_nSignalPoints - 1; - + m_nSignalPoints = nSignalPoints; m_signalInc = signalIncrement; - m_filterMin = -signalIncrement * (m_nSignalPoints - 1); - m_filterMax = signalIncrement * (m_nSignalPoints - 1); - m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1); - m_numIntegral = numint; m_filterParam = param; - m_vecFilter = new double[ m_nFilterPoints ]; + if (m_idFilterMethod == FILTER_METHOD_FOURIER) { - int nFourier = n * n + 1; - double angleIncrement = (2. * PI) / n; + int nFourier = m_nSignalPoints * m_nSignalPoints + 1; + double angleIncrement = (2. * PI) / m_nSignalPoints; m_vecFourierCosTable = new double[ nFourier ]; m_vecFourierSinTable = new double[ nFourier ]; for (int i = 0; i < nFourier; i++) { m_vecFourierCosTable[i] = cos (angleIncrement * i); m_vecFourierSinTable[i] = sin (angleIncrement * i); } + m_nFilterPoints = m_nSignalPoints; + m_filterMin = 0; + m_filterMax = m_nSignalPoints * m_signalInc; + m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1); + m_vecFilter = new double [m_nFilterPoints]; + int halfFilter = m_nFilterPoints / 2; + for (int i = 0; i < halfFilter; i++) + m_vecFilter[i] = static_cast(i) / (halfFilter - 1) / (2 * m_signalInc); + for (int i = 0; i < halfFilter; i++) + m_vecFilter[m_nFilterPoints - i - 1] = static_cast(i) / (halfFilter - 1) / (2 * m_signalInc); + if (halfFilter % 2) // odd + m_vecFilter[halfFilter] = 1; + } else if (m_idFilterMethod == FILTER_METHOD_FFT || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) { + m_nFilterPoints = m_nSignalPoints; + if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_2 || m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) { + double logBase2 = log(m_nSignalPoints) / log(2); + int nextPowerOf2 = static_cast(floor(logBase2)) + 1; + if (m_idFilterMethod == FILTER_METHOD_FFT_ZEROPAD_4) + nextPowerOf2++; + if (logBase2 != floor(logBase2)) + nextPowerOf2++; + m_nFilterPoints = 1 << nextPowerOf2; + cout << "nFilterPoints = " << m_nFilterPoints << endl; + } + m_filterMin = 0; + m_filterMax = m_nSignalPoints * m_signalInc; + m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1); + m_vecFilter = new double [m_nFilterPoints]; + int halfFilter = m_nFilterPoints / 2; + for (int i = 0; i < halfFilter; i++) + m_vecFilter[i] = static_cast(i) / (halfFilter - 1) / (2 * m_signalInc); + for (int i = 0; i < halfFilter; i++) + m_vecFilter[m_nFilterPoints - i - 1] = static_cast(i) / (halfFilter - 1) / (2 * m_signalInc); + if (halfFilter % 2) // odd + m_vecFilter[halfFilter] = 1; + +#if HAVE_FFTW + m_planForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); + m_planBackward = fftw_create_plan (m_nFilterPoints, FFTW_BACKWARD, FFTW_ESTIMATE); +#endif } - if (m_idFilter == FILTER_SHEPP) { - double a = 2 * m_bw; - double c = - 4. / (a * a); - int center = (m_nFilterPoints - 1) / 2; - int sidelen = center; - m_vecFilter[center] = 4. / (a * a); - - for (int i = 1; i <= sidelen; i++ ) - m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1); - } else if (m_idDomain == DOMAIN_FREQUENCY) { - double x; - int i; - for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) - m_vecFilter[i] = frequencyResponse (x, param); - } else if (m_idDomain == DOMAIN_SPATIAL) { - double x; - int i; - for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) - if (numint == 0) - m_vecFilter[i] = spatialResponseAnalytic (x, param); - else - m_vecFilter[i] = spatialResponseCalc (x, param, numint); - } else { + if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) { + m_nFilterPoints = 2 * m_nSignalPoints - 1; + m_filterMin = -m_signalInc * (m_nSignalPoints - 1); + m_filterMax = m_signalInc * (m_nSignalPoints - 1); + m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1); + m_numIntegral = numint; + m_vecFilter = new double[ m_nFilterPoints ]; + + if (m_idFilter == FILTER_SHEPP) { + double a = 2 * m_bw; + double c = - 4. / (a * a); + int center = (m_nFilterPoints - 1) / 2; + int sidelen = center; + m_vecFilter[center] = 4. / (a * a); + + for (int i = 1; i <= sidelen; i++ ) + m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1); + } else if (m_idDomain == DOMAIN_FREQUENCY) { + double x; + int i; + for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) + m_vecFilter[i] = frequencyResponse (x, param); + } else if (m_idDomain == DOMAIN_SPATIAL) { + double x; + int i; + for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) + if (numint == 0) + m_vecFilter[i] = spatialResponseAnalytic (x, param); + else + m_vecFilter[i] = spatialResponseCalc (x, param, numint); + } else { m_failMessage = "Illegal domain name "; m_failMessage += m_idDomain; m_fail = true; + } } } @@ -174,6 +217,12 @@ SignalFilter::~SignalFilter (void) delete m_vecFilter; delete m_vecFourierSinTable; delete m_vecFourierCosTable; +#if HAVE_FFTW + if (m_idFilterMethod == FILTER_METHOD_FFT) { + fftw_destroy_plan(m_planForward); + fftw_destroy_plan(m_planBackward); + } +#endif } @@ -312,11 +361,31 @@ SignalFilter::filterSignal (const float input[], double output[]) const output[i] = convolve (input, m_signalInc, i, m_nSignalPoints); } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) { complex fftSignal[m_nSignalPoints]; - complex complexOutput; - finiteFourierTransform (input, fftSignal, 1); - // finiteFourierTransform (fftSignal, complexOutput, -1); - // for (int i = 0; i < m_nSignalPoints; i++) - // output[i] = complexOutput[i].hypot(); + complex complexOutput[m_nSignalPoints]; + complex filteredSignal[m_nSignalPoints]; + finiteFourierTransform (input, fftSignal, m_nSignalPoints, -1); + dotProduct (m_vecFilter, fftSignal, filteredSignal, m_nSignalPoints); + finiteFourierTransform (filteredSignal, complexOutput, m_nSignalPoints, 1); + for (int i = 0; i < m_nSignalPoints; i++) { + output[i] = abs( complexOutput[i] ); + } + } else if (m_idFilterMethod == FILTER_METHOD_FFT || FILTER_METHOD_FFT_ZEROPAD_2 || FILTER_METHOD_FFT_ZEROPAD_4) { + fftw_complex in[m_nFilterPoints], out[m_nFilterPoints]; + for (int i = 0; i < m_nSignalPoints; i++) { + in[i].re = input[i]; + in[i].im = 0; + } + for (int i = m_nSignalPoints; i < m_nFilterPoints; i++) { + in[i].re = in[i].im = 0; // ZeroPad + } + fftw_one(m_planForward, in, out); + for (int i = 0; i < m_nFilterPoints; i++) { + out[i].re = m_vecFilter[i] * out[i].re / m_nSignalPoints; + out[i].im = m_vecFilter[i] * out[i].im / m_nSignalPoints; + } + fftw_one(m_planBackward, out, in); + for (int i = 0; i < m_nSignalPoints; i++) + output[i] = sqrt (in[i].re * in[i].re + in[i].im * in[i].im); } } @@ -660,10 +729,10 @@ SignalFilter::finiteFourierTransform (const float input[], complex outpu double sumImag = 0; for (int j = 0; j < n; j++) { double angle = i * j * angleIncrement * direction; - sumReal += input[i] * cos(angle); - sumImag += input[i] * sin(angle); + sumReal += input[j] * cos(angle); + sumImag += input[j] * sin(angle); } - if (direction > 0) { + if (direction < 0) { sumReal /= n; sumImag /= n; } @@ -671,6 +740,30 @@ SignalFilter::finiteFourierTransform (const float input[], complex outpu } } + +void +SignalFilter::finiteFourierTransform (const complex input[], complex output[], const int n, int direction) +{ + if (direction < 0) + direction = -1; + else + direction = 1; + + double angleIncrement = 2 * PI / n; + for (int i = 0; i < n; i++) { + complex sum (0,0); + for (int j = 0; j < n; j++) { + double angle = i * j * angleIncrement * direction; + complex exponentTerm (cos(angle), sin(angle)); + sum += input[j] * exponentTerm; + } + if (direction < 0) { + sum /= n; + } + output[i] = sum; + } +} + void SignalFilter::finiteFourierTransform (const float input[], complex output[], int direction) const { @@ -679,7 +772,36 @@ SignalFilter::finiteFourierTransform (const float input[], complex outpu else direction = 1; - double angleIncrement = 2 * PI / m_nSignalPoints; + for (int i = 0; i < m_nSignalPoints; i++) { + double sumReal = 0, sumImag = 0; + for (int j = 0; j < m_nSignalPoints; j++) { + int tableIndex = i * j; + if (direction > 0) { + sumReal += input[i] * m_vecFourierCosTable[tableIndex]; + sumImag += input[i] * m_vecFourierSinTable[tableIndex]; + } else { + sumReal += input[i] * m_vecFourierCosTable[tableIndex]; + sumImag -= input[i] * m_vecFourierSinTable[tableIndex]; + } + } + if (direction < 0) { + sumReal /= m_nSignalPoints; + sumImag /= m_nSignalPoints; + } + output[i] = complex (sumReal, sumImag); + } +} + +// (a+bi) * (c + di) = (ac - db) + (bc + da)i +#if 0 +void +SignalFilter::finiteFourierTransform (const complex input[], complex output[], int direction) const +{ + if (direction < 0) + direction = -1; + else + direction = 1; + for (int i = 0; i < m_nSignalPoints; i++) { double sumReal = 0, sumImag = 0; for (int j = 0; j < m_nSignalPoints; j++) { @@ -699,3 +821,11 @@ SignalFilter::finiteFourierTransform (const float input[], complex outpu output[i] = complex (sumReal, sumImag); } } +#endif + +void +SignalFilter::dotProduct (const double v1[], const complex v2[], complex output[], const int n) +{ + for (int i = 0; i < n; i++) + output[i] = v1[i] * v2[i]; +} diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index b541864..1d90776 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.cpp,v 1.9 2000/07/03 11:02:06 kevin Exp $ +** $Id: projections.cpp,v 1.10 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -516,20 +516,22 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi #if HAVE_SGP SGP_ID gid; - int n_vec_filter = filter.getNFilterPoints(); - double plot_xaxis [n_vec_filter]; // array for plotting + int nVecFilter = filter.getNFilterPoints(); + double plot_xaxis [nVecFilter]; // array for plotting - if (trace > TRACE_TEXT) { + if (trace > TRACE_TEXT && nVecFilter > 0) { int i; double f; double filterInc = filter.getFilterIncrement(); - for (i = 0, f = filter.getFilterMin(); i < n_vec_filter; i++, f += filterInc) + for (i = 0, f = filter.getFilterMin(); i < nVecFilter; i++, f += filterInc) plot_xaxis[i] = f; - - gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter); - cio_put_str ("Press any key to continue"); - cio_kb_getc (); - sgp2_close (gid); + + if (filter.getFilter()) { + gid = ezplot (plot_xaxis, filter.getFilter(), nVecFilter); + cio_put_str ("Press any key to continue"); + cio_kb_getc (); + sgp2_close (gid); + } } if (trace >= TRACE_TEXT) { printf ("nview=%d, ndet=%d, det_start=%.4f, detInc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc); @@ -551,9 +553,6 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi filter.filterSignal (detval, filteredProj); - // for (int j = 0; j < m_nDet; j++) - // filteredProj[j] = filter.convolve (detval, detInc, j, m_nDet); - #ifdef HAVE_SGP if (trace >= TRACE_PLOT) { ezset ("clear."); diff --git a/src/pjrec.cpp b/src/pjrec.cpp index c14312b..0a2b2c8 100644 --- a/src/pjrec.cpp +++ b/src/pjrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: pjrec.cpp,v 1.4 2000/06/29 13:21:14 kevin Exp $ +** $Id: pjrec.cpp,v 1.5 2000/07/04 18:33:35 kevin Exp $ ** ** This program is free software; you can redistribute it and/or modify ** it under the terms of the GNU General Public License (version 2) as @@ -29,12 +29,13 @@ #include "timer.h" -enum {O_INTERP, O_FILTER, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION}; +enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION}; static struct option my_options[] = { {"interp", 1, 0, O_INTERP}, {"filter", 1, 0, O_FILTER}, + {"filter-method", 1, 0, O_FILTER_METHOD}, {"filter-param", 1, 0, O_FILTER_PARAM}, {"backproj", 1, 0, O_BACKPROJ}, {"trace", 1, 0, O_TRACE}, @@ -73,6 +74,9 @@ pjrec_usage (const char *program) cout << " cos Cosine" << endl; cout << " triangle Triangle" << endl; cout << " hamming Hamming" << endl; + cout << " --filter-method Filter method before backprojections\n";; + cout << " convolution Spatial filtering (default)\n"; + cout << " fourier Frequency filtering with discete fourier\n"; cout << " --backproj Backprojection Method" << endl; cout << " trig Trigometric functions at every point" << endl; cout << " table Trigometric functions with precalculated table" << endl; @@ -147,6 +151,9 @@ pjrec_main (int argc, char * argv[]) case O_FILTER: optFilterName = optarg; break; + case O_FILTER_METHOD: + optFilterMethodName = optarg; + break; case O_BACKPROJ: optBackprojName = optarg; break; -- 2.34.1