-2.0.0-b3 - 7/09/00
+2.0.0-b3 - 7/10/00
Added highly optimized backprojection method idiff3
Moved comparative stats to if-2 program from ifinfo
-
-2.0.0-b2 - 7/07/00
+ Fixed image display in if2img using G2 library
+ Added column-plot and row-plot options to if-2
+ Added autoselection of analytic/calculated spatial responses to SignalFilter
+ Added frequency-based preinterpolation to SignalFilter and idiff3
+ backprojection method. Currently, this technique is still under
+ development and debugging
+ Moved graphic file writing to ImageFile class from if2img program
+
+2.0.0-b2 - 7/07/00
Cleaned up SignalFilter class
Added zeropad option to pjrec
Added zeropad options to html and cgi files
fi
-
if test -n "$GCC"; then
CFLAGS="$CFLAGS -Wall"
fi
echo $ac_n "checking whether to enable verbose warnings""... $ac_c" 1>&6
-echo "configure:2813: checking whether to enable verbose warnings" >&5
+echo "configure:2812: checking whether to enable verbose warnings" >&5
# Check whether --enable-verbose-warnings or --disable-verbose-warnings was given.
if test "${enable_verbose_warnings+set}" = set; then
enableval="$enable_verbose_warnings"
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"
fi
if test "$debug" = "true" ; then
- CFLAGS="-g -DDEBUG"
+
+ if test -n "$GCC"; then
+ CFLAGS="$CFLAGS -g -DDEBUG"
+ fi
+
cat >> confdefs.h <<\EOF
#define DEBUG 1
EOF
else
- CFLAGS="-g -O3 -DNDEBUG"
+
+ if test -n "$GCC"; then
+ CFLAGS="$CFLAGS -g -O3 -DNDEBUG --fast-math"
+ fi
+
cat >> confdefs.h <<\EOF
#define NDEBUG 1
EOF
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"
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"
# 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
# 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"
-{ (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_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
-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.
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
-#line 3281 "configure"
+#line 3288 "configure"
#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
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
-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
# 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_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
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
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_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
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
# 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
-#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. */
; 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
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_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
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
# -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
-#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. */
; 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
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_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
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
# 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
-#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. */
; 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
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_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
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
# 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
-#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. */
; 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
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_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
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
# 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_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
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
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_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
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
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_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
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
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_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
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
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
-#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"
-{ (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*
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
AC_PROG_INSTALL
AC_PROG_MAKE_SET
AC_PROG_RANLIB
-AC_PROG_CC([-g])
+AC_PROG_CC
AC_PROG_CXX
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,
AM_CONDITIONAL(DEBUG, test "$debug" = "true")
if test "$debug" = "true" ; then
- CFLAGS="-g -DDEBUG"
+ AC_ADD_GCC_CFLAGS([-g -DDEBUG])
AC_DEFINE(DEBUG)
else
- CFLAGS="-g -O3 -DNDEBUG"
+ AC_ADD_GCC_CFLAGS([-g -O3 -DNDEBUG --fast-math])
AC_DEFINE(NDEBUG)
fi
** 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
I_3BSPLINE,
#endif
INTERP_LINEAR, // Linear interpolation
- INTERP_FREQ_PREINTERPOLATE,
+ INTERP_FREQ_PREINTERPOLATION,
} InterpolationID;
static const char BPROJ_TRIG_STR[]= "trig";
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);
static const BackprojectID convertBackprojectNameToID (const char* const bprojName);
static const char* convertBackprojectIDToName (const BackprojectID bprojID);
- bool initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName);
+ bool initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor);
};
class Backproject
{
public:
- Backproject (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType);
+ Backproject (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor);
virtual ~Backproject ();
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);
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);
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);
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);
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);
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);
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);
** 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
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);
double response (double x);
- static double spatialResponse (FilterID fType, double bw, double x, double param, int nIntegral = 0);
+ static double spatialResponse (FilterID fType, double bw, double x, double param);
static double frequencyResponse (FilterID fType, double bw, double u, double param);
- static double spatialResponseCalc (FilterID fType, double bw, double x, double param, int n);
-
static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param);
+ static double spatialResponseCalc (FilterID fType, double bw, double x, double param, int nIntegral);
+
+ static void setNumIntegral(int nIntegral) {N_INTEGRAL = nIntegral;}
+
private:
double m_bw;
int m_nFilterPoints;
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;
- fftw_complex* m_vecComplexFftInput;
+ fftw_complex* m_vecComplexFftInput, *m_vecComplexFftSignal;
fftw_plan m_complexPlanForward, m_complexPlanBackward;
#endif
FilterMethodID m_idFilterMethod;
DomainID m_idDomain;
double m_filterParam;
- int m_numIntegral;
int m_traceLevel;
int m_zeropad;
+ int m_nOutputPoints;
+ int m_preinterpolationFactor;
+ static int N_INTEGRAL;
+
+ static const bool haveAnalyticSpatial (const FilterID filterID);
static const FilterID convertFilterNameToID (const char* filterName);
static const char* convertFilterIDToName (const FilterID filterID);
static const FilterMethodID convertFilterMethodNameToID (const char* filterMethodName);
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;
};
-
#endif
** 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
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);
};
** 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
sgp2_close (SGP_ID gid)
{
#if HAVE_G2_H
+ if (gid)
g2_close (gid->g2_id);
#endif
- if (gid == _sgp2_cwin)
- _sgp2_cwin = NULL;
+ if (gid == _sgp2_cwin)
+ _sgp2_cwin = NULL;
- delete gid;
+ delete gid;
}
void
** 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
#include "ct.h"
-Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
{
m_fail = false;
m_pBackprojectImplem = NULL;
- initBackprojector (proj, im, backprojName, interpName);
+ initBackprojector (proj, im, backprojName, interpName, interpFactor);
}
void
Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
{
- if (m_pBackprojectImplem)
+ if (m_pBackprojectImplem != NULL)
m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
}
// 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;
}
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)
- 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)
- 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)
- 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)
- 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)
- 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]";
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;
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);
// 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();
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)
// 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);
// 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;
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
// 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) {
- int iDetPos = (curDetPos + halfScale) >> 16;
+ const int iDetPos = (curDetPos + halfScale) >> 16;
+ assert(iDetPos >= 0 && iDetPos < nDet);
+ *pImCol++ += filteredProj[iDetPos];
+ } // end for iy
+ } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
assert(iDetPos >= 0 && iDetPos < nDet);
*pImCol++ += filteredProj[iDetPos];
} // end for iy
} else if (interpType == Backprojector::INTERP_LINEAR) {
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
- kint32 iDetPos = curDetPos >> scaleShift;
- kint32 detRemainder = curDetPos & scaleBitmask;
-#if 0
- double frac = detRemainder * (1. / 65536.);
+ const kint32 iDetPos = curDetPos >> scaleShift;
+ const kint32 detRemainder = curDetPos & scaleBitmask;
assert(iDetPos >= 0 && iDetPos < nDet - 1);
- *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]);
-#else
- assert(iDetPos >= 0 && iDetPos < nDet - 1);
- const double* const detPointer = &filteredProj[iDetPos];
- *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale;
-#endif
+ *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
} // end for iy
} //end linear
} // end for ix
** 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
#include "ct.h"
+int SignalFilter::N_INTEGRAL=500; //static member
+
/* 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 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_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_vecFourierCosTable = NULL;
m_vecFourierSinTable = NULL;
m_filterParam = param;
- m_numIntegral = numIntegral;
m_idFilter = convertFilterNameToID (filterName);
if (m_idFilter == FILTER_INVALID) {
m_fail = true;
}
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_fail = false;
m_nSignalPoints = nSignalPoints;
m_signalInc = signalIncrement;
- m_filterParam = param;
+ m_filterParam = filterParam;
m_zeropad = zeropad;
+ m_preinterpolationFactor = preinterpolationFactor;
m_vecFourierCosTable = NULL;
m_vecFourierSinTable = NULL;
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;
// 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 (m_idFilterMethod == FILTER_METHOD_RFFTW) {
- m_complexPlanForward = m_complexPlanBackward = NULL;
m_realPlanForward = rfftw_create_plan (m_nFilterPoints, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE);
- m_realPlanBackward = rfftw_create_plan (m_nFilterPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
+ m_realPlanBackward = rfftw_create_plan (m_nOutputPoints, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE);
m_vecRealFftInput = new fftw_real [ m_nFilterPoints ];
+ m_vecRealFftSignal = new fftw_real [ m_nOutputPoints ];
for (int i = 0; i < m_nFilterPoints; i++)
m_vecRealFftInput[i] = 0;
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
- m_realPlanForward = m_realPlanBackward = NULL;
- m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
- m_complexPlanBackward = fftw_create_plan (m_nFilterPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
+ m_complexPlanForward = fftw_create_plan (m_nFilterPoints, FFTW_FORWARD, FFTW_ESTIMATE);
+ m_complexPlanBackward = fftw_create_plan (m_nOutputPoints, FFTW_BACKWARD, FFTW_ESTIMATE);
m_vecComplexFftInput = new fftw_complex [ m_nFilterPoints ];
+ m_vecComplexFftSignal = new fftw_complex [ m_nOutputPoints ];
for (int i = 0; i < m_nFilterPoints; i++)
m_vecComplexFftInput[i].re = m_vecComplexFftInput[i].im = 0;
+ for (int i = 0; i < m_nOutputPoints; i++)
+ m_vecComplexFftSignal[i].re = m_vecComplexFftSignal[i].im = 0;
}
#endif
m_filterMin = -m_signalInc * (m_nSignalPoints - 1);
m_filterMax = m_signalInc * (m_nSignalPoints - 1);
m_filterInc = (m_filterMax - m_filterMin) / (m_nFilterPoints - 1);
- m_numIntegral = numint;
m_vecFilter = new double[ m_nFilterPoints ];
if (m_idFilter == FILTER_SHEPP) {
double x;
int i;
for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++)
- m_vecFilter[i] = frequencyResponse (x, param);
+ m_vecFilter[i] = frequencyResponse (x, m_filterParam);
} else if (m_idDomain == DOMAIN_SPATIAL) {
double x;
int i;
for (x = m_filterMin, i = 0; i < m_nFilterPoints; x += m_filterInc, i++)
- if (numint == 0)
- m_vecFilter[i] = spatialResponseAnalytic (x, param);
+ if (haveAnalyticSpatial(m_idFilter))
+ m_vecFilter[i] = spatialResponseAnalytic (x, m_filterParam);
else
- m_vecFilter[i] = spatialResponseCalc (x, param, numint);
+ m_vecFilter[i] = spatialResponseCalc (x, m_filterParam);
} else {
m_failMessage = "Illegal domain name ";
m_failMessage += m_idDomain;
fftw_destroy_plan(m_complexPlanForward);
fftw_destroy_plan(m_complexPlanBackward);
delete [] m_vecComplexFftInput;
+ delete [] m_vecComplexFftSignal;
}
if (m_idFilterMethod == FILTER_METHOD_RFFTW) {
rfftw_destroy_plan(m_realPlanForward);
rfftw_destroy_plan(m_realPlanBackward);
delete [] m_vecRealFftInput;
+ delete [] m_vecRealFftSignal;
}
#endif
}
return (name);
}
-
void
SignalFilter::filterSignal (const float input[], double output[]) const
{
for (int i = 0; i < m_nSignalPoints; i++)
m_vecRealFftInput[i] = input[i];
- fftw_real out[m_nFilterPoints];
- rfftw_one (m_realPlanForward, m_vecRealFftInput, out);
- for (int i = 0; i < m_nFilterPoints; i++) {
- out[i] *= m_vecFilter[i];
- }
- fftw_real outFiltered[m_nFilterPoints];
- rfftw_one(m_realPlanBackward, out, outFiltered);
- for (int i = 0; i < m_nSignalPoints; i++)
- output[i] = outFiltered[i];
+ fftw_real fftOutput [ m_nFilterPoints ];
+ rfftw_one (m_realPlanForward, m_vecRealFftInput, fftOutput);
+ for (int i = 0; i < m_nFilterPoints; i++)
+ m_vecRealFftSignal[i] = m_vecFilter[i] * fftOutput[i];
+ for (int i = m_nFilterPoints; i < m_nOutputPoints; i++)
+ m_vecRealFftSignal[i] = 0;
+
+ fftw_real ifftOutput [ m_nOutputPoints ];
+ rfftw_one(m_realPlanBackward, m_vecRealFftSignal, ifftOutput);
+ for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
+ output[i] = ifftOutput[i];
} else if (m_idFilterMethod == FILTER_METHOD_FFTW) {
for (int i = 0; i < m_nSignalPoints; i++)
m_vecComplexFftInput[i].re = input[i];
- fftw_complex out[m_nFilterPoints];
- fftw_one(m_complexPlanForward, m_vecComplexFftInput, out);
+ fftw_complex fftOutput [ m_nFilterPoints ];
+ fftw_one(m_complexPlanForward, m_vecComplexFftInput, fftOutput);
for (int i = 0; i < m_nFilterPoints; i++) {
- out[i].re *= m_vecFilter[i];
- out[i].im *= m_vecFilter[i];
+ m_vecComplexFftSignal[i].re = m_vecFilter[i] * fftOutput[i].re;
+ m_vecComplexFftSignal[i].im = m_vecFilter[i] * fftOutput[i].im;
}
- fftw_complex outFiltered[m_nFilterPoints];
- fftw_one(m_complexPlanBackward, out, outFiltered);
- for (int i = 0; i < m_nSignalPoints; i++)
- output[i] = outFiltered[i].re;
+ fftw_complex ifftOutput [ m_nOutputPoints ];
+ fftw_one(m_complexPlanBackward, m_vecComplexFftSignal, ifftOutput);
+ for (int i = 0; i < m_nSignalPoints * m_preinterpolationFactor; i++)
+ output[i] = ifftOutput[i].re;
}
#endif
}
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);
double
-SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param, int nIntegral = 0)
+SignalFilter::spatialResponse (FilterID filterID, double bw, double x, double param)
{
- if (nIntegral == 0)
+ if (haveAnalyticSpatial(filterID))
return spatialResponseAnalytic (filterID, bw, x, param);
else
- return spatialResponseCalc (filterID, bw, x, param, nIntegral);
+ return spatialResponseCalc (filterID, bw, x, param, N_INTEGRAL);
}
/* NAME
*/
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
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)
{
** 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
}
+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
+
** 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
*/
bool
-Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
+Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const filterMethodName, const int zeropad, const char* const interpName, int interpFactor, const char* const backprojectName, const int trace)
{
int nview = m_nView;
double detInc = m_detInc;
- int n_filteredProj = m_nDet;
+ int n_filteredProj = m_nDet * interpFactor;
double filteredProj [n_filteredProj]; // filtered projections
#ifdef HAVE_BSPLINE_INTERP
int spline_order = 0, zoom_factor = 0;
if (interp_type == I_BSPLINE) {
- zoom_factor = interp_param;
+ zoom_factor = interpFactor;
spline_order = 3;
zoom_factor = 3;
n_filteredProj = (m_nDet - 1) * (zoom_factor + 1) + 1;
#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()) {
}
#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;
ezset ("xporigin .5.");
ezset ("xlength .5.");
ezset ("box");
+
ezset ("grid");
- gid = ezplot (filteredProj, plot_xaxis, m_nDet);
+ gid = ezplot (filteredProj, plot_xaxis, n_filteredProj);
}
#endif //HAVE_SGP
-bin_PROGRAMS = ctsim pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img
+SUBDIRS = ctsim
+bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img
bin_SCRIPTS = sample-ctsim.sh
EXTRA_PROGRAMS = pjrec-lam phm2pj-lam phm2if-lam
INCLUDES=@my_includes@
EXTRA_DIST=Makefile.nt mpiworld.cpp
SOURCE_DEPEND=../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a
-ctsim_SOURCES = ctsim.cpp $
-ctsim_LDADD = @ctlibs@
-ctsim_DEPENDENCIES=$(SOURCE_DEPEND)
pjrec_SOURCES = pjrec.cpp
pjrec_LDADD=@ctlibs@
pjrec_DEPENDENCIES=$(SOURCE_DEPEND)
+++ /dev/null
-/*****************************************************************************
-** 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);
-}
** 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
#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 };
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 << " 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;
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;
im.printLabels(cout);
if (opt_stats || (! (opt_set_max && opt_set_min))) {
- double minfound = HUGE_VAL, maxfound = -HUGE_VAL;
- double mode = 0, mean = 0, stddev = 0, window = 0;
- double spread;
- int hist[256];
- int ibin, nbin = 256;
- int max_bin, max_binindex;
- double maxbin;
- int nx = im.nx();
- int ny = im.ny();
- ImageFileArray v = im.getArray();
-
- maxbin = nbin - 1;
- for (int ix = 0; ix < nx; ix++) {
- for (int iy = 0; iy < ny; iy++) {
- if (v[ix][iy] > maxfound)
- maxfound = v[ix][iy];
- if (v[ix][iy] < minfound)
- minfound = v[ix][iy];
- mean += v[ix][iy];
- }
- }
- spread = maxfound - minfound;
- if (spread == 0)
- mode = minfound;
- else {
- for (ibin = 0; ibin < nbin; ibin++)
- hist[ibin] = 0;
- for (int ix = 0; ix < nx; ix++) {
- for (int iy = 0; iy < ny; iy++) {
- int b = (int) ((((v[ix][iy] - minfound) / spread) * (double) maxbin) + 0.5);
- hist[b]++;
- }
- }
- max_binindex = 0;
- max_bin = -1;
- for (ibin = 0; ibin < nbin; ibin++) {
- if (hist[ibin] > max_bin) {
- max_bin = hist[ibin];
- max_binindex = ibin;
- }
- }
- mode = (((double) max_binindex) * spread / ((double) maxbin)) + minfound;
- }
+ double min, max, mean, mode, median, stddev;
+ double window;
+ im.statistics(min, max, mean, mode, median, stddev);
- if (opt_auto == O_AUTO_FULL) {
+ if (opt_auto == O_AUTO_FULL) {
if (! opt_set_max)
- densmax = maxfound;
+ densmax = max;
if (! opt_set_min)
- densmin = minfound;
- }
- if (opt_stats || opt_auto != O_AUTO_FULL) {
- mean /= nx * ny;
- for (int ix = 0; ix < nx; ix++) {
- for (int iy = 0; iy < ny; iy++) {
- double diff = (v[ix][iy] - mean);
- stddev += diff * diff;
- }
+ densmin = min;
}
- stddev = sqrt(stddev / (nx * ny));
- if (opt_auto == O_AUTO_FULL)
- ;
- else if (opt_auto == O_AUTO_STD1)
- window = stddev;
- else if (opt_auto == O_AUTO_STD0_1)
- window = stddev * 0.1;
- else if (opt_auto == O_AUTO_STD0_5)
- window = stddev * 0.5;
+ if (opt_stats || opt_auto != O_AUTO_FULL) {
+
+ if (opt_auto == O_AUTO_FULL)
+ ;
+ else if (opt_auto == O_AUTO_STD1)
+ window = stddev;
+ else if (opt_auto == O_AUTO_STD0_1)
+ window = stddev * 0.1;
+ else if (opt_auto == O_AUTO_STD0_5)
+ window = stddev * 0.5;
else if (opt_auto == O_AUTO_STD2)
- window = stddev * 2;
- else if (opt_auto == O_AUTO_STD3)
- window = stddev * 3;
- else {
- sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto);
- return (1);
+ window = stddev * 2;
+ else if (opt_auto == O_AUTO_STD3)
+ window = stddev * 3;
+ else {
+ sys_error (ERR_SEVERE, "Internal Error: Invalid auto mode %d", opt_auto);
+ return (1);
+ }
}
- }
- if (opt_stats) {
- cout <<"nx: " << nx << endl;
- cout <<"ny: " << ny << endl;
- cout <<"min: " << minfound << endl;
- cout <<"max: " << maxfound << endl;
- cout <<"mean: " << mean << endl;
- cout <<"mode: " << mode << endl;
- cout <<"stddev: " << stddev << endl;
- }
- if (opt_auto != O_AUTO_FULL) {
- double center;
-
- if (opt_center == O_CENTER_MODE)
- center = mode;
- else if (opt_center == O_CENTER_MEAN)
- center = mean;
- else {
- sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center);
- return (1);
+ if (opt_stats) {
+ cout <<"nx: " << im.nx() << endl;
+ cout <<"ny: " << im.ny() << endl;
+ cout <<"min: " << min << endl;
+ cout <<"max: " << max << endl;
+ cout <<"mean: " << mean << endl;
+ cout <<"mode: " << mode << endl;
+ cout <<"stddev: " << stddev << endl;
}
- if (! opt_set_max)
- densmax = center + window;
- if (! opt_set_min)
- densmin = center - window;
+ if (opt_auto != O_AUTO_FULL) {
+ double center;
+
+ if (opt_center == O_CENTER_MODE)
+ center = mode;
+ else if (opt_center == O_CENTER_MEAN)
+ center = mean;
+ else {
+ sys_error (ERR_SEVERE, "Internal Error: Invalid center mode %d", opt_center);
+ return (1);
+ }
+ if (! opt_set_max)
+ densmax = center + window;
+ if (! opt_set_min)
+ densmin = center - window;
}
}
}
if (opt_format == O_FORMAT_PGM)
- sdf2d_to_pgm (im, out_file, opt_scale, opt_scale, densmin, densmax);
+ im.writeImagePGM (out_file, opt_scale, opt_scale, densmin, densmax);
else if (opt_format == O_FORMAT_PGMASC)
- sdf2d_to_pgmasc (im, out_file, opt_scale, opt_scale, densmin, densmax);
+ im.writeImagePGMASCII (out_file, opt_scale, opt_scale, densmin, densmax);
#if HAVE_PNG
else if (opt_format == O_FORMAT_PNG)
- sdf2d_to_png (im, out_file, 8, opt_scale, opt_scale, densmin, densmax);
+ im.writeImagePNG (out_file, 8, opt_scale, opt_scale, densmin, densmax);
else if (opt_format == O_FORMAT_PNG16)
- sdf2d_to_png (im, out_file, 16, opt_scale, opt_scale, densmin, densmax);
+ im.writeImagePNG (out_file, 16, opt_scale, opt_scale, densmin, densmax);
#endif
#if HAVE_GIF
else if (opt_format == O_FORMAT_GIF)
- sdf2d_to_gif (im, out_file, opt_scale, opt_scale, densmin, densmax);
+ im.writeImageGIF (out_file, opt_scale, opt_scale, densmin, densmax);
#endif
else if (opt_format == O_FORMAT_DISP) {
#if HAVE_SGP
im.displayScaling (opt_scale, densmin, densmax);
+ cout << "Press enter to continue\n";
cio_kb_getc();
- // sgp2_close(sgp2_get_active_win());
+ sgp2_close(sgp2_get_active_win());
#endif
}
else
}
-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[])
** 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
#include "timer.h"
-enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+enum {O_INTERP, O_FILTER, O_FILTER_METHOD, O_ZEROPAD, O_FILTER_PARAM, O_BACKPROJ, O_PREINTERPOLATION_FACTOR, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
static struct option my_options[] =
{
{"interp", 1, 0, O_INTERP},
+ {"preinterpolation-factor", 1, 0, O_PREINTERPOLATION_FACTOR},
{"filter", 1, 0, O_FILTER},
{"filter-method", 1, 0, O_FILTER_METHOD},
{"zeropad", 1, 0, O_ZEROPAD},
#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;
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;
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;
optFilterParam = strtod(optarg, &endptr);
if (endptr != optarg + strlen(optarg)) {
pjrec_usage(argv[0]);
+ return(1);
}
break;
case O_ZEROPAD:
optZeroPad = strtol(optarg, &endptr, 10);
if (endptr != optarg + strlen(optarg)) {
pjrec_usage(argv[0]);
+ return(1);
}
break;
case O_VERBOSE:
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)
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);
#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)
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
projects / ctsim.git / commitdiff