From 30e455abcd8cac05ce7afe43216ec9e26342e1cf Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Tue, 11 Jul 2000 10:32:44 +0000 Subject: [PATCH] r142: *** empty log message *** --- ChangeLog | 13 +- configure | 133 ++++++------- configure.in | 7 +- include/backprojectors.h | 33 ++-- include/filter.h | 31 +-- include/imagefile.h | 10 +- libctgraphics/sgp.cpp | 9 +- libctsim/backprojectors.cpp | 75 ++++---- libctsim/filter.cpp | 137 ++++++++------ libctsim/imagefile.cpp | 193 ++++++++++++++++++- libctsim/projections.cpp | 15 +- src/Makefile.am | 6 +- src/ctsim.cpp | 33 ---- src/if2img.cpp | 363 ++++++------------------------------ src/pjrec.cpp | 26 ++- 15 files changed, 530 insertions(+), 554 deletions(-) delete mode 100644 src/ctsim.cpp diff --git a/ChangeLog b/ChangeLog index 055fcff..66a5fff 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,15 @@ -2.0.0-b3 - 7/09/00 +2.0.0-b3 - 7/10/00 Added highly optimized backprojection method idiff3 Moved comparative stats to if-2 program from ifinfo - -2.0.0-b2 - 7/07/00 + Fixed image display in if2img using G2 library + Added column-plot and row-plot options to if-2 + Added autoselection of analytic/calculated spatial responses to SignalFilter + Added frequency-based preinterpolation to SignalFilter and idiff3 + backprojection method. Currently, this technique is still under + development and debugging + Moved graphic file writing to ImageFile class from if2img program + +2.0.0-b2 - 7/07/00 Cleaned up SignalFilter class Added zeropad option to pjrec Added zeropad options to html and cgi files diff --git a/configure b/configure index 2b9a3c6..145a261 100755 --- a/configure +++ b/configure @@ -2803,13 +2803,12 @@ fi fi - if test -n "$GCC"; then CFLAGS="$CFLAGS -Wall" fi echo $ac_n "checking whether to enable verbose warnings""... $ac_c" 1>&6 -echo "configure:2813: checking whether to enable verbose warnings" >&5 +echo "configure:2812: 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" @@ -2859,7 +2858,7 @@ fi echo $ac_n "checking debug""... $ac_c" 1>&6 -echo "configure:2863: checking debug" >&5 +echo "configure:2862: checking debug" >&5 # Check whether --enable-debug or --disable-debug was given. if test "${enable_debug+set}" = set; then enableval="$enable_debug" @@ -2888,13 +2887,21 @@ else fi if test "$debug" = "true" ; then - CFLAGS="-g -DDEBUG" + + if test -n "$GCC"; then + CFLAGS="$CFLAGS -g -DDEBUG" + fi + cat >> confdefs.h <<\EOF #define DEBUG 1 EOF else - CFLAGS="-g -O3 -DNDEBUG" + + if test -n "$GCC"; then + CFLAGS="$CFLAGS -g -O3 -DNDEBUG --fast-math" + fi + cat >> confdefs.h <<\EOF #define NDEBUG 1 EOF @@ -2917,7 +2924,7 @@ if test "$withval" != "no" ; then fi echo $ac_n "checking for LAM MPI installation""... $ac_c" 1>&6 -echo "configure:2921: checking for LAM MPI installation" >&5 +echo "configure:2928: 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" @@ -2935,7 +2942,7 @@ else fi echo $ac_n "checking for web access""... $ac_c" 1>&6 -echo "configure:2939: checking for web access" >&5 +echo "configure:2946: 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" @@ -3025,7 +3032,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:3029: checking for X" >&5 +echo "configure:3036: checking for X" >&5 # Check whether --with-x or --without-x was given. if test "${with_x+set}" = set; then @@ -3087,12 +3094,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:3096: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:3103: \"$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* @@ -3161,14 +3168,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:3179: \"$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. @@ -3274,17 +3281,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:3278: checking whether -R must be followed by a space" >&5 +echo "configure:3285: 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:3295: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_nospace=yes else @@ -3300,14 +3307,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:3318: \"$ac_link\") 1>&5; (eval $ac_link) 2>&5; } && test -s conftest${ac_exeext}; then rm -rf conftest* ac_R_space=yes else @@ -3339,7 +3346,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:3343: checking for dnet_ntoa in -ldnet" >&5 +echo "configure:3350: 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 @@ -3347,7 +3354,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:3369: \"$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 @@ -3380,7 +3387,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:3384: checking for dnet_ntoa in -ldnet_stub" >&5 +echo "configure:3391: 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 @@ -3388,7 +3395,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:3410: \"$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 @@ -3428,12 +3435,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:3432: checking for gethostbyname" >&5 +echo "configure:3439: 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:3467: \"$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 @@ -3477,7 +3484,7 @@ fi if test $ac_cv_func_gethostbyname = no; then echo $ac_n "checking for gethostbyname in -lnsl""... $ac_c" 1>&6 -echo "configure:3481: checking for gethostbyname in -lnsl" >&5 +echo "configure:3488: 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 @@ -3485,7 +3492,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:3507: \"$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 +3533,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:3530: checking for connect" >&5 +echo "configure:3537: 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:3565: \"$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 @@ -3575,7 +3582,7 @@ fi if test $ac_cv_func_connect = no; then echo $ac_n "checking for connect in -lsocket""... $ac_c" 1>&6 -echo "configure:3579: checking for connect in -lsocket" >&5 +echo "configure:3586: 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 @@ -3583,7 +3590,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:3605: \"$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 +3625,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:3622: checking for remove" >&5 +echo "configure:3629: 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:3657: \"$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 @@ -3667,7 +3674,7 @@ fi if test $ac_cv_func_remove = no; then echo $ac_n "checking for remove in -lposix""... $ac_c" 1>&6 -echo "configure:3671: checking for remove in -lposix" >&5 +echo "configure:3678: 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 @@ -3675,7 +3682,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:3697: \"$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 @@ -3710,12 +3717,12 @@ fi # BSDI BSD/OS 2.1 needs -lipc for XOpenDisplay. echo $ac_n "checking for shmat""... $ac_c" 1>&6 -echo "configure:3714: checking for shmat" >&5 +echo "configure:3721: 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:3749: \"$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 @@ -3759,7 +3766,7 @@ fi if test $ac_cv_func_shmat = no; then echo $ac_n "checking for shmat in -lipc""... $ac_c" 1>&6 -echo "configure:3763: checking for shmat in -lipc" >&5 +echo "configure:3770: 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 @@ -3767,7 +3774,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:3789: \"$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 @@ -3811,7 +3818,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:3815: checking for IceConnectionNumber in -lICE" >&5 +echo "configure:3822: 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 @@ -3819,7 +3826,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:3841: \"$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 @@ -3867,7 +3874,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:3871: checking for XtToolkitThreadInitialize in -lXt" >&5 +echo "configure:3878: 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 @@ -3875,7 +3882,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:3897: \"$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,7 +3921,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:3918: checking for XdbeQueryExtension in -lXext" >&5 +echo "configure:3925: 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 @@ -3922,7 +3929,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:3944: \"$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 @@ -3957,7 +3964,7 @@ else fi echo $ac_n "checking for XmbufQueryExtension in -lXext""... $ac_c" 1>&6 -echo "configure:3961: checking for XmbufQueryExtension in -lXext" >&5 +echo "configure:3968: 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 @@ -3965,7 +3972,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:3987: \"$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 @@ -4006,17 +4013,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:4010: checking for $ac_hdr" >&5 +echo "configure:4017: 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:4020: \"$ac_try\") 1>&5; (eval $ac_try) 2>&5; } +{ (eval echo configure:4027: \"$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* @@ -4058,7 +4065,7 @@ my_includes="$my_includes -I../include -I.." echo $ac_n "checking interactive graphics""... $ac_c" 1>&6 -echo "configure:4062: checking interactive graphics" >&5 +echo "configure:4069: 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 20cc064..e3ea013 100644 --- a/configure.in +++ b/configure.in @@ -12,7 +12,7 @@ AC_PROG_AWK AC_PROG_INSTALL AC_PROG_MAKE_SET AC_PROG_RANLIB -AC_PROG_CC([-g]) +AC_PROG_CC AC_PROG_CXX @@ -114,7 +114,6 @@ if test "${getopt_long}" = "false" ; then AM_CONDITIONAL(INCLUDED_GETOPT_LONG, test 1==1) fi - AC_ADD_GCC_CFLAGS([-Wall]) AC_MSG_CHECKING(whether to enable verbose warnings) AC_ARG_ENABLE(verbose-warnings, @@ -155,10 +154,10 @@ esac], AM_CONDITIONAL(DEBUG, test "$debug" = "true") if test "$debug" = "true" ; then - CFLAGS="-g -DDEBUG" + AC_ADD_GCC_CFLAGS([-g -DDEBUG]) AC_DEFINE(DEBUG) else - CFLAGS="-g -O3 -DNDEBUG" + AC_ADD_GCC_CFLAGS([-g -O3 -DNDEBUG --fast-math]) AC_DEFINE(NDEBUG) fi diff --git a/include/backprojectors.h b/include/backprojectors.h index f0cc19a..ab15930 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.7 2000/07/09 08:16:17 kevin Exp $ +** $Id: backprojectors.h,v 1.8 2000/07/11 10:32:44 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 @@ -53,7 +53,7 @@ class Backprojector I_3BSPLINE, #endif INTERP_LINEAR, // Linear interpolation - INTERP_FREQ_PREINTERPOLATE, + INTERP_FREQ_PREINTERPOLATION, } InterpolationID; static const char BPROJ_TRIG_STR[]= "trig"; @@ -66,9 +66,9 @@ class Backprojector static const char INTERP_NEAREST_STR[]= "nearest"; static const char INTERP_LINEAR_STR[]= "linear"; static const char INTERP_BSPLINE_STR[]= "bspline"; - static const char INTERP_FREQ_PREINTERPOLATE_STR[]= "freq_preinterpolate"; + static const char INTERP_FREQ_PREINTERPOLATION_STR[]= "freq_preinterpolation"; - Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName); + Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor); ~Backprojector (void); @@ -91,14 +91,14 @@ class Backprojector static const BackprojectID convertBackprojectNameToID (const char* const bprojName); static const char* convertBackprojectIDToName (const BackprojectID bprojID); - bool initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName); + bool initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor); }; class Backproject { public: - Backproject (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType); + Backproject (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor); virtual ~Backproject (); @@ -121,6 +121,7 @@ class Backproject int nDet; double xMin, xMax, yMin, yMax; // Retangular coords of phantom double xInc, yInc; // size of cells + int m_interpFactor; private: Backproject (const Backproject& rhs); @@ -131,8 +132,8 @@ class Backproject class BackprojectTrig : public Backproject { public: - BackprojectTrig (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : Backproject::Backproject (proj, im, interpType) + BackprojectTrig (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : Backproject::Backproject (proj, im, interpType, interpFactor) {} void BackprojectView (const double* const t, double view_angle); @@ -142,7 +143,7 @@ class BackprojectTrig : public Backproject class BackprojectTable : public Backproject { public: - BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType); + BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor); ~BackprojectTable (); void BackprojectView (const double* const t, double view_angle); @@ -158,7 +159,7 @@ class BackprojectTable : public Backproject class BackprojectDiff : public Backproject { public: - BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType); + BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor); ~BackprojectDiff (); void BackprojectView (const double* const t, double view_angle); @@ -172,8 +173,8 @@ class BackprojectDiff : public Backproject class BackprojectDiff2 : public BackprojectDiff { public: - BackprojectDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : BackprojectDiff::BackprojectDiff (proj, im, interpType) + BackprojectDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : BackprojectDiff::BackprojectDiff (proj, im, interpType, interpFactor) {} void BackprojectView (const double* const t, double view_angle); @@ -182,8 +183,8 @@ class BackprojectDiff2 : public BackprojectDiff class BackprojectIntDiff2 : public BackprojectDiff { public: - BackprojectIntDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : BackprojectDiff::BackprojectDiff (proj, im, interpType) + BackprojectIntDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : BackprojectDiff::BackprojectDiff (proj, im, interpType, interpFactor) {} void BackprojectView (const double* const t, double view_angle); @@ -193,8 +194,8 @@ class BackprojectIntDiff2 : public BackprojectDiff class BackprojectIntDiff3 : public BackprojectDiff { public: - BackprojectIntDiff3 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : BackprojectDiff::BackprojectDiff (proj, im, interpType) + BackprojectIntDiff3 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : BackprojectDiff::BackprojectDiff (proj, im, interpType, interpFactor) {} void BackprojectView (const double* const t, double view_angle); diff --git a/include/filter.h b/include/filter.h index 9ee7390..4a25b5b 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.12 2000/07/07 15:30:59 kevin Exp $ +** $Id: filter.h,v 1.13 2000/07/11 10:32:44 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,11 +88,11 @@ class SignalFilter { static const char DOMAIN_SPATIAL_STR[]="spatial"; - SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const int zeropad = 0, const int numIntegral = 0); + SignalFilter (const char* filterName, const char* filterMethodName,double bw, double signalIncrement, int n, double param, const char* domainName, const int zeropad = 0, const int preinterpolationFactor = 1); - SignalFilter (const FilterID filt_type, FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad = 0, const int numIntegral = 0); + SignalFilter (const FilterID filt_type, FilterMethodID filterMethodID, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad = 0, const int preinterpolationFactor = 1); - SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0); + SignalFilter (const char* filterName, const char* domainName, double bw, double param); ~SignalFilter (void); @@ -132,14 +132,16 @@ class SignalFilter { double response (double x); - static double spatialResponse (FilterID fType, double bw, double x, double param, int nIntegral = 0); + static double spatialResponse (FilterID fType, double bw, double x, double param); 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 spatialResponseCalc (FilterID fType, double bw, double x, double param, int nIntegral); + + static void setNumIntegral(int nIntegral) {N_INTEGRAL = nIntegral;} + private: double m_bw; int m_nFilterPoints; @@ -153,9 +155,9 @@ class SignalFilter { double* m_vecFourierSinTable; complex* m_complexVecFilter; #ifdef HAVE_FFTW - fftw_real* m_vecRealFftInput; + fftw_real* m_vecRealFftInput, *m_vecRealFftSignal; rfftw_plan m_realPlanForward, m_realPlanBackward; - fftw_complex* m_vecComplexFftInput; + fftw_complex* m_vecComplexFftInput, *m_vecComplexFftSignal; fftw_plan m_complexPlanForward, m_complexPlanBackward; #endif @@ -168,10 +170,14 @@ class SignalFilter { FilterMethodID m_idFilterMethod; DomainID m_idDomain; double m_filterParam; - int m_numIntegral; int m_traceLevel; int m_zeropad; + int m_nOutputPoints; + int m_preinterpolationFactor; + static int N_INTEGRAL; + + static const bool haveAnalyticSpatial (const FilterID filterID); static const FilterID convertFilterNameToID (const char* filterName); static const char* convertFilterIDToName (const FilterID filterID); static const FilterMethodID convertFilterMethodNameToID (const char* filterMethodName); @@ -179,9 +185,9 @@ class SignalFilter { static const DomainID convertDomainNameToID (const char* domainName); static const char* convertDomainIDToName (const DomainID domainID); - void init (const FilterID filt_type, const FilterMethodID filterMethod, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad, const int numInt); + void init (const FilterID filt_type, const FilterMethodID filterMethod, double bw, double signalIncrement, int n, double param, const DomainID domain, const int zeropad, const int preInterpScale); -double spatialResponseCalc (double x, double param, int n) const; + double spatialResponseCalc (double x, double param) const; double spatialResponseAnalytic (double x, double param) const; @@ -194,5 +200,4 @@ double spatialResponseCalc (double x, double param, int n) const; }; - #endif diff --git a/include/imagefile.h b/include/imagefile.h index ef60d3d..5f8c324 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: imagefile.h,v 1.16 2000/06/28 15:25:34 kevin Exp $ +** $Id: imagefile.h,v 1.17 2000/07/11 10:32:44 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 @@ -147,6 +147,14 @@ class ImageFile : public ImageFileBase int displayScaling (const int scaleFactor, ImageFileValue pmin, ImageFileValue pmax); +#if HAVE_PNG + void writeImagePNG (const char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax); +#endif +#if HAVE_GD + void writeImageGIF (const char *outfile, int nxcell, int nycell, double densmin, double densmax); +#endif + void writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax); + void writeImagePGMASCII (const char *outfile, int nxcell, int nycell, double densmin, double densmax); }; diff --git a/libctgraphics/sgp.cpp b/libctgraphics/sgp.cpp index cc87671..8a5400f 100644 --- a/libctgraphics/sgp.cpp +++ b/libctgraphics/sgp.cpp @@ -7,7 +7,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: sgp.cpp,v 1.3 2000/06/19 19:16:17 kevin Exp $ +** $Id: sgp.cpp,v 1.4 2000/07/11 10:32:44 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 @@ -84,12 +84,13 @@ void sgp2_close (SGP_ID gid) { #if HAVE_G2_H + if (gid) g2_close (gid->g2_id); #endif - if (gid == _sgp2_cwin) - _sgp2_cwin = NULL; + if (gid == _sgp2_cwin) + _sgp2_cwin = NULL; - delete gid; + delete gid; } void diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index d3cb203..ed388ad 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.6 2000/07/09 08:16:17 kevin Exp $ +** $Id: backprojectors.cpp,v 1.7 2000/07/11 10:32:44 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 @@ -26,18 +26,18 @@ #include "ct.h" -Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) +Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) { m_fail = false; m_pBackprojectImplem = NULL; - initBackprojector (proj, im, backprojName, interpName); + initBackprojector (proj, im, backprojName, interpName, interpFactor); } void Backprojector::BackprojectView (const double* const viewData, const double viewAngle) { - if (m_pBackprojectImplem) + if (m_pBackprojectImplem != NULL) m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } @@ -54,7 +54,7 @@ Backprojector::~Backprojector (void) // and initializes the backprojector bool -Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName) +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) { m_nameBackproject = backprojName; m_nameInterpolation = interpName; @@ -78,17 +78,17 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const } if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_TABLE) - m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_DIFF) - m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_DIFF2) - m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_IDIFF2) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); else if (m_idBackproject == BPROJ_IDIFF3) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation)); + m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -152,8 +152,8 @@ Backprojector::convertInterpolationNameToID (const char* const interpName) interpID = INTERP_NEAREST; else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0) interpID = INTERP_LINEAR; - else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATE_STR) == 0) - interpID = INTERP_FREQ_PREINTERPOLATE; + else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0) + interpID = INTERP_FREQ_PREINTERPOLATION; #if HAVE_BSPLINE_INTERP else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0) interpID = INTERP_BSPLINE; @@ -182,8 +182,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID) return (INTERP_NEAREST_STR); else if (interpID == INTERP_LINEAR) return (INTERP_LINEAR_STR); - else if (interpID == INTERP_FREQ_PREINTERPOLATE) - return (INTERP_FREQ_PREINTERPOLATE_STR); + else if (interpID == INTERP_FREQ_PREINTERPOLATION) + return (INTERP_FREQ_PREINTERPOLATION_STR); #if HAVE_BSPLINE_INTERP else if (interpID == INTERP_BSPLINE) return (INTERP_BSPLINE_STR); @@ -199,8 +199,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID) // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType) - : proj(proj), im(im), interpType(interpType) +Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor) + : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) { detInc = proj.detInc(); nDet = proj.nDet(); @@ -219,9 +219,6 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Backproj xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; - - if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR) - sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType); } Backproject::~Backproject (void) @@ -298,8 +295,8 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // PURPOSE // Precalculates trigometric function value for each point in image for backprojection. -BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : Backproject::Backproject (proj, im, interpType) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : Backproject::Backproject (proj, im, interpType, interpFactor) { arrayR.initSetSize (nx, ny); arrayPhi.initSetSize (nx, ny); @@ -360,8 +357,8 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Backprojects by precalculating the change in L position for each x & y step in the image. // Iterates in x & y direction by adding difference in L position -BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType) - : Backproject::Backproject (proj, im, interpType) +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor) + : Backproject::Backproject (proj, im, interpType, interpFactor) { // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; @@ -458,7 +455,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl if (iDetPos < 0 || iDetPos >= nDet - 1) errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); else - *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); + *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); } } // end for y } // end for x @@ -539,29 +536,35 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do // calculate L for first point in image (0, 0) kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + // precalculate scaled difference for linear interpolation + double deltaFilteredProj [nDet - 1]; + if (interpType == Backprojector::INTERP_LINEAR) { + for (int i = 0; i < nDet - 1; i++) + deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale; + } + for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { kint32 curDetPos = detPosColStart; ImageFileColumn pImCol = v[ix]; if (interpType == Backprojector::INTERP_NEAREST) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - int iDetPos = (curDetPos + halfScale) >> 16; + const int iDetPos = (curDetPos + halfScale) >> 16; + assert(iDetPos >= 0 && iDetPos < nDet); + *pImCol++ += filteredProj[iDetPos]; + } // end for iy + } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) { + for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { + const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor; assert(iDetPos >= 0 && iDetPos < nDet); *pImCol++ += filteredProj[iDetPos]; } // end for iy } else if (interpType == Backprojector::INTERP_LINEAR) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - kint32 iDetPos = curDetPos >> scaleShift; - kint32 detRemainder = curDetPos & scaleBitmask; -#if 0 - double frac = detRemainder * (1. / 65536.); + const kint32 iDetPos = curDetPos >> scaleShift; + const kint32 detRemainder = curDetPos & scaleBitmask; assert(iDetPos >= 0 && iDetPos < nDet - 1); - *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]); -#else - assert(iDetPos >= 0 && iDetPos < nDet - 1); - const double* const detPointer = &filteredProj[iDetPos]; - *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale; -#endif + *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy } //end linear } // end for ix diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index 65163f3..0f334b9 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.15 2000/07/07 15:30:59 kevin Exp $ +** $Id: filter.cpp,v 1.16 2000/07/11 10:32:44 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 @@ -28,6 +28,8 @@ #include "ct.h" +int SignalFilter::N_INTEGRAL=500; //static member + /* NAME * SignalFilter::SignalFilter Construct a signal * @@ -40,11 +42,9 @@ * 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 nSignalPoints, double param, const char* domainName, int zeropad = 0, int numIntegral = 0) +SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName, double bw, double signalIncrement, int nSignalPoints, double param, const char* domainName, int zeropad = 0, int preinterpolationFactor = 1) { m_vecFilter = NULL; m_vecFourierCosTable = NULL; @@ -70,15 +70,15 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName m_failMessage += domainName; return; } - init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, zeropad, numIntegral); + init (m_idFilter, m_idFilterMethod, bw, signalIncrement, nSignalPoints, param, m_idDomain, zeropad, preinterpolationFactor); } -SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int zeropad = 0, int numIntegral = 0) +SignalFilter::SignalFilter (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int zeropad = 0, int preinterpolationFactor = 1) { - init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, zeropad, numIntegral); + init (filterID, filterMethodID, bw, signalIncrement, nSignalPoints, param, domainID, zeropad, preinterpolationFactor); } -SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param, int numIntegral = 0) +SignalFilter::SignalFilter (const char* filterName, const char* domainName, double bw, double param) { m_bw = bw; m_nSignalPoints = 0; @@ -87,7 +87,6 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub m_vecFourierCosTable = NULL; m_vecFourierSinTable = NULL; m_filterParam = param; - m_numIntegral = numIntegral; m_idFilter = convertFilterNameToID (filterName); if (m_idFilter == FILTER_INVALID) { m_fail = true; @@ -105,7 +104,7 @@ SignalFilter::SignalFilter (const char* filterName, const char* domainName, doub } void -SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double param, const DomainID domainID, int zeropad, int numint) +SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID, double bw, double signalIncrement, int nSignalPoints, double filterParam, const DomainID domainID, int zeropad, int preinterpolationFactor) { m_bw = bw; m_idFilter = filterID; @@ -122,8 +121,9 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID m_fail = false; m_nSignalPoints = nSignalPoints; m_signalInc = signalIncrement; - m_filterParam = param; + m_filterParam = filterParam; m_zeropad = zeropad; + m_preinterpolationFactor = preinterpolationFactor; m_vecFourierCosTable = NULL; m_vecFourierSinTable = NULL; @@ -154,6 +154,7 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID m_nFilterPoints = 1 << nextPowerOf2; cout << "nFilterPoints = " << m_nFilterPoints << endl; } + m_nOutputPoints = m_nFilterPoints * m_preinterpolationFactor; m_filterMin = -1. / (2 * m_signalInc); m_filterMax = 1. / (2 * m_signalInc); m_filterInc = (m_filterMax - m_filterMin) / m_nFilterPoints; @@ -167,16 +168,16 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID // precalculate sin and cosine tables for fourier transform if (m_idFilterMethod == FILTER_METHOD_FOURIER_TABLE) { - int nFourier = m_nFilterPoints * m_nFilterPoints + 1; - double angleIncrement = (2. * PI) / m_nFilterPoints; - m_vecFourierCosTable = new double[ nFourier ]; - m_vecFourierSinTable = new double[ nFourier ]; - double angle = 0; - for (int i = 0; i < nFourier; i++) { - m_vecFourierCosTable[i] = cos (angle); - m_vecFourierSinTable[i] = sin (angle); - angle += angleIncrement; - } + int nFourier = max(m_nFilterPoints,m_nOutputPoints) * max(m_nFilterPoints, m_nOutputPoints) + 1; + double angleIncrement = (2. * PI) / m_nFilterPoints; + m_vecFourierCosTable = new double[ nFourier ]; + m_vecFourierSinTable = new double[ nFourier ]; + double angle = 0; + for (int i = 0; i < nFourier; i++) { + m_vecFourierCosTable[i] = cos (angle); + m_vecFourierSinTable[i] = sin (angle); + angle += angleIncrement; + } } #if HAVE_FFTW @@ -186,19 +187,21 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID } if (m_idFilterMethod == FILTER_METHOD_RFFTW) { - m_complexPlanForward = m_complexPlanBackward = NULL; m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE); - m_realPlanBackward = rfftw_create_plan (m_nFilterPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); + m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE); m_vecRealFftInput = new fftw_real [ m_nFilterPoints ]; + m_vecRealFftSignal = new fftw_real [ m_nOutputPoints ]; for (int i = 0; i < m_nFilterPoints; i++) m_vecRealFftInput[i] = 0; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { - m_realPlanForward = m_realPlanBackward = NULL; - m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); - m_complexPlanBackward = fftw_create_plan (m_nFilterPoints, FFTW_BACKWARD, FFTW_ESTIMATE); + m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE); + m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE); m_vecComplexFftInput = new fftw_complex [ m_nFilterPoints ]; + m_vecComplexFftSignal = new fftw_complex [ m_nOutputPoints ]; for (int i = 0; i < m_nFilterPoints; i++) m_vecComplexFftInput[i].re = m_vecComplexFftInput[i].im = 0; + for (int i = 0; i < m_nOutputPoints; i++) + m_vecComplexFftSignal[i].re = m_vecComplexFftSignal[i].im = 0; } #endif @@ -207,7 +210,6 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID 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) { @@ -223,15 +225,15 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID double x; int i; for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++) - m_vecFilter[i] = frequencyResponse (x, param); + m_vecFilter[i] = frequencyResponse (x, m_filterParam); } 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); + if (haveAnalyticSpatial(m_idFilter)) + m_vecFilter[i] = spatialResponseAnalytic (x, m_filterParam); else - m_vecFilter[i] = spatialResponseCalc (x, param, numint); + m_vecFilter[i] = spatialResponseCalc (x, m_filterParam); } else { m_failMessage = "Illegal domain name "; m_failMessage += m_idDomain; @@ -251,11 +253,13 @@ SignalFilter::~SignalFilter (void) fftw_destroy_plan(m_complexPlanForward); fftw_destroy_plan(m_complexPlanBackward); delete [] m_vecComplexFftInput; + delete [] m_vecComplexFftSignal; } if (m_idFilterMethod == FILTER_METHOD_RFFTW) { rfftw_destroy_plan(m_realPlanForward); rfftw_destroy_plan(m_realPlanBackward); delete [] m_vecRealFftInput; + delete [] m_vecRealFftSignal; } #endif } @@ -391,7 +395,6 @@ SignalFilter::convertDomainIDToName (const DomainID domain) return (name); } - void SignalFilter::filterSignal (const float input[], double output[]) const { @@ -432,29 +435,31 @@ SignalFilter::filterSignal (const float input[], double output[]) const for (int i = 0; i < m_nSignalPoints; i++) m_vecRealFftInput[i] = input[i]; - fftw_real out[m_nFilterPoints]; - rfftw_one (m_realPlanForward, m_vecRealFftInput, out); - for (int i = 0; i < m_nFilterPoints; i++) { - out[i] *= m_vecFilter[i]; - } - fftw_real outFiltered[m_nFilterPoints]; - rfftw_one(m_realPlanBackward, out, outFiltered); - for (int i = 0; i < m_nSignalPoints; i++) - output[i] = outFiltered[i]; + fftw_real fftOutput [ m_nFilterPoints ]; + rfftw_one (m_realPlanForward, m_vecRealFftInput, fftOutput); + for (int i = 0; i < m_nFilterPoints; i++) + m_vecRealFftSignal[i] = m_vecFilter[i] * fftOutput[i]; + for (int i = m_nFilterPoints; i < m_nOutputPoints; i++) + m_vecRealFftSignal[i] = 0; + + fftw_real ifftOutput [ m_nOutputPoints ]; + rfftw_one(m_realPlanBackward, m_vecRealFftSignal, ifftOutput); + for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++) + output[i] = ifftOutput[i]; } else if (m_idFilterMethod == FILTER_METHOD_FFTW) { for (int i = 0; i < m_nSignalPoints; i++) m_vecComplexFftInput[i].re = input[i]; - fftw_complex out[m_nFilterPoints]; - fftw_one(m_complexPlanForward, m_vecComplexFftInput, out); + fftw_complex fftOutput [ m_nFilterPoints ]; + fftw_one(m_complexPlanForward, m_vecComplexFftInput, fftOutput); for (int i = 0; i < m_nFilterPoints; i++) { - out[i].re *= m_vecFilter[i]; - out[i].im *= m_vecFilter[i]; + m_vecComplexFftSignal[i].re = m_vecFilter[i] * fftOutput[i].re; + m_vecComplexFftSignal[i].im = m_vecFilter[i] * fftOutput[i].im; } - fftw_complex outFiltered[m_nFilterPoints]; - fftw_one(m_complexPlanBackward, out, outFiltered); - for (int i = 0; i < m_nSignalPoints; i++) - output[i] = outFiltered[i].re; + fftw_complex ifftOutput [ m_nOutputPoints ]; + fftw_one(m_complexPlanBackward, m_vecComplexFftSignal, ifftOutput); + for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++) + output[i] = ifftOutput[i].re; } #endif } @@ -465,7 +470,7 @@ SignalFilter::response (double x) double response = 0; if (m_idDomain == DOMAIN_SPATIAL) - response = spatialResponse (m_idFilter, m_bw, x, m_filterParam, m_numIntegral); + response = spatialResponse (m_idFilter, m_bw, x, m_filterParam); else if (m_idDomain == DOMAIN_FREQUENCY) response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam); @@ -474,12 +479,12 @@ SignalFilter::response (double x) double -SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param, int nIntegral = 0) +SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param) { - if (nIntegral == 0) + if (haveAnalyticSpatial(filterID)) return spatialResponseAnalytic (filterID, bw, x, param); else - return spatialResponseCalc (filterID, bw, x, param, nIntegral); + return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL); } /* NAME @@ -498,9 +503,9 @@ SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double pa */ double -SignalFilter::spatialResponseCalc (double x, double param, int nIntegral) const +SignalFilter::spatialResponseCalc (double x, double param) const { - return (spatialResponseCalc (m_idFilter, m_bw, x, param, nIntegral)); + return (spatialResponseCalc (m_idFilter, m_bw, x, param, N_INTEGRAL)); } double @@ -632,6 +637,28 @@ SignalFilter::spatialResponseAnalytic (double x, double param) const return spatialResponseAnalytic (m_idFilter, m_bw, x, param); } +const bool +SignalFilter::haveAnalyticSpatial (FilterID filterID) +{ + bool haveAnalytic = false; + + switch (filterID) { + case FILTER_BANDLIMIT: + case FILTER_TRIANGLE: + case FILTER_COSINE: + case FILTER_G_HAMMING: + case FILTER_ABS_BANDLIMIT: + case FILTER_ABS_COSINE: + case FILTER_ABS_G_HAMMING: + case FILTER_SHEPP: + case FILTER_SINC: + haveAnalytic = true; + break; + } + + return (haveAnalytic); +} + double SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param) { diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index c3d2f67..9ae326f 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.6 2000/06/28 15:25:34 kevin Exp $ +** $Id: imagefile.cpp,v 1.7 2000/07/11 10:32:44 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 @@ -240,3 +240,194 @@ ImageFile::statistics (double& min, double& max, double& mean, double& mode, dou } +void +ImageFile::writeImagePGM (const char *outfile, int nxcell, int nycell, double densmin, double densmax) +{ + FILE *fp; + int nx = m_nx; + int ny = m_ny; + ImageFileArray v = getArray(); + + unsigned char rowp [nx * nxcell]; + + if ((fp = fopen (outfile, "wb")) == NULL) + return; + + fprintf(fp, "P5\n"); + fprintf(fp, "%d %d\n", nx, ny); + fprintf(fp, "255\n"); + + for (int irow = ny - 1; irow >= 0; irow--) { + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + for (int p = pos; p < pos + nxcell; p++) { + rowp[p] = static_cast (dens * 255.); + } + } + for (int ir = 0; ir < nycell; ir++) { + for (int ic = 0; ic < nx * nxcell; ic++) + fprintf(fp, "%c ", rowp[ic]); + } + } + + fclose(fp); +} + +void +ImageFile::writeImagePGMASCII (const char *outfile, int nxcell, int nycell, double densmin, double densmax) +{ + FILE *fp; + int nx = m_nx; + int ny = m_ny; + ImageFileArray v = getArray(); + + unsigned char rowp [nx * nxcell]; + if (rowp == NULL) + return; + + if ((fp = fopen (outfile, "wb")) == NULL) + return; + + fprintf(fp, "P2\n"); + fprintf(fp, "%d %d\n", nx, ny); + fprintf(fp, "255\n"); + + for (int irow = ny - 1; irow >= 0; irow--) { + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + for (int p = pos; p < pos + nxcell; p++) { + rowp[p] = static_cast (dens * 255.); + } + } + for (int ir = 0; ir < nycell; ir++) { + for (int ic = 0; ic < nx * nxcell; ic++) + fprintf(fp, "%d ", rowp[ic]); + fprintf(fp, "\n"); + } + } + + fclose(fp); +} + + +#ifdef HAVE_PNG +void +ImageFile::writeImagePNG (const char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax) +{ + FILE *fp; + png_structp png_ptr; + png_infop info_ptr; + double max_out_level = (1 << bitdepth) - 1; + int nx = m_nx; + int ny = m_ny; + ImageFileArray v = getArray(); + + unsigned char rowp [nx * nxcell * (bitdepth / 8)]; + + if ((fp = fopen (outfile, "wb")) == NULL) + return; + + png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); + if (! png_ptr) + return; + + info_ptr = png_create_info_struct(png_ptr); + if (! info_ptr) { + png_destroy_write_struct(&png_ptr, (png_infopp) NULL); + fclose(fp); + return; + } + + if (setjmp(png_ptr->jmpbuf)) { + png_destroy_write_struct(&png_ptr, &info_ptr); + fclose(fp); + return; + } + + png_init_io(png_ptr, fp); + + png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT); + + png_write_info(png_ptr, info_ptr); + for (int irow = ny - 1; irow >= 0; irow--) { + png_bytep row_pointer = rowp; + + for (int icol = 0; icol < nx; icol++) { + int pos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp (dens, 0., 1.); + unsigned int outval = static_cast (dens * max_out_level); + + for (int p = pos; p < pos + nxcell; p++) { + if (bitdepth == 8) + rowp[p] = outval; + else { + int rowpos = p * 2; + rowp[rowpos] = (outval >> 8) & 0xFF; + rowp[rowpos+1] = (outval & 0xFF); + } + } + } + for (int ir = 0; ir < nycell; ir++) + png_write_rows (png_ptr, &row_pointer, 1); + } + + png_write_end(png_ptr, info_ptr); + png_destroy_write_struct(&png_ptr, &info_ptr); + + fclose(fp); +} +#endif + +#ifdef HAVE_GD +#include "gd.h" +static const int N_GRAYSCALE=256; + +void +ImageFile::writeImageGIF (const char *outfile, int nxcell, int nycell, double densmin, double densmax) +{ + gdImagePtr gif; + FILE *out; + int gs_indices[N_GRAYSCALE]; + int nx = m_nx; + int ny = m_ny; + ImageFileArray v = getArray(); + + unsigned char rowp [nx * nxcell]; + if (rowp == NULL) + return; + + gif = gdImageCreate(nx * nxcell, ny * nycell); + for (int i = 0; i < N_GRAYSCALE; i++) + gs_indices[i] = gdImageColorAllocate(gif, i, i, i); + + int lastrow = ny * nycell - 1; + for (int irow = 0; irow < ny; irow++) { + int rpos = irow * nycell; + for (int ir = rpos; ir < rpos + nycell; ir++) { + for (int icol = 0; icol < nx; icol++) { + int cpos = icol * nxcell; + double dens = (v[icol][irow] - densmin) / (densmax - densmin); + dens = clamp(dens, 0., 1.); + for (int ic = cpos; ic < cpos + nxcell; ic++) { + rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1)); + gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]); + } + } + } + } + + if ((out = fopen(outfile,"w")) == NULL) { + sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile); + return (1); + } + gdImageGif(gif,out); + fclose(out); + gdImageDestroy(gif); +} +#endif + diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index db64314..8f1edf7 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.12 2000/07/06 08:30:30 kevin Exp $ +** $Id: projections.cpp,v 1.13 2000/07/11 10:32:44 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 @@ -487,17 +487,17 @@ Projections::printScanInfo (void) const */ bool -Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interp_param, const char* const backprojectName, const int trace) +Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace) { int nview = m_nView; double detInc = m_detInc; - int n_filteredProj = m_nDet; + int n_filteredProj = m_nDet * interpFactor; double filteredProj [n_filteredProj]; // filtered projections #ifdef HAVE_BSPLINE_INTERP int spline_order = 0, zoom_factor = 0; if (interp_type == I_BSPLINE) { - zoom_factor = interp_param; + zoom_factor = interpFactor; spline_order = 3; zoom_factor = 3; n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1; @@ -505,7 +505,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi #endif double filterBW = 1. / detInc; - SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", zeropad); + SignalFilter filter (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", zeropad, interpFactor); filter.setTraceLevel(trace); if (filter.fail()) { @@ -540,7 +540,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi } #endif //HAVE_SGP - Backprojector bj (*this, im, backprojectName, interpName); + Backprojector bj (*this, im, backprojectName, interpName, interpFactor); if (bj.fail()) { sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", bj.failMessage().c_str()); return false; @@ -574,8 +574,9 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi ezset ("xporigin .5."); ezset ("xlength .5."); ezset ("box"); + ezset ("grid"); - gid = ezplot (filteredProj, plot_xaxis, m_nDet); + gid = ezplot (filteredProj, plot_xaxis, n_filteredProj); } #endif //HAVE_SGP diff --git a/src/Makefile.am b/src/Makefile.am index 9b48e9c..203ab18 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,13 +1,11 @@ -bin_PROGRAMS = ctsim pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img +SUBDIRS = ctsim +bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img bin_SCRIPTS = sample-ctsim.sh EXTRA_PROGRAMS = pjrec-lam phm2pj-lam phm2if-lam INCLUDES=@my_includes@ EXTRA_DIST=Makefile.nt mpiworld.cpp SOURCE_DEPEND=../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a -ctsim_SOURCES = ctsim.cpp $ -ctsim_LDADD = @ctlibs@ -ctsim_DEPENDENCIES=$(SOURCE_DEPEND) pjrec_SOURCES = pjrec.cpp pjrec_LDADD=@ctlibs@ pjrec_DEPENDENCIES=$(SOURCE_DEPEND) diff --git a/src/ctsim.cpp b/src/ctsim.cpp deleted file mode 100644 index 57b5549..0000000 --- a/src/ctsim.cpp +++ /dev/null @@ -1,33 +0,0 @@ -/***************************************************************************** -** FILE IDENTIFICATION -** -** Name: ctsim.cpp -** Purpose: GUI shell for CT Simulation -** Programmer: Kevin Rosenberg -** Date Started: June 2000 -** -** -** 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 -** published by the Free Software Foundation. -** -** This program is distributed in the hope that it will be useful, -** but WITHOUT ANY WARRANTY; without even the implied warranty of -** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -** GNU General Public License for more details. -** -** You should have received a copy of the GNU General Public License -** along with this program; if not, write to the Free Software -** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA -******************************************************************************/ - -#include -#include - -int -main (int argc, char **argv) -{ - cout << "Hello, CTSim!" << endl; - - return (0); -} diff --git a/src/if2img.cpp b/src/if2img.cpp index 8211d6d..e6ec073 100644 --- a/src/if2img.cpp +++ b/src/if2img.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if2img.cpp,v 1.12 2000/07/09 08:47:41 kevin Exp $ +** $Id: if2img.cpp,v 1.13 2000/07/11 10:32:44 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 @@ -28,15 +28,6 @@ #include "ct.h" -#if HAVE_PNG -void sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax); -#endif -#if HAVE_GIF -void sdf2d_to_gif (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax); -#endif -void sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax); -void sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax); - enum { O_SCALE, O_MIN, O_MAX, O_AUTO, O_CENTER, O_STATS, O_FORMAT, O_LABELS, O_HELP, O_VERBOSE, O_VERSION, O_DEBUG }; @@ -83,7 +74,7 @@ if2img_usage (const char *program) cout << "usage: " << fileBasename(program) << " ifname outfile [OPTIONS]" << endl; cout << "Convert IF file to an image file" << endl; cout << endl; - cout << " sdfname Name of input SDF file" << endl; + cout << " ifname Name of input file" << endl; cout << " outfile Name of output file" << endl; cout << " --format Output format" << endl; cout << " pgm PGM (portable graymap) format (default)" << endl; @@ -92,7 +83,7 @@ if2img_usage (const char *program) cout << " png PNG (8-bit) format" << endl; cout << " png16 PNG (16-bit) format" << endl; #endif -#ifdef HAVE_GIF +#if HAVE_G2 cout << " gif GIF format" << endl; #endif cout << " disp Display on screen" << endl; @@ -122,7 +113,6 @@ int if2img_main (int argc, char *const argv[]) { ImageFile* pim = NULL; - double densmin = HUGE_VAL, densmax = -HUGE_VAL; char *in_file, *out_file; int opt_verbose = 0; @@ -287,106 +277,59 @@ if2img_main (int argc, char *const argv[]) im.printLabels(cout); if (opt_stats || (! (opt_set_max && opt_set_min))) { - double minfound = HUGE_VAL, maxfound = -HUGE_VAL; - double mode = 0, mean = 0, stddev = 0, window = 0; - double spread; - int hist[256]; - int ibin, nbin = 256; - int max_bin, max_binindex; - double maxbin; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); - - maxbin = nbin - 1; - for (int ix = 0; ix < nx; ix++) { - for (int 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 (int ix = 0; ix < nx; ix++) { - for (int 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; - } + double min, max, mean, mode, median, stddev; + double window; + im.statistics(min, max, mean, mode, median, stddev); - if (opt_auto == O_AUTO_FULL) { + if (opt_auto == O_AUTO_FULL) { if (! opt_set_max) - densmax = maxfound; + densmax = max; if (! opt_set_min) - densmin = minfound; - } - if (opt_stats || opt_auto != O_AUTO_FULL) { - mean /= nx * ny; - for (int ix = 0; ix < nx; ix++) { - for (int iy = 0; iy < ny; iy++) { - double diff = (v[ix][iy] - mean); - stddev += diff * diff; - } + densmin = min; } - stddev = sqrt(stddev / (nx * ny)); - if (opt_auto == O_AUTO_FULL) - ; - else if (opt_auto == O_AUTO_STD1) - window = stddev; - else if (opt_auto == O_AUTO_STD0_1) - window = stddev * 0.1; - else if (opt_auto == O_AUTO_STD0_5) - window = stddev * 0.5; + if (opt_stats || opt_auto != O_AUTO_FULL) { + + if (opt_auto == O_AUTO_FULL) + ; + else if (opt_auto == O_AUTO_STD1) + window = stddev; + else if (opt_auto == O_AUTO_STD0_1) + window = stddev * 0.1; + else if (opt_auto == O_AUTO_STD0_5) + window = stddev * 0.5; else if (opt_auto == O_AUTO_STD2) - window = stddev * 2; - else if (opt_auto == O_AUTO_STD3) - window = stddev * 3; - else { - sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto); - return (1); + window = stddev * 2; + else if (opt_auto == O_AUTO_STD3) + window = stddev * 3; + else { + sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto); + return (1); + } } - } - if (opt_stats) { - cout <<"nx: " << nx << endl; - cout <<"ny: " << ny << endl; - cout <<"min: " << minfound << endl; - cout <<"max: " << maxfound << endl; - cout <<"mean: " << mean << endl; - cout <<"mode: " << mode << endl; - cout <<"stddev: " << stddev << endl; - } - if (opt_auto != O_AUTO_FULL) { - double center; - - if (opt_center == O_CENTER_MODE) - center = mode; - else if (opt_center == O_CENTER_MEAN) - center = mean; - else { - sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center); - return (1); + if (opt_stats) { + cout <<"nx: " << im.nx() << endl; + cout <<"ny: " << im.ny() << endl; + cout <<"min: " << min << endl; + cout <<"max: " << max << endl; + cout <<"mean: " << mean << endl; + cout <<"mode: " << mode << endl; + cout <<"stddev: " << stddev << endl; } - if (! opt_set_max) - densmax = center + window; - if (! opt_set_min) - densmin = center - window; + if (opt_auto != O_AUTO_FULL) { + double center; + + if (opt_center == O_CENTER_MODE) + center = mode; + else if (opt_center == O_CENTER_MEAN) + center = mean; + else { + sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center); + return (1); + } + if (! opt_set_max) + densmax = center + window; + if (! opt_set_min) + densmin = center - window; } } @@ -396,24 +339,25 @@ if2img_main (int argc, char *const argv[]) } if (opt_format == O_FORMAT_PGM) - sdf2d_to_pgm (im, out_file, opt_scale, opt_scale, densmin, densmax); + im.writeImagePGM (out_file, opt_scale, opt_scale, densmin, densmax); else if (opt_format == O_FORMAT_PGMASC) - sdf2d_to_pgmasc (im, out_file, opt_scale, opt_scale, densmin, densmax); + im.writeImagePGMASCII (out_file, opt_scale, opt_scale, densmin, densmax); #if HAVE_PNG else if (opt_format == O_FORMAT_PNG) - sdf2d_to_png (im, out_file, 8, opt_scale, opt_scale, densmin, densmax); + im.writeImagePNG (out_file, 8, opt_scale, opt_scale, densmin, densmax); else if (opt_format == O_FORMAT_PNG16) - sdf2d_to_png (im, out_file, 16, opt_scale, opt_scale, densmin, densmax); + im.writeImagePNG (out_file, 16, opt_scale, opt_scale, densmin, densmax); #endif #if HAVE_GIF else if (opt_format == O_FORMAT_GIF) - sdf2d_to_gif (im, out_file, opt_scale, opt_scale, densmin, densmax); + im.writeImageGIF (out_file, opt_scale, opt_scale, densmin, densmax); #endif else if (opt_format == O_FORMAT_DISP) { #if HAVE_SGP im.displayScaling (opt_scale, densmin, densmax); + cout << "Press enter to continue\n"; cio_kb_getc(); - // sgp2_close(sgp2_get_active_win()); + sgp2_close(sgp2_get_active_win()); #endif } else @@ -425,201 +369,6 @@ if2img_main (int argc, char *const argv[]) } -void -sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax) -{ - FILE *fp; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); - - unsigned char rowp [nx * nxcell]; - - if ((fp = fopen (outfile, "wb")) == NULL) - return; - - fprintf(fp, "P5\n"); - fprintf(fp, "%d %d\n", nx, ny); - fprintf(fp, "255\n"); - - for (int irow = ny - 1; irow >= 0; irow--) { - for (int icol = 0; icol < nx; icol++) { - int pos = icol * nxcell; - double dens = (v[icol][irow] - densmin) / (densmax - densmin); - dens = clamp (dens, 0., 1.); - for (int p = pos; p < pos + nxcell; p++) { - rowp[p] = static_cast (dens * 255.); - } - } - for (int ir = 0; ir < nycell; ir++) { - for (int ic = 0; ic < nx * nxcell; ic++) - fprintf(fp, "%c ", rowp[ic]); - } - } - - fclose(fp); -} - -void -sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax) -{ - FILE *fp; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); - - unsigned char rowp [nx * nxcell]; - if (rowp == NULL) - return; - - if ((fp = fopen (outfile, "wb")) == NULL) - return; - - fprintf(fp, "P2\n"); - fprintf(fp, "%d %d\n", nx, ny); - fprintf(fp, "255\n"); - - for (int irow = ny - 1; irow >= 0; irow--) { - for (int icol = 0; icol < nx; icol++) { - int pos = icol * nxcell; - double dens = (v[icol][irow] - densmin) / (densmax - densmin); - dens = clamp (dens, 0., 1.); - for (int p = pos; p < pos + nxcell; p++) { - rowp[p] = static_cast (dens * 255.); - } - } - for (int ir = 0; ir < nycell; ir++) { - for (int ic = 0; ic < nx * nxcell; ic++) - fprintf(fp, "%d ", rowp[ic]); - fprintf(fp, "\n"); - } - } - - fclose(fp); -} - - -#ifdef HAVE_PNG -void -sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell, double densmin, double densmax) -{ - FILE *fp; - png_structp png_ptr; - png_infop info_ptr; - double max_out_level = (1 << bitdepth) - 1; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); - - unsigned char rowp [nx * nxcell * (bitdepth / 8)]; - - if ((fp = fopen (outfile, "wb")) == NULL) - return; - - png_ptr = png_create_write_struct (PNG_LIBPNG_VER_STRING, NULL, NULL, NULL); - if (! png_ptr) - return; - - info_ptr = png_create_info_struct(png_ptr); - if (! info_ptr) { - png_destroy_write_struct(&png_ptr, (png_infopp) NULL); - fclose(fp); - return; - } - - if (setjmp(png_ptr->jmpbuf)) { - png_destroy_write_struct(&png_ptr, &info_ptr); - fclose(fp); - return; - } - - png_init_io(png_ptr, fp); - - png_set_IHDR(png_ptr, info_ptr, nx * nxcell, ny * nycell, bitdepth, PNG_COLOR_TYPE_GRAY, PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_BASE, PNG_FILTER_TYPE_DEFAULT); - - png_write_info(png_ptr, info_ptr); - for (int irow = ny - 1; irow >= 0; irow--) { - png_bytep row_pointer = rowp; - - for (int icol = 0; icol < nx; icol++) { - int pos = icol * nxcell; - double dens = (v[icol][irow] - densmin) / (densmax - densmin); - dens = clamp (dens, 0., 1.); - unsigned int outval = static_cast (dens * max_out_level); - - for (int p = pos; p < pos + nxcell; p++) { - if (bitdepth == 8) - rowp[p] = outval; - else { - int rowpos = p * 2; - rowp[rowpos] = (outval >> 8) & 0xFF; - rowp[rowpos+1] = (outval & 0xFF); - } - } - } - for (int ir = 0; ir < nycell; ir++) - png_write_rows (png_ptr, &row_pointer, 1); - } - - png_write_end(png_ptr, info_ptr); - png_destroy_write_struct(&png_ptr, &info_ptr); - - fclose(fp); -} -#endif - -#ifdef HAVE_GD -#include "gd.h" -static const int N_GRAYSCALE=256; -#endif - -void -sdf2d_to_gif (ImageFile& im, char *outfile, int nxcell, int nycell, double densmin, double densmax) -{ -#ifdef HAVE_GD - gdImagePtr gif; - FILE *out; - int gs_indices[N_GRAYSCALE]; - int nx = im.nx(); - int ny = im.ny(); - ImageFileArray v = im.getArray(); - - unsigned char rowp [nx * nxcell]; - if (rowp == NULL) - return; - - gif = gdImageCreate(nx * nxcell, ny * nycell); - for (int i = 0; i < N_GRAYSCALE; i++) - gs_indices[i] = gdImageColorAllocate(gif, i, i, i); - - int lastrow = ny * nycell - 1; - for (int irow = 0; irow < ny; irow++) { - int rpos = irow * nycell; - for (int ir = rpos; ir < rpos + nycell; ir++) { - for (int icol = 0; icol < nx; icol++) { - int cpos = icol * nxcell; - double dens = (v[icol][irow] - densmin) / (densmax - densmin); - dens = clamp(dens, 0., 1.); - for (int ic = cpos; ic < cpos + nxcell; ic++) { - rowp[ic] = (unsigned int) (dens * (double) (N_GRAYSCALE - 1)); - gdImageSetPixel(gif, ic, lastrow - ir, gs_indices[rowp[ic]]); - } - } - } - } - - if ((out = fopen(outfile,"w")) == NULL) { - sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile); - return (1); - } - gdImageGif(gif,out); - fclose(out); - gdImageDestroy(gif); -#else - cout << "This version does not support GIF" << endl; -#endif -} - #ifndef NO_MAIN int main (int argc, char *const argv[]) diff --git a/src/pjrec.cpp b/src/pjrec.cpp index 40a8cba..d6a229c 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.9 2000/07/09 08:16:18 kevin Exp $ +** $Id: pjrec.cpp,v 1.10 2000/07/11 10:32:44 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,11 +29,12 @@ #include "timer.h" -enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION}; +enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION}; static struct option my_options[] = { {"interp", 1, 0, O_INTERP}, + {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR}, {"filter", 1, 0, O_FILTER}, {"filter-method", 1, 0, O_FILTER_METHOD}, {"zeropad", 1, 0, O_ZEROPAD}, @@ -64,6 +65,8 @@ pjrec_usage (const char *program) #if HAVE_BSPLINE_INTERP cout << " bspline B-spline interpolation" << endl; #endif + cout << " --preinterpolate Preinterpolation factor (default = 1)\n"; + cout << " Used only with frequency-based filtering\n"; cout << " --filter Filter name" << endl; cout << " abs_bandlimit Abs * Bandlimiting (default)" << endl; cout << " abs_sinc Abs * Sinc" << endl; @@ -133,7 +136,7 @@ pjrec_main (int argc, char * argv[]) string optInterpName = Backprojector::INTERP_LINEAR_STR; string optBackprojName = Backprojector::BPROJ_IDIFF2_STR; // string optFilterMethodName = SignalFilter::FILTER_METHOD_CONVOLUTION_STR; - int optInterpParam = 1; + int optPreinterpolationFactor = 1; int nx, ny; #ifdef HAVE_MPI ImageFile* imLocal; @@ -159,6 +162,13 @@ pjrec_main (int argc, char * argv[]) case O_INTERP: optInterpName = optarg; break; + case O_PREINTERPOLATION_FACTOR: + optPreinterpolationFactor = strtol(optarg, &endptr, 10); + if (endptr != optarg + strlen(optarg)) { + pjrec_usage(argv[0]); + return(1); + } + break; case O_FILTER: optFilterName = optarg; break; @@ -172,12 +182,14 @@ pjrec_main (int argc, char * argv[]) optFilterParam = strtod(optarg, &endptr); if (endptr != optarg + strlen(optarg)) { pjrec_usage(argv[0]); + return(1); } break; case O_ZEROPAD: optZeroPad = strtol(optarg, &endptr, 10); if (endptr != optarg + strlen(optarg)) { pjrec_usage(argv[0]); + return(1); } break; case O_VERBOSE: @@ -228,7 +240,7 @@ pjrec_main (int argc, char * argv[]) filterDesc << optFilterName; ostringstream label; - label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName; + label << "pjrec: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", preinterpolation=" << optPreinterpolationFactor << ", " << optBackprojName; remark = label.str(); if (optVerbose) @@ -260,7 +272,7 @@ pjrec_main (int argc, char * argv[]) mpiWorld.getComm().Bcast (&optTrace, 1, MPI::INT, 0); mpiWorld.getComm().Bcast (&optFilterParam, 1, MPI::DOUBLE, 0); mpiWorld.getComm().Bcast (&optZeroPad, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&optInterpParam, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&optPreinterpolationFactor, 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); @@ -298,7 +310,7 @@ pjrec_main (int argc, char * argv[]) #ifdef HAVE_MPI TimerCollectiveMPI timerReconstruct (mpiWorld.getComm()); - projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace); + projLocal.reconstruct (*imLocal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace); if (optVerbose) timerReconstruct.timerEndAndReport ("Time to reconstruct"); @@ -307,7 +319,7 @@ pjrec_main (int argc, char * argv[]) if (optVerbose) timerReduce.timerEndAndReport ("Time to reduce image"); #else - projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optInterpParam, optBackprojName.c_str(), optTrace); + projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), optFilterParam, optFilterMethodName.c_str(), optZeroPad, optInterpName.c_str(), optPreinterpolationFactor, optBackprojName.c_str(), optTrace); #endif #ifdef HAVE_MPI -- 2.34.1