+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
# 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 = "";
$error .= "Phantom NX and NY must be between 5 and 1024<br>" if ($Phantom_Nx < 5 || $Phantom_Nx > 1024 || $Phantom_Ny < 5 || $Phantom_Ny > 1024);
$error .= "Phantom NSample must be between 1 and 10<br>" 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<br>" if ($RS_NDet < 5 || $RS_NDet > 1800);
-$error .= "Raysum NRot must be between 5 and 2048<br>" if ($RS_NRot < 5 || $RS_NRot > 2048);
-$error .= "Raysum RotAngle must be between 0.1 and 2<br>" 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<br>" if ($PJ_NDet < 5 || $PJ_NDet > 1800);
+$error .= "Projection NRot must be between 5 and 2048<br>" if ($PJ_NRot < 5 || $PJ_NRot > 2048);
+$error .= "Projection RotAngle must be between 0.1 and 2<br>" if ($PJ_RotAngle < 0.1 || $PJ_RotAngle > 2);
#my $IR_Nx = FilterToNumber($in{'IR_Nx'});
#my $IR_Ny = FilterToNumber($in{'IR_Ny'});
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";
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";
if ($opt_d) {
$out .= "<H2>Commands</H2>\n";
- $out .= "$gp_cmd<br>\n$rs_cmd<br>\n$rs_if_cmd<br>\n$ir_cmd<br>\n$diff_cmd<br>\n$png1_cmd<br>\n$png2_cmd<br>\n" .
+ $out .= "$gp_cmd<br>\n$pj_cmd<br>\n$pj_if_cmd<br>\n$ir_cmd<br>\n$diff_cmd<br>\n$png1_cmd<br>\n$png2_cmd<br>\n" .
"$png3_cmd<br>\n$png4_cmd<br>\n";
}
$out .= "<i>$error</i>";
} 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`;
}
}
- $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);
}
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/<br>/gms;
$png_ir_out_html =~ s/\n/<br>/gms;
- $png_rs_out_html =~ s/\n/<br>/gms;
+ $png_pj_out_html =~ s/\n/<br>/gms;
$png_diff_out_html =~ s/\n/<br>/gms;
$out .= "<TABLE><TR><TD>Phantom Image</TD><TD>Reconstructed Image</TD></TR>\n";
$out .= "<TR><TD><IMG SRC=\"$phantom_png_url\"><br><FONT SIZE=1>$png_gp_out</FONT></TD>\n";
$out .= "<TD><IMG SRC=\"$ir_png_url\"><br><FONT SIZE=1>$png_ir_out</FONT></TD></TR>\n";
- $out .= "<TR><TD>Raysum Sinusoid</TD><TD>Phantom/Reconst Error</TD></TR>\n";
- $out .= "<TR><TD><IMG SRC=\"$rs_png_url\"><br><FONT SIZE=1>$png_rs_out</FONT></TD>\n";
+ $out .= "<TR><TD>Projection Sinusoid</TD><TD>Phantom/Reconst Error</TD></TR>\n";
+ $out .= "<TR><TD><IMG SRC=\"$pj_png_url\"><br><FONT SIZE=1>$png_pj_out</FONT></TD>\n";
$out .= "<TD><IMG SRC=\"$diff_png_url\"><br><FONT SIZE=2>$diff_out</FONT><br><FONT SIZE=1>$png_diff_out</FONT></TD></TR>\n";
$out .= "</TABLE>";
$out .= "Execution time: $execution_time seconds\n";
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";
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;
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.
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
<INPUT type="radio" name="MPI" value="no" checked>No (Single CPU)<br>
</td><td valign="top">
<h3>Simulate X-Ray acquistion</h3>
-Number of detectors: <input type="text" name="RS_NDet" size="4" value="367"><p>
-Number of Rotations: <input type="text" name="RS_NRot" size="4" value="320"><p>
-Number of Rays<br>(samples) per detector: <input type="text" name="RS_NRay" size="2" value="1"><p>
-Rotation Angle<br>as a multiple of PI: <input type="text" name="RS_RotAngle" size="3" value="1.0"><br>
+Number of detectors: <input type="text" name="PJ_NDet" size="4" value="367"><p>
+Number of Rotations: <input type="text" name="PJ_NRot" size="4" value="320"><p>
+Number of Rays<br>(samples) per detector: <input type="text" name="PJ_NRay" size="2" value="1"><p>
+Rotation Angle<br>as a multiple of PI: <input type="text" name="PJ_RotAngle" size="3" value="1.0"><br>
</td><td valign="top">
<H3>Image Reconstruction</H3>
Interpolation:<br>
-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
+
** 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
**
class Backproject
{
public:
- Backproject (const RAYSUM* rs, ImageFile& im, InterpolationType interpType);
+ Backproject (const Projections& proj, ImageFile& im, InterpolationType interpType);
virtual ~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;
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);
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);
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);
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);
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);
/* 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);
** 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
#include "array2d.h"
#include "imagefile.h"
-#include "ir.h"
#include "phantom.h"
+#include "projections.h"
+#include "scanner.h"
+#include "ir.h"
#include "backprojectors.h"
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;
~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;
};
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);
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);
}
{
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<const void*>(label.label_str.c_str()), strlength);
+ write (file_id, static_cast<const void*>(label.m_strLabel.c_str()), strlength);
num_labels++;
** 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
#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 */
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";
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);
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
** 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
inline T clamp (T value, T lowerBounds, T upperBounds)
{ return (value >= upperBounds ? upperBounds : (value <= lowerBounds ? lowerBounds : value )); }
+template<class T>
+inline T lineLength (T x1, T y1, T x2, T y2)
+{ return static_cast<T>( 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]);
--- /dev/null
+/*****************************************************************************
+** 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
--- /dev/null
+/*****************************************************************************
+** 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
-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
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
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
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
** 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
#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
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;
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);
return (1);
}
- rs_name = argv[optind];
+ pj_name = argv[optind];
im_filename = argv[optind + 1];
#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());
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");
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();
#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");
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();
//////////////////////////////////////////////////////////////////////////////////////
#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
** 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
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;
}
}
}
** 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
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;
}
}
}
}
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);
--- /dev/null
+/*****************************************************************************
+** 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
+
+++ /dev/null
-/*****************************************************************************
-** 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
-
--- /dev/null
+/*****************************************************************************
+** 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
+++ /dev/null
-/*****************************************************************************
-** 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
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