From: Kevin M. Rosenberg Date: Sat, 17 Jun 2000 20:12:15 +0000 (+0000) Subject: r97: Converted Scanner and Projections to C++ X-Git-Tag: debian-4.5.3-3~920 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=f4a23943110823118f35756dd41fbd6707f04511 r97: Converted Scanner and Projections to C++ --- diff --git a/ChangeLog b/ChangeLog index 784dfb0..1fcd9b2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,12 @@ +1.9.1 - 6/20/2000 + Converted all file i/o to C++ library + +1.9.0 - 6/15/2000 + Skip versions to make version 2.0 the first C++ version + Renamed Raysums to Projections throughout files + Renamed Detector to Scanner + Converted Scanner and Projecions to full object-oriented + 0.6.1 - 6/12/2000 Converted Phantom and PhantomElements to Object-Oriented Converted Detector to Object-Oriented diff --git a/cgi-bin/ctsim.cgi.in b/cgi-bin/ctsim.cgi.in index 6b0da48..898aa89 100755 --- a/cgi-bin/ctsim.cgi.in +++ b/cgi-bin/ctsim.cgi.in @@ -22,7 +22,7 @@ CGI::ReadParse(\%in); # Incoming variables # Phantom_Name, Phantom_Nx, Phantom_Ny, Phantom_NSample -# RS_NDet, RS_NRot, RS_NRay, RS_RotAngle, +# PJ_NDet, PJ_NRot, PJ_NRay, PJ_RotAngle, # IR_Nx, IR_Ny, IR_Filter, IR_Filter_Param my $error = ""; @@ -35,13 +35,13 @@ $error .= "Phantom name must not be blank
" if ($Phantom_Name eq ""); $error .= "Phantom NX and NY must be between 5 and 1024
" if ($Phantom_Nx < 5 || $Phantom_Nx > 1024 || $Phantom_Ny < 5 || $Phantom_Ny > 1024); $error .= "Phantom NSample must be between 1 and 10
" if ($Phantom_NSample < 1 || $Phantom_NSample > 10); -my $RS_NDet = FilterToNumber($in{'RS_NDet'}); -my $RS_NRot = FilterToNumber($in{'RS_NRot'}); -my $RS_NRay = FilterToNumber($in{'RS_NRay'}); -my $RS_RotAngle = FilterToNumber($in{'RS_RotAngle'}); -$error .= "Raysum NDet must be between 5 and 1800
" if ($RS_NDet < 5 || $RS_NDet > 1800); -$error .= "Raysum NRot must be between 5 and 2048
" if ($RS_NRot < 5 || $RS_NRot > 2048); -$error .= "Raysum RotAngle must be between 0.1 and 2
" if ($RS_RotAngle < 0.1 || $RS_RotAngle > 2); +my $PJ_NDet = FilterToNumber($in{'PJ_NDet'}); +my $PJ_NRot = FilterToNumber($in{'PJ_NRot'}); +my $PJ_NRay = FilterToNumber($in{'PJ_NRay'}); +my $PJ_RotAngle = FilterToNumber($in{'PJ_RotAngle'}); +$error .= "Projection NDet must be between 5 and 1800
" if ($PJ_NDet < 5 || $PJ_NDet > 1800); +$error .= "Projection NRot must be between 5 and 2048
" if ($PJ_NRot < 5 || $PJ_NRot > 2048); +$error .= "Projection RotAngle must be between 0.1 and 2
" if ($PJ_RotAngle < 0.1 || $PJ_RotAngle > 2); #my $IR_Nx = FilterToNumber($in{'IR_Nx'}); #my $IR_Ny = FilterToNumber($in{'IR_Ny'}); @@ -74,38 +74,38 @@ $error .= "IR Filter Parameter must be between 0 and 1
" if ($IR_Filter_Param my $tmpid = $$; my $auto_window_img = "std0.1"; my $auto_window_diff = "std1"; -my $auto_window_rs = "full"; +my $auto_window_pj = "full"; my $logfile = "$::jobdir/ctsim.log"; my $result_fname = "$::datadir/result-$tmpid.html"; my $phantom_fname = "$::datadir/phantom-$tmpid.if"; -my $rs_fname = "$::datadir/rs-$tmpid.rs"; +my $pj_fname = "$::datadir/pj-$tmpid.pj"; my $ir_fname = "$::datadir/ir-$tmpid.if"; -my $rs_if_fname = "$::datadir/rs-$tmpid.if"; +my $pj_if_fname = "$::datadir/pj-$tmpid.if"; my $diff_fname = "$::datadir/diff-$tmpid.if"; my $phantom_png = "$::datadir/phantom-$tmpid.png"; my $ir_png = "$::datadir/ir-$tmpid.png"; -my $rs_png = "$::datadir/rs-$tmpid.png"; +my $pj_png = "$::datadir/pj-$tmpid.png"; my $diff_png = "$::datadir/diff-$tmpid.png"; my $result_url = "$::url_datadir/result-$tmpid.html"; my $phantom_png_url = "$::url_datadir/phantom-$tmpid.png"; my $ir_png_url = "$::url_datadir/ir-$tmpid.png"; -my $rs_png_url = "$::url_datadir/rs-$tmpid.png"; +my $pj_png_url = "$::url_datadir/pj-$tmpid.png"; my $diff_png_url = "$::url_datadir/diff-$tmpid.png"; my $ctrec_ver = "$::bindir/ctrec"; -my $phm2rs_ver = "$::bindir/phm2rs"; +my $phm2pj_ver = "$::bindir/phm2pj"; my $phm2if_ver = "$::bindir/phm2if"; my $diff_ver = "$::bindir/if-2"; $ctrec_ver = "/opt/lam/bin/mpirun N N $::lamrundir/ctrec-lam" if $MPI; -$phm2rs_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2rs-lam" if $MPI; +$phm2pj_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2pj-lam" if $MPI; $phm2if_ver = "/opt/lam/bin/mpirun N N $::lamrundir/phm2if-lam" if $MPI; my $gp_cmd = "$phm2if_ver $phantom_fname $Phantom_Nx $Phantom_Ny --phantom $Phantom_Name --nsample $Phantom_NSample"; -my $rs_cmd = "$phm2rs_ver $rs_fname $RS_NDet $RS_NRot --phantom $Phantom_Name --nray $RS_NRay --rotangle $RS_RotAngle"; -my $rs_if_cmd = "$::bindir/rs2if $rs_fname $rs_if_fname"; -my $ir_cmd = "$ctrec_ver $rs_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj"; +my $pj_cmd = "$phm2pj_ver $pj_fname $PJ_NDet $PJ_NRot --phantom $Phantom_Name --nray $PJ_NRay --rotangle $PJ_RotAngle"; +my $pj_if_cmd = "$::bindir/pj2if $pj_fname $pj_if_fname"; +my $ir_cmd = "$ctrec_ver $pj_fname $ir_fname $IR_Nx $IR_Ny --filter $IR_Filter --filter-param $IR_Filter_Param --interp $IR_Interp --backproj $IR_Backproj"; my $diff_cmd = "$diff_ver $phantom_fname $ir_fname $diff_fname --comp"; my $window_options = "--auto $auto_window_img"; @@ -118,7 +118,7 @@ if ($Disp_Max ne 'auto') { my $png1_cmd = "$::bindir/if2img $phantom_fname $phantom_png $window_options --stats --format png"; my $png2_cmd = "$::bindir/if2img $ir_fname $ir_png $window_options --stats --format png"; -my $png3_cmd = "$::bindir/if2img $rs_if_fname $rs_png --auto $auto_window_rs --stats --format png"; +my $png3_cmd = "$::bindir/if2img $pj_if_fname $pj_png --auto $auto_window_pj --stats --format png"; my $png4_cmd = "$::bindir/if2img $diff_fname $diff_png --auto $auto_window_diff --stats --format png"; my $title = "CT Simulation Results"; @@ -130,7 +130,7 @@ $out .= "

$title


\n"; if ($opt_d) { $out .= "

Commands

\n"; - $out .= "$gp_cmd
\n$rs_cmd
\n$rs_if_cmd
\n$ir_cmd
\n$diff_cmd
\n$png1_cmd
\n$png2_cmd
\n" . + $out .= "$gp_cmd
\n$pj_cmd
\n$pj_if_cmd
\n$ir_cmd
\n$diff_cmd
\n$png1_cmd
\n$png2_cmd
\n" . "$png3_cmd
\n$png4_cmd
\n"; } @@ -141,21 +141,21 @@ if ($error ne "") { $out .= "$error"; } else { my $gp_out; - my $rs_out; - my $rs_if_out; + my $pj_out; + my $pj_if_out; my $ir_out; my $diff_out; my $png_gp_out; my $png_ir_out; - my $png_rs_out; + my $png_pj_out; my $png_diff_out; $gp_out = `$gp_cmd`; if (-s $phantom_fname) { - $rs_out .= `$rs_cmd`; + $pj_out .= `$pj_cmd`; $png_gp_out .= `$png1_cmd`; - if (-s $rs_fname) { - $rs_if_out .= `$rs_if_cmd`; - $png_rs_out .= `$png3_cmd`; + if (-s $pj_fname) { + $pj_if_out .= `$pj_if_cmd`; + $png_pj_out .= `$png3_cmd`; $ir_out .= `$ir_cmd`; if (-s $ir_fname) { $png_ir_out .= `$png2_cmd`; @@ -165,7 +165,7 @@ if ($error ne "") { } } - $cmdout = "$gp_cmd\n $gp_out $rs_cmd\n $rs_out $rs_if_cmd\n $rs_if_out $ir_cmd\n $ir_out $diff_cmd\n $diff_out $png1_cmd\n $png_gp_out $png2_cmd\n $png_ir_out $png3_cmd\n $png_rs_out $png4_cmd\n $png_diff_out"; + $cmdout = "$gp_cmd\n $gp_out $pj_cmd\n $pj_out $pj_if_cmd\n $pj_if_out $ir_cmd\n $ir_out $diff_cmd\n $diff_out $png1_cmd\n $png_gp_out $png2_cmd\n $png_ir_out $png3_cmd\n $png_pj_out $png4_cmd\n $png_diff_out"; if (open(LOGFILE,">> $logfile")) { flock(LOGFILE,LOCK_EX); seek(LOGFILE, 0, 2); @@ -182,17 +182,17 @@ if ($error ne "") { } my $png_gp_out_html = $png_gp_out; my $png_ir_out_html = $png_ir_out; - my $png_rs_out_html = $png_rs_out; + my $png_pj_out_html = $png_pj_out; my $png_diff_out_html = $png_diff_out; $png_gp_out_html =~ s/\n/
/gms; $png_ir_out_html =~ s/\n/
/gms; - $png_rs_out_html =~ s/\n/
/gms; + $png_pj_out_html =~ s/\n/
/gms; $png_diff_out_html =~ s/\n/
/gms; $out .= "\n"; $out .= "\n"; $out .= "\n"; - $out .= "\n"; - $out .= "\n"; + $out .= "\n"; + $out .= "\n"; $out .= "\n"; $out .= "
Phantom ImageReconstructed Image

$png_gp_out

$png_ir_out
Raysum SinusoidPhantom/Reconst Error

$png_rs_out
Projection SinusoidPhantom/Reconst Error

$png_pj_out

$diff_out
$png_diff_out
"; $out .= "Execution time: $execution_time seconds\n"; @@ -227,10 +227,10 @@ if (open(JOBFILES,"> $::jobdir/$tmpid")) print JOBFILES "Phantom_Nx=$Phantom_Nx\n"; print JOBFILES "Phantom_Ny=$Phantom_Nx\n"; print JOBFILES "Phantom_NSample=$Phantom_NSample\n"; - print JOBFILES "RS_NDet=$RS_NDet\n"; - print JOBFILES "RS_NRot=$RS_NRot\n"; - print JOBFILES "RS_NRay=$RS_NRay\n"; - print JOBFILES "RS_RotAngle=$RS_RotAngle\n"; + print JOBFILES "PJ_NDet=$PJ_NDet\n"; + print JOBFILES "PJ_NRot=$PJ_NRot\n"; + print JOBFILES "PJ_NRay=$PJ_NRay\n"; + print JOBFILES "PJ_RotAngle=$PJ_RotAngle\n"; print JOBFILES "IR_Nx=$IR_Nx\n"; print JOBFILES "IR_Ny=$IR_Ny\n"; print JOBFILES "IR_Interp=$IR_Interp\n"; @@ -240,7 +240,7 @@ if (open(JOBFILES,"> $::jobdir/$tmpid")) print JOBFILES "Disp_Min=$Disp_Min\n"; print JOBFILES "Disp_Max=$Disp_Max\n"; print JOBFILES "MPI=$MPI\n"; - print JOBFILES "Files=$result_fname,$phantom_fname,$rs_fname,$ir_fname,$phantom_png,$ir_png,$rs_if_fname,$rs_png\n" if ($error eq ""); + print JOBFILES "Files=$result_fname,$phantom_fname,$pj_fname,$ir_fname,$phantom_png,$ir_png,$pj_if_fname,$pj_png\n" if ($error eq ""); printf JOBFILES "cmdout=$cmdout\n"; flock(JOBFILES,LOCK_UN); close JOBFILES; diff --git a/configure.in b/configure.in index 9cc9819..7814d8e 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/ctrec.cpp) -AM_INIT_AUTOMAKE(ctsim,0.6.1) +AM_INIT_AUTOMAKE(ctsim,1.9.0) AM_CONFIG_HEADER(config.h) dnl Checks for programs. @@ -310,7 +310,7 @@ ctlibs="$ctlibs_base -lir $ctlibs_graphics -lkmath -lk -lcio $ctlibs_tools" AC_SUBST(ctlibs) if test -n "$lamdir" ; then - lamprograms="ctrec-lam phm2if-lam phm2rs-lam" + lamprograms="ctrec-lam phm2if-lam phm2pj-lam" AC_SUBST(lamprograms) fi diff --git a/html/simulate.html.in b/html/simulate.html.in index 5424cbf..8b855a5 100644 --- a/html/simulate.html.in +++ b/html/simulate.html.in @@ -28,10 +28,10 @@ MPI Supercomputing:
No (Single CPU)

Simulate X-Ray acquistion

-Number of detectors:

-Number of Rotations:

-Number of Rays
(samples) per detector:

-Rotation Angle
as a multiple of PI:
+Number of detectors:

+Number of Rotations:

+Number of Rays
(samples) per detector:

+Rotation Angle
as a multiple of PI:

Image Reconstruction

Interpolation:
diff --git a/include/Makefile.am b/include/Makefile.am index 74f7f81..901cad4 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -1,4 +1,5 @@ -noinst_HEADERS=ascii.h cio.h ct.h ezplot.h ir.h keyboard.h kmath.h kstddef.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h byteorder.h phantom.h timer.h sstream +noinst_HEADERS=ascii.h cio.h ct.h ezplot.h ir.h keyboard.h kmath.h kstddef.h pol.h sgp.h array2d.h imagefile.h backprojectors.h mpiworld.h byteorder.h phantom.h timer.h sstream scanner.h projections.h + diff --git a/include/backprojectors.h b/include/backprojectors.h index ea0e85c..c28a11b 100644 --- a/include/backprojectors.h +++ b/include/backprojectors.h @@ -9,8 +9,11 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: backprojectors.h,v 1.2 2000/06/13 16:20:31 kevin Exp $ +** $Id: backprojectors.h,v 1.3 2000/06/17 20:12:14 kevin Exp $ ** $Log: backprojectors.h,v $ +** Revision 1.3 2000/06/17 20:12:14 kevin +** Converted Scanner and Projections to C++ +** ** Revision 1.2 2000/06/13 16:20:31 kevin ** finished c++ conversions ** @@ -35,7 +38,7 @@ class Backproject { public: - Backproject (const RAYSUM* rs, ImageFile& im, InterpolationType interpType); + Backproject (const Projections& proj, ImageFile& im, InterpolationType interpType); virtual ~Backproject (); @@ -47,7 +50,7 @@ class Backproject void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int ni); InterpolationType interpType; - const RAYSUM* rs; + const Projections& proj; ImageFile& im; ImageFileArray v; kuint32 nx; @@ -68,8 +71,8 @@ class Backproject class BackprojectTrig : public Backproject { public: - BackprojectTrig (const RAYSUM* rs, ImageFile& im, InterpolationType interpType) - : Backproject::Backproject (rs, im, interpType) + BackprojectTrig (const Projections& proj, ImageFile& im, InterpolationType interpType) + : Backproject::Backproject (proj, im, interpType) {} void BackprojectView (const double* const t, double view_angle); @@ -79,7 +82,7 @@ class BackprojectTrig : public Backproject class BackprojectTable : public Backproject { public: - BackprojectTable (const RAYSUM* rs, ImageFile& im, InterpolationType interpType); + BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType); ~BackprojectTable (); void BackprojectView (const double* const t, double view_angle); @@ -95,7 +98,7 @@ class BackprojectTable : public Backproject class BackprojectDiff : public Backproject { public: - BackprojectDiff (const RAYSUM* rs, ImageFile& im, InterpolationType interpType); + BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType); ~BackprojectDiff (); void BackprojectView (const double* const t, double view_angle); @@ -109,8 +112,8 @@ class BackprojectDiff : public Backproject class BackprojectDiff2 : public BackprojectDiff { public: - BackprojectDiff2 (const RAYSUM* rs, ImageFile& im, InterpolationType interpType) - : BackprojectDiff::BackprojectDiff (rs, im, interpType) + BackprojectDiff2 (const Projections& proj, ImageFile& im, InterpolationType interpType) + : BackprojectDiff::BackprojectDiff (proj, im, interpType) {} void BackprojectView (const double* const t, double view_angle); @@ -119,12 +122,12 @@ class BackprojectDiff2 : public BackprojectDiff class BackprojectIntDiff2 : public BackprojectDiff { public: - BackprojectIntDiff2 (const RAYSUM* rs, ImageFile& im, InterpolationType interpType) - : BackprojectDiff::BackprojectDiff (rs, im, interpType) + BackprojectIntDiff2 (const Projections& proj, ImageFile& im, InterpolationType interpType) + : BackprojectDiff::BackprojectDiff (proj, im, interpType) {} void BackprojectView (const double* const t, double view_angle); }; -Backproject* selectBackprojector (BackprojType type, const RAYSUM*, ImageFile& im, InterpolationType interpType); +Backproject* selectBackprojector (BackprojType type, const Projections& proj, ImageFile& im, InterpolationType interpType); diff --git a/include/byteorder.h b/include/byteorder.h index a3c51e8..187e4fb 100644 --- a/include/byteorder.h +++ b/include/byteorder.h @@ -4,14 +4,14 @@ /* netorder.cpp */ void *strreverse (void *dest, const void *src, size_t n); -int read_nint16 (kuint16 *n, int fd); -int write_nint16 (kuint16 const *n, int fd); -int read_nint32 (kuint32 *n, int fd); -int write_nint32 (kuint32 const *n, int fd); -int read_nfloat32 (float *f, int fd); -int write_nfloat32 (float const *f, int fd); -int read_nfloat64 (double *d, int fd); -int write_nfloat64 (double const *d, int fd); +bool read_nint16 (kuint16 *n, int fd); +bool write_nint16 (kuint16 const *n, int fd); +bool read_nint32 (kuint32 *n, int fd); +bool write_nint32 (kuint32 const *n, int fd); +bool read_nfloat32 (float *f, int fd); +bool write_nfloat32 (float const *f, int fd); +bool read_nfloat64 (double *d, int fd); +bool write_nfloat64 (double const *d, int fd); void ConvertNetworkOrder (void* buffer, size_t bytes); void ConvertReverseNetworkOrder (void* buffer, size_t bytes); diff --git a/include/ct.h b/include/ct.h index 0289653..4817fb3 100644 --- a/include/ct.h +++ b/include/ct.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ct.h,v 1.14 2000/06/15 19:07:10 kevin Exp $ +** $Id: ct.h,v 1.15 2000/06/17 20:12:14 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 @@ -135,8 +135,10 @@ using namespace std; #include "array2d.h" #include "imagefile.h" -#include "ir.h" #include "phantom.h" +#include "projections.h" +#include "scanner.h" +#include "ir.h" #include "backprojectors.h" diff --git a/include/imagefile.h b/include/imagefile.h index 70c3fe4..ccf3769 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -17,15 +17,16 @@ private: Array2dFileLabel& operator= (const Array2dFileLabel&); public: - kfloat64 calc_time; - kuint16 label_type; - kuint16 year; - kuint16 month; - kuint16 day; - kuint16 hour; - kuint16 minute; - kuint16 second; - string label_str; + kfloat64 m_calcTime; + kuint16 m_labelType; + kuint16 m_year; + kuint16 m_month; + kuint16 m_day; + kuint16 m_hour; + kuint16 m_minute; + kuint16 m_second; + string m_strLabel; + mutable string m_strDate; static const int L_EMPTY = 0; static const int L_HISTORY = 1; @@ -39,22 +40,22 @@ public: ~Array2dFileLabel(); - string getLabelString (void) const - { return label_str; } + const string& getLabelString (void) const + { return m_strLabel; } kfloat64 getCalcTime (void) const - { return calc_time; } + { return m_calcTime; } kfloat64 getLabelType (void) const - { return label_type; } + { return m_labelType; } string& setLabelString (const char* const str) - { label_str = str; return (label_str); } + { m_strLabel = str; return (m_strLabel); } string& setLabelString (const string& str) - { label_str = str; return (label_str); } + { m_strLabel = str; return (m_strLabel); } - void getDateString (string& str) const; + const string& getDateString () const; }; @@ -118,9 +119,9 @@ public: void labelAdd (const Array2dFileLabel& label); - void labelAdd (const char* const label_str, double calc_time=0.); + void labelAdd (const char* const m_strLabel, double calc_time=0.); - void labelAdd (int type, const char* const label_str, double calc_time=0.); + void labelAdd (int type, const char* const m_strLabel, double calc_time=0.); void labelsCopy (Array2dFile& file, const char* const idStr = NULL); @@ -524,20 +525,21 @@ Array2dFile::labelRead (Array2dFileLabel& label, unsigned int label_num) return (false); } - read_nint16 (&label.label_type, file_id); - read_nint16 (&label.year, file_id); - read_nint16 (&label.month, file_id); - read_nint16 (&label.day, file_id); - read_nint16 (&label.hour, file_id); - read_nint16 (&label.minute, file_id); - read_nint16 (&label.second, file_id); - read_nfloat64 (&label.calc_time, file_id); + read_nint16 (&label.m_labelType, file_id); + read_nint16 (&label.m_year, file_id); + read_nint16 (&label.m_month, file_id); + read_nint16 (&label.m_day, file_id); + read_nint16 (&label.m_hour, file_id); + read_nint16 (&label.m_minute, file_id); + read_nint16 (&label.m_second, file_id); + read_nfloat64 (&label.m_calcTime, file_id); kuint16 strlength; read_nint16 (&strlength, file_id); - char *str = new char [strlength+1]; + char str [strlength+1]; read (file_id, str, strlength); - label.label_str = str; + str[strlength] = 0; + label.m_strLabel = str; return (true); } @@ -565,17 +567,17 @@ Array2dFile::labelAdd (const Array2dFileLabel& label) { labelSeek (num_labels); - write_nint16 (&label.label_type, file_id); - write_nint16 (&label.year, file_id); - write_nint16 (&label.month, file_id); - write_nint16 (&label.day, file_id); - write_nint16 (&label.hour, file_id); - write_nint16 (&label.minute, file_id); - write_nint16 (&label.second, file_id); - write_nfloat64 (&label.calc_time, file_id); - kuint16 strlength = label.label_str.length(); + write_nint16 (&label.m_labelType, file_id); + write_nint16 (&label.m_year, file_id); + write_nint16 (&label.m_month, file_id); + write_nint16 (&label.m_day, file_id); + write_nint16 (&label.m_hour, file_id); + write_nint16 (&label.m_minute, file_id); + write_nint16 (&label.m_second, file_id); + write_nfloat64 (&label.m_calcTime, file_id); + kuint16 strlength = label.m_strLabel.length(); write_nint16 (&strlength, file_id); - write (file_id, static_cast(label.label_str.c_str()), strlength); + write (file_id, static_cast(label.m_strLabel.c_str()), strlength); num_labels++; diff --git a/include/ir.h b/include/ir.h index c7b2f7b..8af8572 100644 --- a/include/ir.h +++ b/include/ir.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ir.h,v 1.27 2000/06/15 19:07:10 kevin Exp $ +** $Id: ir.h,v 1.28 2000/06/17 20:12:14 kevin Exp $ ** ** ** This program is free software; you can redistribute it and/or modify @@ -31,82 +31,8 @@ #include "phantom.h" - - -/*----------------------------------------------------------------------*/ -/* RAYSUM SYMBOLS */ -/*----------------------------------------------------------------------*/ - -/* Ray sums are collected along an array of ndet detectors. The data - * for these detectors is stored in the structure DETECTARRAY - */ - -typedef float DETECT_TYPE; - -struct DetectorArray -{ - DETECT_TYPE *detval; /* Pointer to array of values recorded by detector */ - int ndet; /* Number of detectors in array */ - double view_angle; /* View angle in radians */ -}; -typedef struct DetectorArray DETARRAY; - - -typedef enum { - DETECTOR_PARALLEL, - DETECTOR_EQUIANGLE, - DETECTOR_EQUILINEAR -} ScannerGeometry; - - -class Detector -{ - public: - ScannerGeometry geometry; /* Geometry of detectory */ - int ndet; /* Number of detectors in array */ - int nview; /* Number of rotated views */ - int nsample; /* Number of rays per detector */ - double detlen; /* Total length of detector array */ - double rotlen; /* Rotation angle length in radians (norm 2PI) */ - double det_inc; /* Increment between centers of detectors */ - double rot_inc; /* Increment in rotation angle between views */ - double radius; /* Radius of rotation. Distance from */ - /* center of phm to center of det */ - double phmlen; /* Maximum Length of phantom or area of interest */ - struct { - double xd1,yd1,xd2,yd2; /* Coordinates of detector endpoints */ - double xs1,ys1,xs2,ys2; /* Coordinates of source endpoints */ - double angle; /* Starting angle */ - } init; - - Detector (const Phantom& phm, const ScannerGeometry geometry, int ndet, int nview, int nsample, const double rot_anglen); - ~Detector(); - -}; - - -class Projections -{ - public: - int fd; /* Disk file descriptor */ - int file_mode; /* Current file mode (read or write) */ - int header_size; /* Size of disk file header */ - int geometry; /* Geometry of scanner */ - struct DetectorArray **view; /* Pointer to array of detarray_st pointers */ - - string remark; /* description of raysum data */ - double calctime; /* time required to calculate raysums */ - - int ndet; /* number of detectors in array */ - int nview; /* number of rotated views */ - double rot_start; /* starting view rotation */ - double rot_inc; /* angle between rotations */ - double det_start; /* distance of beginning detector to center */ - /* of PHANTOM */ - double det_inc; /* increment between detectors */ - double phmlen; /* Length of PHANTOM edge (phm is square) */ -}; -typedef class Projections RAYSUM; +#include "scanner.h" +#include "projections.h" /*----------------------------------------------------------------------*/ /* USER SYMBOLS */ @@ -156,8 +82,6 @@ typedef enum { /* Interpolation methods */ I_LINEAR /* Linear interpolation */ } InterpolationType; -static const int N_EXTRA_DETECTORS=4; /* Number of extra detectors widths when calculating detlen */ - static const char O_TRACE_NONE_STR[]= "none"; static const char O_TRACE_TEXT_STR[]= "text"; static const char O_TRACE_PHM_STR[]= "phm"; @@ -263,32 +187,6 @@ const char *name_of_filter_domain(const DomainType domain); BackprojType opt_set_backproj(const char *optarg); const char *name_of_backproj(const BackprojType backproj); -/* raycollect.cpp */ -int raysum_collect (RAYSUM *rs, const Detector& det, const Phantom& phm, const int start_view, const int trace, const int unit_pulse); -void rayview (const Phantom& phm, DETARRAY *darray, const Detector& det, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, const int unit_pulse); -double phm_ray_attenuation (const Phantom& phm, const double x1, const double y1, const double x2, const double y2); -double pelem_ray_attenuation (PhantomElement& pelem, const double x1, const double y1, const double x2, const double y2); -void raysum_trace_show_param (const char *label, const char *fmt, int row, int color, ...); - -/* rayio.cpp */ -RAYSUM *raysum_create(const char *fname, const int nview, const int ndet); -RAYSUM *raysum_create_from_det(const char *fname, const Detector& det); -RAYSUM *raysum_open(const char *filename); -void raysum_alloc_views(RAYSUM *rs); -void raysum_free(RAYSUM *rs); -int raysum_is_open(RAYSUM *rs); -int raysum_close(RAYSUM *rs); -int raysum_read_header(RAYSUM *rs); -int raysum_write_header(RAYSUM *rs); -int raysum_read(RAYSUM *rs); -int raysum_write(RAYSUM *rs); -void raysum_print(const RAYSUM *rs); -void raysum_print_info(const RAYSUM *rs); -DETARRAY *detarray_alloc(const int n); -void detarray_free(DETARRAY *darray); -int detarray_read(RAYSUM *rs, DETARRAY *darray, const int view_num); -int detarray_write(RAYSUM *rs, const DETARRAY *darray, const int view_num); - /* From phm2image.cpp */ void phm_to_imagefile (const Phantom& phm, ImageFile& im, const int col_start, const int col_count, const int nsample, const int trace); @@ -298,36 +196,6 @@ int image_display (const ImageFile& im); int image_display_scale (const ImageFile& im, const int scale, const double pmin, const double pmax); /* From reconstr.cpp */ -ImageFile& proj_reconst (ImageFile& im, RAYSUM *rs, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, int const ir_trace); - -/* From bproj.cpp */ -void backproj_init (const RAYSUM *rs, ImageFile& im, const BackprojType bproj_method); -int backproj_calc (const RAYSUM *rs, ImageFile& im, const double *t, const double view_angle, const int interp_type, const int bproj_method); -void backproj_term (const RAYSUM *rs, ImageFile& im, const int bproj_method); - -void backproj_init_trig (const RAYSUM *rs, ImageFile& im); -int backproj_calc_trig (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_trig (const RAYSUM *rs, ImageFile& im); -void backproj_init_table (const RAYSUM *rs, ImageFile& im); -int backproj_calc_table (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_table (const RAYSUM *rs, ImageFile& im); -void backproj_init_d (const RAYSUM *rs, ImageFile& im); -int backproj_calc_d (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_d (const RAYSUM *rs, ImageFile& im); -void backproj_init_d2 (const RAYSUM *rs, ImageFile& im); -int backproj_calc_d2 (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_d2 (const RAYSUM *rs, ImageFile& im); -void backproj_init_id (const RAYSUM *rs, ImageFile& im); -int backproj_calc_id (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_id (const RAYSUM *rs, ImageFile& im); -void backproj_init_id2 (const RAYSUM *rs, ImageFile& im); -int backproj_calc_id2 (const RAYSUM *rs, ImageFile& im, const double *t, - const double view_angle, const int interp_type); -void backproj_term_id2 (const RAYSUM *rs, ImageFile& im); +ImageFile& proj_reconst (ImageFile& im, Projections& rs, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, int const ir_trace); #endif diff --git a/include/kmath.h b/include/kmath.h index 4837026..bb30f05 100644 --- a/include/kmath.h +++ b/include/kmath.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: kmath.h,v 1.13 2000/06/15 19:07:10 kevin Exp $ +** $Id: kmath.h,v 1.14 2000/06/17 20:12:14 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 @@ -64,6 +64,9 @@ template inline T clamp (T value, T lowerBounds, T upperBounds) { return (value >= upperBounds ? upperBounds : (value <= lowerBounds ? lowerBounds : value )); } +template +inline T lineLength (T x1, T y1, T x2, T y2) +{ return static_cast( sqrt ((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) ); } /* clip.cpp */ int clip_rect(double *x1, double *y1, double *x2, double *y2, const double rect[4]); diff --git a/include/projections.h b/include/projections.h new file mode 100644 index 0000000..bbf4305 --- /dev/null +++ b/include/projections.h @@ -0,0 +1,95 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: projections.h +** Purpose: Header file for Projections class +** Programmer: Kevin Rosenberg +** Date Started: July 1, 1984 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: projections.h,v 1.1 2000/06/17 20:12:14 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 +** 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 +******************************************************************************/ + +#ifndef PROJECTIONS_H +#define PROJECTIONS_H + +class Scanner; +class DetectorArray; + +// Projections +class Projections +{ + public: + Projections (const int nView, const int nDet); + Projections (const Scanner& scanner); + Projections (void); + ~Projections (void); + + void init (const int nView, const int nDet); + + void printProjectionData (void); + void printScanInfo (void) const; + + bool read (const char* fname); + bool write (const char* fname); + bool detarrayRead (DetectorArray& darray, const int view_num); + bool detarrayWrite (const DetectorArray& darray, const int view_num); + + void setRotInc (double rotInc) { m_rotInc = rotInc;} + void setDetInc (double detInc) { m_detInc = detInc;} + void setPhmLen (double phmLen) { m_phmLen = phmLen;} + void setCalcTime (double calcTime) {m_calcTime = calcTime;} + void setRemark (const char* remark) {m_remark = remark;} + void setRemark (const string remark) {m_remark = remark;} + + const double detStart(void) const {return m_detStart;} + const double rotStart(void) const {return m_rotStart;} + const double calcTime(void) const {return m_calcTime;} + const double phmLen(void) const {return m_phmLen;} + const char* remark(void) const {return m_remark.c_str();} + const double detInc(void) const {return m_detInc;} + const double rotInc(void) const {return m_rotInc;} + const int nDet(void) const {return m_nDet;} + const int nView(void) const {return m_nView;} + DetectorArray& getDetectorArray (const int iview) + { return (*m_projData[iview]); } + + private: + int m_fd; // Disk file descriptor + int m_headerSize; // Size of disk file header + int m_geometry; // Geometry of scanner + struct DetectorArray **m_projData; // Pointer to array of detarray_st pointers + string m_remark; // description of raysum data + int m_nDet; // number of detectors in array + int m_nView; // number of rotated views + double m_calcTime; // time required to calculate raysums + double m_rotStart; // starting view rotation + double m_rotInc; // angle between rotations + double m_detStart; // distance of beginning detector to center phantom + double m_detInc; // increment between detectors + double m_phmLen; // Length of phantom edge (phm is square) + + bool headerRead (void); + bool headerWrite (void); + void newProjData (void); + void deleteProjData (void); +}; + + +#endif diff --git a/include/scanner.h b/include/scanner.h new file mode 100644 index 0000000..8c5ab45 --- /dev/null +++ b/include/scanner.h @@ -0,0 +1,116 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: scanner.h +** Purpose: Scanner header file +** Programmer: Kevin Rosenberg +** Date Started: July 1, 1984 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: scanner.h,v 1.1 2000/06/17 20:12:14 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 +** 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 +******************************************************************************/ + +#ifndef SCANNER_H +#define SCANNER_H + +#include "projections.h" + +// Projections are collected along an array of ndet detectors. The data +// for these detectors is stored in the class DetectorArray + +typedef float DetectorValue; + +class DetectorArray +{ + public: + DetectorArray (const int ndet); + ~DetectorArray (void); + + const int nDet(void) const {return m_nDet;} + const double viewAngle(void) const {return m_viewAngle;} + DetectorValue* detValues(void) {return m_detValues;} + const DetectorValue* detValues(void) const {return m_detValues;} + + void setViewAngle (double viewAngle) + { m_viewAngle = viewAngle; } + + private: + DetectorValue* m_detValues; /* Pointer to array of values recorded by detector */ + int m_nDet; /* Number of detectors in array */ + double m_viewAngle; /* View angle in radians */ +}; + +typedef enum { + DETECTOR_PARALLEL, + DETECTOR_EQUIANGLE, + DETECTOR_EQUILINEAR +} ScannerGeometry; + + +class Scanner +{ + public: + Scanner (const Phantom& phm, const ScannerGeometry geometry, int nDet, int nView, int nSample, const double rot_anglen); + ~Scanner(); + + void collectProjections (Projections& proj, const Phantom& phm, const int start_view, const int trace); + + const unsigned int nDet(void) const {return m_nDet;} + const unsigned int nView(void) const {return m_nView;} + const double phmLen(void) const {return m_phmLen;} + const double rotInc(void) const {return m_rotInc;} + const double detInc(void) const {return m_detInc;} + const double radius(void) const {return m_radius;} + + + private: + void projectSingleView (const Phantom& phm, DetectorArray& darray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2); + + double projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2); + + double projectLineAgainstPElem (const PhantomElement& pelem, const double x1, const double y1, const double x2, const double y2); + + void traceShowParam (const char *label, const char *fmt, int row, int color, ...); + + + ScannerGeometry m_geometry; /* Geometry of detectory */ + unsigned int m_nDet; /* Number of detectors in array */ + unsigned int m_nView; /* Number of rotated views */ + unsigned int m_nSample; /* Number of rays per detector */ + double m_detLen; /* Total length of detector array */ + double m_rotLen; /* Rotation angle length in radians (norm 2PI) */ + double m_detInc; /* Increment between centers of detectors */ + double m_rotInc; /* Increment in rotation angle between views */ + double m_radius; /* Radius of rotation. Distance from */ + /* center of phm to center of det */ + double m_phmLen; /* Maximum Length of phantom or area of interest */ + + struct { + double xd1,yd1,xd2,yd2; /* Coordinates of detector endpoints */ + double xs1,ys1,xs2,ys2; /* Coordinates of source endpoints */ + double angle; /* Starting angle */ + } m_initPos; + + static const int N_EXTRA_DETECTORS=4; /* Number of extra detectors widths when calculating detlen */ + + int m_trace; + +}; + +#endif diff --git a/src/Makefile.am b/src/Makefile.am index 7fda6a4..ec47a60 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ -bin_PROGRAMS = ctsim ctrec phm2rs rs2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img +bin_PROGRAMS = ctsim ctrec phm2pj pj2if @lamprograms@ ifinfo phm2if if-1 if-2 if2img bin_SCRIPTS = sample-ctrec.sh -EXTRA_PROGRAMS = ctrec-lam phm2rs-lam phm2if-lam +EXTRA_PROGRAMS = ctrec-lam phm2pj-lam phm2if-lam INCLUDES=@my_includes@ EXTRA_DIST=Makefile.nt mpiworld.cpp @@ -8,14 +8,14 @@ ctsim_SOURCES = ctsim.cpp ctsim_LDADD = @ctlibs@ ctrec_SOURCES = ctrec.cpp ctrec_LDADD=@ctlibs@ -phm2rs_SOURCES=phm2rs.cpp -phm2rs_LDADD=@ctlibs@ +phm2pj_SOURCES=phm2pj.cpp +phm2pj_LDADD=@ctlibs@ phm2if_SOURCES = phm2if.cpp phm2if_LDADD=@ctlibs@ if2img_SOURCES = if2img.cpp if2img_LDADD=@ctlibs@ -rs2if_SOURCES = rs2if.cpp -rs2if_LDADD=@ctlibs@ +pj2if_SOURCES = pj2if.cpp +pj2if_LDADD=@ctlibs@ if_1_SOURCES=if-1.cpp if_1_LDADD=@ctlibs@ if_2_SOURCES=if-2.cpp @@ -27,11 +27,11 @@ ctrec_lam_SOURCES=ctrec.cpp ctrec_lam_LDADD=@ctlamlibs@ phm2if_lam_SOURCES=phm2if.cpp phm2if_lam_LDADD=@ctlamlibs@ -phm2rs_lam_SOURCES=phm2rs.cpp -phm2rs_lam_LDADD=@ctlamlibs@ +phm2pj_lam_SOURCES=phm2pj.cpp +phm2pj_lam_LDADD=@ctlamlibs@ realclean: - rm -f *.pgm *.if *~ *.rs + rm -f *.pgm *.if *~ *.pj if USE_LAM CC_LAM = $(lamdir)/bin/balky @@ -40,16 +40,16 @@ LAM_EXTRA_SRC = mpiworld.cpp ctrec-lam: ctrec.cpp mpiworld.cpp $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI ctrec.cpp -o ctrec-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ -phm2rs-lam: phm2rs.cpp mpiworld.cpp - $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2rs.cpp -o phm2rs-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ +phm2pj-lam: phm2pj.cpp mpiworld.cpp + $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2pj.cpp -o phm2pj-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ phm2if-lam: phm2if.cpp mpiworld.cpp $(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2if.cpp -o phm2if-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@ endif -shared: ctrec.cpp phm2rs.cpp - $(CC) @DEFS@ $(CFLAGS) $(INCLUDES) -DNO_MAIN -shared phm2rs.cpp ctrec.cpp @ctlibs@ $(LDFLAGS) -o ir.so +shared: ctrec.cpp phm2pj.cpp + $(CC) @DEFS@ $(CFLAGS) $(INCLUDES) -DNO_MAIN -shared phm2pj.cpp ctrec.cpp @ctlibs@ $(LDFLAGS) -o ir.so diff --git a/src/ctrec.cpp b/src/ctrec.cpp index 586f2b6..0ec8402 100644 --- a/src/ctrec.cpp +++ b/src/ctrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctrec.cpp,v 1.9 2000/06/15 19:07:10 kevin Exp $ +** $Id: ctrec.cpp,v 1.10 2000/06/17 20:12:15 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 @@ -95,7 +95,7 @@ ctrec_usage (const char *program) #ifdef HAVE_MPI -static void ScatterProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int debug); +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int debug); static void ReduceImageMPI (MPIWorld& mpiWorld, ImageFile* imLocal, ImageFile* imGlobal); #endif @@ -104,8 +104,8 @@ int ctrec_main (int argc, char * argv[]) { ImageFile *imGlobal = NULL; - RAYSUM *rsGlobal = NULL; - char *rs_name, *im_filename = NULL; + Projections projGlobal; + char *pj_name, *im_filename = NULL; string remark; char filt_name[80]; char *endptr; @@ -120,7 +120,6 @@ ctrec_main (int argc, char * argv[]) int nx, ny; #ifdef HAVE_MPI ImageFile *imLocal; - RAYSUM *rsLocal; int mpi_nview, mpi_ndet; double mpi_detinc, mpi_rotinc, mpi_phmlen; MPIWorld mpiWorld (argc, argv); @@ -198,7 +197,7 @@ ctrec_main (int argc, char * argv[]) return (1); } - rs_name = argv[optind]; + pj_name = argv[optind]; im_filename = argv[optind + 1]; @@ -223,16 +222,15 @@ ctrec_main (int argc, char * argv[]) #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) { - rsGlobal = raysum_open (rs_name); - raysum_read (rsGlobal); + projGlobal.read (pj_name); if (opt_verbose) - raysum_print_info(rsGlobal); + projGlobal.printScanInfo(); - mpi_ndet = rsGlobal->ndet; - mpi_nview = rsGlobal->nview; - mpi_detinc = rsGlobal->det_inc; - mpi_phmlen = rsGlobal->phmlen; - mpi_rotinc = rsGlobal->rot_inc; + mpi_ndet = projGlobal.nDet(); + mpi_nview = projGlobal.nView(); + mpi_detinc = projGlobal.detInc(); + mpi_phmlen = projGlobal.phmLen(); + mpi_rotinc = projGlobal.rotInc(); } TimerCollectiveMPI timerBcast (mpiWorld.getComm()); @@ -256,15 +254,13 @@ ctrec_main (int argc, char * argv[]) mpiWorld.setTotalWorkUnits (mpi_nview); - rsLocal = raysum_create (NULL, mpiWorld.getMyLocalWorkUnits(), mpi_ndet); - rsLocal->ndet = mpi_ndet; - rsLocal->nview = mpi_nview; - rsLocal->det_inc = mpi_detinc; - rsLocal->phmlen = mpi_phmlen; - rsLocal->rot_inc = mpi_rotinc; + Projections projLocal (mpiWorld.getMyLocalWorkUnits(), mpi_ndet); + projLocal.setDetInc (mpi_detinc); + projLocal.setPhmLen (mpi_phmlen); + projLocal.setRotInc (mpi_rotinc); TimerCollectiveMPI timerScatter (mpiWorld.getComm()); - ScatterProjectionsMPI (mpiWorld, rsGlobal, rsLocal, opt_debug); + ScatterProjectionsMPI (mpiWorld, projGlobal, projLocal, opt_debug); if (opt_verbose) timerScatter.timerEndAndReport ("Time to scatter projections"); @@ -275,10 +271,9 @@ ctrec_main (int argc, char * argv[]) imLocal = new ImageFile (nx, ny); #else - rsGlobal = raysum_open (rs_name); - raysum_read (rsGlobal); + projGlobal.read (pj_name); if (opt_verbose) - raysum_print_info(rsGlobal); + projGlobal.printScanInfo(); imGlobal = new ImageFile (im_filename, nx, ny); imGlobal->fileCreate(); @@ -286,7 +281,7 @@ ctrec_main (int argc, char * argv[]) #ifdef HAVE_MPI TimerCollectiveMPI timerReconstruct (mpiWorld.getComm()); - proj_reconst (*imLocal, rsLocal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); + proj_reconst (*imLocal, projLocal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); if (opt_verbose) timerReconstruct.timerEndAndReport ("Time to reconstruct"); @@ -295,21 +290,20 @@ ctrec_main (int argc, char * argv[]) if (opt_verbose) timerReduce.timerEndAndReport ("Time to reduce image"); #else - proj_reconst (*imGlobal, rsGlobal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); + proj_reconst (*imGlobal, projGlobal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace); #endif #ifdef HAVE_MPI if (mpiWorld.getRank() == 0) #endif { - raysum_close (rsGlobal); - double calctime = timerProgram.timerEnd(); + double calcTime = timerProgram.timerEnd(); imGlobal->arrayDataWrite (); - imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, rsGlobal->remark.c_str(), rsGlobal->calctime); - imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calctime); + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, projGlobal.remark(), projGlobal.calcTime()); + imGlobal->labelAdd (Array2dFileLabel::L_HISTORY, remark.c_str(), calcTime); imGlobal->fileClose (); if (opt_verbose) - cout << "Run time: " << calctime << " seconds" << endl; + cout << "Run time: " << calcTime << " seconds" << endl; } #ifdef HAVE_MPI MPI::Finalize(); @@ -325,25 +319,34 @@ ctrec_main (int argc, char * argv[]) ////////////////////////////////////////////////////////////////////////////////////// #ifdef HAVE_MPI -static void ScatterProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug) +static void ScatterProjectionsMPI (MPIWorld& mpiWorld, Projections& projGlobal, Projections& projLocal, const int opt_debug) { if (mpiWorld.getRank() == 0) { for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) { for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) { - mpiWorld.getComm().Send(&rsGlobal->view[iw]->ndet, 1, MPI::INT, iProc, 0); - mpiWorld.getComm().Send(&rsGlobal->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0); - mpiWorld.getComm().Send(rsGlobal->view[iw]->detval, rsGlobal->ndet, MPI::FLOAT, iProc, 0); + DetectorArray& detarray = projGlobal.getDetectorArray( iw ); + int nDet = detarray.nDet(); + DetectorValue* detval = detarray.detValues(); + + double viewAngle = detarray.viewAngle(); + mpiWorld.getComm().Send(&nDet, 1, MPI::INT, iProc, 0); + mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, iProc, 0); + mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, iProc, 0); } } } for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) { MPI::Status status; - mpiWorld.getComm().Recv(&rsLocal->view[iw]->ndet, 1, MPI::INT, 0, 0, status); - mpiWorld.getComm().Recv(&rsLocal->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0, status); - mpiWorld.getComm().Recv(rsLocal->view[iw]->detval, rsLocal->ndet, MPI::FLOAT, 0, 0, status); + int nDet; + double viewAngle; + DetectorValue* detval = projLocal.getDetectorArray(iw).detValues(); + + mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, 0, 0, status); + mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, 0, 0, status); + mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, 0, 0, status); + projLocal.getDetectorArray(iw).setViewAngle( viewAngle ); } - rsLocal->nview = mpiWorld.getMyLocalWorkUnits(); } static void diff --git a/src/if2img.cpp b/src/if2img.cpp index 6971fbb..85237a7 100644 --- a/src/if2img.cpp +++ b/src/if2img.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if2img.cpp,v 1.4 2000/06/13 16:20:31 kevin Exp $ +** $Id: if2img.cpp,v 1.5 2000/06/17 20:12:15 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 @@ -292,16 +292,14 @@ if2img_main (int argc, char *const argv[]) for (int i = 0; i < nlabels; i++) { Array2dFileLabel label; im.labelRead (label, i); - string str; - label.getDateString (str); if (label.getLabelType() == Array2dFileLabel::L_HISTORY) { cout << "History: " << label.getLabelString() << endl; cout << " calc time = " << label.getCalcTime() << " secs" << endl; - cout << " Timestamp = " << str << endl; + cout << " Timestamp = " << label.getDateString() << endl; } else if (label.getLabelType() == Array2dFileLabel::L_USER) { cout << "Note: " << label.getLabelString() << endl; - cout << " Timestamp = %s" << str << endl; + cout << " Timestamp = %s" << label.getDateString() << endl; } } } diff --git a/src/ifinfo.cpp b/src/ifinfo.cpp index 7c771f7..a054f2a 100644 --- a/src/ifinfo.cpp +++ b/src/ifinfo.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ifinfo.cpp,v 1.6 2000/06/13 16:20:31 kevin Exp $ +** $Id: ifinfo.cpp,v 1.7 2000/06/17 20:12:15 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 @@ -142,17 +142,16 @@ ifinfo_main (int argc, char *const argv[]) Array2dFileLabel label; im->labelRead (label, i); - string str; - label.getDateString (str); - if (label.getLabelType() == Array2dFileLabel::L_HISTORY) { - cout << "History: " << label.getLabelString() << endl; + cout << "History: " << endl; + cout << " " << label.getLabelString() << endl; cout << " calc time = " << label.getCalcTime() << " secs" << endl; - cout << " Timestamp = " << str << endl; + cout << " Timestamp = " << label.getDateString() << endl; } else if (label.getLabelType() == Array2dFileLabel::L_USER) { cout << "Note: " << label.getLabelString() << endl; - cout << " Timestamp = %s" << str << endl; - } + cout << " Timestamp = %s" << label.getDateString() << endl; + } + cout << endl; } } @@ -221,13 +220,13 @@ ifinfo_main (int argc, char *const argv[]) } } stddev = sqrt(stddev / (nx * ny)); - cout << "nx=" << nx << endl; - cout << "nx=" << ny << endl; - cout << "min=" << minfound << endl; - cout << "max=" << maxfound << endl; - cout << "mean=" << mean << endl; - cout << "mode=" << mode << endl; - cout << "stddef=" << stddev << endl; + cout << " nx: " << nx << endl; + cout << " nx: " << ny << endl; + cout << " min: " << minfound << endl; + cout << " max: " << maxfound << endl; + cout << " mean: " << mean << endl; + cout << " mode: " << mode << endl; + cout << "stddef: " << stddev << endl; } return (0); diff --git a/src/phm2pj.cpp b/src/phm2pj.cpp new file mode 100644 index 0000000..fe83154 --- /dev/null +++ b/src/phm2pj.cpp @@ -0,0 +1,342 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: phm2pj.cpp +** Purpose: Take projections of a phantom object +** Programmer: Kevin Rosenberg +** Date Started: 1984 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: phm2pj.cpp,v 1.1 2000/06/17 20:12:15 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 +** 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 "timer.h" + + +enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, + O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION }; + +static struct option phm2pj_options[] = +{ + {"phantom", 1, 0, O_PHANTOM}, + {"phmfile", 1, 0, O_PHMFILE}, + {"desc", 1, 0, O_DESC}, + {"nray", 1, 0, O_NRAY}, + {"rotangle", 1, 0, O_ROTANGLE}, + {"trace", 1, 0, O_TRACE}, + {"verbose", 0, 0, O_VERBOSE}, + {"help", 0, 0, O_HELP}, + {"debug", 0, 0, O_DEBUG}, + {"version", 0, 0, O_VERSION}, + {0, 0, 0, 0} +}; + + +void +phm2pj_usage (const char *program) +{ + cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]" << endl; + cout << "Calculate (projections) through phantom object, either a predefined --phantom or a --phmfile" << endl; + cout << "" << endl; + cout << " outfile Name of output file for raysums" << endl; + cout << " ndet Number of detectors" << endl; + cout << " nview Number of rotated views" << endl; + cout << " --phantom Phantom to use for projection" << endl; + cout << " herman Herman head phantom" << endl; + cout << " rowland Rowland head phantom" << endl; + cout << " browland Bordered Rowland head phantom" << endl; + cout << " unitpulse Unit pulse phantom" << endl; + cout << " --phmfile Get Phantom from phantom file" << endl; + cout << " --desc Description of raysum" << endl; + cout << " --nray Number of rays per detector (default = 1)" << endl; + cout << " --rotangle Degrees to rotate view through, multiple of PI (default = 1)" << endl; + cout << " --trace Trace level to use" << endl; + cout << " none No tracing (default)" << endl; + cout << " text Trace text level" << endl; + cout << " phm Trace phantom image" << endl; + cout << " rays Trace rays" << endl; + cout << " plot Trace plot" << endl; + cout << " clipping Trace clipping" << endl; + cout << " --verbose Verbose mode" << endl; + cout << " --debug Debug mode" << endl; + cout << " --version Print version" << endl; + cout << " --help Print this help message" << endl; +} + +#ifdef HAVE_MPI +void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug); +#endif + +int +phm2pj_main (int argc, char* argv[]) +{ + Phantom phm; + ScannerGeometry opt_geometry = DETECTOR_PARALLEL; + char *opt_outfile = NULL; + string opt_desc; + string opt_phmfilename; + int opt_ndet, opt_nview; + int opt_nray = 1; + int opt_trace = 0; + int opt_phmnum = -1; + int opt_verbose = 0; + int opt_debug = 0; + double opt_rotangle = 1; + char* endptr = NULL; + char* endstr; + +#ifdef HAVE_MPI + MPIWorld mpiWorld (argc, argv); +#endif + + Timer timerProgram; + +#ifdef HAVE_MPI + if (mpiWorld.getRank() == 0) { +#endif + while (1) { + int c = getopt_long(argc, argv, "", phm2pj_options, NULL); + + if (c == -1) + break; + + switch (c) { + case O_PHANTOM: + if ((opt_phmnum = opt_set_phantom (optarg)) < 0) { + phm2pj_usage(argv[0]); + return (1); + } + phm.create (opt_phmnum); + break; + case O_PHMFILE: +#ifdef HAVE_MPI + if (mpiWorld.getRank() == 0) + cerr << "Can not read phantom from file in MPI mode" << endl; + return (1); +#endif + opt_phmfilename = optarg; + phm.createFromFile (opt_phmfilename.c_str()); + break; + case O_VERBOSE: + opt_verbose = 1; + break; + case O_DEBUG: + opt_debug = 1; + break; + break; + case O_TRACE: + if ((opt_trace = opt_set_trace(optarg)) < 0) { + phm2pj_usage(argv[0]); + return (1); + } + break; + case O_DESC: + opt_desc = optarg; + break; + case O_ROTANGLE: + opt_rotangle = strtod(optarg, &endptr); + endstr = optarg + strlen(optarg); + if (endptr != endstr) { + cerr << "Error setting --rotangle to " << optarg << endl; + phm2pj_usage(argv[0]); + return (1); + } + break; + case O_NRAY: + opt_nray = strtol(optarg, &endptr, 10); + endstr = optarg + strlen(optarg); + if (endptr != endstr) { + cerr << "Error setting --nray to %s" << optarg << endl; + phm2pj_usage(argv[0]); + return (1); + } + break; + case O_VERSION: +#ifdef VERSION + cout << "Version: " << VERSION << endl; +#else + cout << "Unknown version number" << endl; +#endif + exit(0); + case O_HELP: + case '?': + phm2pj_usage(argv[0]); + return (0); + default: + phm2pj_usage(argv[0]); + return (1); + } + } + + if (phm.nPElem() == 0) { + cerr << "No phantom defined" << endl; + phm2pj_usage(argv[0]); + return (1); + } + if (optind + 3 != argc) { + phm2pj_usage(argv[0]); + return (1); + } + + opt_outfile = argv[optind]; + opt_ndet = strtol(argv[optind+1], &endptr, 10); + endstr = argv[optind+1] + strlen(argv[optind+1]); + if (endptr != endstr) { + cerr << "Error setting --ndet to " << argv[optind+1] << endl; + phm2pj_usage(argv[0]); + return (1); + } + opt_nview = strtol(argv[optind+2], &endptr, 10); + endstr = argv[optind+2] + strlen(argv[optind+2]); + if (endptr != endstr) { + cerr << "Error setting --nview to " << argv[optind+2] << endl; + phm2pj_usage(argv[0]); + return (1); + } + + ostringstream desc; + desc << "Raysum_Collect: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", "; + if (opt_phmfilename.length()) { + desc << "PhantomFile=" << opt_phmfilename; + } else if (opt_phmnum != -1) { + desc << "Phantom=" << name_of_phantom(opt_phmnum); + } + if (opt_desc.length()) { + desc << ": " << opt_desc; + } + opt_desc = desc.str(); +#ifdef HAVE_MPI + } +#endif + +#ifdef HAVE_MPI + mpiWorld.getComm().Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0); + mpiWorld.getComm().Bcast (&opt_nview, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_ndet, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_nray, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_phmnum, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0); + mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0); + + if (mpiWorld.getRank() > 0 && opt_phmnum >= 0) + phm.create (opt_phmnum); +#endif + + opt_rotangle *= PI; + Scanner scanner (phm, opt_geometry, opt_ndet, opt_nview, opt_nray, opt_rotangle); + +#ifdef HAVE_MPI + mpiWorld.setTotalWorkUnits (opt_nview); + + Projections pjGlobal; + if (mpiWorld.getRank() == 0) { + pjGlobal = Projections (scanner); + } + + Scanner localScanner (phm, opt_geometry, opt_ndet, mpiWorld.getMyLocalWorkUnits(), opt_nray, opt_rotangle); + Projections pjLocal (localScanner); + if (opt_debug) + cout << "pjLocal->nview = " << pjLocal.nView() << " (process " << mpiWorld.getRank() << ")" << endl;; + + TimerMPI timerProject (mpiWorld.getComm()); + localScanner.collectProjections (pjLocal, phm, mpiWorld.getMyStartWorkUnit(), opt_trace); + if (opt_verbose) + timerProject.timerEndAndReport ("Time to collect projections"); + + TimerMPI timerGather (mpiWorld.getComm()); + GatherProjectionsMPI (mpiWorld, pjGlobal, pjLocal, opt_debug); + if (opt_verbose) + timerGather.timerEndAndReport ("Time to gather projections"); +#else + Projections pjGlobal (scanner); + scanner.collectProjections (pjGlobal, phm, 0, opt_trace); +#endif + +#ifdef HAVE_MPI + if (mpiWorld.getRank() == 0) +#endif + { + pjGlobal.setCalcTime (timerProgram.timerEnd()); + pjGlobal.setRemark (opt_desc); + pjGlobal.write (opt_outfile); + if (opt_verbose) { + phm.print(); + cout << endl; + pjGlobal.printScanInfo(); + cout << endl; + cout << "Remark: " << pjGlobal.remark() << endl; + cout << "Run time: " << pjGlobal.calcTime() << " seconds" << endl; + } + } + + return (0); +} + + +/* FUNCTION + * GatherProjectionsMPI + * + * SYNOPSIS + * Gather's raysums from all processes in pjLocal in pjGlobal in process 0 + */ + +#ifdef HAVE_MPI +void GatherProjectionsMPI (MPIWorld& mpiWorld, Projections& pjGlobal, Projections& pjLocal, const int opt_debug) +{ + for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) { + DetectorArray& detArray = pjLocal.getDetectorArray(iw); + double viewAngle = detArray.viewAngle(); + int nDet = detArray.nDet(); + DetectorValue* detval = detArray.detValues(); + + mpiWorld.getComm().Send(&viewAngle, 1, MPI::DOUBLE, 0, 0); + mpiWorld.getComm().Send(&nDet, 1, MPI::INT, 0, 0); + mpiWorld.getComm().Send(detval, nDet, MPI::FLOAT, 0, 0); + } + + if (mpiWorld.getRank() == 0) { + for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) { + for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) { + MPI::Status status; + double viewAngle; + int nDet; + DetectorArray detArray = pjGlobal.getDetectorArray(iw); + DetectorValue* detval = detArray.detValues(); + + mpiWorld.getComm().Recv(&viewAngle, 1, MPI::DOUBLE, iProc, 0, status); + mpiWorld.getComm().Recv(&nDet, 1, MPI::INT, iProc, 0, status); + mpiWorld.getComm().Recv(detval, nDet, MPI::FLOAT, iProc, 0, status); + detArray.setViewAngle (viewAngle); + } + } + } + +} +#endif + + +#ifndef NO_MAIN +int +main (int argc, char* argv[]) +{ + return (phm2pj_main(argc, argv)); +} +#endif + diff --git a/src/phm2rs.cpp b/src/phm2rs.cpp deleted file mode 100644 index 0459349..0000000 --- a/src/phm2rs.cpp +++ /dev/null @@ -1,335 +0,0 @@ -/***************************************************************************** -** FILE IDENTIFICATION -** -** Name: phm2rs.cpp -** Purpose: Take projections of a phantom object -** Programmer: Kevin Rosenberg -** Date Started: 1984 -** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** $Id: phm2rs.cpp,v 1.6 2000/06/15 19:07:10 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 -** 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 "timer.h" - - -enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, - O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION }; - -static struct option phm2rs_options[] = -{ - {"phantom", 1, 0, O_PHANTOM}, - {"phmfile", 1, 0, O_PHMFILE}, - {"desc", 1, 0, O_DESC}, - {"nray", 1, 0, O_NRAY}, - {"rotangle", 1, 0, O_ROTANGLE}, - {"trace", 1, 0, O_TRACE}, - {"verbose", 0, 0, O_VERBOSE}, - {"help", 0, 0, O_HELP}, - {"debug", 0, 0, O_DEBUG}, - {"version", 0, 0, O_VERSION}, - {0, 0, 0, 0} -}; - - -void -phm2rs_usage (const char *program) -{ - cout << "usage: " << fileBasename(program) << " outfile ndet nview [--phantom phantom-name] [--phmfile filename] [OPTIONS]" << endl; - cout << "Calculate raysums (projections) through phantom object, either" << endl; - cout << "a predefined --phantom or a --phmfile" << endl; - cout << "" << endl; - cout << " outfile Name of output file for raysums" << endl; - cout << " ndet Number of detectors" << endl; - cout << " nview Number of rotated views" << endl; - cout << " --phantom Phantom to use for projection" << endl; - cout << " herman Herman head phantom" << endl; - cout << " rowland Rowland head phantom" << endl; - cout << " browland Bordered Rowland head phantom" << endl; - cout << " unitpulse Unit pulse phantom" << endl; - cout << " --phmfile Get Phantom from phantom file" << endl; - cout << " --desc Description of raysum" << endl; - cout << " --nray Number of rays per detector (default = 1)" << endl; - cout << " --rotangle Degrees to rotate view through, multiple of PI (default = 1)" << endl; - cout << " --trace Trace level to use" << endl; - cout << " none No tracing (default)" << endl; - cout << " text Trace text level" << endl; - cout << " phm Trace phantom image" << endl; - cout << " rays Trace rays" << endl; - cout << " plot Trace plot" << endl; - cout << " clipping Trace clipping" << endl; - cout << " --verbose Verbose mode" << endl; - cout << " --debug Debug mode" << endl; - cout << " --version Print version" << endl; - cout << " --help Print this help message" << endl; -} - -#ifdef HAVE_MPI -void GatherProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug); -#endif - -int -phm2rs_main (int argc, char* argv[]) -{ - Detector *det; - Phantom phm; - RAYSUM *rsGlobal = NULL; - char *opt_outfile = NULL; - string opt_desc; - string opt_phmfilename; - char *endptr, *endstr; - int opt_ndet, opt_nview; - int opt_nray = 1; - int opt_trace = 0; - int opt_phmnum = -1; - int opt_verbose = 0; - int opt_debug = 0; - double opt_rotangle = 1; - -#ifdef HAVE_MPI - RAYSUM *rsLocal; - MPIWorld mpiWorld (argc, argv); -#endif - - Timer timerProgram; - -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) { -#endif - while (1) { - int c = getopt_long(argc, argv, "", phm2rs_options, NULL); - char *endptr = NULL; - char *endstr; - - if (c == -1) - break; - - switch (c) { - case O_PHANTOM: - if ((opt_phmnum = opt_set_phantom (optarg)) < 0) { - phm2rs_usage(argv[0]); - return (1); - } - phm.create (opt_phmnum); - break; - case O_PHMFILE: -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) - cerr << "Can not read phantom from file in MPI mode" << endl; - return (1); -#endif - opt_phmfilename = optarg; - phm.createFromFile (opt_phmfilename.c_str()); - break; - case O_VERBOSE: - opt_verbose = 1; - break; - case O_DEBUG: - opt_debug = 1; - break; - break; - case O_TRACE: - if ((opt_trace = opt_set_trace(optarg)) < 0) { - phm2rs_usage(argv[0]); - return (1); - } - break; - case O_DESC: - opt_desc = optarg; - break; - case O_ROTANGLE: - opt_rotangle = strtod(optarg, &endptr); - endstr = optarg + strlen(optarg); - if (endptr != endstr) { - cerr << "Error setting --rotangle to " << optarg << endl; - phm2rs_usage(argv[0]); - return (1); - } - break; - case O_NRAY: - opt_nray = strtol(optarg, &endptr, 10); - endstr = optarg + strlen(optarg); - if (endptr != endstr) { - cerr << "Error setting --nray to %s" << optarg << endl; - phm2rs_usage(argv[0]); - return (1); - } - break; - case O_VERSION: -#ifdef VERSION - cout << "Version: " << VERSION << endl; -#else - cout << "Unknown version number" << endl; -#endif - exit(0); - case O_HELP: - case '?': - phm2rs_usage(argv[0]); - return (0); - default: - phm2rs_usage(argv[0]); - return (1); - } - } - - if (phm.nPElem() == 0) { - cerr << "No phantom defined" << endl; - phm2rs_usage(argv[0]); - return (1); - } - if (optind + 3 != argc) { - phm2rs_usage(argv[0]); - return (1); - } - - opt_outfile = argv[optind]; - opt_ndet = strtol(argv[optind+1], &endptr, 10); - endstr = argv[optind+1] + strlen(argv[optind+1]); - if (endptr != endstr) { - cerr << "Error setting --ndet to " << argv[optind+1] << endl; - phm2rs_usage(argv[0]); - return (1); - } - opt_nview = strtol(argv[optind+2], &endptr, 10); - endstr = argv[optind+2] + strlen(argv[optind+2]); - if (endptr != endstr) { - cerr << "Error setting --nview to " << argv[optind+2] << endl; - phm2rs_usage(argv[0]); - return (1); - } - - ostringstream desc; - desc << "Raysum_Collect: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", "; - if (opt_phmfilename.length()) { - desc << "PhantomFile=" << opt_phmfilename; - } else if (opt_phmnum != -1) { - desc << "Phantom=" << name_of_phantom(opt_phmnum); - } - if (opt_desc.length()) { - desc << ": " << opt_desc; - } - opt_desc = desc.str(); -#ifdef HAVE_MPI - } -#endif - -#ifdef HAVE_MPI - mpiWorld.getComm().Bcast (&opt_rotangle, 1, MPI::DOUBLE, 0); - mpiWorld.getComm().Bcast (&opt_nview, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_ndet, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_nray, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_phmnum, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_verbose, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_debug, 1, MPI::INT, 0); - mpiWorld.getComm().Bcast (&opt_trace, 1, MPI::INT, 0); - - if (mpiWorld.getRank() > 0 && opt_phmnum >= 0) - phm.create (opt_phmnum); -#endif - - opt_rotangle *= PI; - det = new Detector (phm, DETECTOR_PARALLEL, opt_ndet, opt_nview, opt_nray, opt_rotangle); - -#ifdef HAVE_MPI - mpiWorld.setTotalWorkUnits (opt_nview); - - if (mpiWorld.getRank() == 0) { - rsGlobal = raysum_create_from_det (opt_outfile, *det); - raysum_alloc_views(rsGlobal); - } - - rsLocal = raysum_create_from_det (NULL, *det); - rsLocal->nview = mpiWorld.getMyLocalWorkUnits(); - if (opt_debug) - cout << "rsLocal->nview = " << rsLocal->nview << " (process " << mpiWorld.getRank() << ")" << endl;; - - TimerMPI timerProject (mpiWorld.getComm()); - raysum_collect (rsLocal, *det, phm, mpiWorld.getMyStartWorkUnit(), opt_trace, FALSE); - if (opt_verbose) - timerProject.timerEndAndReport ("Time to collect projections"); - - TimerMPI timerGather (mpiWorld.getComm()); - GatherProjectionsMPI (mpiWorld, rsGlobal, rsLocal, opt_debug); - if (opt_verbose) - timerGather.timerEndAndReport ("Time to gather projections"); -#else - rsGlobal = raysum_create_from_det (opt_outfile, *det); - raysum_collect (rsGlobal, *det, phm, 0, opt_trace, FALSE); -#endif - -#ifdef HAVE_MPI - if (mpiWorld.getRank() == 0) -#endif - { - rsGlobal->calctime = timerProgram.timerEnd(); - rsGlobal->remark = opt_desc; - raysum_write (rsGlobal); - raysum_close (rsGlobal); - if (opt_verbose) { - cout << "Remark: " << rsGlobal->remark << endl; - cout << "Run time: " << rsGlobal->calctime << " seconds" << endl; - } - } - - delete det; - - return (0); -} - - -/* FUNCTION - * GatherProjectionsMPI - * - * SYNOPSIS - * Gather's raysums from all processes in rsLocal in rsGlobal in process 0 - */ - -#ifdef HAVE_MPI -void GatherProjectionsMPI (MPIWorld& mpiWorld, RAYSUM *rsGlobal, RAYSUM *rsLocal, const int opt_debug) -{ - for (int iw = 0; iw < mpiWorld.getMyLocalWorkUnits(); iw++) { - mpiWorld.getComm().Send(&rsLocal->view[iw]->view_angle, 1, MPI::DOUBLE, 0, 0); - mpiWorld.getComm().Send(&rsLocal->view[iw]->ndet, 1, MPI::INT, 0, 0); - mpiWorld.getComm().Send(rsLocal->view[iw]->detval, rsLocal->ndet, MPI::FLOAT, 0, 0); - } - - if (mpiWorld.getRank() == 0) { - for (int iProc = 0; iProc < mpiWorld.getNumProcessors(); iProc++) { - for (int iw = mpiWorld.getStartWorkUnit(iProc); iw <= mpiWorld.getEndWorkUnit(iProc); iw++) { - MPI::Status status; - - mpiWorld.getComm().Recv(&rsGlobal->view[iw]->view_angle, 1, MPI::DOUBLE, iProc, 0, status); - mpiWorld.getComm().Recv(&rsGlobal->view[iw]->ndet, 1, MPI::INT, iProc, 0, status); - mpiWorld.getComm().Recv(rsGlobal->view[iw]->detval, rsGlobal->ndet, MPI::FLOAT, iProc, 0, status); - } - } - } - -} -#endif - - -#ifndef NO_MAIN -int -main (int argc, char* argv[]) -{ - return (phm2rs_main(argc, argv)); -} -#endif - diff --git a/src/pj2if.cpp b/src/pj2if.cpp new file mode 100644 index 0000000..c77f901 --- /dev/null +++ b/src/pj2if.cpp @@ -0,0 +1,145 @@ +/***************************************************************************** +** FILE IDENTIFICATION +** +** Name: proj2if.cpp +** Purpose: Convert an projection data file to an image file +** Programmer: Kevin Rosenberg +** Date Started: April 2000 +** +** This is part of the CTSim program +** Copyright (C) 1983-2000 Kevin Rosenberg +** +** $Id: pj2if.cpp,v 1.1 2000/06/17 20:12:15 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 +** 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 +******************************************************************************/ + +/* FILE + * proj2if.c Convert Raysum to image + * + * DATE + * Apr 1999 + */ + +#include "ct.h" + + +enum { O_VERBOSE, O_HELP, O_VERSION }; + +static struct option my_options[] = +{ + {"verbose", 0, 0, O_VERBOSE}, + {"help", 0, 0, O_HELP}, + {"version", 0, 0, O_VERSION}, + {0, 0, 0, 0} +}; + + +void +proj2if_usage (const char *program) +{ + fprintf(stdout, "usage: %s in-proj-file out-if-file [OPTIONS]\n", fileBasename(program)); + fprintf(stdout, "This program converts a projection file to a IF file\n"); + fprintf(stdout, "\n"); + fprintf(stdout, " --verbose Verbose mode\n"); + fprintf(stdout, " --version Print version\n"); + fprintf(stdout, " --help Print this help message\n"); +} + + + +int +proj2if_main (const int argc, char *const argv[]) +{ + char *proj_name, *im_name; + int ix, iy; + bool opt_verbose = false; + extern int optind; + + while (1) + { + int c = getopt_long (argc, argv, "", my_options, NULL); + if (c == -1) + break; + + switch (c) + { + case O_VERBOSE: + opt_verbose = true; + break; + case O_VERSION: +#ifdef VERSION + fprintf(stdout, "Version %s\n", VERSION); +#else + fprintf(stderr, "Unknown version number"); +#endif + exit(0); + case O_HELP: + case '?': + proj2if_usage(argv[0]); + return (0); + default: + proj2if_usage(argv[0]); + return (1); + } + } + + if (argc - optind != 2) { + proj2if_usage(argv[0]); + return (1); + } + + proj_name = argv[optind]; + im_name = argv[optind + 1]; + + Projections proj; + if (! proj.read (proj_name)) { + sys_error (ERR_SEVERE, "Can not open projection file %s", proj_name); + return (1); + } + + if (opt_verbose) + proj.printScanInfo(); + + ImageFile im (im_name, proj.nDet(), proj.nView()); + + ImageFileArray v = im.getArray(); + + for (iy = 0; iy < proj.nView(); iy++) + { + DetectorArray& detarray = proj.getDetectorArray (iy); + const DetectorValue* detval = detarray.detValues(); + for (ix = 0; ix < proj.nDet(); ix++) + { + v[ix][iy] = detval[ix]; + } + } + + im.fileCreate (); + im.arrayDataWrite (); + im.labelAdd (Array2dFileLabel::L_HISTORY, proj.remark(), proj.calcTime()); + im.labelAdd (Array2dFileLabel::L_HISTORY, "Conversion from .pj to .if"); + im.fileClose (); + + return(0); +} + + +#ifndef NO_MAIN +int +main (const int argc, char *const argv[]) +{ + return (proj2if_main(argc, argv)); +} +#endif diff --git a/src/rs2if.cpp b/src/rs2if.cpp deleted file mode 100644 index d82ee09..0000000 --- a/src/rs2if.cpp +++ /dev/null @@ -1,155 +0,0 @@ -/***************************************************************************** -** FILE IDENTIFICATION -** -** Name: rs2if.cpp -** Purpose: Convert an projection data file to an image file -** Programmer: Kevin Rosenberg -** Date Started: April 2000 -** -** This is part of the CTSim program -** Copyright (C) 1983-2000 Kevin Rosenberg -** -** $Id: rs2if.cpp,v 1.6 2000/06/15 19:07:10 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 -** 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 -******************************************************************************/ - -/* FILE - * rs2if.c Convert Raysum to image - * - * DATE - * Apr 1999 - */ - -#include "ct.h" - - -enum { O_VERBOSE, O_HELP, O_VERSION }; - -static struct option my_options[] = -{ - {"verbose", 0, 0, O_VERBOSE}, - {"help", 0, 0, O_HELP}, - {"version", 0, 0, O_VERSION}, - {0, 0, 0, 0} -}; - - -void -rs2if_usage (const char *program) -{ - fprintf(stdout, "usage: %s in-rs-file out-if-file [OPTIONS]\n", fileBasename(program)); - fprintf(stdout, "This program converts a raysum file to a IF file\n"); - fprintf(stdout, "\n"); - fprintf(stdout, " --verbose Verbose mode\n"); - fprintf(stdout, " --version Print version\n"); - fprintf(stdout, " --help Print this help message\n"); -} - - - -int -rs2if_main (const int argc, char *const argv[]) -{ - ImageFile *im; - RAYSUM *rs; - DETARRAY *detarray; - char *rs_name, *im_name; - int ix, iy; - int opt_v = 0; - extern int optind; - - while (1) - { - int c = getopt_long (argc, argv, "", my_options, NULL); - if (c == -1) - break; - - switch (c) - { - case O_VERBOSE: - opt_v = 1; - break; - case O_VERSION: -#ifdef VERSION - fprintf(stdout, "Version %s\n", VERSION); -#else - fprintf(stderr, "Unknown version number"); -#endif - exit(0); - case O_HELP: - case '?': - rs2if_usage(argv[0]); - return (0); - default: - rs2if_usage(argv[0]); - return (1); - } - } - - if (argc - optind != 2) { - rs2if_usage(argv[0]); - return (1); - } - - rs_name = argv[optind]; - im_name = argv[optind + 1]; - - if (file_exists(rs_name) == FALSE) { - fprintf (stderr, "Raysum file %s does not exist\n", rs_name); - return (1); - } - - rs = raysum_open (rs_name); - - if (opt_v) - { - printf ("Number of detectors: %d\n", rs->ndet); - printf (" Number of views: %d\n", rs->nview); - printf (" Remark: %s\n", rs->remark.c_str()); - } - - im = new ImageFile (im_name, rs->ndet, rs->nview); - - detarray = detarray_alloc (rs->ndet); - ImageFileArray v = im->getArray(); - - for (iy = 0; iy < rs->nview; iy++) - { - detarray_read (rs, detarray, iy); - for (ix = 0; ix < rs->ndet; ix++) - { - v[ix][iy] = detarray->detval[ix]; - } - } - detarray_free (detarray); - raysum_close (rs); - - im->fileCreate (); - im->arrayDataWrite (); - im->labelAdd (Array2dFileLabel::L_HISTORY, rs->remark.c_str(), rs->calctime); - im->labelAdd (Array2dFileLabel::L_HISTORY, "Conversion from .rs to .idf"); - im->fileClose (); - - return(0); -} - - -#ifndef NO_MAIN -int -main (const int argc, char *const argv[]) -{ - return (rs2if_main(argc, argv)); -} -#endif diff --git a/src/sample-ctrec.sh.in b/src/sample-ctrec.sh.in index 3d2568c..11809a2 100755 --- a/src/sample-ctrec.sh.in +++ b/src/sample-ctrec.sh.in @@ -15,20 +15,20 @@ if [ -f sample-phm.if ] ; then fi # Simulate CT data collection and generate raysum sinugram for display -${bin}phm2rs sample-rs.rs 367 320 --nray 2 --phantom herman -if [ -f sample-rs.rs ]; then - ${bin}rs2if sample-rs.rs sample-rs.if +${bin}phm2pj sample-pj.pj 367 320 --nray 2 --phantom herman +if [ -f sample-pj.pj ]; then + ${bin}pj2if sample-pj.pj sample-pj.if fi -if [ -f sample-rs.if ]; then - ${bin}if2img sample-rs.if sample-rs.png --format png - ${bin}if2img sample-rs.if sample-rs16.png --format png16 +if [ -f sample-pj.if ]; then + ${bin}if2img sample-pj.if sample-pj.png --format png + ${bin}if2img sample-pj.if sample-pj16.png --format png16 fi # Reconstruct raysums and generate image for display -${bin}ctrec sample-rs.rs sample-rec.if 256 256 +${bin}ctrec sample-pj.pj sample-rec.if 256 256 if [ -f sample-rec.if ]; then ${bin}if2img sample-rec.if sample-rec.png --format png ${bin}if2img sample-rec.if sample-rec16.png --format png16 fi -# Files sample-phm.png, sample-rs.png, and sample-rec.png are ready for display +# Files sample-phm.png, sample-pj.png, and sample-rec.png are ready for display