r142: *** empty log message ***
authorKevin M. Rosenberg <kevin@rosenberg.net>
Tue, 11 Jul 2000 10:32:44 +0000 (10:32 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Tue, 11 Jul 2000 10:32:44 +0000 (10:32 +0000)
15 files changed:
ChangeLog
configure
configure.in
include/backprojectors.h
include/filter.h
include/imagefile.h
libctgraphics/sgp.cpp
libctsim/backprojectors.cpp
libctsim/filter.cpp
libctsim/imagefile.cpp
libctsim/projections.cpp
src/Makefile.am
src/ctsim.cpp [deleted file]
src/if2img.cpp
src/pjrec.cpp

index 055fcff298ab21a0cdb08754e7ef657b15d7c12c..66a5fff5cac6de41a8d11bcba29ab153f2cc5c07 100644 (file)
--- 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
    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
    Cleaned up SignalFilter class
    Added zeropad option to pjrec
    Added zeropad options to html and cgi files
index 2b9a3c6ccead6cb165f7845786e03bd6cf53d62d..145a261d1ed925034053beba333c77915a07950a 100755 (executable)
--- a/configure
+++ b/configure
@@ -2803,13 +2803,12 @@ fi
 fi
 
 
 fi
 
 
-
  if test -n "$GCC"; then
     CFLAGS="$CFLAGS -Wall"
  fi
  
 echo $ac_n "checking whether to enable verbose warnings""... $ac_c" 1>&6
  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"
 # 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 $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"
 # 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
 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
   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
   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
   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"
   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
 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"
 # 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
 # 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
 
 # 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
 
   # First, try using that file with no special directory specified.
 cat > conftest.$ac_ext <<EOF
-#line 3091 "configure"
+#line 3098 "configure"
 #include "confdefs.h"
 #include <$x_direct_test_include>
 EOF
 ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out"
 #include "confdefs.h"
 #include <$x_direct_test_include>
 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*
 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 <<EOF
   ac_save_LIBS="$LIBS"
   LIBS="-l$x_direct_test_library $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3165 "configure"
+#line 3172 "configure"
 #include "confdefs.h"
 
 int main() {
 ${x_direct_test_function}()
 ; return 0; }
 EOF
 #include "confdefs.h"
 
 int main() {
 ${x_direct_test_function}()
 ; return 0; }
 EOF
-if { (eval echo configure:3172: \"$ac_link\") 1>&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.
   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
     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 <<EOF
       ac_xsave_LIBS="$LIBS"; LIBS="$LIBS -R$x_libraries"
       cat > conftest.$ac_ext <<EOF
-#line 3281 "configure"
+#line 3288 "configure"
 #include "confdefs.h"
 
 int main() {
 
 ; return 0; }
 EOF
 #include "confdefs.h"
 
 int main() {
 
 ; return 0; }
 EOF
-if { (eval echo configure:3288: \"$ac_link\") 1>&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
   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 <<EOF
       else
        LIBS="$ac_xsave_LIBS -R $x_libraries"
        cat > conftest.$ac_ext <<EOF
-#line 3304 "configure"
+#line 3311 "configure"
 #include "confdefs.h"
 
 int main() {
 
 ; return 0; }
 EOF
 #include "confdefs.h"
 
 int main() {
 
 ; return 0; }
 EOF
-if { (eval echo configure:3311: \"$ac_link\") 1>&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
   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
     # 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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-ldnet  $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3351 "configure"
+#line 3358 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3358,7 +3365,7 @@ int main() {
 dnet_ntoa()
 ; return 0; }
 EOF
 dnet_ntoa()
 ; return 0; }
 EOF
-if { (eval echo configure:3362: \"$ac_link\") 1>&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
   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
 
     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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-ldnet_stub  $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3392 "configure"
+#line 3399 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3399,7 +3406,7 @@ int main() {
 dnet_ntoa()
 ; return 0; }
 EOF
 dnet_ntoa()
 ; return 0; }
 EOF
-if { (eval echo configure:3403: \"$ac_link\") 1>&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
   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
     # 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 <<EOF
 if eval "test \"`echo '$''{'ac_cv_func_gethostbyname'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 3437 "configure"
+#line 3444 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char gethostbyname(); below.  */
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char gethostbyname(); below.  */
@@ -3456,7 +3463,7 @@ gethostbyname();
 
 ; return 0; }
 EOF
 
 ; return 0; }
 EOF
-if { (eval echo configure:3460: \"$ac_link\") 1>&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
   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
 
     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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lnsl  $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3489 "configure"
+#line 3496 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3496,7 +3503,7 @@ int main() {
 gethostbyname()
 ; return 0; }
 EOF
 gethostbyname()
 ; return 0; }
 EOF
-if { (eval echo configure:3500: \"$ac_link\") 1>&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
   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
     # -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 <<EOF
 if eval "test \"`echo '$''{'ac_cv_func_connect'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 3535 "configure"
+#line 3542 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char connect(); below.  */
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char connect(); below.  */
@@ -3554,7 +3561,7 @@ connect();
 
 ; return 0; }
 EOF
 
 ; return 0; }
 EOF
-if { (eval echo configure:3558: \"$ac_link\") 1>&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
   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
 
     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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lsocket $X_EXTRA_LIBS $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3587 "configure"
+#line 3594 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3594,7 +3601,7 @@ int main() {
 connect()
 ; return 0; }
 EOF
 connect()
 ; return 0; }
 EOF
-if { (eval echo configure:3598: \"$ac_link\") 1>&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
   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
 
     # 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 <<EOF
 if eval "test \"`echo '$''{'ac_cv_func_remove'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 3627 "configure"
+#line 3634 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char remove(); below.  */
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char remove(); below.  */
@@ -3646,7 +3653,7 @@ remove();
 
 ; return 0; }
 EOF
 
 ; return 0; }
 EOF
-if { (eval echo configure:3650: \"$ac_link\") 1>&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
   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
 
     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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lposix  $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3679 "configure"
+#line 3686 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3686,7 +3693,7 @@ int main() {
 remove()
 ; return 0; }
 EOF
 remove()
 ; return 0; }
 EOF
-if { (eval echo configure:3690: \"$ac_link\") 1>&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
   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
 
     # 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 <<EOF
 if eval "test \"`echo '$''{'ac_cv_func_shmat'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   cat > conftest.$ac_ext <<EOF
-#line 3719 "configure"
+#line 3726 "configure"
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char shmat(); below.  */
 #include "confdefs.h"
 /* System header to define __stub macros and hopefully few prototypes,
     which can conflict with char shmat(); below.  */
@@ -3738,7 +3745,7 @@ shmat();
 
 ; return 0; }
 EOF
 
 ; return 0; }
 EOF
-if { (eval echo configure:3742: \"$ac_link\") 1>&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
   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
 
     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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lipc  $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3771 "configure"
+#line 3778 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3778,7 +3785,7 @@ int main() {
 shmat()
 ; return 0; }
 EOF
 shmat()
 ; return 0; }
 EOF
-if { (eval echo configure:3782: \"$ac_link\") 1>&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
   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
   # 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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lICE $X_EXTRA_LIBS $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3823 "configure"
+#line 3830 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3830,7 +3837,7 @@ int main() {
 IceConnectionNumber()
 ; return 0; }
 EOF
 IceConnectionNumber()
 ; return 0; }
 EOF
-if { (eval echo configure:3834: \"$ac_link\") 1>&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
   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
 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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lXt $X_PRE_LIBS $X_BASIC_LIBS $X_EXTRA_LIBS $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3879 "configure"
+#line 3886 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3886,7 +3893,7 @@ int main() {
 XtToolkitThreadInitialize()
 ; return 0; }
 EOF
 XtToolkitThreadInitialize()
 ; return 0; }
 EOF
-if { (eval echo configure:3890: \"$ac_link\") 1>&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
   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
 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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lXext -lX11 "$X_EXTRA_LIBS" $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3926 "configure"
+#line 3933 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3933,7 +3940,7 @@ int main() {
 XdbeQueryExtension()
 ; return 0; }
 EOF
 XdbeQueryExtension()
 ; return 0; }
 EOF
-if { (eval echo configure:3937: \"$ac_link\") 1>&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
   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
 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
 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 <<EOF
   ac_save_LIBS="$LIBS"
 LIBS="-lXext -lX11 "$X_EXTRA_LIBS" $LIBS"
 cat > conftest.$ac_ext <<EOF
-#line 3969 "configure"
+#line 3976 "configure"
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
 #include "confdefs.h"
 /* Override any gcc2 internal prototype to avoid an error.  */
 /* We use char because int might match the return type of a gcc2
@@ -3976,7 +3983,7 @@ int main() {
 XmbufQueryExtension()
 ; return 0; }
 EOF
 XmbufQueryExtension()
 ; return 0; }
 EOF
-if { (eval echo configure:3980: \"$ac_link\") 1>&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
   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
 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
 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
-#line 4015 "configure"
+#line 4022 "configure"
 #include "confdefs.h"
 #include <$ac_hdr>
 EOF
 ac_try="$ac_cpp conftest.$ac_ext >/dev/null 2>conftest.out"
 #include "confdefs.h"
 #include <$ac_hdr>
 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*
 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 $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
 if test "$no_x" != "yes" ; then
   cat >> confdefs.h <<\EOF
 #define HAVE_X11 1
index 20cc06498eb751008de0e1040a56d0a579d5f1fb..e3ea013d03db779b9419bfdfa1b05f46779242a3 100644 (file)
@@ -12,7 +12,7 @@ AC_PROG_AWK
 AC_PROG_INSTALL
 AC_PROG_MAKE_SET
 AC_PROG_RANLIB
 AC_PROG_INSTALL
 AC_PROG_MAKE_SET
 AC_PROG_RANLIB
-AC_PROG_CC([-g])
+AC_PROG_CC
 AC_PROG_CXX
 
 
 AC_PROG_CXX
 
 
@@ -114,7 +114,6 @@ if test "${getopt_long}" = "false" ; then
   AM_CONDITIONAL(INCLUDED_GETOPT_LONG, test 1==1)
 fi
 
   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,
 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
 AM_CONDITIONAL(DEBUG, test "$debug" = "true")
 
 if test "$debug" = "true" ; then
-  CFLAGS="-g -DDEBUG"
+  AC_ADD_GCC_CFLAGS([-g -DDEBUG])
   AC_DEFINE(DEBUG)
 else
   AC_DEFINE(DEBUG)
 else
-  CFLAGS="-g -O3 -DNDEBUG"
+  AC_ADD_GCC_CFLAGS([-g -O3 -DNDEBUG --fast-math])
   AC_DEFINE(NDEBUG)
 fi
 
   AC_DEFINE(NDEBUG)
 fi
 
index f0cc19a1d0147ce6da2dab7169f44735418cc2df..ab15930d600eb6ca695f3adb2e407b0cf736977a 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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 
     I_3BSPLINE,
 #endif
     INTERP_LINEAR,        // Linear interpolation 
-    INTERP_FREQ_PREINTERPOLATE,
+    INTERP_FREQ_PREINTERPOLATION,
   } InterpolationID;
 
   static const char BPROJ_TRIG_STR[]=     "trig";
   } 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_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);
                 
 
   ~Backprojector (void);
                 
@@ -91,14 +91,14 @@ class Backprojector
   static const BackprojectID convertBackprojectNameToID (const char* const bprojName);
   static const char* convertBackprojectIDToName (const BackprojectID bprojID);
 
   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:
 };
 
 
 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 ();
 
 
     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 nDet;
     double xMin, xMax, yMin, yMax;     // Retangular coords of phantom
     double xInc, yInc; // size of cells
+    int m_interpFactor;
 
  private:
     Backproject (const Backproject& rhs);
 
  private:
     Backproject (const Backproject& rhs);
@@ -131,8 +132,8 @@ class Backproject
 class BackprojectTrig : public Backproject
 {
  public:
 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);
       {}
 
   void BackprojectView (const double* const t, double view_angle);
@@ -142,7 +143,7 @@ class BackprojectTrig : public Backproject
 class BackprojectTable : public Backproject
 {
  public:
 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);
   ~BackprojectTable ();
 
   void BackprojectView (const double* const t, double view_angle);
@@ -158,7 +159,7 @@ class BackprojectTable : public Backproject
 class BackprojectDiff : public Backproject
 {
  public:
 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);
   ~BackprojectDiff ();
 
   void BackprojectView (const double* const t, double view_angle);
@@ -172,8 +173,8 @@ class BackprojectDiff : public Backproject
 class BackprojectDiff2 : public BackprojectDiff
 {
  public:
 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);
     {}
 
   void BackprojectView (const double* const t, double view_angle);
@@ -182,8 +183,8 @@ class BackprojectDiff2 : public BackprojectDiff
 class BackprojectIntDiff2 : public BackprojectDiff
 {
  public:
 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);
     {}
   
   void BackprojectView (const double* const t, double view_angle);
@@ -193,8 +194,8 @@ class BackprojectIntDiff2 : public BackprojectDiff
 class BackprojectIntDiff3 : public BackprojectDiff
 {
  public:
 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);
     {}
   
   void BackprojectView (const double* const t, double view_angle);
index 9ee7390a4acd83c96bd9a72b7f9990fbfe7bede9..4a25b5bcd6eee6a91eca50a4407dcf33af30cb1b 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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";
 
 
     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);
 
 
     ~SignalFilter (void);
 
@@ -132,14 +132,16 @@ class SignalFilter {
 
     double response (double x);
 
 
     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 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 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;
  private:
     double m_bw;
     int m_nFilterPoints;
@@ -153,9 +155,9 @@ class SignalFilter {
     double* m_vecFourierSinTable;
     complex<double>* m_complexVecFilter;
 #ifdef HAVE_FFTW
     double* m_vecFourierSinTable;
     complex<double>* m_complexVecFilter;
 #ifdef HAVE_FFTW
-    fftw_real* m_vecRealFftInput;
+    fftw_real* m_vecRealFftInput, *m_vecRealFftSignal;
     rfftw_plan m_realPlanForward, m_realPlanBackward;
     rfftw_plan m_realPlanForward, m_realPlanBackward;
-    fftw_complex* m_vecComplexFftInput;
+    fftw_complex* m_vecComplexFftInput, *m_vecComplexFftSignal;
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #endif
 
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #endif
 
@@ -168,10 +170,14 @@ class SignalFilter {
     FilterMethodID m_idFilterMethod;
     DomainID m_idDomain;
     double m_filterParam;
     FilterMethodID m_idFilterMethod;
     DomainID m_idDomain;
     double m_filterParam;
-    int m_numIntegral;
     int m_traceLevel;
     int m_zeropad;
     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);
     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);
 
     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;
 
 
     double spatialResponseAnalytic (double x, double param) const;
 
@@ -194,5 +200,4 @@ double spatialResponseCalc (double x, double param, int n) const;
 
 };
 
 
 };
 
-
 #endif
 #endif
index ef60d3ddebf39f26a92ae2352c6ae9a6828b4009..5f8c3248e2f7e7709adce82672aad04a85115524 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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);
 
 
   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);
 };
 
 
 };
 
 
index cc87671d0c45388327cba801bc1c69b2a0cce518..8a5400fc00e620ffbbc802a558a1ed2dd47afe98 100644 (file)
@@ -7,7 +7,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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    
 sgp2_close (SGP_ID gid)
 {
 #if HAVE_G2_H    
+  if (gid)
     g2_close (gid->g2_id);
 #endif
     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
 }
 
 void
index d3cb2035dc125af3415006a2d204cc61b5cc7907..ed388addbfb292718db205f1e5e13b2561a25c5f 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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
 
 #include "ct.h"
 
 
 #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;
 
 {
   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)
 {
 }
 
 void 
 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
 {
-  if (m_pBackprojectImplem)
+  if (m_pBackprojectImplem != NULL)
     m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
 }
 
     m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
 }
 
@@ -54,7 +54,7 @@ Backprojector::~Backprojector (void)
 //     and initializes the backprojector
 
 bool
 //     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;
 {
   m_nameBackproject = backprojName;
   m_nameInterpolation = interpName;
@@ -78,17 +78,17 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
   }
 
   if (m_idBackproject == BPROJ_TRIG)
   }
 
   if (m_idBackproject == BPROJ_TRIG)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_TABLE)
   else if (m_idBackproject == BPROJ_TABLE)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_DIFF)
   else if (m_idBackproject == BPROJ_DIFF)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_DIFF2)
   else if (m_idBackproject == BPROJ_DIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_IDIFF2)
   else if (m_idBackproject == BPROJ_IDIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_IDIFF3)
   else if (m_idBackproject == BPROJ_IDIFF3)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
   else {
     m_fail = true;
     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
   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;
     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;
 #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);
     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);
 #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.
 
 // 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();
 {
   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;
 
   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)
 }
 
 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.
 
 // 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);
 {
   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
 
 //   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;
 {
   // 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
        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
       }
     }  // 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<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
        
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> ((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) {
   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) {
        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);
        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
       }        // end for iy
     } //end linear
   } // end for ix
index 65163f3185b0be2de8beb7342ad643f0701ff875..0f334b9177d3ae7c9816ecb539f5f7f728ecfea5 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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"
 
 
 #include "ct.h"
 
 
+int SignalFilter::N_INTEGRAL=500;  //static member
+
 /* NAME
  *   SignalFilter::SignalFilter     Construct a signal
  *
 /* NAME
  *   SignalFilter::SignalFilter     Construct a signal
  *
  *   int nSignalPoints Number of points in signal
  *   double param      General input parameter to filters
  *   int domain                FREQUENCY or SPATIAL domain wanted
  *   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;
 {
   m_vecFilter = NULL;
   m_vecFourierCosTable = NULL;
@@ -70,15 +70,15 @@ SignalFilter::SignalFilter (const char* filterName, const char* filterMethodName
     m_failMessage += domainName;
     return;
   }
     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;
 {
   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_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;
   m_filterParam = param;  
-  m_numIntegral = numIntegral;
   m_idFilter = convertFilterNameToID (filterName);
   if (m_idFilter == FILTER_INVALID) {
     m_fail = true;
   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
 }
 
 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;
 {
   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_fail = false;
   m_nSignalPoints = nSignalPoints;
   m_signalInc = signalIncrement;
-  m_filterParam = param;  
+  m_filterParam = filterParam;  
   m_zeropad = zeropad;
   m_zeropad = zeropad;
+  m_preinterpolationFactor = preinterpolationFactor;
 
   m_vecFourierCosTable = NULL;
   m_vecFourierSinTable = NULL;
 
   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_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;
     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) {
 
   // 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
   }
 
 #if HAVE_FFTW
@@ -186,19 +187,21 @@ SignalFilter::init (const FilterID filterID, const FilterMethodID filterMethodID
   }
 
   if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
   }
 
   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_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_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) {
       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_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_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
 
   }
 #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_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) {
     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++)
       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++)
     } 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
        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;
     } 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;
        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;
     }
     if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
        rfftw_destroy_plan(m_realPlanForward);
        rfftw_destroy_plan(m_realPlanBackward);
        delete [] m_vecRealFftInput;
+       delete [] m_vecRealFftSignal;
     }
 #endif
 }
     }
 #endif
 }
@@ -391,7 +395,6 @@ SignalFilter::convertDomainIDToName (const DomainID domain)
   return (name);
 }
 
   return (name);
 }
 
-
 void
 SignalFilter::filterSignal (const float input[], double output[]) const
 {
 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];
 
       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];
 
   } 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++) {
       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
 }
   }
 #endif
 }
@@ -465,7 +470,7 @@ SignalFilter::response (double x)
   double response = 0;
 
   if (m_idDomain == DOMAIN_SPATIAL)
   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);
 
   else if (m_idDomain == DOMAIN_FREQUENCY)
     response = frequencyResponse (m_idFilter, m_bw, x, m_filterParam);
 
@@ -474,12 +479,12 @@ SignalFilter::response (double x)
 
 
 double 
 
 
 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 spatialResponseAnalytic (filterID, bw, x, param);
   else
-    return spatialResponseCalc (filterID, bw, x, param, nIntegral);
+    return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
 }
 
 /* NAME
 }
 
 /* NAME
@@ -498,9 +503,9 @@ SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double pa
  */
 
 double 
  */
 
 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 
 }
 
 double 
@@ -632,6 +637,28 @@ SignalFilter::spatialResponseAnalytic (double x, double param) const
   return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
 }
 
   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)
 {
 double 
 SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param)
 {
index c3d2f67896820bd65fec117fad0b3958cb4b0c11..9ae326f73ed6de8e3bd1038ef736c645d7c7b5b5 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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<unsigned int> (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<unsigned int> (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<unsigned int> (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
+
index db6431411cff6a435d78a2dfb31a9407a7a1f95b..8f1edf7b4e7f8ed32757b0c522cae49264784751 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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
  */
 
 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 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) {
   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;
     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;
 #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()) {
   filter.setTraceLevel(trace);
 
   if (filter.fail()) {
@@ -540,7 +540,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
   }
 #endif  //HAVE_SGP
 
   }
 #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;
   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  ("xporigin .5.");
       ezset  ("xlength .5.");
       ezset ("box");
+
       ezset ("grid");
       ezset ("grid");
-      gid = ezplot (filteredProj, plot_xaxis, m_nDet);
+      gid = ezplot (filteredProj, plot_xaxis, n_filteredProj);
     }
 #endif  //HAVE_SGP
 
     }
 #endif  //HAVE_SGP
 
index 9b48e9c22d45f1da8313fad2af0cd42cf96934e7..203ab187997d24c6e1607e0de169d7d98dd6988a 100644 (file)
@@ -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
 
 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)
 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 (file)
index 57b5549..0000000
+++ /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 <ct.h>
-#include <iostream>
-
-int
-main (int argc, char **argv)
-{
-  cout << "Hello, CTSim!" << endl;
-
-  return (0);
-}
index 8211d6d3ea7284cbe0932d126d841c2769ef1ad2..e6ec07396fc96efb7bd099962558442b589d84b8 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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
 #include "ct.h"
 
 
 #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 };
 
 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 << "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;
   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
   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;
   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;
 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;
   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))) {
     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)
        if (! opt_set_max)
-         densmax = maxfound;
+           densmax = max;
        if (! opt_set_min)
        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)
          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)
   }
   
   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)
   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)
 #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)
   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) 
 #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);
 #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();
     cio_kb_getc();
-    //    sgp2_close(sgp2_get_active_win());
+    sgp2_close(sgp2_get_active_win());
 #endif
   }
   else
 #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<unsigned int> (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<unsigned int> (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<unsigned int> (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[])
 #ifndef NO_MAIN
 int 
 main (int argc, char *const argv[])
index 40a8cba5633e0636aab73888f9498fcf6eafcc6c..d6a229caeb35cdc7ffb0ff3ccc24716a9dc1159c 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
 **  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
 **
 **  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
 #include "timer.h"
 
 
 #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},
 
 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},
   {"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
 #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;
   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;
   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;
   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_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;
        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]);
          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]);
          }
          break;
        case O_ZEROPAD:
          optZeroPad = strtol(optarg, &endptr, 10);
          if (endptr != optarg + strlen(optarg)) {
            pjrec_usage(argv[0]);
+           return(1);
          }
          break;
        case O_VERBOSE:
          }
          break;
        case O_VERBOSE:
@@ -228,7 +240,7 @@ pjrec_main (int argc, char * argv[])
       filterDesc << optFilterName;
 
     ostringstream label;
       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)
     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 (&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);
   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());
 
 #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");
 
   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
   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
 #endif
 
 #ifdef HAVE_MPI