From 08f34bf3ba14d4f436f4d2ef0ee5af1d6eb266ac Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Sun, 25 Jun 2000 17:32:24 +0000 Subject: [PATCH] r120: *** empty log message *** --- ChangeLog | 12 ++ cgi-bin/ctsim.cgi.in | 7 + configure | 253 +++++++++++++++++++++--------------- configure.in | 3 +- include/backprojectors.h | 27 ++-- include/filter.h | 47 ++++--- include/imagefile.h | 102 ++++++++++++--- include/phantom.h | 25 ++-- include/projections.h | 6 +- include/scanner.h | 72 +++++----- libctsim/backprojectors.cpp | 19 ++- libctsim/filter.cpp | 87 +++++++++---- libctsim/imagefile.cpp | 176 ++++++++++++++++++++++--- libctsim/phantom.cpp | 19 +-- libctsim/projections.cpp | 12 +- libctsim/scanner.cpp | 30 ++++- src/ctrec.cpp | 84 ++++++------ src/if-2.cpp | 8 +- src/ifinfo.cpp | 143 ++++++-------------- src/phm2if.cpp | 55 ++++---- src/phm2pj.cpp | 54 ++++---- 21 files changed, 782 insertions(+), 459 deletions(-) diff --git a/ChangeLog b/ChangeLog index 58f4e14..332a837 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,15 @@ +1.9.7 - 6/25/2000 + Standardized string option processing by classes. All classes use character strings + to select options rather than numeric constants. Added fail() and failMessage() + methods to verify that objects are created correctly by character strings. + Hid C++ assignment and copy constructors in classes that should not have assignment/copy + Rewrote ImageFile class + Started support for dmallocxx library, not finished + Added G.T. Herman image comparision statistics to ifinfo (see imagefile.cpp) + Updated ifinfo to show comparative statistics + Added printLabels() to Array2dFile class + Added printStatistics() to ImageFile class + 1.9.6 - 6/22/2000 Moved conversion filter name/id to Filter class Moved conversion backprojection name/id to Backproj class diff --git a/cgi-bin/ctsim.cgi.in b/cgi-bin/ctsim.cgi.in index 898aa89..5e27ccb 100755 --- a/cgi-bin/ctsim.cgi.in +++ b/cgi-bin/ctsim.cgi.in @@ -98,6 +98,8 @@ my $ctrec_ver = "$::bindir/ctrec"; my $phm2pj_ver = "$::bindir/phm2pj"; my $phm2if_ver = "$::bindir/phm2if"; my $diff_ver = "$::bindir/if-2"; +my $ifinfo_ver = "$::bindir/ifinfo"; + $ctrec_ver = "/opt/lam/bin/mpirun N N $::lamrundir/ctrec-lam" if $MPI; $phm2pj_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2pj-lam" if $MPI; $phm2if_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2if-lam" if $MPI; @@ -107,6 +109,7 @@ my $pj_cmd = "$phm2pj_ver $pj_fname $PJ_NDet $PJ_NRot --phantom $Phantom_Name -- my $pj_if_cmd = "$::bindir/pj2if $pj_fname $pj_if_fname"; my $ir_cmd = "$ctrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj"; my $diff_cmd = "$diff_ver $phantom_fname $ir_fname $diff_fname --comp"; +my $compare_cmd = "$ifinfo_ver $phantom_fname $ir_fname"; my $window_options = "--auto $auto_window_img"; if ($Disp_Min ne 'auto') { @@ -149,6 +152,7 @@ if ($error ne "") { my $png_ir_out; my $png_pj_out; my $png_diff_out; + my $compare_out; $gp_out = `$gp_cmd`; if (-s $phantom_fname) { $pj_out .= `$pj_cmd`; @@ -161,6 +165,7 @@ if ($error ne "") { $png_ir_out .= `$png2_cmd`; $diff_out .= `$diff_cmd`; $png_diff_out .= `$png4_cmd`; + $compare_out = `$compare_cmd`; } } } @@ -196,6 +201,8 @@ if ($error ne "") { $out .= "
$diff_out
$png_diff_out\n"; $out .= ""; $out .= "Execution time: $execution_time seconds\n"; + $out .= "

\nStatistics
"; + $out .= "$compare_out"; } $out .= "

\n"; diff --git a/configure b/configure index 5402f8e..0143ad3 100755 --- a/configure +++ b/configure @@ -711,7 +711,7 @@ fi PACKAGE=ctsim -VERSION=1.9.6 +VERSION=1.9.7 if test "`cd $srcdir && pwd`" != "`pwd`" && test -f $srcdir/config.status; then { echo "configure: error: source directory already configured; run "make distclean" there first" 1>&2; exit 1; } @@ -1953,10 +1953,47 @@ else wxwin=false fi +echo $ac_n "checking for main in -ldmallocxx""... $ac_c" 1>&6 +echo "configure:1958: checking for main in -ldmallocxx" >&5 +ac_lib_var=`echo dmallocxx'_'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="-ldmallocxx $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 + dmalloc=true +else + echo "$ac_t""no" 1>&6 +dmalloc=false +fi + if test "$zlib" = "true" ; then echo $ac_n "checking for main in -lpng""... $ac_c" 1>&6 -echo "configure:1960: checking for main in -lpng" >&5 +echo "configure:1997: 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 @@ -1964,14 +2001,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:2012: \"$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 @@ -1998,7 +2035,7 @@ fi fi echo $ac_n "checking how to run the C preprocessor""... $ac_c" 1>&6 -echo "configure:2002: checking how to run the C preprocessor" >&5 +echo "configure:2039: checking how to run the C preprocessor" >&5 # On Suns, sometimes $CPP names a directory. if test -n "$CPP" && test -d "$CPP"; then CPP= @@ -2013,13 +2050,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:2023: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2060: \"$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 : @@ -2030,13 +2067,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:2040: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2077: \"$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 : @@ -2047,13 +2084,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:2057: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2094: \"$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 : @@ -2078,12 +2115,12 @@ fi echo "$ac_t""$CPP" 1>&6 echo $ac_n "checking for ANSI C header files""... $ac_c" 1>&6 -echo "configure:2082: checking for ANSI C header files" >&5 +echo "configure:2119: 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 @@ -2091,7 +2128,7 @@ else #include EOF ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out" -{ (eval echo configure:2095: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2132: \"$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* @@ -2108,7 +2145,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 @@ -2126,7 +2163,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 @@ -2147,7 +2184,7 @@ if test "$cross_compiling" = yes; then : else cat > conftest.$ac_ext < #define ISLOWER(c) ('a' <= (c) && (c) <= 'z') @@ -2158,7 +2195,7 @@ if (XOR (islower (i), ISLOWER (i)) || toupper (i) != TOUPPER (i)) exit(2); exit (0); } EOF -if { (eval echo configure:2162: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext} && (./conftest; exit) 2>/dev/null +if { (eval echo configure:2199: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext} && (./conftest; exit) 2>/dev/null then : else @@ -2185,17 +2222,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:2189: checking for $ac_hdr" >&5 +echo "configure:2226: 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:2199: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:2236: \"$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* @@ -2223,12 +2260,12 @@ done echo $ac_n "checking for working const""... $ac_c" 1>&6 -echo "configure:2227: checking for working const" >&5 +echo "configure:2264: 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:2318: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then rm -rf conftest* ac_cv_c_const=yes else @@ -2298,12 +2335,12 @@ EOF fi echo $ac_n "checking for off_t""... $ac_c" 1>&6 -echo "configure:2302: checking for off_t" >&5 +echo "configure:2339: 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 @@ -2331,12 +2368,12 @@ EOF fi echo $ac_n "checking for size_t""... $ac_c" 1>&6 -echo "configure:2335: checking for size_t" >&5 +echo "configure:2372: 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 @@ -2364,12 +2401,12 @@ EOF fi echo $ac_n "checking whether struct tm is in sys/time.h or time.h""... $ac_c" 1>&6 -echo "configure:2368: checking whether struct tm is in sys/time.h or time.h" >&5 +echo "configure:2405: 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 @@ -2377,7 +2414,7 @@ int main() { struct tm *tp; tp->tm_sec; ; return 0; } EOF -if { (eval echo configure:2381: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then +if { (eval echo configure:2418: \"$ac_compile\") 1>&5; (eval $ac_compile) 2>&5; }; then rm -rf conftest* ac_cv_struct_tm=time.h else @@ -2399,12 +2436,12 @@ fi echo $ac_n "checking for vprintf""... $ac_c" 1>&6 -echo "configure:2403: checking for vprintf" >&5 +echo "configure:2440: 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:2468: \"$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 @@ -2451,12 +2488,12 @@ fi if test "$ac_cv_func_vprintf" != yes; then echo $ac_n "checking for _doprnt""... $ac_c" 1>&6 -echo "configure:2455: checking for _doprnt" >&5 +echo "configure:2492: 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:2520: \"$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 @@ -2506,12 +2543,12 @@ fi for ac_func in strtod strtol snprintf htonl do echo $ac_n "checking for $ac_func""... $ac_c" 1>&6 -echo "configure:2510: checking for $ac_func" >&5 +echo "configure:2547: 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:2575: \"$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 @@ -2559,12 +2596,12 @@ fi done echo $ac_n "checking for basename""... $ac_c" 1>&6 -echo "configure:2563: checking for basename" >&5 +echo "configure:2600: 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:2628: \"$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 @@ -2607,12 +2644,12 @@ else fi echo $ac_n "checking for setjmp""... $ac_c" 1>&6 -echo "configure:2611: checking for setjmp" >&5 +echo "configure:2648: 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:2676: \"$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 @@ -2658,12 +2695,12 @@ if test "${OSTYPE}" = "cygwin" ; then getopt_long=false else echo $ac_n "checking for getopt_long""... $ac_c" 1>&6 -echo "configure:2662: checking for getopt_long" >&5 +echo "configure:2699: 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:2727: \"$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 @@ -2733,7 +2770,7 @@ fi fi echo $ac_n "checking whether to enable verbose warnings""... $ac_c" 1>&6 -echo "configure:2737: checking whether to enable verbose warnings" >&5 +echo "configure:2774: 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" @@ -2825,7 +2862,7 @@ if test "$withval" != "no" ; then fi echo $ac_n "checking for LAM MPI installation""... $ac_c" 1>&6 -echo "configure:2829: checking for LAM MPI installation" >&5 +echo "configure:2866: 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" @@ -2843,7 +2880,7 @@ else fi echo $ac_n "checking for web access""... $ac_c" 1>&6 -echo "configure:2847: checking for web access" >&5 +echo "configure:2884: 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" @@ -2933,7 +2970,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:2937: checking for X" >&5 +echo "configure:2974: checking for X" >&5 # Check whether --with-x or --without-x was given. if test "${with_x+set}" = set; then @@ -2995,12 +3032,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:3004: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:3041: \"$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* @@ -3069,14 +3106,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:3117: \"$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. @@ -3182,17 +3219,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:3186: checking whether -R must be followed by a space" >&5 +echo "configure:3223: 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:3233: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_nospace=yes else @@ -3208,14 +3245,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:3256: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_space=yes else @@ -3247,7 +3284,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:3251: checking for dnet_ntoa in -ldnet" >&5 +echo "configure:3288: 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 @@ -3255,7 +3292,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:3307: \"$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 @@ -3288,7 +3325,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:3292: checking for dnet_ntoa in -ldnet_stub" >&5 +echo "configure:3329: 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 @@ -3296,7 +3333,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:3348: \"$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 @@ -3336,12 +3373,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:3340: checking for gethostbyname" >&5 +echo "configure:3377: 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:3405: \"$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 @@ -3385,7 +3422,7 @@ fi if test $ac_cv_func_gethostbyname = no; then echo $ac_n "checking for gethostbyname in -lnsl""... $ac_c" 1>&6 -echo "configure:3389: checking for gethostbyname in -lnsl" >&5 +echo "configure:3426: 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 @@ -3393,7 +3430,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:3445: \"$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 @@ -3434,12 +3471,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:3438: checking for connect" >&5 +echo "configure:3475: 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:3503: \"$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 @@ -3483,7 +3520,7 @@ fi if test $ac_cv_func_connect = no; then echo $ac_n "checking for connect in -lsocket""... $ac_c" 1>&6 -echo "configure:3487: checking for connect in -lsocket" >&5 +echo "configure:3524: 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 @@ -3491,7 +3528,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:3543: \"$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 @@ -3526,12 +3563,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:3530: checking for remove" >&5 +echo "configure:3567: 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:3595: \"$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 @@ -3575,7 +3612,7 @@ fi if test $ac_cv_func_remove = no; then echo $ac_n "checking for remove in -lposix""... $ac_c" 1>&6 -echo "configure:3579: checking for remove in -lposix" >&5 +echo "configure:3616: 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 @@ -3583,7 +3620,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:3635: \"$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 @@ -3618,12 +3655,12 @@ fi # BSDI BSD/OS 2.1 needs -lipc for XOpenDisplay. echo $ac_n "checking for shmat""... $ac_c" 1>&6 -echo "configure:3622: checking for shmat" >&5 +echo "configure:3659: 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:3687: \"$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 @@ -3667,7 +3704,7 @@ fi if test $ac_cv_func_shmat = no; then echo $ac_n "checking for shmat in -lipc""... $ac_c" 1>&6 -echo "configure:3671: checking for shmat in -lipc" >&5 +echo "configure:3708: 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 @@ -3675,7 +3712,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:3727: \"$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 @@ -3719,7 +3756,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:3723: checking for IceConnectionNumber in -lICE" >&5 +echo "configure:3760: 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 @@ -3727,7 +3764,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:3779: \"$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 @@ -3775,7 +3812,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:3779: checking for XtToolkitThreadInitialize in -lXt" >&5 +echo "configure:3816: 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 @@ -3783,7 +3820,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:3835: \"$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 @@ -3822,7 +3859,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:3826: checking for XdbeQueryExtension in -lXext" >&5 +echo "configure:3863: 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 @@ -3830,7 +3867,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:3882: \"$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 @@ -3865,7 +3902,7 @@ else fi echo $ac_n "checking for XmbufQueryExtension in -lXext""... $ac_c" 1>&6 -echo "configure:3869: checking for XmbufQueryExtension in -lXext" >&5 +echo "configure:3906: 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 @@ -3873,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 @@ -3914,17 +3951,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:3918: checking for $ac_hdr" >&5 +echo "configure:3955: 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:3928: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:3965: \"$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* @@ -3966,7 +4003,7 @@ my_includes="$my_includes -I../include -I.." echo $ac_n "checking interactive graphics""... $ac_c" 1>&6 -echo "configure:3970: checking interactive graphics" >&5 +echo "configure:4007: checking interactive graphics" >&5 if test "$no_x" != "yes" ; then cat >> confdefs.h <<\EOF #define HAVE_X11 1 diff --git a/configure.in b/configure.in index 084ed5a..43e4a12 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/ctrec.cpp) -AM_INIT_AUTOMAKE(ctsim,1.9.6) +AM_INIT_AUTOMAKE(ctsim,1.9.7) AM_CONFIG_HEADER(config.h) dnl Checks for programs. @@ -79,6 +79,7 @@ AC_CHECK_LIB(ncurses, main, [ncurses=true], [ncurses=false]) 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]) if test "$zlib" = "true" ; then AC_CHECK_LIB(png, main, [png=true ; AC_DEFINE(HAVE_PNG)], [png=false]) diff --git a/include/backprojectors.h b/include/backprojectors.h index 08b72d4..6b8b9c4 100644 --- a/include/backprojectors.h +++ b/include/backprojectors.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.h,v 1.4 2000/06/22 10:17:28 kevin Exp $ +** $Id: backprojectors.h,v 1.5 2000/06/25 17:32:24 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 @@ -54,6 +54,17 @@ class Backprojector INTERP_LINEAR // Linear interpolation } InterpolationID; + static const char BPROJ_TRIG_STR[]= "trig"; + static const char BPROJ_TABLE_STR[]= "table"; + static const char BPROJ_DIFF_STR[]= "diff"; + static const char BPROJ_DIFF2_STR[]= "diff2"; + static const char BPROJ_IDIFF2_STR[]= "idiff2"; + + static const char INTERP_NEAREST_STR[]= "nearest"; + static const char INTERP_LINEAR_STR[]= "linear"; + static const char INTERP_BSPLINE_STR[]= "bspline"; + + Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName); ~Backprojector (void); @@ -61,24 +72,16 @@ class Backprojector void BackprojectView (const double* const viewData, const double viewAngle); bool fail(void) const {return m_fail;} + const string& failMessage(void) const {return m_failMessage;} private: string m_nameBackproject; string m_nameInterpolation; BackprojectID m_idBackproject; InterpolationID m_idInterpolation; - bool m_fail; Backproject* m_pBackprojectImplem; - - static const char BPROJ_TRIG_STR[]= "trig"; - static const char BPROJ_TABLE_STR[]= "table"; - static const char BPROJ_DIFF_STR[]= "diff"; - static const char BPROJ_DIFF2_STR[]= "diff2"; - static const char BPROJ_IDIFF2_STR[]= "idiff2"; - - static const char INTERP_NEAREST_STR[]= "nearest"; - static const char INTERP_LINEAR_STR[]= "linear"; - static const char INTERP_BSPLINE_STR[]= "bspline"; + bool m_fail; + string m_failMessage; static const InterpolationID convertInterpolationNameToID (const char* const interpName); static const char* convertInterpolationIDToName (const InterpolationID interpID); diff --git a/include/filter.h b/include/filter.h index b544a1f..9bbf554 100644 --- a/include/filter.h +++ b/include/filter.h @@ -25,10 +25,26 @@ class SignalFilter { DOMAIN_SPATIAL } DomainID; + static const char FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit"; + static const char FILTER_ABS_SINC_STR[]= "abs_sinc"; + static const char FILTER_ABS_COS_STR[]= "abs_cos"; + static const char FILTER_ABS_HAMMING_STR[]= "abs_hamming"; + static const char FILTER_SHEPP_STR[]= "shepp"; + static const char FILTER_BANDLIMIT_STR[]= "bandlimit"; + static const char FILTER_SINC_STR[]= "sinc"; + static const char FILTER_COS_STR[]= "cos"; + static const char FILTER_HAMMING_STR[]= "hamming"; + static const char FILTER_TRIANGLE_STR[]= "triangle"; + + static const char DOMAIN_FREQ_STR[]= "freq"; + static const char DOMAIN_SPATIAL_STR[]= "spatial"; + - SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, const int numInt); + SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, const int numIntegral = 0); - SignalFilter (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numInt); + SignalFilter (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numIntegral = 0); + + SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0); ~SignalFilter (void); @@ -40,19 +56,23 @@ class SignalFilter { double convolve (const float f[], const double dx, const int n, const int np) const; bool fail(void) const {return m_fail;} + const string& failMessage(void) const {return m_failMessage;} + const string& nameFilter(void) const { return m_nameFilter;} const string& nameDomain(void) const { return m_nameDomain;} const FilterID idFilter(void) const { return m_idFilter;} const DomainID idDomain(void) const { return m_idDomain;} - static double response (const char* const FilterName, const char* const DomainName, double bw, double x, double param); + double response (double x); + + static double spatialResponse (FilterID fType, double bw, double x, double param, int nIntegral = 0); + + static double frequencyResponse (FilterID fType, double bw, double u, double param); static double spatialResponseCalc (FilterID fType, double bw, double x, double param, int n); static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param); - static double frequencyResponse (FilterID fType, double bw, double u, double param); - private: double m_bw; int m_nPoints; @@ -60,24 +80,13 @@ class SignalFilter { double m_xmax; double* m_vecFilter; bool m_fail; + string m_failMessage; string m_nameFilter; string m_nameDomain; FilterID m_idFilter; DomainID m_idDomain; - - static const char FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit"; - static const char FILTER_ABS_SINC_STR[]= "abs_sinc"; - static const char FILTER_ABS_COS_STR[]= "abs_cos"; - static const char FILTER_ABS_HAMMING_STR[]= "abs_hamming"; - static const char FILTER_SHEPP_STR[]= "shepp"; - static const char FILTER_BANDLIMIT_STR[]= "bandlimit"; - static const char FILTER_SINC_STR[]= "sinc"; - static const char FILTER_COS_STR[]= "cos"; - static const char FILTER_HAMMING_STR[]= "hamming"; - static const char FILTER_TRIANGLE_STR[]= "triangle"; - - static const char DOMAIN_FREQ_STR[]= "freq"; - static const char DOMAIN_SPATIAL_STR[]= "spatial"; + double m_filterParam; + int m_numIntegral; static FilterID convertFilterNameToID (const char* filterName); static const char* convertFilterIDToName (const FilterID filterID); diff --git a/include/imagefile.h b/include/imagefile.h index 3b1882f..7c6c84a 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -15,11 +15,6 @@ using namespace std; class Array2dFileLabel { -private: - void init (void); - Array2dFileLabel (const Array2dFileLabel&); - Array2dFileLabel& operator= (const Array2dFileLabel&); - public: kfloat64 m_calcTime; kuint16 m_labelType; @@ -60,6 +55,11 @@ public: { m_strLabel = str; return (m_strLabel); } const string& getDateString () const; + +private: + void init (void); + Array2dFileLabel (const Array2dFileLabel&); + Array2dFileLabel& operator= (const Array2dFileLabel&); }; @@ -68,14 +68,6 @@ class Array2dFile { private: void init (void); - kuint16 signature; - kuint16 num_labels; - kuint16 headersize; - string filename; - int file_id; - frnetorderstream *m_pFS; - bool bHeaderWritten; - bool bDataWritten; bool headerWrite (void); @@ -85,6 +77,19 @@ private: bool labelSeek (unsigned int label_num); + Array2dFile (const Array2dFile& rhs); // copy constructor + Array2dFile& operator= (const Array2dFile&); // assignment operator + + protected: + kuint16 signature; + kuint16 num_labels; + kuint16 headersize; + string filename; + int file_id; + frnetorderstream *m_pFS; + bool bHeaderWritten; + bool bDataWritten; + kuint16 mPixelSize; kuint16 axis_increment_known; kfloat64 mIncX, mIncY; @@ -95,8 +100,6 @@ private: kuint32 mNX; kuint32 mNY; - Array2dFile (const Array2dFile&); - Array2dFile& operator= (const Array2dFile&); public: @@ -148,6 +151,8 @@ public: void doPixelOffsetScale (double offset, double scale); + void printLabels (ostream& os); + T** getArray (void) const { return (array.getArray()); } @@ -618,6 +623,28 @@ Array2dFile::arrayDataClear (void) } } +template +void +Array2dFile::printLabels (ostream& os) +{ + int nlabels = getNumLabels(); + + for (int i = 0; i < nlabels; i++) { + Array2dFileLabel label; + labelRead (label, i); + + if (label.getLabelType() == Array2dFileLabel::L_HISTORY) { + os << "History: " << endl; + os << " " << label.getLabelString() << endl; + os << " calc time = " << label.getCalcTime() << " secs" << endl; + os << " Timestamp = " << label.getDateString() << endl; + } else if (label.getLabelType() == Array2dFileLabel::L_USER) { + os << "Note: " << label.getLabelString() << endl; + os << " Timestamp = %s" << label.getDateString() << endl; + } + os << endl; + } +} #ifdef HAVE_MPI #include @@ -648,6 +675,10 @@ public: MPI::Datatype getMPIDataType (void) const { return MPI::FLOAT; } #endif + + private: + F32Image (const F32Image& rhs); //copy constructor + F32Image& operator= (const F32Image& rhs); // assignment operator }; @@ -677,26 +708,55 @@ class F64Image : public Array2dFile MPI::Datatype getMPIDataType (void) const { return MPI::DOUBLE; } #endif + private: + F64Image (const F64Image& rhs); //copy constructor + F64Image& operator= (const F64Image& rhs); // assignment operator }; #undef IMAGEFILE_64_BITS #ifdef IMAGEFILE_64_BITS -typedef F64Image ImageFile; +typedef F64Image ImageFileBase; typedef kfloat64 ImageFileValue; typedef kfloat64* ImageFileColumn; typedef kfloat64** ImageFileArray; #else -typedef F32Image ImageFile; +typedef F32Image ImageFileBase; typedef kfloat32 ImageFileValue; typedef kfloat32* ImageFileColumn; typedef kfloat32** ImageFileArray; #endif -// imagefile.cpp -void image_filter_response (ImageFile& im, const char* const domainName, double bw, const char* const filterName, double filt_param, const int opt_trace); -int image_display (const ImageFile& im); -int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax); +class ImageFile : public ImageFileBase +{ + public: + ImageFile (const char* const fname, unsigned int nx, unsigned int ny) + : ImageFileBase (fname, nx, ny) + {} + + ImageFile (unsigned int nx, unsigned int ny) + : ImageFileBase (nx, ny) + {} + + ImageFile (const char* const fname) + : ImageFileBase (fname) + {} + + void filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param); + + void statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const; + + void printStatistics (ostream& os) const; + + bool comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const; + + bool printComparativeStatistics (const ImageFile& imComp, ostream& os) const; + + int display (void); + + int displayScaling (const int scaleFactor, ImageFileValue pmin, ImageFileValue pmax); + +}; #endif diff --git a/include/phantom.h b/include/phantom.h index b3765c6..4a2680b 100644 --- a/include/phantom.h +++ b/include/phantom.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phantom.h,v 1.3 2000/06/22 10:17:28 kevin Exp $ +** $Id: phantom.h,v 1.4 2000/06/25 17:32:24 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 @@ -132,6 +132,11 @@ class Phantom PHM_UNITPULSE /* Unit pulse phantom */ } PhantomID; + static const char PHM_HERMAN_STR[]= "herman"; + static const char PHM_BHERMAN_STR[]= "bherman"; + static const char PHM_ROWLAND_STR[]= "rowland"; + static const char PHM_BROWLAND_STR[]= "browland"; + static const char PHM_UNITPULSE_STR[]= "unitpulse"; Phantom (void); Phantom (const char* const phmName); @@ -158,14 +163,10 @@ class Phantom void convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const; - bool fail(void) const - {return m_fail;} - - const string& name(void) const - {return m_name;} - - const PhantomID id(void) const - {return m_id;} + bool fail(void) const {return m_fail;} + const string& failMessage(void) const {return m_failMessage;} + const string& name(void) const {return m_name;} + const PhantomID id(void) const {return m_id;} #if HAVE_SGP void show (void) const; @@ -199,13 +200,9 @@ class Phantom string m_name; PhantomID m_id; bool m_fail; + string m_failMessage; // Standard phantomsa - static const char PHM_HERMAN_STR[]= "herman"; - static const char PHM_BHERMAN_STR[]= "bherman"; - static const char PHM_ROWLAND_STR[]= "rowland"; - static const char PHM_BROWLAND_STR[]= "browland"; - static const char PHM_UNITPULSE_STR[]= "unitpulse"; static PhantomID convertNameToPhantomID (const char* const phmName); static const char* convertPhantomIDToName (const PhantomID phmID); diff --git a/include/projections.h b/include/projections.h index a9970b3..6446f57 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: projections.h,v 1.4 2000/06/22 10:17:28 kevin Exp $ +** $Id: projections.h,v 1.5 2000/06/25 17:32:24 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -95,6 +95,10 @@ class Projections void deleteProjData (void); void init (const int nView, const int nDet); + + // prevent default methods + Projections& operator= (const Projections& rhs); // assignment + Projections(const Projections& rhs); // copy }; diff --git a/include/scanner.h b/include/scanner.h index d383498..39a9002 100644 --- a/include/scanner.h +++ b/include/scanner.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: scanner.h,v 1.4 2000/06/22 10:17:28 kevin Exp $ +** $Id: scanner.h,v 1.5 2000/06/25 17:32:24 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -52,30 +52,37 @@ class DetectorArray { m_viewAngle = viewAngle; } private: - DetectorValue* m_detValues; /* Pointer to array of values recorded by detector */ - int m_nDet; /* Number of detectors in array */ - double m_viewAngle; /* View angle in radians */ + DetectorValue* m_detValues; // Pointer to array of values recorded by detector + int m_nDet; // Number of detectors in array */ + double m_viewAngle; // View angle in radians - DetectorArray& operator=(const DetectorArray& rhs); - DetectorArray& operator()(const DetectorArray& rhs); + DetectorArray& operator=(const DetectorArray& rhs); // assignment + DetectorArray (const DetectorArray& rhs); // copy constructor }; -typedef enum { - DETECTOR_PARALLEL, - DETECTOR_EQUIANGLE, - DETECTOR_EQUILINEAR -} ScannerGeometry; - - class Scanner { public: - Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, int nView, int nSample, const double rot_anglen); + typedef enum { + GEOMETRY_INVALID, + GEOMETRY_PARALLEL, + GEOMETRY_EQUILINEAR, + GEOMETRY_EQUIANGLE + } GeometryID; + + static const char GEOMETRY_PARALLEL_STR[] = "parallel"; + static const char GEOMETRY_EQUILINEAR_STR[] = "equilinear"; + static const char GEOMETRY_EQUIANGLE_STR[] = "equiangle"; + + Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen); ~Scanner(); void collectProjections (Projections& proj, const Phantom& phm, const int start_view, const int trace); void setNView (int nView); + + const bool fail(void) const {return m_fail;} + const string& failMessage(void) const {return m_failMessage;} const unsigned int nDet(void) const {return m_nDet;} const unsigned int nView(void) const {return m_nView;} const double phmLen(void) const {return m_phmLen;} @@ -85,27 +92,19 @@ class Scanner private: - void projectSingleView (const Phantom& phm, DetectorArray& darray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2); - - double projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2); - - double projectLineAgainstPElem (const PhantomElement& pelem, const double x1, const double y1, const double x2, const double y2); - - void traceShowParam (const char *label, const char *fmt, int row, int color, ...); - - - ScannerGeometry m_geometry; /* Geometry of detectory */ - unsigned int m_nDet; /* Number of detectors in array */ - unsigned int m_nView; /* Number of rotated views */ - unsigned int m_nSample; /* Number of rays per detector */ + bool m_fail; + string m_failMessage; + GeometryID m_idGeometry; + unsigned int m_nDet; /* Number of detectors in array */ + unsigned int m_nView; /* Number of rotated views */ + unsigned int m_nSample; /* Number of rays per detector */ double m_detLen; /* Total length of detector array */ double m_rotLen; /* Rotation angle length in radians (norm 2PI) */ double m_detInc; /* Increment between centers of detectors */ double m_rotInc; /* Increment in rotation angle between views */ - double m_radius; /* Radius of rotation. Distance from */ - /* center of phm to center of det */ - double m_phmLen; /* Maximum Length of phantom or area of interest */ - + double m_radius; /* Radius of rotation. Distance from center of phm to center of det */ + double m_phmLen; /* Maximum Length of phantom or area of interest */ + int m_trace; struct { double xd1,yd1,xd2,yd2; /* Coordinates of detector endpoints */ double xs1,ys1,xs2,ys2; /* Coordinates of source endpoints */ @@ -114,7 +113,16 @@ class Scanner static const int N_EXTRA_DETECTORS=4; /* Number of extra detectors widths when calculating detlen */ - int m_trace; + static GeometryID convertGeometryNameToID (const char* const geometryName); + + void projectSingleView (const Phantom& phm, DetectorArray& darray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2); + + double projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2); + + double projectLineAgainstPElem (const PhantomElement& pelem, const double x1, const double y1, const double x2, const double y2); + + void traceShowParam (const char *label, const char *fmt, int row, int color, ...); + }; diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 9e35092..1f6e92c 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $ +** $Id: backprojectors.cpp,v 1.4 2000/06/25 17:32:24 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 @@ -58,11 +58,21 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const { m_nameBackproject = backprojName; m_nameInterpolation = interpName; + m_pBackprojectImplem = NULL; m_idBackproject = convertBackprojectNameToID (backprojName); + if (m_idBackproject == BPROJ_INVALID) { + m_fail = true; + m_failMessage = "Invalid backprojection name "; + m_failMessage += backprojName; + } m_idInterpolation = convertInterpolationNameToID (interpName); - m_pBackprojectImplem = NULL; + if (m_idInterpolation == INTERP_INVALID) { + m_fail = true; + m_failMessage = "Invalid interpolation name "; + m_failMessage += interpName; + } - if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) { + if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) { m_fail = true; return false; } @@ -79,6 +89,7 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation)); else { m_fail = true; + m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; return false; } @@ -139,7 +150,7 @@ Backprojector::convertInterpolationNameToID (const char* const interpName) else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) interpID = INTERP_BSPLINE; #endif - + return (interpID); } diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index c62eb24..e8038aa 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.3 2000/06/22 10:17:28 kevin Exp $ +** $Id: filter.cpp,v 1.4 2000/06/25 17:32:24 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 @@ -44,16 +44,52 @@ * for spatial domain filters. For ANALYTIC solutions, use numint = 0 */ -SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numint) +SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numIntegral = 0) { + m_vecFilter = NULL; m_idFilter = convertFilterNameToID (filterName); + if (m_idFilter == FILTER_INVALID) { + m_fail = true; + m_failMessage = "Invalid Filter name "; + m_failMessage += filterName; + return; + } m_idDomain = convertDomainNameToID (domainName); - init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numint); + if (m_idDomain == DOMAIN_INVALID) { + m_fail = true; + m_failMessage = "Invalid domain name "; + m_failMessage += domainName; + return; + } + init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numIntegral); } -SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint) +SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numIntegral = 0) { - init (filterID, bw, xmin, xmax, n, param, domainID, numint); + init (filterID, bw, xmin, xmax, n, param, domainID, numIntegral); +} + +SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0) +{ + m_bw = bw; + m_nPoints = 0; + m_vecFilter = NULL; + m_filterParam = param; + m_numIntegral = numIntegral; + m_idFilter = convertFilterNameToID (filterName); + if (m_idFilter == FILTER_INVALID) { + m_fail = true; + m_failMessage = "Invalid Filter name "; + m_failMessage += filterName; + return; + } + m_idDomain = convertDomainNameToID (domainName); + if (m_idDomain == DOMAIN_INVALID) { + m_fail = true; + m_failMessage = "Invalid domain name "; + m_failMessage += domainName; + return; + } } void @@ -72,6 +108,8 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax m_nPoints = n; m_xmin = xmin; m_xmax = xmax; + m_numIntegral = numint; + m_filterParam = param; m_vecFilter = new double[n]; double xinc = (m_xmax - m_xmin) / (m_nPoints - 1); @@ -99,7 +137,8 @@ SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax else m_vecFilter[i] = spatialResponseCalc (x, param, numint); } else { - sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", m_idDomain); + m_failMessage = "Illegal domain name "; + m_failMessage += m_idDomain; m_fail = true; } } @@ -113,7 +152,7 @@ SignalFilter::~SignalFilter (void) SignalFilter::FilterID SignalFilter::convertFilterNameToID (const char *filterName) { - FilterID filterID; + FilterID filterID = FILTER_INVALID; if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0) filterID = FILTER_BANDLIMIT; @@ -135,10 +174,6 @@ SignalFilter::convertFilterNameToID (const char *filterName) filterID = FILTER_ABS_COSINE; else if (strcasecmp (filterName, FILTER_SHEPP_STR) == 0) filterID = FILTER_SHEPP; - else { - sys_error(ERR_WARNING, "Invalid filter type %s\n", filterName); - filterID = FILTER_INVALID; - } return (filterID); } @@ -175,14 +210,12 @@ SignalFilter::convertFilterIDToName (const FilterID filterID) const SignalFilter::DomainID SignalFilter::convertDomainNameToID (const char* const domainName) { - DomainID dID; + DomainID dID = DOMAIN_INVALID; if (strcasecmp (domainName, DOMAIN_SPATIAL_STR) == 0) dID = DOMAIN_SPATIAL; else if (strcasecmp (domainName, DOMAIN_FREQ_STR) == 0) dID = DOMAIN_FREQ; - else - dID = DOMAIN_INVALID; return (dID); } @@ -202,20 +235,28 @@ SignalFilter::convertDomainIDToName (const DomainID domain) double -SignalFilter::response (const char* filterName, const char* domainName, double bw, double x, double filt_param) +SignalFilter::response (double x) { double response = 0; - FilterID filterID = convertFilterNameToID (filterName); - DomainID domainID = convertDomainNameToID (domainName); - if (domainID == DOMAIN_SPATIAL) - response = spatialResponseAnalytic (filterID, bw, x, filt_param); - else if (domainID == DOMAIN_FREQ) - response = frequencyResponse (filterID, bw, x, filt_param); + if (m_idDomain == DOMAIN_SPATIAL) + response = spatialResponse (m_idFilter, m_bw, x, m_filterParam, m_numIntegral); + else if (m_idDomain == DOMAIN_FREQ) + response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam); return (response); } + +double +SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param, int nIntegral = 0) +{ + if (nIntegral == 0) + return spatialResponseAnalytic (filterID, bw, x, param); + else + return spatialResponseCalc (filterID, bw, x, param, nIntegral); +} + /* NAME * filter_spatial_response_calc Calculate filter by discrete inverse fourier * transform of filters's frequency @@ -232,9 +273,9 @@ SignalFilter::response (const char* filterName, const char* domainName, double b */ double -SignalFilter::spatialResponseCalc (double x, double param, int n) const +SignalFilter::spatialResponseCalc (double x, double param, int nIntegral) const { - return (spatialResponseCalc (m_idFilter, m_bw, x, param, n)); + return (spatialResponseCalc (m_idFilter, m_bw, x, param, nIntegral)); } double diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index bcb07e2..7f1c6a0 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $ +** $Id: imagefile.cpp,v 1.4 2000/06/25 17:32:24 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 @@ -93,39 +93,39 @@ Array2dFileLabel::getDateString (void) const void -image_filter_response (ImageFile& im, const char* const domainName, double bw, const char* const filterName, double filt_param, const int opt_trace) +ImageFile::filterResponse (const char* const domainName, double bw, const char* const filterName, double filt_param) { - int hx = im.nx() / 2; - int hy = im.ny() / 2; - ImageFileArray v = im.getArray(); + int hx = (mNX - 1) / 2; + int hy = (mNY - 1) / 2; + ImageFileArray v = getArray(); + SignalFilter filter (filterName, domainName, bw, filt_param); for (int i = -hx; i <= hx; i++) { for (int j = -hy; j <= hy; j++) { - double r = sqrt(i * i + j * j); + double r = sqrt (i * i + j * j); - v[i+hx][j+hy] = SignalFilter::response (filterName, domainName, bw, r, filt_param); - if (opt_trace >= TRACE_PHM) - printf ("r=%8.4f, v=%8.4f\n", r, v[i+hx][j+hy]); + v[i+hx][j+hy] = filter.response (r); } } } int -image_display (const ImageFile& im) +ImageFile::display (void) { ImageFileValue pmin, pmax; - im.getPixelValueRange (pmin, pmax); + getPixelValueRange (pmin, pmax); - return (image_display_scale (im, 1, pmin, pmax)); + return (displayScaling (1, pmin, pmax)); } -int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax) +int +ImageFile::displayScaling (const int scale, const ImageFileValue pmin, const ImageFileValue pmax) { int grayscale[256]; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); + int nx = mNX; + int ny = mNY; + ImageFileArray v = getArray(); #if HAVE_G2_H int pens [nx * ny * scale * scale ]; @@ -158,4 +158,148 @@ int image_display_scale (const ImageFile& im, const int scale, const double pmin +// ImageFile::comparativeStatistics Calculate comparative stats +// +// OUTPUT +// d Normalized root mean squared distance measure +// r Normalized mean absolute distance measure +// e Worst case distance measure +// +// REFERENCES +// G.T. Herman, Image Reconstruction From Projections, 1980 + +bool +ImageFile::comparativeStatistics (const ImageFile& imComp, double& d, double& r, double& e) const +{ + if (imComp.nx() != mNX && imComp.ny() != mNY) { + sys_error (ERR_WARNING, "Image sizes differ [ImageFile::comparativeStatistics]"); + return false; + } + ImageFileArray v = getArray(); + ImageFileArray vComp = imComp.getArray(); + + double myMean = 0.; + for (int ix = 0; ix < mNX; ix++) { + for (int iy = 0; iy < mNY; iy++) { + myMean += v[ix][iy]; + } + } + myMean /= (mNX * mNY); + + double sqErrorSum = 0.; + double absErrorSum = 0.; + double sqDiffFromMeanSum = 0.; + double absValueSum = 0.; + for (int ix = 0; ix < mNX; ix++) { + for (int iy = 0; iy < mNY; iy++) { + double diff = v[ix][iy] - vComp[ix][iy]; + sqErrorSum += diff * diff; + absErrorSum += fabs(diff); + double diffFromMean = v[ix][iy] - myMean; + sqDiffFromMeanSum += diffFromMean * diffFromMean; + absValueSum += fabs(v[ix][iy]); + } + } + + d = sqrt (sqErrorSum / sqDiffFromMeanSum); + r = absErrorSum / absValueSum; + + int hx = mNX / 2; + int hy = mNY / 2; + double eMax = -1; + for (int ix = 0; ix < hx; ix++) { + for (int iy = 0; iy < hy; iy++) { + double avgPixel = 0.25 * (v[2*ix][2*iy] + v[2*ix+1][2*iy] + v[2*ix][2*iy+1] + v[2*ix+1][2*iy+1]); + double avgPixelComp = 0.25 * (vComp[2*ix][2*iy] + vComp[2*ix+1][2*iy] + vComp[2*ix][2*iy+1] + vComp[2*ix+1][2*iy+1]); + double error = fabs (avgPixel - avgPixelComp); + if (error > eMax) + eMax = error; + } + } + + e = eMax; + + return true; +} + + +bool +ImageFile::printComparativeStatistics (const ImageFile& imComp, ostream& os) const +{ + double d, r, e; + + if (comparativeStatistics (imComp, d, r, e)) { + os << " Normalized root mean squared distance (d): " << d << endl; + os << " Normalized mean absolute distance (r): " << r << endl; + os << "Worst case distance (2x2 pixel average) (e): " << e << endl; + } +} + + +void +ImageFile::printStatistics (ostream& os) const +{ + double min, max, mean, mode, median, stddev; + + statistics (min, max, mean, mode, median, stddev); + + os << " min: " << min << endl; + os << " max: " << max << endl; + os << " mean: " << mean << endl; + os << " mode: " << mode << endl; + os << "median: " << median << endl; + os << "stddev: " << stddev << endl; +} + + +void +ImageFile::statistics (double& min, double& max, double& mean, double& mode, double& median, double& stddev) const +{ + int nx = mNX; + int ny = mNY; + ImageFileArray v = getArray(); + + mean = 0; + min = v[0][0]; + max = v[0][0]; + for (int ix = 0; ix < nx; ix++) { + for (int iy = 0; iy < ny; iy++) { + if (v[ix][iy] > max) + max = v[ix][iy]; + if (v[ix][iy] < min) + min = v[ix][iy]; + mean += v[ix][iy]; + } + } + mean /= (nx * ny); + + static const int nbin = 256; + int hist[ nbin ] = {0}; + double spread = max - min; + mode = 0; + stddev = 0; + for (int ix = 0; ix < nx; ix++) { + for (int iy = 0; iy < ny; iy++) { + int b = static_cast((((v[ix][iy] - min) / spread) * (nbin - 1)) + 0.5); + hist[b]++; + double diff = (v[ix][iy] - mean); + stddev += diff * diff; + } + } + stddev = sqrt(stddev / (nx * ny)); + + int max_binindex = 0; + int max_bin = -1; + for (int ibin = 0; ibin < nbin; ibin++) { + if (hist[ibin] > max_bin) { + max_bin = hist[ibin]; + max_binindex = ibin; + } + } + + mode = (max_binindex * spread / (nbin - 1)) + min; + + median = 0.; +} + diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp index c080493..aff7b5d 100644 --- a/libctsim/phantom.cpp +++ b/libctsim/phantom.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phantom.cpp,v 1.4 2000/06/22 10:17:28 kevin Exp $ +** $Id: phantom.cpp,v 1.5 2000/06/25 17:32:24 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 @@ -88,7 +88,7 @@ Phantom::convertPhantomIDToName (const PhantomID phmID) Phantom::PhantomID Phantom::convertNameToPhantomID (const char* const phmName) { - PhantomID id; + PhantomID id = PHM_INVALID; if (strcasecmp (phmName, PHM_HERMAN_STR) == 0) id = PHM_HERMAN; @@ -100,8 +100,6 @@ Phantom::convertNameToPhantomID (const char* const phmName) id = PHM_BROWLAND; else if (strcasecmp (phmName, PHM_UNITPULSE_STR) == 0) id = PHM_UNITPULSE; - else - id = PHM_INVALID; return (id); } @@ -111,8 +109,14 @@ bool Phantom::createFromPhantom (const char* const phmName) { PhantomID phmid = convertNameToPhantomID (phmName); - m_name = phmName; + if (phmid == PHM_INVALID) { + m_fail = true; + m_failMessage = "Invalid phantom name "; + m_failMessage += phmName; + return false; + } + m_name = phmName; createFromPhantom (phmid); } @@ -139,11 +143,10 @@ Phantom::createFromPhantom (const PhantomID phmid) addPElem ("ellipse", 0., 0., 1., 1., 0., 1.); // pulse break; default: - sys_error (ERR_WARNING, "Illegal phantom id %d\n", phmid); - m_name += " -- INVALID"; m_fail = true; + m_failMessage = "Illegal phantom id "; + m_failMessage += phmid; return false; - break; } m_id = phmid; diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 5ce0a8e..cd151be 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.4 2000/06/22 10:17:28 kevin Exp $ +** $Id: projections.cpp,v 1.5 2000/06/25 17:32:24 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 @@ -474,8 +474,10 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi double filterMax = detlen; SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0); - if (filter.fail()) - return false; + if (filter.fail()) { + sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", filter.failMessage().c_str()); + return false; + } if (trace) cout << "Reconstruct: filter="<= TRACE_TEXT) diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 1d7eaf7..6ebb04e 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: scanner.cpp,v 1.1 2000/06/19 02:59:34 kevin Exp $ +** $Id: scanner.cpp,v 1.2 2000/06/25 17:32:24 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 @@ -60,10 +60,19 @@ DetectorArray::~DetectorArray (void) * int nSample Number of rays per detector */ -Scanner::Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, int nView, int nSample, const double rot_anglen) +Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen) { m_phmLen = phm.maxAxisLength(); // maximal length along an axis + m_fail = false; + m_idGeometry = convertGeometryNameToID (geometryName); + if (m_idGeometry == GEOMETRY_INVALID) { + m_fail = true; + m_failMessage = "Invalid geometry name "; + m_failMessage += geometryName; + return; + } + if (nView < 1) nView = 1; if (nSample < 1) @@ -73,7 +82,6 @@ Scanner::Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, if ((nDet % 2) == 0) ++nDet; // ensure odd number of detectors - m_geometry = geometry; m_nDet = nDet; m_nView = nView; m_nSample = nSample; @@ -101,6 +109,22 @@ Scanner::~Scanner (void) } +Scanner::GeometryID +Scanner::convertGeometryNameToID (const char* const geometryName) +{ + GeometryID geometryID = GEOMETRY_INVALID; + + if (strcasecmp (geometryName, GEOMETRY_PARALLEL_STR) == 0) + geometryID = GEOMETRY_PARALLEL; + else if (strcasecmp (geometryName, GEOMETRY_EQUILINEAR_STR) == 0) + geometryID = GEOMETRY_EQUILINEAR; + else if (strcasecmp (geometryName, GEOMETRY_EQUIANGLE_STR) == 0) + geometryID = GEOMETRY_EQUIANGLE; + + return (geometryID); +} + + /* NAME * raysum_collect Calculate ray sums for a Phantom diff --git a/src/ctrec.cpp b/src/ctrec.cpp index 39688cf..dc945e3 100644 --- a/src/ctrec.cpp +++ b/src/ctrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctrec.cpp,v 1.13 2000/06/22 10:17:28 kevin Exp $ +** $Id: ctrec.cpp,v 1.14 2000/06/25 17:32:24 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 @@ -103,22 +103,24 @@ static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* i int ctrec_main (int argc, char * argv[]) { - ImageFile *imGlobal = NULL; Projections projGlobal; - char *pj_name, *im_filename = NULL; + ImageFile* imGlobal = NULL; + char* filenameProj = NULL; + char* filenameImage = NULL; string remark; char *endptr; - int opt_verbose = 0; - int opt_debug = 0; - int opt_trace = TRACE_NONE; - double opt_filter_param = 1; - string optFilterName = "abs_bandlimit"; - string optInterpName = "linear"; - string optBackprojName = "idiff2"; - int opt_interp_param = 1; + int optVerbose = 0; + int optDebug = 0; + int optTrace = TRACE_NONE; + double optFilterParam = -1; + string optFilterName = SignalFilter::FILTER_ABS_BANDLIMIT_STR; + string optInterpName = Backprojector::INTERP_LINEAR_STR; + string optBackprojName = Backprojector::BPROJ_IDIFF2_STR; + // string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR; + int optInterpParam = 1; int nx, ny; #ifdef HAVE_MPI - ImageFile *imLocal; + ImageFile* imLocal; int mpi_nview, mpi_ndet; double mpi_detinc, mpi_rotinc, mpi_phmlen; MPIWorld mpiWorld (argc, argv); @@ -148,19 +150,19 @@ ctrec_main (int argc, char * argv[]) optBackprojName = optarg; break; case O_FILTER_PARAM: - opt_filter_param = strtod(optarg, &endptr); + optFilterParam = strtod(optarg, &endptr); if (endptr != optarg + strlen(optarg)) { ctrec_usage(argv[0]); } break; case O_VERBOSE: - opt_verbose = 1; + optVerbose = 1; break; case O_DEBUG: - opt_debug = 1; + optDebug = 1; break; case O_TRACE: - if ((opt_trace = opt_set_trace(optarg)) < 0) { + if ((optTrace = opt_set_trace(optarg)) < 0) { ctrec_usage(argv[0]); return (1); } @@ -187,16 +189,16 @@ ctrec_main (int argc, char * argv[]) return (1); } - pj_name = argv[optind]; + filenameProj = argv[optind]; - im_filename = argv[optind + 1]; + filenameImage = argv[optind + 1]; nx = strtol(argv[optind + 2], &endptr, 10); ny = strtol(argv[optind + 3], &endptr, 10); ostringstream filterDesc; - if (opt_filter_param >= 0) - filterDesc << optFilterName << ": alpha=" << opt_filter_param; + if (optFilterParam >= 0) + filterDesc << optFilterName << ": alpha=" << optFilterParam; else filterDesc << optFilterName; @@ -204,7 +206,7 @@ ctrec_main (int argc, char * argv[]) label << "Reconstruct: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName; remark = label.str(); - if (opt_verbose) + if (optVerbose) cout << "Remark: " << remark << endl; #ifdef HAVE_MPI } @@ -212,8 +214,8 @@ ctrec_main (int argc, char * argv[]) #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { - projGlobal.read (pj_name); - if (opt_verbose) + projGlobal.read (filenameProj); + if (optVerbose) projGlobal.printScanInfo(); mpi_ndet = projGlobal.nDet(); @@ -227,11 +229,11 @@ ctrec_main (int argc, char * argv[]) mpiWorld.BcastString (optBackprojName); mpiWorld.BcastString (optFilterName); mpiWorld.BcastString (optInterpName); - mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0); - mpiWorld.getComm().Bcast (&opt_interp_param, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&optVerbose, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&optDebug, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0); + mpiWorld.getComm().Bcast (&optInterpParam, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0); @@ -239,7 +241,7 @@ ctrec_main (int argc, char * argv[]) mpiWorld.getComm().Bcast (&mpi_rotinc, 1, MPI::DOUBLE, 0); mpiWorld.getComm().Bcast (&nx, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&ny, 1, MPI::INT, 0); - if (opt_verbose) + if (optVerbose) timerBcast.timerEndAndReport ("Time to broadcast variables"); mpiWorld.setTotalWorkUnits (mpi_nview); @@ -250,37 +252,37 @@ ctrec_main (int argc, char * argv[]) projLocal.setRotInc (mpi_rotinc); TimerCollectiveMPI timerScatter (mpiWorld.getComm()); - ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, opt_debug); - if (opt_verbose) + ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, optDebug); + if (optVerbose) timerScatter.timerEndAndReport ("Time to scatter projections"); if (mpiWorld.getRank() == 0) { - imGlobal = new ImageFile (im_filename, nx, ny); + imGlobal = new ImageFile (filenameImage, nx, ny); imGlobal->fileCreate(); } imLocal = new ImageFile (nx, ny); #else - projGlobal.read (pj_name); - if (opt_verbose) + projGlobal.read (filenameProj); + if (optVerbose) projGlobal.printScanInfo(); - imGlobal = new ImageFile (im_filename, nx, ny); + imGlobal = new ImageFile (filenameImage, nx, ny); imGlobal->fileCreate(); #endif #ifdef HAVE_MPI TimerCollectiveMPI timerReconstruct (mpiWorld.getComm()); - projLocal.reconstruct (*imLocal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace); - if (opt_verbose) + projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace); + if (optVerbose) timerReconstruct.timerEndAndReport ("Time to reconstruct"); TimerCollectiveMPI timerReduce (mpiWorld.getComm()); ReduceImageMPI (mpiWorld, imLocal, imGlobal); - if (opt_verbose) + if (optVerbose) timerReduce.timerEndAndReport ("Time to reduce image"); #else - projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace); + projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace); #endif #ifdef HAVE_MPI @@ -292,7 +294,7 @@ ctrec_main (int argc, char * argv[]) imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, projGlobal.remark(), projGlobal.calcTime()); imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime); imGlobal->fileClose (); - if (opt_verbose) + if (optVerbose) cout << "Run time: " << calcTime << " seconds" << endl; } #ifdef HAVE_MPI @@ -309,7 +311,7 @@ ctrec_main (int argc, char * argv[]) ////////////////////////////////////////////////////////////////////////////////////// #ifdef HAVE_MPI -static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int opt_debug) +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int optDebug) { if (mpiWorld.getRank() == 0) { for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) { diff --git a/src/if-2.cpp b/src/if-2.cpp index 0d35a12..951b431 100644 --- a/src/if-2.cpp +++ b/src/if-2.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if-2.cpp,v 1.6 2000/06/19 17:58:13 kevin Exp $ +** $Id: if-2.cpp,v 1.7 2000/06/25 17:32:24 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 @@ -195,11 +195,11 @@ if2_main (int argc, char *const argv[]) for (int iy = 0; iy < im_in1.ny(); iy++) { double diff = *in1++ - *in2++; *out++ = diff; - abs_error += fabs(diff); } } - abs_error /= (im_in1.nx() * im_in1.ny()); - cout << "Average Error: " << abs_error << endl; + double d, r, e; + im_in1.comparativeStatistics (im_in2, d, r, e); + cout << "d=" << d << ", r=" << r << ", e=" << e << endl; } im_out.arrayDataWrite(); diff --git a/src/ifinfo.cpp b/src/ifinfo.cpp index c911aa5..f6c8da7 100644 --- a/src/ifinfo.cpp +++ b/src/ifinfo.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ifinfo.cpp,v 1.9 2000/06/19 17:58:13 kevin Exp $ +** $Id: ifinfo.cpp,v 1.10 2000/06/25 17:32:24 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 @@ -49,7 +49,7 @@ static struct option my_options[] = void ifinfo_usage (const char *program) { - cout << "usage: " << fileBasename(program) << " infile [OPTIONS]" << endl; + cout << "usage: " << fileBasename(program) << " image1 [image2] [OPTIONS]" << endl; cout << "Imagefile information" << endl; cout << endl; cout << " infile Name of input IF file" << endl; @@ -67,8 +67,10 @@ ifinfo_usage (const char *program) int ifinfo_main (int argc, char *const argv[]) { - ImageFile *im; - char *in_file; + ImageFile *im = NULL; + ImageFile* im2 = NULL; + string in_file; + string in2_file; int opt_verbose = 0; int opt_stats = 1; int opt_labels = 1; @@ -118,116 +120,53 @@ ifinfo_main (int argc, char *const argv[]) } } - if (optind + 1 != argc) - { - ifinfo_usage(argv[0]); - return (1); - } + if (optind + 2 == argc) { + in2_file = argv [optind+1]; + } else if (optind + 1 != argc) { + ifinfo_usage (argv[0]); + return (1); + } in_file = argv[optind]; - im = new ImageFile (in_file); + im = new ImageFile (in_file.c_str()); if (! im->fileRead ()) { - sys_error (ERR_WARNING, "Unable to read file %s", in_file); + sys_error (ERR_WARNING, "Unable to read file %s", in_file.c_str()); return (1); } + if (in2_file != "") { + im2 = new ImageFile(in2_file.c_str()); + if (! im2->fileRead ()) { + sys_error (ERR_WARNING, "Unable to read file %s", in2_file.c_str()); + return (1); + } + } + + if (opt_stats) + cout << "Image size: (" << im->nx() << "," << im->ny() << ")" << endl << endl; if (opt_labels) - { - int nlabels = im->getNumLabels(); - int i; + im->printLabels (cout); - for (i = 0; i < nlabels; i++) - { - Array2dFileLabel label; - im->labelRead (label, i); - - if (label.getLabelType() == Array2dFileLabel::L_HISTORY) { - cout << "History: " << endl; - cout << " " << label.getLabelString() << endl; - cout << " calc time = " << label.getCalcTime() << " secs" << endl; - cout << " Timestamp = " << label.getDateString() << endl; - } else if (label.getLabelType() == Array2dFileLabel::L_USER) { - cout << "Note: " << label.getLabelString() << endl; - cout << " Timestamp = %s" << label.getDateString() << endl; - } - cout << endl; - } - } + if (opt_stats) { + if (im2) + cout << "Image 1" << endl; + + im->printStatistics (cout); - if (opt_stats) - { - double minfound = HUGE_VAL, maxfound = -HUGE_VAL; - double mode = 0, mean = 0, stddev = 0; - double spread; - int hist[256]; - int ibin, nbin = 256; - int max_bin, max_binindex; - double maxbin; - int ix, iy; - - maxbin = nbin - 1; - ImageFileArray v = im->getArray(); - int nx = im->nx(); - int ny = im->ny(); - - for (ix = 0; ix < nx; ix++) - { - for (iy = 0; iy < ny; iy++) - { - if (v[ix][iy] > maxfound) - maxfound = v[ix][iy]; - if (v[ix][iy] < minfound) - minfound = v[ix][iy]; - mean += v[ix][iy]; - } - } - spread = maxfound - minfound; - if (spread == 0) - mode = minfound; - else - { - for (ibin = 0; ibin < nbin; ibin++) - hist[ibin] = 0; - for (ix = 0; ix < nx; ix++) - { - for (iy = 0; iy < ny; iy++) - { - int b = (int) ((((v[ix][iy] - minfound) / spread) * (double) maxbin) + 0.5); - hist[b]++; - } - } - max_binindex = 0; - max_bin = -1; - for (ibin = 0; ibin < nbin; ibin++) - { - if (hist[ibin] > max_bin) - { - max_bin = hist[ibin]; - max_binindex = ibin; - } - } - mode = (((double) max_binindex) * spread / ((double) maxbin)) + minfound; - } + if (im2) { + cout << endl; - mean /= nx * ny; - for (ix = 0; ix < nx; ix++) - { - for (iy = 0; iy < ny; iy++) - { - double diff = (v[ix][iy] - mean); - stddev += diff * diff; - } - } - stddev = sqrt(stddev / (nx * ny)); - cout << " nx: " << nx << endl; - cout << " nx: " << ny << endl; - cout << " min: " << minfound << endl; - cout << " max: " << maxfound << endl; - cout << " mean: " << mean << endl; - cout << " mode: " << mode << endl; - cout << "stddef: " << stddev << endl; + if (opt_labels) + im2->printLabels(cout); + + cout << "Image 2" << endl; + im2->printStatistics (cout); + cout << endl; + + im->printComparativeStatistics (*im2, cout); } + } return (0); } diff --git a/src/phm2if.cpp b/src/phm2if.cpp index 75119c9..91d0521 100644 --- a/src/phm2if.cpp +++ b/src/phm2if.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phm2if.cpp,v 1.13 2000/06/22 10:17:28 kevin Exp $ +** $Id: phm2if.cpp,v 1.14 2000/06/25 17:32:24 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 @@ -108,13 +108,13 @@ phm2if_main (int argc, char* argv[]) Phantom phm; int opt_nx = 0, opt_ny = 0; int opt_nsample = 1; - string optPhmName; + string optPhmName = Phantom::PHM_HERMAN_STR; string optFilterName; - string optDomainName = "spatial"; + string optDomainName = SignalFilter::DOMAIN_SPATIAL_STR; char *opt_outfile = NULL; int opt_debug = 0; string opt_desc; - string opt_phmfilename; + string opt_phmFileName; char *endstr, *endptr; double opt_filter_param = 1; double opt_filter_bw = 1.; @@ -141,20 +141,9 @@ phm2if_main (int argc, char* argv[]) switch (c) { case O_PHANTOM: optPhmName = optarg; - if (! phm.createFromPhantom (optPhmName.c_str())) { - cout << "Invalid phantom name " << optPhmName << endl << endl; - phm2if_usage(argv[0]); - return (1); - } break; case O_PHMFILE: - opt_phmfilename = optarg; - phm.createFromFile(opt_phmfilename.c_str()); -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) - cerr << "Can't use phantom from file in MPI mode" << endl; - return (1); -#endif + opt_phmFileName = optarg; break; case O_VERBOSE: opt_verbose = true; @@ -220,8 +209,8 @@ phm2if_main (int argc, char* argv[]) } } - if (phm.nPElem() == 0 && optPhmName == "" && optFilterName == "") { - cerr << "No phantom defined" << endl; + if (optPhmName == "" && optFilterName == "" && opt_phmFileName == "") { + cerr << "No phantom defined" << endl << endl; phm2if_usage(argv[0]); return (1); } @@ -248,17 +237,35 @@ phm2if_main (int argc, char* argv[]) ostringstream oss; oss << "nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", "; - if (opt_phmfilename.length()) - oss << "phantom=" << opt_phmfilename; + if (opt_phmFileName != "") + oss << "phantom=" << opt_phmFileName; else if (optPhmName != "") oss << "phantom=" << optPhmName; else if (optFilterName != "") { oss << "filter=" << optFilterName << " - " << optDomainName; } - if (opt_desc.length()) + if (opt_desc != "") oss << ": " << opt_desc; opt_desc = oss.str(); + if (optPhmName != "") { + phm.createFromPhantom (optPhmName.c_str()); + if (phm.fail()) { + cout << phm.failMessage() << endl << endl; + phm2if_usage(argv[0]); + return (1); + } + } + + if (opt_phmFileName != "") { + phm.createFromFile(opt_phmFileName.c_str()); +#ifdef HAVE_MPI + if (mpiWorld.getRank() == 0) + cerr << "Can't use phantom from file in MPI mode" << endl; + return (1); +#endif + } + if (opt_verbose) cout << "Rasterize Phantom to Image" << endl << endl; #ifdef HAVE_MPI @@ -309,7 +316,7 @@ phm2if_main (int argc, char* argv[]) } } else if (optFilterName != "") { if (mpiWorld.getRank() == 0) { - image_filter_response (*imGlobal, optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param, opt_trace); + imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param); } } else { TimerCollectiveMPI timerRasterize (mpiWorld.getComm()); @@ -327,7 +334,7 @@ phm2if_main (int argc, char* argv[]) if (phm.getComposition() == P_UNIT_PULSE) { v[opt_nx/2][opt_ny/2] = 1.; } else if (optFilterName != "") { - image_filter_response (*imGlobal, optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param, opt_trace); + imGlobal->filterResponse (optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param); } else { #if HAVE_SGP if (opt_trace >= TRACE_PHM) @@ -356,7 +363,7 @@ phm2if_main (int argc, char* argv[]) scanf ("%d", &nscale); printf ("Enter minimum and maximum densities (min, max): "); scanf ("%lf %lf", &dmin, &dmax); - image_display_scale (*imGlobal, nscale, dmin, dmax); + imGlobal->displayScaling (nscale, dmin, dmax); } } diff --git a/src/phm2pj.cpp b/src/phm2pj.cpp index a4510c4..4ed549e 100644 --- a/src/phm2pj.cpp +++ b/src/phm2pj.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: phm2pj.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $ +** $Id: phm2pj.cpp,v 1.4 2000/06/25 17:32:24 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 @@ -88,15 +88,15 @@ int phm2pj_main (int argc, char* argv[]) { Phantom phm; - ScannerGeometry opt_geometry = DETECTOR_PARALLEL; + string optGeometryName = Scanner::GEOMETRY_PARALLEL_STR; char *opt_outfile = NULL; string opt_desc; - string opt_phmfilename; + string optPhmFileName; int opt_ndet; int opt_nview; int opt_nray = 1; int opt_trace = 0; - string optPhmName; + string optPhmName = Phantom::PHM_HERMAN_STR; int opt_verbose = 0; int opt_debug = 0; double opt_rotangle = 1; @@ -121,20 +121,9 @@ phm2pj_main (int argc, char* argv[]) switch (c) { case O_PHANTOM: optPhmName = optarg; - if (! phm.createFromPhantom (optPhmName.c_str())) { - cout << "ERROR: Invalid phantom name " << optPhmName << endl << endl; - phm2pj_usage(argv[0]); - return (1); - } break; case O_PHMFILE: -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) - cerr << "Can not read phantom from file in MPI mode" << endl; - return (1); -#endif - opt_phmfilename = optarg; - phm.createFromFile (opt_phmfilename.c_str()); + optPhmFileName = optarg; break; case O_VERBOSE: opt_verbose = 1; @@ -187,8 +176,8 @@ phm2pj_main (int argc, char* argv[]) } } - if (phm.nPElem() == 0) { - cerr << "No phantom defined" << endl; + if (optPhmName == "" && optPhmFileName == "") { + cerr << "No phantom defined" << endl << endl; phm2pj_usage(argv[0]); return (1); } @@ -215,8 +204,8 @@ phm2pj_main (int argc, char* argv[]) ostringstream desc; desc << "Raysum_Collect: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", "; - if (opt_phmfilename.length()) { - desc << "PhantomFile=" << opt_phmfilename; + if (optPhmFileName.length()) { + desc << "PhantomFile=" << optPhmFileName; } else if (optPhmName != "") { desc << "Phantom=" << optPhmName; } @@ -224,6 +213,24 @@ phm2pj_main (int argc, char* argv[]) desc << ": " << opt_desc; } opt_desc = desc.str(); + + if (optPhmName != "") { + phm.createFromPhantom (optPhmName.c_str()); + if (phm.fail()) { + cout << phm.failMessage() << endl << endl; + phm2pj_usage(argv[0]); + return (1); + } + } + + if (optPhmFileName != "") { +#ifdef HAVE_MPI + cerr << "Can not read phantom from file in MPI mode" << endl; + return (1); +#endif + phm.createFromFile (optPhmFileName.c_str()); + } + #ifdef HAVE_MPI } #endif @@ -245,8 +252,11 @@ phm2pj_main (int argc, char* argv[]) #endif opt_rotangle *= PI; - Scanner scanner (phm, opt_geometry, opt_ndet, opt_nview, opt_nray, opt_rotangle); - + Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_nray, opt_rotangle); + if (scanner.fail()) { + cout << "Scanner Creation Error: " << scanner.failMessage() << endl; + return (1); + } #ifdef HAVE_MPI mpiWorld.setTotalWorkUnits (opt_nview); -- 2.34.1