+1.9.6 - 6/22/2000
+ Moved conversion filter name/id to Filter class
+ Moved conversion backprojection name/id to Backproj class
+ Added MPI broadcasting of strings
+
+1.9.5 - 6/21/2000
+ Merged proj_reconstr into class Projections
+ Used auto_ptr in Projections::reconstruct to make sure destructor is always called
+ Code cleanup in projections.cpp
+ Moved conversion of phantom names/id to Phantom class
+
1.9.4 - 6/20/2000
Converted projection files to C++ library with frnetorderstream
Converted image files to C++ library with frnetorderstream
Integrate low-level X11 graphics, replace all low-level driver code in
libgraph. All the high-level graphics routines are in place.
-Integrate audio drivers into libcio -- used for the interactive CT display.
-
PROPOSED MENU STRUCTURE
=======================
PACKAGE=ctsim
-VERSION=1.9.4
+VERSION=1.9.5
if test "`cd $srcdir && pwd`" != "`pwd`" && test -f $srcdir/config.status; then
{ echo "configure: error: source directory already configured; run "make distclean" there first" 1>&2; exit 1; }
dnl CDPATH=
AC_INIT(src/ctrec.cpp)
-AM_INIT_AUTOMAKE(ctsim,1.9.4)
+AM_INIT_AUTOMAKE(ctsim,1.9.5)
AM_CONFIG_HEADER(config.h)
dnl Checks for programs.
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $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
-**
-** Revision 1.1 2000/06/10 23:00:17 kevin
-** *** empty log message ***
-**
+** $Id: backprojectors.h,v 1.4 2000/06/22 10:17:28 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
** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
******************************************************************************/
+
+#undef HAVE_BSPLINE_INTERP
+
+class Backproject;
+
+class Backprojector
+{
+ public:
+ typedef enum {
+ BPROJ_INVALID,
+ BPROJ_TRIG,
+ BPROJ_TABLE,
+ BPROJ_DIFF,
+ BPROJ_DIFF2,
+ BPROJ_IDIFF2
+} BackprojectID;
+
+ typedef enum {
+ INTERP_INVALID,
+ INTERP_NEAREST, // Nearest neighbor
+#if HAVE_BSPLINE_INTERP
+ I_BSPLINE,
+ I_1BSPLINE, // 1st order B-Spline
+ I_2BSPLINE,
+ I_3BSPLINE,
+#endif
+ INTERP_LINEAR // Linear interpolation
+ } InterpolationID;
+
+ Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName);
+
+ ~Backprojector (void);
+
+ void BackprojectView (const double* const viewData, const double viewAngle);
+
+ bool fail(void) const {return m_fail;}
+
+ private:
+ string m_nameBackproject;
+ string m_nameInterpolation;
+ BackprojectID m_idBackproject;
+ InterpolationID m_idInterpolation;
+ bool m_fail;
+ Backproject* m_pBackprojectImplem;
+
+ static const char BPROJ_TRIG_STR[]= "trig";
+ static const char BPROJ_TABLE_STR[]= "table";
+ static const char BPROJ_DIFF_STR[]= "diff";
+ static const char BPROJ_DIFF2_STR[]= "diff2";
+ static const char BPROJ_IDIFF2_STR[]= "idiff2";
+
+ static const char INTERP_NEAREST_STR[]= "nearest";
+ static const char INTERP_LINEAR_STR[]= "linear";
+ static const char INTERP_BSPLINE_STR[]= "bspline";
+
+ static const InterpolationID convertInterpolationNameToID (const char* const interpName);
+ static const char* convertInterpolationIDToName (const InterpolationID interpID);
+ static const BackprojectID convertBackprojectNameToID (const char* const bprojName);
+ static const char* convertBackprojectIDToName (const BackprojectID bprojID);
+
+ bool initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName);
+};
+
+
class Backproject
{
public:
- Backproject (const Projections& proj, ImageFile& im, InterpolationType interpType);
+ Backproject (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType);
virtual ~Backproject ();
- virtual void BackprojectView (const double* const t, const double view_angle) {};
+ virtual void BackprojectView (const double* const viewData, const double viewAngle) {};
protected:
void ScaleImageByRotIncrement (void);
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int ni);
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int ni);
- InterpolationType interpType;
+ Backprojector::InterpolationID interpType;
const Projections& proj;
ImageFile& im;
ImageFileArray v;
class BackprojectTrig : public Backproject
{
public:
- BackprojectTrig (const Projections& proj, ImageFile& im, InterpolationType interpType)
+ BackprojectTrig (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
: Backproject::Backproject (proj, im, interpType)
{}
class BackprojectTable : public Backproject
{
public:
- BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType);
+ BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType);
~BackprojectTable ();
void BackprojectView (const double* const t, double view_angle);
class BackprojectDiff : public Backproject
{
public:
- BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType);
+ BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType);
~BackprojectDiff ();
void BackprojectView (const double* const t, double view_angle);
class BackprojectDiff2 : public BackprojectDiff
{
public:
- BackprojectDiff2 (const Projections& proj, ImageFile& im, InterpolationType interpType)
+ BackprojectDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
: BackprojectDiff::BackprojectDiff (proj, im, interpType)
{}
class BackprojectIntDiff2 : public BackprojectDiff
{
public:
- BackprojectIntDiff2 (const Projections& proj, ImageFile& im, InterpolationType interpType)
+ BackprojectIntDiff2 (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
: BackprojectDiff::BackprojectDiff (proj, im, interpType)
{}
};
-Backproject* selectBackprojector (BackprojType type, const Projections& proj, ImageFile& im, InterpolationType interpType);
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: ct.h,v 1.20 2000/06/20 17:54:51 kevin Exp $
+** $Id: ct.h,v 1.21 2000/06/22 10:17:28 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 <stdint.h> /* Standard ints on Linux */
#endif
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <iterator>
+#include <algorithm>
+#include <exception>
+#include <stdexcept>
+#include <memory>
+
+using namespace std;
+
#ifdef HAVE_MPI
#include "mpi++.h"
#include "mpiworld.h"
#include "sgp.h"
#endif
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <string>
-#include <iterator>
-#include <algorithm>
-#include <exception>
-#include <stdexcept>
-
-using namespace std;
-
#include "array2d.h"
#include "fnetorderstream.h"
#include "imagefile.h"
#include "phantom.h"
-#include "projections.h"
#include "scanner.h"
+#include "backprojectors.h"
+#include "filter.h"
+#include "projections.h"
//----------------------------------------------------------------------//
// USER SYMBOLS //
//----------------------------------------------------------------------//
-// Filter types
-static const char O_FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit";
-static const char O_FILTER_ABS_SINC_STR[]= "abs_sinc";
-static const char O_FILTER_ABS_COS_STR[]= "abs_cos";
-static const char O_FILTER_ABS_HAMMING_STR[]= "abs_hamming";
-static const char O_FILTER_SHEPP_STR[]= "shepp";
-static const char O_FILTER_BANDLIMIT_STR[]= "bandlimit";
-static const char O_FILTER_SINC_STR[]= "sinc";
-static const char O_FILTER_COS_STR[]= "cos";
-static const char O_FILTER_HAMMING_STR[]= "hamming";
-static const char O_FILTER_TRIANGLE_STR[]= "triangle";
-
-typedef enum {
- FILTER_BANDLIMIT,
- FILTER_SINC,
- FILTER_G_HAMMING,
- FILTER_COSINE,
- FILTER_TRIANGLE,
- FILTER_ABS_BANDLIMIT, // filter times |x|
- FILTER_ABS_SINC,
- FILTER_ABS_G_HAMMING,
- FILTER_ABS_COSINE,
- FILTER_SHEPP
-} FilterType;
-
-
-// Function domains
-static const char D_FREQ_STR[]= "freq";
-static const char D_SPATIAL_STR[]= "spatial";
-
-typedef enum {
- D_FREQ = 1,
- D_SPATIAL
-} DomainType;
-
-
-/* interpolation methods */
-static const char O_INTERP_NEAREST_STR[]= "nearest";
-static const char O_INTERP_LINEAR_STR[]= "linear";
-static const char O_INTERP_BSPLINE_STR[]= "bspline";
-
-#undef HAVE_BSPLINE_INTERP
-typedef enum { // Interpolation methods
- I_NEAREST = 1, // Nearest neighbor
-#if HAVE_BSPLINE_INTERP
- I_BSPLINE,
- I_1BSPLINE, // 1st order B-Spline
- I_2BSPLINE,
- I_3BSPLINE,
-#endif
- I_LINEAR // Linear interpolation
-} InterpolationType;
-
-
// Trace levels
static const char O_TRACE_NONE_STR[]= "none";
static const char O_TRACE_TEXT_STR[]= "text";
TRACE_CLIPPING /* Plot clipping */
};
-// Standard phantomsa
-static const char O_PHM_HERMAN_STR[]= "herman";
-static const char O_PHM_ROWLAND_STR[]= "rowland";
-static const char O_PHM_BROWLAND_STR[]= "browland";
-static const char O_PHM_UNITPULSE_STR[]= "unitpulse";
-typedef enum {
- O_PHM_HERMAN, /* Herman head phantom */
- O_PHM_ROWLAND, /* Rowland head phantom */
- O_PHM_BROWLAND, /* Bordered Rowland head phantom */
- O_PHM_UNITPULSE /* Unit pulse phantom */
-} PhantomType;
-
-// Backproject types
-static const char O_BPROJ_TRIG_STR[]= "trig";
-static const char O_BPROJ_TABLE_STR[]= "table";
-static const char O_BPROJ_DIFF_STR[]= "diff";
-static const char O_BPROJ_DIFF2_STR[]= "diff2";
-static const char O_BPROJ_IDIFF2_STR[]= "idiff2";
-
-typedef enum {
- O_BPROJ_TRIG,
- O_BPROJ_TABLE,
- O_BPROJ_DIFF,
- O_BPROJ_DIFF2,
- O_BPROJ_IDIFF2
-} BackprojType;
-
-// Convolution symmetries
-typedef enum {
- FUNC_EVEN = 1, // function types, f[-n] = f[n]
- FUNC_ODD, // f[-n] = -f[n]
- FUNC_BOTH // function has both odd & even components
-} FunctionSymmetry;
-
-
-
-#include "backprojectors.h"
-#include "filter.h"
-
/*************************************************************************
* FUNCTION DECLARATIONS
// options.cpp
int opt_set_trace(const char *optarg);
-const char *name_of_phantom(const int phmid);
-int opt_set_phantom(const char *optarg);
-InterpolationType opt_set_interpolation(const char *optarg);
-const char *name_of_interpolation(int interp_type);
-FilterType opt_set_filter(const char *optarg);
-const char *name_of_filter(const int filter);
-DomainType opt_set_filter_domain(const char *optarg);
-const char *name_of_filter_domain(const DomainType domain);
-BackprojType opt_set_backproj(const char *optarg);
-const char *name_of_backproj(const BackprojType backproj);
-
-// 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);
-// image.cpp
-void image_filter_response(ImageFile& im, const DomainType domain, double bw, const FilterType filt_type, double filt_param, const int opt_trace);
+// imagefile.cpp
+void image_filter_response (ImageFile& im, const char* const domainName, double bw, const char* const filterName, double filt_param, const int opt_trace);
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, 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: ctsupport.h,v 1.3 2000/06/19 19:10:08 kevin Exp $
+** $Id: ctsupport.h,v 1.4 2000/06/22 10:17:28 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
/* filefuncs.cpp */
-bool file_exists(const char* fname);
-const char* fileBasename(const char* filename);
-FILE *sys_fopen(const char *filename, const char *mode, const char *progname);
+bool fileExists (const char* fname);
+const char* fileBasename (const char* filename);
/* strfuncs.cpp */
char* str_skip_head(const char* str, const char* const charlist);
int clip_circle (double& x1, double& y1, double& x2, double& y2, const double cx, const double cy, const double radius, double t1, double t2);
int clip_triangle (double& x1, double& y1, double& x2, double& y2, const double u, const double v, const int clip_xaxis);
-// norm_ang.cpp
-double norm_ang (double theta);
// xform.cpp
void indent_mtx2 (GRFMTX_2D m);
void xlat2d (double x[], double y[], int pts, double xoffset, double yoffset);
void scale2d (double x[], double y[], int pts, double xfact, double yfact);
-// simpson.cpp
+// mathfuncs.cpp
+double normalizeAngle (double theta);
double integrateSimpson (const double xmin, const double xmax, const double *y, const int np);
+#ifndef FILTER_H
+#define FILTER_H
+
+
class SignalFilter {
public:
- SignalFilter (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, const int numInt);
+ typedef enum {
+ FILTER_INVALID,
+ FILTER_BANDLIMIT,
+ FILTER_SINC,
+ FILTER_G_HAMMING,
+ FILTER_COSINE,
+ FILTER_TRIANGLE,
+ FILTER_ABS_BANDLIMIT, // filter times |x|
+ FILTER_ABS_SINC,
+ FILTER_ABS_G_HAMMING,
+ FILTER_ABS_COSINE,
+ FILTER_SHEPP
+ } FilterID;
+
+ typedef enum {
+ DOMAIN_INVALID,
+ DOMAIN_FREQ,
+ DOMAIN_SPATIAL
+ } DomainID;
+
+
+ SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, const int numInt);
+
+ SignalFilter (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numInt);
~SignalFilter (void);
double* getFilter (void) const
{ return m_vecFilter; }
- double convolve (const double f[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const;
+ double convolve (const double f[], const double dx, const int n, const int np) const;
+
+ double convolve (const float f[], const double dx, const int n, const int np) const;
- double convolve (const float f[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const;
+ bool fail(void) const {return m_fail;}
+ const string& nameFilter(void) const { return m_nameFilter;}
+ const string& nameDomain(void) const { return m_nameDomain;}
+ const FilterID idFilter(void) const { return m_idFilter;}
+ const DomainID idDomain(void) const { return m_idDomain;}
- static double spatialResponseCalc (FilterType fType, double bw, double x, double param, int n);
+ static double response (const char* const FilterName, const char* const DomainName, double bw, double x, double param);
- static double spatialResponseAnalytic (FilterType fType, double bw, double x, double param);
+ static double spatialResponseCalc (FilterID fType, double bw, double x, double param, int n);
- static double frequencyResponse (FilterType fType, double bw, double u, double param);
+ static double spatialResponseAnalytic (FilterID fType, double bw, double x, double param);
+
+ static double frequencyResponse (FilterID fType, double bw, double u, double param);
private:
- double *m_vecFilter;
double m_bw;
- FilterType m_filterType;
-
- double spatialResponseCalc (double x, double param, int n) const;
+ int m_nPoints;
+ double m_xmin;
+ double m_xmax;
+ double* m_vecFilter;
+ bool m_fail;
+ string m_nameFilter;
+ string m_nameDomain;
+ FilterID m_idFilter;
+ DomainID m_idDomain;
+
+ static const char FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit";
+ static const char FILTER_ABS_SINC_STR[]= "abs_sinc";
+ static const char FILTER_ABS_COS_STR[]= "abs_cos";
+ static const char FILTER_ABS_HAMMING_STR[]= "abs_hamming";
+ static const char FILTER_SHEPP_STR[]= "shepp";
+ static const char FILTER_BANDLIMIT_STR[]= "bandlimit";
+ static const char FILTER_SINC_STR[]= "sinc";
+ static const char FILTER_COS_STR[]= "cos";
+ static const char FILTER_HAMMING_STR[]= "hamming";
+ static const char FILTER_TRIANGLE_STR[]= "triangle";
+
+ static const char DOMAIN_FREQ_STR[]= "freq";
+ static const char DOMAIN_SPATIAL_STR[]= "spatial";
+
+ static FilterID convertFilterNameToID (const char* filterName);
+ static const char* convertFilterIDToName (const FilterID filterID);
+ static const DomainID convertDomainNameToID (const char* domainName);
+ static const char* convertDomainIDToName (const DomainID domainID);
+
+ void init (const FilterID filt_type, double bw, double xmin, double xmax, int n, double param, const DomainID domain, const int numInt);
+
+double spatialResponseCalc (double x, double param, int n) const;
double spatialResponseAnalytic (double x, double param) const;
};
+
+#endif
#include <mpi++.h>
#include <vector.h>
-
+#include <string>
class MPIWorld
{
MPI::Intracomm& getComm()
{ return m_comm; }
-
+
+ void BcastString (string& str);
+
private:
int m_myRank;
int m_nProcessors;
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: phantom.h,v 1.2 2000/06/19 20:08:09 kevin Exp $
+** $Id: phantom.h,v 1.3 2000/06/22 10:17:28 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
double m_rectLimits[4];
static const int POINTS_PER_CIRCLE = 360;
+ static const double SCALE_PELEM_EXTENT=0.005; // increase pelem limits by 0.5%
+
+ static PhmElemType PhantomElement::convertNameToType (const char* const typeName);
void makeTransformMatrices (void);
void calcEllipsePoints (double x[], double y[], const int pts, const double u, const double v);
- int numCirclePoints (double theta) const;
+ static int numCirclePoints (double theta);
};
class Phantom
{
public:
+ typedef enum {
+ PHM_INVALID,
+ PHM_HERMAN, /* Herman head phantom */
+ PHM_BHERMAN, /* Bordered Herman head phantom */
+ PHM_ROWLAND, /* Rowland head phantom */
+ PHM_BROWLAND, /* Bordered Rowland head phantom */
+ PHM_UNITPULSE /* Unit pulse phantom */
+ } PhantomID;
+
+
Phantom (void);
+ Phantom (const char* const phmName);
+
~Phantom (void);
void setComposition (PhantomComposition composition)
const PhantomComposition getComposition (void) const
{ return m_composition; }
- void create (const int phmid);
+ bool createFromPhantom (const char* const phmName);
+
+ bool createFromPhantom (const PhantomID phmid);
bool createFromFile (const char* const fname);
void convertToImagefile (ImageFile& im, const int in_nsample, const int trace) const;
+ bool fail(void) const
+ {return m_fail;}
+
+ const string& name(void) const
+ {return m_name;}
+
+ const PhantomID id(void) const
+ {return m_id;}
+
#if HAVE_SGP
void show (void) const;
void draw (void) const;
#endif
void addStdHerman (void);
+ void addStdHermanBordered (void);
void addStdRowland (void);
void addStdRowlandBordered (void);
double m_diameter; // diameter of object
mutable slist<PhantomElement*> m_listPElem; // pelem lists
string m_name;
+ PhantomID m_id;
+ bool m_fail;
+
+ // Standard phantomsa
+ static const char PHM_HERMAN_STR[]= "herman";
+ static const char PHM_BHERMAN_STR[]= "bherman";
+ static const char PHM_ROWLAND_STR[]= "rowland";
+ static const char PHM_BROWLAND_STR[]= "browland";
+ static const char PHM_UNITPULSE_STR[]= "unitpulse";
+ static PhantomID convertNameToPhantomID (const char* const phmName);
+ static const char* convertPhantomIDToName (const PhantomID phmID);
+
+ void init(void);
};
typedef slist<PhantomElement*>::iterator PElemIterator;
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: projections.h,v 1.3 2000/06/20 17:54:51 kevin Exp $
+** $Id: projections.h,v 1.4 2000/06/22 10:17:28 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
bool detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int view_num);
bool detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int view_num);
+ bool reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const interpName, int interp_param, const char* const backprojName, const int trace);
+
void setNView (int nView); // used in MPI to restrict # of views
void setRotInc (double rotInc) { m_rotInc = rotInc;}
void setDetInc (double detInc) { m_detInc = detInc;}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: scanner.h,v 1.3 2000/06/19 17:58:20 kevin Exp $
+** $Id: scanner.h,v 1.4 2000/06/22 10:17:28 kevin Exp $
**
**
** This program is free software; you can redistribute it and/or modify
#ifndef SCANNER_H
#define SCANNER_H
-#include "projections.h"
+class Projections;
+
// Projections are collected along an array of ndet detectors. The data
// for these detectors is stored in the class DetectorArray
--- /dev/null
+#ifndef TRACE_H
+#define TRACE_H
+
+// Trace levels
+static const char TRACE_NONE_STR[]= "none";
+static const char TRACE_TEXT_STR[]= "text";
+static const char TRACE_PHM_STR[]= "phm";
+static const char TRACE_RAYS_STR[]= "rays";
+static const char TRACE_PLOT_STR[]= "plot";
+static const char TRACE_CLIPPING_STR[]= "clipping";
+
+enum {
+ TRACE_NONE, /* No tracing */
+ TRACE_TEXT, /* Minimal status */
+ TRACE_PHM, /* Show phantom */
+ TRACE_RAYS, /* Show all rays */
+ TRACE_PLOT, /* Plot raysums */
+ TRACE_CLIPPING /* Plot clipping */
+};
+
+inline int
+opt_set_trace (const char *optarg)
+{
+ int opt;
+
+ if (strcmp(optarg, TRACE_NONE_STR) == 0)
+ opt = TRACE_NONE;
+ else if (strcmp(optarg, TRACE_TEXT_STR) == 0)
+ opt = TRACE_TEXT;
+ else if (strcmp(optarg, TRACE_PHM_STR) == 0)
+ opt = TRACE_PHM;
+ else if (strcmp(optarg, TRACE_PLOT_STR) == 0)
+ opt = TRACE_PLOT;
+ else if (strcmp(optarg, TRACE_CLIPPING_STR) == 0)
+ opt = TRACE_CLIPPING;
+ else if (strcmp(optarg, TRACE_RAYS_STR) == 0)
+ opt = TRACE_RAYS;
+ else {
+ sys_error(ERR_WARNING,"Invalid trace option %s\n", optarg);
+ opt = -1;
+ }
+
+ return (opt);
+}
+
+
+#endif
+
noinst_LIBRARIES = libctsim.a
-libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp options.cpp imagefile.cpp reconstr.cpp backprojectors.cpp
+libctsim_a_SOURCES = filter.cpp scanner.cpp projections.cpp phantom.cpp imagefile.cpp backprojectors.cpp
INCLUDES=@my_includes@
EXTRA_DIST=Makefile.nt
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.2 2000/06/19 19:07:33 kevin Exp $
+** $Id: backprojectors.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
#include "ct.h"
+Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+{
+ m_fail = false;
+ m_pBackprojectImplem = NULL;
+
+ initBackprojector (proj, im, backprojName, interpName);
+}
+
+void
+Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
+{
+ if (m_pBackprojectImplem)
+ m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
+}
+
+Backprojector::~Backprojector (void)
+{
+ delete m_pBackprojectImplem;
+}
// FUNCTION IDENTIFICATION
// Backproject* projector = selectBackprojector (...)
// Selects a backprojector based on BackprojType
// and initializes the backprojector
-Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
+bool
+Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
{
- Backproject* bj = NULL;
-
- if (bjType == O_BPROJ_TRIG)
- bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
- else if (bjType == O_BPROJ_TABLE)
- bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
- else if (bjType == O_BPROJ_DIFF)
- bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
- else if (bjType == O_BPROJ_DIFF2)
- bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
- else if (bjType == O_BPROJ_IDIFF2)
- bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
- else
- sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
-
- return (bj);
+ m_nameBackproject = backprojName;
+ m_nameInterpolation = interpName;
+ m_idBackproject = convertBackprojectNameToID (backprojName);
+ m_idInterpolation = convertInterpolationNameToID (interpName);
+ m_pBackprojectImplem = NULL;
+
+ if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
+ m_fail = true;
+ return false;
+ }
+
+ if (m_idBackproject == BPROJ_TRIG)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
+ else if (m_idBackproject == BPROJ_TABLE)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
+ else if (m_idBackproject == BPROJ_DIFF)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
+ else if (m_idBackproject == BPROJ_DIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
+ else if (m_idBackproject == BPROJ_IDIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+ else {
+ m_fail = true;
+ return false;
+ }
+
+ return true;
+}
+
+
+const Backprojector::BackprojectID
+Backprojector::convertBackprojectNameToID (const char* const backprojName)
+{
+ BackprojectID backprojID = BPROJ_INVALID;
+
+ if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
+ backprojID = BPROJ_TRIG;
+ else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
+ backprojID = BPROJ_TABLE;
+ else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
+ backprojID = BPROJ_DIFF;
+ else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
+ backprojID = BPROJ_DIFF2;
+ else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
+ backprojID = BPROJ_IDIFF2;
+
+ return (backprojID);
+}
+
+const char*
+Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
+{
+ const char *bprojName = "";
+
+ if (bprojID == BPROJ_TRIG)
+ bprojName = BPROJ_TRIG_STR;
+ else if (bprojID == BPROJ_TABLE)
+ bprojName = BPROJ_TABLE_STR;
+ else if (bprojID == BPROJ_DIFF)
+ bprojName = BPROJ_DIFF_STR;
+ else if (bprojID == BPROJ_DIFF2)
+ bprojName = BPROJ_DIFF2_STR;
+ else if (bprojID == BPROJ_IDIFF2)
+ bprojName = BPROJ_IDIFF2_STR;
+
+ return (bprojName);
+}
+
+
+
+const Backprojector::InterpolationID
+Backprojector::convertInterpolationNameToID (const char* const interpName)
+{
+ InterpolationID interpID = INTERP_INVALID;
+
+ if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
+ interpID = INTERP_NEAREST;
+ else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
+ interpID = INTERP_LINEAR;
+#if HAVE_BSPLINE_INTERP
+ else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
+ interpID = INTERP_BSPLINE;
+#endif
+
+ return (interpID);
+}
+
+
+/* NAME
+ * name_of_interp Return name of interpolation method
+ *
+ * SYNOPSIS
+ * name = name_of_interp (interp_type)
+ * char *name Name of interpolation method
+ * int interp_type Method of interpolation
+ *
+ * NOTES
+ * Returns NULL if interp_type is invalid
+ */
+
+const char*
+Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
+{
+ if (interpID == INTERP_NEAREST)
+ return (INTERP_NEAREST_STR);
+ else if (interpID == INTERP_LINEAR)
+ return (INTERP_LINEAR_STR);
+#if HAVE_BSPLINE_INTERP
+ else if (interpID == INTERP_BSPLINE)
+ return (INTERP_BSPLINE_STR);
+#endif
+ else
+ return ("");
}
// PURPOSE
// Pure virtual base class for all backprojectors.
-Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType)
+Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType)
: proj(proj), im(im), interpType(interpType)
{
detInc = proj.detInc();
xInc = (xMax - xMin) / nx; // size of cells
yInc = (yMax - yMin) / ny;
- if (interpType != I_NEAREST && interpType != I_LINEAR)
+ if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR)
sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
}
double phi = atan2 (y, x); // angle of cell from center
double L = r * cos (theta - phi); // position on detector
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
else
v[ix][iy] += rotInc * filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double p = L / detInc; // position along detector
double pFloor = floor (p);
int iDetPos = iDetCenter + static_cast<int>(pFloor);
// PURPOSE
// Precalculates trigometric function value for each point in image for backprojection.
-BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType)
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
: Backproject::Backproject (proj, im, interpType)
{
arrayR.initSetSize (nx, ny);
for (int iy = 0; iy < ny; iy++) {
double L = r[ix][iy] * cos (theta - phi[ix][iy]);
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
else
pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double dPos = L / detInc; // position along detector
double dPosFloor = floor (dPos);
int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
// Backprojects by precalculating the change in L position for each x & y step in the image.
// Iterates in x & y direction by adding difference in L position
-BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType)
+BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
: Backproject::Backproject (proj, im, interpType)
{
// calculate center of first pixel v[0][0]
#ifdef DEBUG
printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
#endif
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double detPos = curDetPos / detInc; // position along detector
double detPosFloor = floor (detPos);
int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
#ifdef DEBUG
printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
#endif
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double detPosFloor = floor (curDetPos);
int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
double frac = curDetPos - detPosFloor; // fraction distance from det
ImageFileColumn pImCol = v[ix];
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
kint32 detPosFloor = curDetPos / scale;
kint32 detPosRemainder = curDetPos % scale;
if (detPosRemainder < 0) {
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: filter.cpp,v 1.2 2000/06/20 17:54:51 kevin Exp $
+** $Id: filter.cpp,v 1.3 2000/06/22 10:17:28 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
/* NAME
- * filter_generate Generate a filter
+ * SignalFilter::SignalFilter Construct a signal
*
* SYNOPSIS
- * f = filter_generate (filt_type, bw, xmin, xmax, n, param, domain, analytic)
- * double f Generated filter vector
- * int filt_type Type of filter wanted
- * double bw Bandwidth of filter
- * double xmin, xmax Filter limits
- * int n Number of points in filter
- * double param General input parameter to filters
- * int domain FREQ or SPATIAL domain wanted
- * int numint Number if intervals for calculating
- * discrete inverse fourier xform
- * for spatial domain filters. For
- * ANALYTIC solutions, use numint = 0
+ * f = SignalFilter (filt_type, bw, xmin, xmax, n, param, domain, analytic)
+ * double f Generated filter vector
+ * int filt_type Type of filter wanted
+ * double bw Bandwidth of filter
+ * double xmin, xmax Filter limits
+ * int n Number of points in filter
+ * double param General input parameter to filters
+ * int domain FREQ or SPATIAL domain wanted
+ * int numint Number if intervals for calculating discrete inverse fourier xform
+ * for spatial domain filters. For ANALYTIC solutions, use numint = 0
*/
-SignalFilter::SignalFilter (const FilterType filt_type, double bw, double xmin, double xmax, int n, double param, const DomainType domain, int numint)
+SignalFilter::SignalFilter (const char* filterName, double bw, double xmin, double xmax, int n, double param, const char* domainName, int numint)
+{
+ m_idFilter = convertFilterNameToID (filterName);
+ m_idDomain = convertDomainNameToID (domainName);
+ init (m_idFilter, bw, xmin, xmax, n, param, m_idDomain, numint);
+}
+
+SignalFilter::SignalFilter (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
+{
+ init (filterID, bw, xmin, xmax, n, param, domainID, numint);
+}
+
+void
+SignalFilter::init (const FilterID filterID, double bw, double xmin, double xmax, int n, double param, const DomainID domainID, int numint)
{
- m_vecFilter = new double [n];
- m_filterType = filt_type;
m_bw = bw;
+ m_idFilter = filterID;
+ m_idDomain = domainID;
+ if (m_idFilter == FILTER_INVALID || m_idDomain == DOMAIN_INVALID) {
+ m_fail = true;
+ return;
+ }
+ m_nameFilter = convertFilterIDToName (m_idFilter);
+ m_nameDomain = convertDomainIDToName (m_idDomain);
+ m_fail = false;
+ m_nPoints = n;
+ m_xmin = xmin;
+ m_xmax = xmax;
+ m_vecFilter = new double[n];
- double xinc = (xmax - xmin) / (n - 1);
+ double xinc = (m_xmax - m_xmin) / (m_nPoints - 1);
- if (m_filterType == FILTER_SHEPP) {
+ if (m_idFilter == FILTER_SHEPP) {
double a = 2 * m_bw;
double c = - 4. / (a * a);
- int center = (n - 1) / 2;
+ int center = (m_nPoints - 1) / 2;
int sidelen = center;
m_vecFilter[center] = 4. / (a * a);
-
+
for (int i = 1; i <= sidelen; i++ )
m_vecFilter [center + i] = m_vecFilter [center - i] = c / (4 * (i * i) - 1);
- } else if (domain == D_FREQ) {
+ } else if (m_idDomain == DOMAIN_FREQ) {
double x;
int i;
- for (x = xmin, i = 0; i < n; x += xinc, i++)
+ for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
m_vecFilter[i] = frequencyResponse (x, param);
- } else if (domain == D_SPATIAL) {
+ } else if (m_idDomain == DOMAIN_SPATIAL) {
double x;
int i;
- for (x = xmin, i = 0; i < n; x += xinc, i++)
+ for (x = m_xmin, i = 0; i < m_nPoints; x += xinc, i++)
if (numint == 0)
m_vecFilter[i] = spatialResponseAnalytic (x, param);
else
m_vecFilter[i] = spatialResponseCalc (x, param, numint);
} else {
- sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", domain);
+ sys_error (ERR_WARNING, "Illegal domain %d [filt_generate]", m_idDomain);
+ m_fail = true;
}
}
SignalFilter::~SignalFilter (void)
{
- delete m_vecFilter;
+ delete m_vecFilter;
+}
+
+
+SignalFilter::FilterID
+SignalFilter::convertFilterNameToID (const char *filterName)
+{
+ FilterID filterID;
+
+ if (strcasecmp (filterName, FILTER_BANDLIMIT_STR) == 0)
+ filterID = FILTER_BANDLIMIT;
+ else if (strcasecmp (filterName, FILTER_HAMMING_STR) == 0)
+ filterID = FILTER_G_HAMMING;
+ else if (strcasecmp (filterName, FILTER_SINC_STR) == 0)
+ filterID = FILTER_SINC;
+ else if (strcasecmp (filterName, FILTER_COS_STR) == 0)
+ filterID = FILTER_COSINE;
+ else if (strcasecmp (filterName, FILTER_TRIANGLE_STR) == 0)
+ filterID = FILTER_TRIANGLE;
+ else if (strcasecmp (filterName, FILTER_ABS_BANDLIMIT_STR) == 0)
+ filterID = FILTER_ABS_BANDLIMIT;
+ else if (strcasecmp (filterName, FILTER_ABS_HAMMING_STR) == 0)
+ filterID = FILTER_ABS_G_HAMMING;
+ else if (strcasecmp (filterName, FILTER_ABS_SINC_STR) == 0)
+ filterID = FILTER_ABS_SINC;
+ else if (strcasecmp (filterName, FILTER_ABS_COS_STR) == 0)
+ filterID = FILTER_ABS_COSINE;
+ else if (strcasecmp (filterName, FILTER_SHEPP_STR) == 0)
+ filterID = FILTER_SHEPP;
+ else {
+ sys_error(ERR_WARNING, "Invalid filter type %s\n", filterName);
+ filterID = FILTER_INVALID;
+ }
+
+ return (filterID);
+}
+
+const char *
+SignalFilter::convertFilterIDToName (const FilterID filterID)
+{
+ const char *name = "";
+
+ if (filterID == FILTER_SHEPP)
+ name = FILTER_SHEPP_STR;
+ else if (filterID == FILTER_ABS_COSINE)
+ name = FILTER_ABS_COS_STR;
+ else if (filterID == FILTER_ABS_SINC)
+ name = FILTER_ABS_SINC_STR;
+ else if (filterID == FILTER_ABS_G_HAMMING)
+ name = FILTER_ABS_HAMMING_STR;
+ else if (filterID == FILTER_ABS_BANDLIMIT)
+ name = FILTER_ABS_BANDLIMIT_STR;
+ else if (filterID == FILTER_COSINE)
+ name = FILTER_COS_STR;
+ else if (filterID == FILTER_SINC)
+ name = FILTER_SINC_STR;
+ else if (filterID == FILTER_G_HAMMING)
+ name = FILTER_HAMMING_STR;
+ else if (filterID == FILTER_BANDLIMIT)
+ name = FILTER_BANDLIMIT_STR;
+ else if (filterID == FILTER_TRIANGLE)
+ name = FILTER_TRIANGLE_STR;
+
+ return (name);
+}
+
+const SignalFilter::DomainID
+SignalFilter::convertDomainNameToID (const char* const domainName)
+{
+ DomainID dID;
+
+ if (strcasecmp (domainName, DOMAIN_SPATIAL_STR) == 0)
+ dID = DOMAIN_SPATIAL;
+ else if (strcasecmp (domainName, DOMAIN_FREQ_STR) == 0)
+ dID = DOMAIN_FREQ;
+ else
+ dID = DOMAIN_INVALID;
+
+ return (dID);
}
+const char *
+SignalFilter::convertDomainIDToName (const DomainID domain)
+{
+ const char *name = "";
+
+ if (domain == DOMAIN_SPATIAL)
+ return (DOMAIN_SPATIAL_STR);
+ else if (domain == DOMAIN_FREQ)
+ return (DOMAIN_FREQ_STR);
+
+ return (name);
+}
+
+
+double
+SignalFilter::response (const char* filterName, const char* domainName, double bw, double x, double filt_param)
+{
+ double response = 0;
+ FilterID filterID = convertFilterNameToID (filterName);
+ DomainID domainID = convertDomainNameToID (domainName);
+
+ if (domainID == DOMAIN_SPATIAL)
+ response = spatialResponseAnalytic (filterID, bw, x, filt_param);
+ else if (domainID == DOMAIN_FREQ)
+ response = frequencyResponse (filterID, bw, x, filt_param);
+
+ return (response);
+}
/* NAME
* filter_spatial_response_calc Calculate filter by discrete inverse fourier
double
SignalFilter::spatialResponseCalc (double x, double param, int n) const
{
- return (spatialResponseCalc (m_filterType, m_bw, x, param, n));
+ return (spatialResponseCalc (m_idFilter, m_bw, x, param, n));
}
double
-SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double param, int n)
+SignalFilter::spatialResponseCalc (FilterID filterID, double bw, double x, double param, int n)
{
double zmin, zmax;
- if (fType == FILTER_TRIANGLE) {
+ if (filterID == FILTER_TRIANGLE) {
zmin = 0;
zmax = bw;
} else {
double z = zmin;
double q [n];
for (int i = 0; i < n; i++, z += zinc)
- q[i] = frequencyResponse (fType, bw, z, param) * cos (TWOPI * z * x);
+ q[i] = frequencyResponse (filterID, bw, z, param) * cos (TWOPI * z * x);
double y = 2 * integrateSimpson (zmin, zmax, q, n);
double
SignalFilter::frequencyResponse (double u, double param) const
{
- return frequencyResponse (m_filterType, m_bw, u, param);
+ return frequencyResponse (m_idFilter, m_bw, u, param);
}
double
-SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double param)
+SignalFilter::frequencyResponse (FilterID filterID, double bw, double u, double param)
{
double q;
double au = fabs (u);
- switch (fType) {
+ switch (filterID) {
case FILTER_BANDLIMIT:
if (au >= bw / 2)
q = 0.;
break;
default:
q = 0;
- sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", fType);
+ sys_error (ERR_WARNING, "Frequency response for filter %d not implemented [filter_frequency_response]", filterID);
break;
}
return (q);
double
SignalFilter::spatialResponseAnalytic (double x, double param) const
{
- return spatialResponseAnalytic (m_filterType, m_bw, x, param);
+ return spatialResponseAnalytic (m_idFilter, m_bw, x, param);
}
double
-SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, double param)
+SignalFilter::spatialResponseAnalytic (FilterID filterID, double bw, double x, double param)
{
double q, temp;
double u = TWOPI * x;
double b = PI / bw;
double b2 = TWOPI / bw;
- switch (fType) {
+ switch (filterID) {
case FILTER_BANDLIMIT:
q = bw * sinc(u * w, 1.0);
break;
break;
case FILTER_ABS_SINC:
default:
- sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", fType);
+ sys_error (ERR_WARNING, "Analytic filter type %d not implemented [filter_spatial_response_analytic]", filterID);
q = 0;
break;
}
*/
double
-SignalFilter::convolve (const double func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+SignalFilter::convolve (const double func[], const double dx, const int n, const int np) const
{
double sum = 0.0;
- if (func_type == FUNC_BOTH) {
#if UNOPTIMIZED_CONVOLUTION
- for (int i = 0; i < np; i++)
- sum += func[i] * m_vecFilter[n - i + (np - 1)];
+ for (int i = 0; i < np; i++)
+ sum += func[i] * m_vecFilter[n - i + (np - 1)];
#else
- double* f2 = m_vecFilter + n + (np - 1);
- for (int i = 0; i < np; i++)
- sum += *func++ * *f2--;
+ double* f2 = m_vecFilter + n + (np - 1);
+ for (int i = 0; i < np; i++)
+ sum += *func++ * *f2--;
#endif
- } else if (func_type == FUNC_EVEN) {
- for (int i = 0; i < np; i++) {
- int k = abs (n - i);
- sum += func[i] * m_vecFilter[k];
- }
- } else if (func_type == FUNC_ODD) {
- for (int i = 0; i < np; i++) {
- int k = n - i;
- if (k < 0)
- sum -= func[i] * m_vecFilter[k];
- else
- sum += func[i] * m_vecFilter[k];
- }
- } else
- sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
return (sum * dx);
}
double
-SignalFilter::convolve (const float func[], const double dx, const int n, const int np, const FunctionSymmetry func_type) const
+SignalFilter::convolve (const float func[], const double dx, const int n, const int np) const
{
double sum = 0.0;
- if (func_type == FUNC_BOTH) {
#if UNOPTIMIZED_CONVOLUTION
- for (int i = 0; i < np; i++)
- sum += func[i] * m_vecFilter[n - i + (np - 1)];
+for (int i = 0; i < np; i++)
+ sum += func[i] * m_vecFilter[n - i + (np - 1)];
#else
- double* f2 = m_vecFilter + n + (np - 1);
- for (int i = 0; i < np; i++)
- sum += *func++ * *f2--;
+double* f2 = m_vecFilter + n + (np - 1);
+for (int i = 0; i < np; i++)
+ sum += *func++ * *f2--;
#endif
- } else if (func_type == FUNC_EVEN) {
- for (int i = 0; i < np; i++) {
- int k = abs (n - i);
- sum += func[i] * m_vecFilter[k];
- }
- } else if (func_type == FUNC_ODD) {
- for (int i = 0; i < np; i++) {
- int k = n - i;
- if (k < 0)
- sum -= func[i] * m_vecFilter[k];
- else
- sum += func[i] * m_vecFilter[k];
- }
- } else
- sys_error (ERR_WARNING, "Illegal function type %d [convolve]", func_type);
return (sum * dx);
}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: imagefile.cpp,v 1.2 2000/06/20 17:54:51 kevin Exp $
+** $Id: imagefile.cpp,v 1.3 2000/06/22 10:17:28 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
/* FILE
* image.c Routines for managing images
- *
- * PROGRAMMER: Kevin Rosenberg
- * DATE: Aug 1984
- *
- * FUNCTION
- * Provides a set of routines for handling image files
*/
void
-image_filter_response (ImageFile& im, const DomainType domain, double bw, const FilterType filt_type, double filt_param, const int opt_trace)
+image_filter_response (ImageFile& im, const char* const domainName, double bw, const char* const filterName, double filt_param, const int opt_trace)
{
int hx = im.nx() / 2;
int hy = im.ny() / 2;
for (int j = -hy; j <= hy; j++) {
double r = sqrt(i * i + j * j);
- if (domain == D_SPATIAL)
- v[i+hx][j+hy] = SignalFilter::spatialResponseAnalytic (filt_type, bw, r, filt_param);
- else
- v[i+hx][j+hy] = SignalFilter::frequencyResponse (filt_type, bw, r, filt_param);
+ v[i+hx][j+hy] = SignalFilter::response (filterName, domainName, bw, r, filt_param);
if (opt_trace >= TRACE_PHM)
printf ("r=%8.4f, v=%8.4f\n", r, v[i+hx][j+hy]);
}
ImageFileArray v = im.getArray();
#if HAVE_G2_H
- int* pens = new int [nx * ny * scale * scale ];
+ int pens [nx * ny * scale * scale ];
double view_scale = 255 / (pmax - pmin);
int id_X11 = g2_open_X11 (nx * scale, ny * scale);
for (int i= 0, iy = ny - 1; iy >= 0; iy--) {
for (int ix = 0; ix < nx; ix++) {
- int cval = (int) ((v[ix][iy] - pmin) * view_scale);
+ int cval = static_cast<int>((v[ix][iy] - pmin) * view_scale);
if (cval < 0)
cval = 0;
else if (cval > 255)
g2_image (id_X11, 0., 0., nx, ny, pens);
- delete [] pens;
return (id_X11);
#endif
+++ /dev/null
-/*****************************************************************************
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: options.cpp,v 1.1 2000/06/19 02:59:34 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"
-
-int
-opt_set_trace (const char *optarg)
-{
- int opt;
-
- if (strcmp(optarg, O_TRACE_NONE_STR) == 0)
- opt = TRACE_NONE;
- else if (strcmp(optarg, O_TRACE_TEXT_STR) == 0)
- opt = TRACE_TEXT;
- else if (strcmp(optarg, O_TRACE_PHM_STR) == 0)
- opt = TRACE_PHM;
- else if (strcmp(optarg, O_TRACE_PLOT_STR) == 0)
- opt = TRACE_PLOT;
- else if (strcmp(optarg, O_TRACE_CLIPPING_STR) == 0)
- opt = TRACE_CLIPPING;
- else if (strcmp(optarg, O_TRACE_RAYS_STR) == 0)
- opt = TRACE_RAYS;
- else {
- sys_error(ERR_WARNING,"Invalid trace option %s\n", optarg);
- opt = -1;
- }
-
- return (opt);
-}
-
-const char *
-name_of_phantom (const int phmnum)
-{
- const char *name = "Unknown phantom";
-
- if (phmnum == O_PHM_HERMAN)
- name = O_PHM_HERMAN_STR;
- else if (phmnum == O_PHM_ROWLAND)
- name = O_PHM_ROWLAND_STR;
- else if (phmnum == O_PHM_BROWLAND)
- name = O_PHM_BROWLAND_STR;
- else if (phmnum == O_PHM_UNITPULSE)
- name = O_PHM_UNITPULSE_STR;
-
- return (name);
-}
-
-int
-opt_set_phantom (const char *optarg)
-{
- int opt;
-
- if (strcmp(optarg, O_PHM_HERMAN_STR) == 0)
- opt = O_PHM_HERMAN;
- else if (strcmp(optarg, O_PHM_ROWLAND_STR) == 0)
- opt = O_PHM_ROWLAND;
- else if (strcmp(optarg, O_PHM_BROWLAND_STR) == 0)
- opt = O_PHM_BROWLAND;
- else if (strcmp(optarg, O_PHM_UNITPULSE_STR) == 0)
- opt = O_PHM_UNITPULSE;
- else {
- sys_error(ERR_WARNING,"Invalid phantom option %s\n", optarg);
- opt = -1;
- }
-
- return (opt);
-
-}
-
-InterpolationType
-opt_set_interpolation (const char *optarg)
-{
- InterpolationType opt;
-
- if (strcmp(optarg, O_INTERP_NEAREST_STR) == 0)
- opt = I_NEAREST;
- else if (strcmp(optarg, O_INTERP_LINEAR_STR) == 0)
- opt = I_LINEAR;
-#if HAVE_BSPLINE_INTERP
- else if (strcmp(optarg, O_INTERP_BSPLINE_STR) == 0)
- opt = I_BSPLINE;
-#endif
- else {
- sys_error(ERR_WARNING, "Invalid interpolation type %s\n", optarg);
- opt = static_cast<InterpolationType>(-1);
- }
-
- return (opt);
-}
-
-FilterType
-opt_set_filter (const char *optarg)
-{
- FilterType opt;
-
- if (strcmp(optarg, O_FILTER_BANDLIMIT_STR) == 0)
- opt = FILTER_BANDLIMIT;
- else if (strcmp(optarg, O_FILTER_HAMMING_STR) == 0)
- opt = FILTER_G_HAMMING;
- else if (strcmp(optarg, O_FILTER_SINC_STR) == 0)
- opt = FILTER_SINC;
- else if (strcmp(optarg, O_FILTER_COS_STR) == 0)
- opt = FILTER_COSINE;
- else if (strcmp(optarg, O_FILTER_TRIANGLE_STR) == 0)
- opt = FILTER_TRIANGLE;
- else if (strcmp(optarg, O_FILTER_ABS_BANDLIMIT_STR) == 0)
- opt = FILTER_ABS_BANDLIMIT;
- else if (strcmp(optarg, O_FILTER_ABS_HAMMING_STR) == 0)
- opt = FILTER_ABS_G_HAMMING;
- else if (strcmp(optarg, O_FILTER_ABS_SINC_STR) == 0)
- opt = FILTER_ABS_SINC;
- else if (strcmp(optarg, O_FILTER_ABS_COS_STR) == 0)
- opt = FILTER_ABS_COSINE;
- else if (strcmp(optarg, O_FILTER_SHEPP_STR) == 0)
- opt = FILTER_SHEPP;
- else {
- sys_error(ERR_WARNING, "Invalid filter type %s\n", optarg);
- opt = static_cast<FilterType>(-1);
- }
-
- return (opt);
-}
-
-const char *
-name_of_filter (const int filter)
-{
- const char *name = "Unknown filter";
-
- if (filter == FILTER_SHEPP)
- name = O_FILTER_SHEPP_STR;
- else if (filter == FILTER_ABS_COSINE)
- name = O_FILTER_ABS_COS_STR;
- else if (filter == FILTER_ABS_SINC)
- name = O_FILTER_ABS_SINC_STR;
- else if (filter == FILTER_ABS_G_HAMMING)
- name = O_FILTER_ABS_HAMMING_STR;
- else if (filter == FILTER_ABS_BANDLIMIT)
- name = O_FILTER_ABS_BANDLIMIT_STR;
- else if (filter == FILTER_COSINE)
- name = O_FILTER_COS_STR;
- else if (filter == FILTER_SINC)
- name = O_FILTER_SINC_STR;
- else if (filter == FILTER_G_HAMMING)
- name = O_FILTER_HAMMING_STR;
- else if (filter == FILTER_BANDLIMIT)
- name = O_FILTER_BANDLIMIT_STR;
- else if (filter == FILTER_TRIANGLE)
- name = O_FILTER_TRIANGLE_STR;
-
- return (name);
-}
-
-DomainType
-opt_set_filter_domain (const char *optarg)
-{
- DomainType opt;
-
- if (strcmp(optarg, D_SPATIAL_STR) == 0)
- opt = D_SPATIAL;
- else if (strcmp(optarg, D_FREQ_STR) == 0)
- opt = D_FREQ;
- else {
- sys_error(ERR_WARNING, "Invalid filter domain %s\n", optarg);
- opt = static_cast<DomainType>(-1);
- }
-
- return (opt);
-}
-
-const char *
-name_of_filter_domain (const DomainType domain)
-{
- const char *name = "Unknown domain";
-
- if (domain == D_SPATIAL)
- return(D_SPATIAL_STR);
- else if (domain == D_FREQ)
- return(D_FREQ_STR);
-
- return (name);
-}
-
-
-BackprojType
-opt_set_backproj (const char *optarg)
-{
- BackprojType opt;
-
- if (strcmp(optarg, O_BPROJ_TRIG_STR) == 0)
- opt = O_BPROJ_TRIG;
- else if (strcmp(optarg, O_BPROJ_TABLE_STR) == 0)
- opt = O_BPROJ_TABLE;
- else if (strcmp(optarg, O_BPROJ_DIFF_STR) == 0)
- opt = O_BPROJ_DIFF;
- else if (strcmp(optarg, O_BPROJ_DIFF2_STR) == 0)
- opt = O_BPROJ_DIFF2;
- else if (strcmp(optarg, O_BPROJ_IDIFF2_STR) == 0)
- opt = O_BPROJ_IDIFF2;
- else {
- sys_error(ERR_WARNING, "Invalid backprojection method %s\n", optarg);
- opt = static_cast<BackprojType>(-1);
- }
-
- return (opt);
-}
-
-const char *
-name_of_backproj(const BackprojType bproj)
-{
- const char *name = "Unknown backprojection method";
-
- if (bproj == O_BPROJ_TRIG)
- name = O_BPROJ_TRIG_STR;
- else if (bproj == O_BPROJ_TABLE)
- name = O_BPROJ_TABLE_STR;
- else if (bproj == O_BPROJ_DIFF)
- name = O_BPROJ_DIFF_STR;
- else if (bproj == O_BPROJ_DIFF2)
- name = O_BPROJ_DIFF2_STR;
- else if (bproj == O_BPROJ_IDIFF2)
- name = O_BPROJ_IDIFF2_STR;
-
- return (name);
-}
-
-
-
-/* NAME
- * name_of_interp Return name of interpolation method
- *
- * SYNOPSIS
- * name = name_of_interp (interp_type)
- * char *name Name of interpolation method
- * int interp_type Method of interpolation
- *
- * NOTES
- * Returns NULL if interp_type is invalid
- */
-
-const char *
-name_of_interpolation (int interp_type)
-{
- if (interp_type == I_NEAREST)
- return (O_INTERP_NEAREST_STR);
- else if (interp_type == I_LINEAR)
- return (O_INTERP_LINEAR_STR);
-#if HAVE_BSPLINE_INTERP
- else if (interp_type == I_BSPLINE)
- return (O_INTERP_BSPLINE_STR);
-#endif
- else
- return ("Unknown interpolation method");
-}
-
-
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: phantom.cpp,v 1.3 2000/06/19 20:08:09 kevin Exp $
+** $Id: phantom.cpp,v 1.4 2000/06/22 10:17:28 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
//
Phantom::Phantom (void)
+{
+ init ();
+}
+
+
+Phantom::Phantom (const char* const phmName)
+{
+ init ();
+ createFromPhantom (phmName);
+}
+
+void
+Phantom::init (void)
{
m_nPElem = 0;
m_xmin = 1E30;
m_ymax = -1E30;
m_diameter = 0;
m_composition = P_PELEMS;
+ m_fail = false;
+ m_id = PHM_INVALID;
}
-
Phantom::~Phantom (void)
{
for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++) {
}
-void
-Phantom::create (const int phmid)
+const char*
+Phantom::convertPhantomIDToName (const PhantomID phmID)
+{
+ const char *name = "";
+
+ if (phmID == PHM_HERMAN)
+ name = PHM_HERMAN_STR;
+ else if (phmID == PHM_BHERMAN)
+ name = PHM_BHERMAN_STR;
+ else if (phmID == PHM_ROWLAND)
+ name = PHM_ROWLAND_STR;
+ else if (phmID == PHM_BROWLAND)
+ name = PHM_BROWLAND_STR;
+ else if (phmID == PHM_UNITPULSE)
+ name = PHM_UNITPULSE_STR;
+
+ return (name);
+}
+
+Phantom::PhantomID
+Phantom::convertNameToPhantomID (const char* const phmName)
+{
+ PhantomID id;
+
+ if (strcasecmp (phmName, PHM_HERMAN_STR) == 0)
+ id = PHM_HERMAN;
+ else if (strcasecmp (phmName, PHM_BHERMAN_STR) == 0)
+ id = PHM_BHERMAN;
+ else if (strcasecmp (phmName, PHM_ROWLAND_STR) == 0)
+ id = PHM_ROWLAND;
+ else if (strcasecmp (phmName, PHM_BROWLAND_STR) == 0)
+ id = PHM_BROWLAND;
+ else if (strcasecmp (phmName, PHM_UNITPULSE_STR) == 0)
+ id = PHM_UNITPULSE;
+ else
+ id = PHM_INVALID;
+
+ return (id);
+}
+
+
+bool
+Phantom::createFromPhantom (const char* const phmName)
+{
+ PhantomID phmid = convertNameToPhantomID (phmName);
+ m_name = phmName;
+
+ createFromPhantom (phmid);
+}
+
+bool
+Phantom::createFromPhantom (const PhantomID phmid)
{
switch (phmid)
{
- case O_PHM_HERMAN:
+ case PHM_HERMAN:
addStdHerman();
break;
- case O_PHM_ROWLAND:
+ case PHM_BHERMAN:
+ addStdHermanBordered();
+ break;
+ case PHM_ROWLAND:
addStdRowland();
break;
- case O_PHM_BROWLAND:
+ case PHM_BROWLAND:
addStdRowlandBordered ();
break;
- case O_PHM_UNITPULSE:
+ case PHM_UNITPULSE:
m_composition = P_UNIT_PULSE;
addPElem ("rectangle", 0., 0., 100., 100., 0., 0.); // outline
addPElem ("ellipse", 0., 0., 1., 1., 0., 1.); // pulse
break;
default:
sys_error (ERR_WARNING, "Illegal phantom id %d\n", phmid);
+ m_name += " -- INVALID";
+ m_fail = true;
+ return false;
break;
}
+
+ m_id = phmid;
+
+ return true;
}
* draw ()
*/
+#ifdef HAVE_SGP
void
Phantom::draw (void) const
{
-#ifdef HAVE_SGP
for (PElemIterator i = m_listPElem.begin(); i != m_listPElem.end(); i++)
sgp2_polyline_abs ((*i)->xOutline(), (*i)->yOutline(), (*i)->nOutlinePoints());
-#endif
}
+#endif
/* NAME
addPElem("ellipse", 0.000, 0.00, 7.875, 5.7187, 90.00, -0.206);
}
+void
+Phantom::addStdHermanBordered (void)
+{
+ addPElem("ellipse", 0., 0., 6.6, 5.9, 90., 0.);
+}
+
/* NAME
* convertToImagefile Make image array from Phantom
{
m_rot = convertDegreesToRadians (rot); // convert angle to radians
- if (strcasecmp (type, "rectangle") == 0)
- m_type = PELEM_RECTANGLE;
- else if (strcasecmp (type, "triangle") == 0)
- m_type = PELEM_TRIANGLE;
- else if (strcasecmp (type, "ellipse") == 0)
- m_type = PELEM_ELLIPSE;
- else if (strcasecmp (type, "sector") == 0)
- m_type = PELEM_SECTOR;
- else if (strcasecmp (type, "segment") == 0)
- m_type = PELEM_SEGMENT;
- else {
- sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
- m_type = PELEM_INVALID;
- }
+ m_type = convertNameToType (type);
makeTransformMatrices (); // calc transform matrices between phantom and normalized phantomelement
makeVectorOutline (); // calculate vector outline of pelem
delete m_yOutline;
}
+PhmElemType
+PhantomElement::convertNameToType (const char* const typeName)
+{
+ PhmElemType type = PELEM_INVALID;
+
+ if (strcasecmp (typeName, "rectangle") == 0)
+ type = PELEM_RECTANGLE;
+ else if (strcasecmp (typeName, "triangle") == 0)
+ type = PELEM_TRIANGLE;
+ else if (strcasecmp (typeName, "ellipse") == 0)
+ type = PELEM_ELLIPSE;
+ else if (strcasecmp (typeName, "sector") == 0)
+ type = PELEM_SECTOR;
+ else if (strcasecmp (typeName, "segment") == 0)
+ type = PELEM_SEGMENT;
+ else
+ sys_error (ERR_WARNING, "Unknown PhantomElement type %s [PhantomElement::PhantomElement]", type);
+
+ return (type);
+}
+
void
PhantomElement::makeTransformMatrices (void)
{
* Called by phm_add_pelem()
*/
-static const double SCALE_PELEM_EXTENT=0.005; // increase pelem limits by 0.5%
-
void
PhantomElement::makeVectorOutline (void)
{
}
-
/* NAME
* calc_arc Calculate outline of a arc of a circle
*
*/
int
-PhantomElement::numCirclePoints (double theta) const
+PhantomElement::numCirclePoints (double theta)
{
- if (theta < 0.0 || theta > TWOPI)
- sys_error(ERR_WARNING, "illegal values sent to circle_pts");
+ theta = clamp (theta, 0., TWOPI);
return static_cast<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.3 2000/06/20 17:54:51 kevin Exp $
+** $Id: projections.cpp,v 1.4 2000/06/22 10:17:28 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
}
+
+/* NAME
+ * Projections::reconstruct Reconstruct Image from Projections
+ *
+ * SYNOPSIS
+ * im = proj.reconstruct (im, filt_type, filt_param, interp_type)
+ * IMAGE *im Output image
+ * int filt_type Type of convolution filter to use
+ * double filt_param Filter specific parameter
+ * Currently, used only with Hamming filters
+ * int interp_type Type of interpolation method to use
+ *
+ * ALGORITHM
+ *
+ * Calculate one-dimensional filter in spatial domain
+ * Allocate & clear (zero) the 2d output image array
+ * For each projection view
+ * Convolve raysum array with filter
+ * Backproject raysums and add (summate) to image array
+ * end
+ */
+
+bool
+Projections::reconstruct (ImageFile& im, const char* const filterName, double filt_param, const char* const interpName, int interp_param, const char* const backprojectName, const int trace)
+{
+ int nview = m_nView;
+ double det_inc = m_detInc;
+ double detlen = (m_nDet - 1) * det_inc;
+ int n_filtered_proj = m_nDet;
+ double filtered_proj [n_filtered_proj]; // convolved result
+
+#ifdef HAVE_BSPLINE_INTERP
+ int spline_order = 0, zoom_factor = 0;
+ if (interp_type == I_BSPLINE) {
+ zoom_factor = interp_param;
+ spline_order = 3;
+ zoom_factor = 3;
+ n_filtered_proj = (m_nDet - 1) * (zoom_factor + 1) + 1;
+ }
+#endif
+
+ int n_vec_filter = 2 * m_nDet - 1;
+ double filterBW = 1. / det_inc;
+ double filterMin = -detlen;
+ double filterMax = detlen;
+
+ SignalFilter filter (filterName, filterBW, filterMin, filterMax, n_vec_filter, filt_param, "spatial", 0);
+ if (filter.fail())
+ return false;
+
+ if (trace)
+ cout << "Reconstruct: filter="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
+
+#if HAVE_SGP
+ SGP_ID gid;
+ double plot_xaxis [n_vec_filter]; // array for plotting
+
+ if (trace > TRACE_TEXT) {
+ int i;
+ double f;
+ double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
+ for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
+ plot_xaxis[i] = f;
+
+ gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
+ cio_put_str ("Press any key to continue");
+ cio_kb_getc ();
+ sgp2_close (gid);
+ }
+ if (trace >= TRACE_TEXT) {
+ printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
+ }
+#endif //HAVE_SGP
+
+ Backprojector bj (*this, im, backprojectName, interpName);
+ if (bj.fail())
+ return false;
+
+ for (int iview = 0; iview < m_nView; iview++) {
+ if (trace >= TRACE_TEXT)
+ printf ("Reconstructing view %d (last = %d)\n", iview, m_nView - 1);
+
+ DetectorArray& darray = getDetectorArray (iview);
+ DetectorValue* detval = darray.detValues();
+
+ for (int j = 0; j < m_nDet; j++)
+ filtered_proj[j] = filter.convolve (detval, det_inc, j, m_nDet);
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT) {
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel ");
+ ezset ("ylabel ");
+ ezset ("xlength .5.");
+ ezset ("box.");
+ ezset ("grid.");
+ ezset ("ufinish yes.");
+ ezplot (detval, plot_xaxis, m_nDet);
+ ezset ("clear.");
+ ezset ("xticks major 5.");
+ ezset ("xlabel ");
+ ezset ("ylabel ");
+ ezset ("ustart yes.");
+ ezset ("xporigin .5.");
+ ezset ("xlength .5.");
+ ezset ("box");
+ ezset ("grid");
+ gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
+ }
+#endif //HAVE_SGP
+
+#ifdef HAVE_BSPLINE_INTERP
+ if (interp_type == I_BSPLINE)
+ bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
+ bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+ ezplot_1d (filtered_proj, n_filtered_proj);
+ }
+#endif
+#endif
+
+ bj.BackprojectView (filtered_proj, darray.viewAngle());
+
+#ifdef HAVE_SGP
+ if (trace >= TRACE_PLOT) {
+ char str[256];
+ printf ("Do you want to exit with current pic (y/n) -- ");
+ fgets(str, sizeof(str), stdin);
+ sgp2_close (sgp2_get_active_win());
+ if (tolower(str[0]) == 'y') {
+ break;
+ }
+ }
+#endif //HAVE_SGP
+ }
+
+ return true;
+}
+
+++ /dev/null
-/*****************************************************************************
-** FILE IDENTIFICATION
-**
-** Name: rec_t.c Image Reconstruction Procedures
-** Programmer: Kevin Rosenberg
-** Date Started: Aug 84
-**
-** GLOBAL FUNCTIONS
-** proj_reconst() Reconstruct Image from Projections
-**
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: reconstr.cpp,v 1.2 2000/06/20 17:54:51 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"
-
-
-/* NAME
- * proj_reconst Reconstruct Image from Projections
- *
- * SYNOPSIS
- * im = proj_reconst (im, proj, filt_type, filt_param, interp_type)
- * IMAGE *im Output image
- * RAYSUM *proj Raysum data if in memory
- * int filt_type Type of convolution filter to use
- * double filt_param Filter specific parameter
- * Currently, used only with Hamming filters
- * int interp_type Type of interpolation method to use
- *
- * ALGORITHM
- *
- * Calculate one-dimensional filter in spatial domain
- * Allocate & clear (zero) the 2d output image array
- * For each projection view
- * Convolve raysum array with filter
- * Backproject raysums and add (summate) to image array
- * end
- */
-
-ImageFile&
-proj_reconst (ImageFile& im, Projections& proj, const FilterType filt_type, double filt_param, InterpolationType interp_type, int interp_param, const BackprojType backproj_type, const int trace)
-{
- int ndet = proj.nDet();
- int nview = proj.nView();
- double det_inc = proj.detInc();
- double detlen = (ndet - 1) * det_inc;
- int n_filtered_proj = ndet;
- double filtered_proj [n_filtered_proj]; // convolved result
-
-#ifdef HAVE_BSPLINE_INTERP
- int spline_order = 0, zoom_factor = 0;
- if (interp_type == I_BSPLINE) {
- zoom_factor = interp_param;
- spline_order = 3;
- zoom_factor = 3;
- n_filtered_proj = (ndet - 1) * (zoom_factor + 1) + 1;
- }
-#endif
-
- int n_vec_filter = 2 * ndet - 1;
- double filterBW = 1. / det_inc;
- double filterMin = -detlen;
- double filterMax = detlen;
-
- SignalFilter filter (filt_type, filterBW, filterMin, filterMax, n_vec_filter, filt_param, D_SPATIAL, 0);
-
-#if HAVE_SGP
- SGP_ID gid;
- double plot_xaxis [n_vec_filter]; // array for plotting
-
- if (trace > TRACE_TEXT) {
- int i;
- double f;
- double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
- for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
- plot_xaxis[i] = f;
-
- gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
- cio_put_str ("Press any key to continue");
- cio_kb_getc ();
- sgp2_close (gid);
- }
- if (trace >= TRACE_TEXT) {
- printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", proj.nView(), proj.nDet(), proj.detStart(), proj.detInc());
- }
-#endif //HAVE_SGP
-
- Backproject* bj = selectBackprojector (backproj_type, proj, im, interp_type);
-
- for (int iview = 0; iview < nview; iview++) {
- if (trace >= TRACE_TEXT)
- printf ("Reconstructing view %d (last = %d)\n", iview, nview - 1);
-
- DetectorArray& darray = proj.getDetectorArray (iview);
- DetectorValue* detval = darray.detValues();
-
- for (int j = 0; j < ndet; j++)
- filtered_proj[j] = filter.convolve (detval, det_inc, j, ndet, FUNC_BOTH);
-
-#ifdef HAVE_SGP
- if (trace >= TRACE_PLOT) {
- ezset ("clear.");
- ezset ("xticks major 5.");
- ezset ("xlabel ");
- ezset ("ylabel ");
- ezset ("xlength .5.");
- ezset ("box.");
- ezset ("grid.");
- ezset ("ufinish yes.");
- ezplot (detval, plot_xaxis, proj.nDet());
- ezset ("clear.");
- ezset ("xticks major 5.");
- ezset ("xlabel ");
- ezset ("ylabel ");
- ezset ("ustart yes.");
- ezset ("xporigin .5.");
- ezset ("xlength .5.");
- ezset ("box");
- ezset ("grid");
- gid = ezplot (filtered_proj, plot_xaxis, proj.nDet());
- }
-#endif //HAVE_SGP
-
-#ifdef HAVE_BSPLINE_INTERP
- if (interp_type == I_BSPLINE)
- bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
-
-#ifdef HAVE_SGP
- if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
- bspline (proj.nDet(), zoom_factor, spline_order, filtered_proj, filtered_proj);
- ezplot_1d (filtered_proj, n_filtered_proj);
- }
-#endif
-#endif
-
- bj->BackprojectView (filtered_proj, darray.viewAngle());
-
-#ifdef HAVE_SGP
- if (trace >= TRACE_PLOT) {
- char str[256];
- printf ("Do you want to exit with current pic (y/n) -- ");
- fgets(str, sizeof(str), stdin);
- sgp2_close (sgp2_get_active_win());
- if (tolower(str[0]) == 'y') {
- break;
- }
- }
-#endif //HAVE_SGP
- }
-
- delete bj;
-
- return (im);
-}
-
noinst_LIBRARIES = libctsupport.a
INCLUDES=@my_includes@
-libctsupport_a_SOURCES= filefuncs.cpp strfuncs.cpp syserror.cpp timedate.cpp fnetorderstream.cpp audio.cpp consoleio.cpp normangle.cpp simpson.cpp xform.cpp clip.cpp
+libctsupport_a_SOURCES= filefuncs.cpp strfuncs.cpp syserror.cpp timedate.cpp fnetorderstream.cpp consoleio.cpp mathfuncs.cpp xform.cpp clip.cpp
EXTRA_DIST=Makefile.nt
+++ /dev/null
-/*****************************************************************************
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: audio.cpp,v 1.2 2000/06/19 19:04:05 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
-******************************************************************************/
-
-/* NAME
- * beep sound a beep to user
- *
- * SYNOPSIS
- * beep()
- */
-
-#include "ctsupport.h"
-
-void cio_beep (void)
-{
- cio_tone (2000.0, 0.3);
-}
-
-/* NAME
- * tone play a frequency sound for some duration
- *
- * SYNOPSIS
- * tone (freq, length)
- * double freq frequency to play in Hertz
- * double length duration to play note in seconds
- */
-
-void
-cio_tone (double freq, double length)
-{
-#if 1
- fprintf(stdout, "\007");
-#else
- cio_spkr_freq (freq); /* Set frequency of tone */
- cio_spkr_on (); /* Turn on speaker */
- pause (length); /* Pause for length seconds */
- cio_spkr_off (); /* Turn off speaker */
-#endif
-}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: clip.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $
+** $Id: clip.cpp,v 1.3 2000/06/22 10:17:28 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
xform_mtx2 (rotmtx, ccx, ccy);
t1 += theta; // rotate start and stop angles
t2 += theta;
- t1 = norm_ang (t1);
- t2 = norm_ang (t2);
+ t1 = normalizeAngle (t1);
+ t2 = normalizeAngle (t2);
if (xc2 < -D_EPSILON || fabs(yc2) > F_EPSILON) {
sys_error (ERR_SEVERE, "Internal error in clip_circle\n x1=%6.2f, y1=%6.2f, x2=%6.2f, y2=%6.2f, xc2=%6.2f, yc2=%6.2f, theta=%6.2f", x1, y1, x2, y2, xc2, yc2, theta);
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: consoleio.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $
+** $Id: consoleio.cpp,v 1.3 2000/06/22 10:17:28 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
}
+/* NAME
+ * beep sound a beep to user
+ *
+ * SYNOPSIS
+ * beep()
+ */
+
+void cio_beep (void)
+{
+ cio_tone (2000.0, 0.3);
+}
+
+/* NAME
+ * tone play a frequency sound for some duration
+ *
+ * SYNOPSIS
+ * tone (freq, length)
+ * double freq frequency to play in Hertz
+ * double length duration to play note in seconds
+ */
+
+void
+cio_tone (double freq, double length)
+{
+#if 1
+ fprintf(stdout, "\007");
+#else
+ cio_spkr_freq (freq); /* Set frequency of tone */
+ cio_spkr_on (); /* Turn on speaker */
+ pause (length); /* Pause for length seconds */
+ cio_spkr_off (); /* Turn off speaker */
+#endif
+}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: filefuncs.cpp,v 1.2 2000/06/19 19:04:05 kevin Exp $
+** $Id: filefuncs.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $
**
** Revision 1.1.1.1 2000/04/28 13:02:44 kevin
** Initial CVS import for first public release
/* NAME
- * file_exists Checks if a specified disk fie exists
+ * fileExists Checks if a specified disk fie exists
*
* SYNOPSIS
- * exist = file_exists (fname)
+ * exist = fileExists (fname)
* bool exist TRUE if specified file exists
*/
bool
-file_exists (const char *fname)
+fileExists (const char *fname)
{
FILE *fp;
bool exist;
return (exist);
}
-/*-----------------------------------------------------------------------------
- *
- * FUNCTION IDENTIFICATION
- *
- * Name: sys_fopen Open a file for user
- * Date: 12-17-84
- * Programmer: Kevin Rosenberg
- *
- * SYNOPSIS
- * fp = sys_fopen (filename, mode, progname)
- * FILE *fp Standard pointer to 'C' file
- * char *filename Name of file to open
- * If user enters a new name, it goes here
- * char *mode Mode to open file (std. 'C')
- * char *progname Name of program calling this routine
- *
- * DESCRIPTION
- * This routine opens a file using the standard C fopen() routine. If
- * the file is not found, the user is given to option to:
- *
- * 1 - Retry opening file with same name
- * 2 - Enter new file name
- * 3 - Abort and return to DOS
- *
- * CAUTIONS
- * If the the requested file is not found, the name of the file given
- * entered at keyboard will be returned in filename. So, make sure there
- * is room for a maximum length filename (MAXFULLNAME)
- *
- *---------------------------------------------------------------------------*/
-
-FILE *
-sys_fopen (const char *filename, const char *mode, const char *progname)
-{
- FILE *fp;
- char fname[256]; /* name used for call to fopen() */
- char c; /* keyboard response */
-
- strncpy (fname, filename, sizeof(fname));
-
- do {
- if ((fp = fopen (fname, mode)) == NULL) {
- cerr << endl;
- cerr << "Can't open file " << fname << " [" << progname << "]" << endl;
- cerr << "Enter: <R> - Retry | <N> - New name | <A> - Abort program --> ";
- c = cio_kb_waitc ("RrNnAa", TRUE);
- c = tolower(c);
- cerr << c << endl;
-
- if (c == 'r') // Retry -- Nothing to do here
- ;
- else if (c == 'a') // Abort -- Exit to OS
- exit (1);
- else if (c == 'n') { // New name - get from console
- cerr << "Enter new file name -- ";
- fgets (fname, sizeof(fname), stdin);
- str_wrm_tail (fname);
- cerr << endl;
- }
- }
- } while (fp == NULL);
-
- return (fp);
-}
--- /dev/null
+/*****************************************************************************
+** This is part of the CTSim program
+** Copyright (C) 1983-2000 Kevin Rosenberg
+**
+** $Id: mathfuncs.cpp,v 1.1 2000/06/22 10:17:28 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 "ctsupport.h"
+
+
+/* NAME
+ * integrateSimpson Integrate array of data by Simpson's rule
+ *
+ * SYNOPSIS
+ * double integrateSimpson (xmin, xmax, y, np)
+ * double xmin, xmax Extent of integration
+ * double y[] Function values to be integrated
+ * int np number of data points
+ * (must be an odd number and at least 3)
+ *
+ * RETURNS
+ * integrand of function
+ */
+
+double
+integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
+{
+ if (np < 2)
+ return (0.);
+ else if (np == 2)
+ return ((xmax - xmin) * (y[0] + y[1]) / 2);
+
+ double area = 0;
+ int nDiv = (np - 1) / 2; // number of divisions
+ double width = (xmax - xmin) / (double) (np - 1); // width of cells
+
+ for (int i = 1; i <= nDiv; i++) {
+ int xr = 2 * i;
+ int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
+ int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
+
+ area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
+ }
+
+ if ((np & 1) == 0) /* do last trapazoid */
+ area += width * (y[np-2] + y[np-1]) / 2;
+
+ return (area);
+}
+
+
+/* NAME
+ * norm_angle Normalize angle to 0 to 2 * PI range
+ *
+ * SYNOPSIS
+ * t = norm_angle (theta)
+ * double t Normalized angle
+ * double theta Input angle
+ */
+
+double
+normalizeAngle (double theta)
+{
+ while (theta < 0.)
+ theta += TWOPI;
+ while (theta >= TWOPI)
+ theta -= TWOPI;
+
+ return (theta);
+}
+++ /dev/null
-/*****************************************************************************
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: normangle.cpp,v 1.2 2000/06/19 19:04:05 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 "ctsupport.h"
-
-
-/* NAME
- * norm_ang Normalize angle to 0 to 2 * PI range
- *
- * SYNOPSIS
- * t = norm_ang (theta)
- * double t Normalized angle
- * double theta Input angle
- */
-
-double
-norm_ang (double theta)
-{
- while (theta < 0.)
- theta += TWOPI;
- while (theta >= TWOPI)
- theta -= TWOPI;
-
- return (theta);
-}
+++ /dev/null
-/*****************************************************************************
-** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
-**
-** $Id: simpson.cpp,v 1.2 2000/06/19 19:04:05 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 "ctsupport.h"
-
-
-/* NAME
- * integrateSimpson Integrate array of data by Simpson's rule
- *
- * SYNOPSIS
- * double integrateSimpson (xmin, xmax, y, np)
- * double xmin, xmax Extent of integration
- * double y[] Function values to be integrated
- * int np number of data points
- * (must be an odd number and at least 3)
- *
- * RETURNS
- * integrand of function
- */
-
-double
-integrateSimpson (const double xmin, const double xmax, const double *y, const int np)
-{
- if (np < 2)
- return (0.);
- else if (np == 2)
- return ((xmax - xmin) * (y[0] + y[1]) / 2);
-
- double area = 0;
- int nDiv = (np - 1) / 2; // number of divisions
- double width = (xmax - xmin) / (double) (np - 1); // width of cells
-
- for (int i = 1; i <= nDiv; i++) {
- int xr = 2 * i;
- int xl = xr - 2; // 2 * (i - 1) == 2 * i - 2 == xr - 2
- int xm = xr - 1; // (xl+xr)/2 == (xr+xr-2)/2 == (2*xr-2)/2 = xr-1
-
- area += (width / 3.0) * (y[xl] + 4.0 * y[xm] + y[xr]);
- }
-
- if ((np & 1) == 0) /* do last trapazoid */
- area += width * (y[np-2] + y[np-1]) / 2;
-
- return (area);
-}
ctrec-lam: ctrec.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
$(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI ctrec.cpp -o ctrec-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
-phm2pj-lam: phm2pj.cpp mpiworld.cpp ../include/ct.h
+phm2pj-lam: phm2pj.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
$(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2pj.cpp -o phm2pj-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
-phm2if-lam: phm2if.cpp mpiworld.cpp ../include/ct.h
+phm2if-lam: phm2if.cpp mpiworld.cpp ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
$(CC_LAM) @DEFS@ $(CFLAGS) $(INCLUDES) -DHAVE_MPI phm2if.cpp -o phm2if-lam $(LDFLAGS) $(LAM_EXTRA_SRC) @ctlibs@
endif
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: ctrec.cpp,v 1.12 2000/06/19 17:58:13 kevin Exp $
+** $Id: ctrec.cpp,v 1.13 2000/06/22 10:17:28 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
int opt_debug = 0;
int opt_trace = TRACE_NONE;
double opt_filter_param = 1;
- FilterType opt_filter = FILTER_ABS_BANDLIMIT;
- InterpolationType opt_interp = I_LINEAR;
+ string optFilterName = "abs_bandlimit";
+ string optInterpName = "linear";
+ string optBackprojName = "idiff2";
int opt_interp_param = 1;
- BackprojType opt_backproj = O_BPROJ_DIFF2;
int nx, ny;
#ifdef HAVE_MPI
ImageFile *imLocal;
switch (c)
{
case O_INTERP:
- if ((opt_interp = opt_set_interpolation(optarg)) < 0) {
- ctrec_usage(argv[0]);
- return (1);
- }
+ optInterpName = optarg;
break;
case O_FILTER:
- if ((opt_filter = opt_set_filter(optarg)) < 0) {
- ctrec_usage(argv[0]);
- return (1);
- }
+ optFilterName = optarg;
break;
case O_BACKPROJ:
- if ((opt_backproj = opt_set_backproj(optarg)) < 0) {
- ctrec_usage(argv[0]);
- return (1);
- }
+ optBackprojName = optarg;
break;
case O_FILTER_PARAM:
opt_filter_param = strtod(optarg, &endptr);
nx = strtol(argv[optind + 2], &endptr, 10);
ny = strtol(argv[optind + 3], &endptr, 10);
- ostringstream filt_name;
- if (opt_filter == FILTER_G_HAMMING || opt_filter == FILTER_ABS_G_HAMMING)
- filt_name << name_of_filter (opt_filter) << ": alpha=" << opt_filter_param;
+ ostringstream filterDesc;
+ if (opt_filter_param >= 0)
+ filterDesc << optFilterName << ": alpha=" << opt_filter_param;
else
- filt_name << name_of_filter (opt_filter);
-
+ filterDesc << optFilterName;
+
ostringstream label;
- label << "Reconstruct: " << nx << "x" << ny << ", " << filt_name.str() << ", " << name_of_interpolation (opt_interp) << ", " << name_of_backproj(opt_backproj);
+ label << "Reconstruct: " << nx << "x" << ny << ", " << filterDesc.str() << ", " << optInterpName << ", " << optBackprojName;
remark = label.str();
if (opt_verbose)
}
TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+ mpiWorld.BcastString (optBackprojName);
+ mpiWorld.BcastString (optFilterName);
+ mpiWorld.BcastString (optInterpName);
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);
- mpiWorld.getComm().Bcast (&opt_filter, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_interp, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
mpiWorld.getComm().Bcast (&opt_interp_param, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_backproj, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&mpi_ndet, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&mpi_nview, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&mpi_detinc, 1, MPI::DOUBLE, 0);
#ifdef HAVE_MPI
TimerCollectiveMPI timerReconstruct (mpiWorld.getComm());
- proj_reconst (*imLocal, projLocal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace);
+ projLocal.reconstruct (*imLocal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace);
if (opt_verbose)
timerReconstruct.timerEndAndReport ("Time to reconstruct");
if (opt_verbose)
timerReduce.timerEndAndReport ("Time to reduce image");
#else
- proj_reconst (*imGlobal, projGlobal, opt_filter, opt_filter_param, opt_interp, opt_interp_param, opt_backproj, opt_trace);
+ projGlobal.reconstruct (*imGlobal, optFilterName.c_str(), opt_filter_param, optInterpName.c_str(), opt_interp_param, optBackprojName.c_str(), opt_trace);
#endif
#ifdef HAVE_MPI
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: if2img.cpp,v 1.7 2000/06/19 17:58:13 kevin Exp $
+** $Id: if2img.cpp,v 1.8 2000/06/22 10:17:28 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
int ny = im.ny();
ImageFileArray v = im.getArray();
- unsigned char* rowp = new unsigned char [nx * nxcell];
+ unsigned char rowp [nx * nxcell];
if ((fp = fopen (outfile, "wb")) == NULL)
return;
fprintf(fp, "%c ", rowp[ic]);
}
}
- delete rowp;
fclose(fp);
}
int ny = im.ny();
ImageFileArray v = im.getArray();
- unsigned char* rowp = new unsigned char [nx * nxcell];
+ unsigned char rowp [nx * nxcell];
if (rowp == NULL)
return;
fprintf(fp, "\n");
}
}
- delete rowp;
fclose(fp);
}
int ny = im.ny();
ImageFileArray v = im.getArray();
- unsigned char* rowp = new unsigned char [nx * nxcell * (bitdepth / 8)];
+ unsigned char rowp [nx * nxcell * (bitdepth / 8)];
if ((fp = fopen (outfile, "wb")) == NULL)
return;
for (int ir = 0; ir < nycell; ir++)
png_write_rows (png_ptr, &row_pointer, 1);
}
- delete rowp;
png_write_end(png_ptr, info_ptr);
png_destroy_write_struct(&png_ptr, &info_ptr);
int ny = im.ny();
ImageFileArray v = im.getArray();
- usnigned char* rowp = new unsigned char [nx * nxcell];
+ unsigned char rowp [nx * nxcell];
if (rowp == NULL)
return;
}
}
}
- delete rowp;
if ((out = fopen(outfile,"w")) == NULL) {
sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile);
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: mpiworld.cpp,v 1.1 2000/06/09 01:35:33 kevin Exp $
+** $Id: mpiworld.cpp,v 1.2 2000/06/22 10:17:28 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
}
+void
+MPIWorld::BcastString (string& str)
+{
+ int len;
+
+ if (m_myRank == 0)
+ len = str.length();
+ m_comm.Bcast (&len, 1, MPI::INT, 0);
+ char buf [ len + 1];
+
+ if (m_myRank == 0)
+ strcpy (buf, str.c_str());
+
+ m_comm.Bcast (buf, len + 1, MPI::CHAR, 0);
+
+ if (m_myRank > 0)
+ str = buf;
+}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: phm2if.cpp,v 1.12 2000/06/19 20:08:09 kevin Exp $
+** $Id: phm2if.cpp,v 1.13 2000/06/22 10:17:28 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
cout << " ny Number of pixels Y-axis" << endl;
cout << " --phantom Phantom to use for projection" << endl;
cout << " herman Herman head phantom" << endl;
+ cout << " bherman Bordered Herman head phantom" << endl;
cout << " rowland Rowland head phantom" << endl;
cout << " browland Bordered Rowland head phantom" << endl;
cout << " unitpulse Unit pulse phantom" << endl;
Phantom phm;
int opt_nx = 0, opt_ny = 0;
int opt_nsample = 1;
- int opt_phmnum = -1;
- FilterType opt_filter = static_cast<FilterType>(-1);
- DomainType opt_filter_domain = D_SPATIAL;
+ string optPhmName;
+ string optFilterName;
+ string optDomainName = "spatial";
char *opt_outfile = NULL;
int opt_debug = 0;
string opt_desc;
switch (c) {
case O_PHANTOM:
- if ((opt_phmnum = opt_set_phantom(optarg)) < 0) {
+ optPhmName = optarg;
+ if (! phm.createFromPhantom (optPhmName.c_str())) {
+ cout << "Invalid phantom name " << optPhmName << endl << endl;
phm2if_usage(argv[0]);
return (1);
}
}
break;
case O_FILTER:
- if ((opt_filter = opt_set_filter(optarg)) < 0) {
- phm2if_usage (argv[0]);
- return (1);
- }
+ optFilterName = optarg;
break;
case O_FILTER_DOMAIN:
- if ((opt_filter_domain = opt_set_filter_domain(optarg)) < 0) {
- phm2if_usage (argv[0]);
- return (1);
- }
+ optDomainName = optarg;
break;
case O_DESC:
opt_desc = optarg;
}
}
- if (phm.nPElem() == 0 && opt_phmnum == -1 && opt_filter == -1) {
+ if (phm.nPElem() == 0 && optPhmName == "" && optFilterName == "") {
cerr << "No phantom defined" << endl;
phm2if_usage(argv[0]);
return (1);
oss << "nx=" << opt_nx << ", ny=" << opt_ny << ", nsample=" << opt_nsample << ", ";
if (opt_phmfilename.length())
oss << "phantom=" << opt_phmfilename;
- else if (opt_phmnum != -1)
- oss << "phantom=" << name_of_phantom(opt_phmnum);
- else if (opt_filter != -1) {
- oss << "filter=" << name_of_filter(opt_filter);
- oss << " - " << name_of_filter_domain(opt_filter_domain);
+ else if (optPhmName != "")
+ oss << "phantom=" << optPhmName;
+ else if (optFilterName != "") {
+ oss << "filter=" << optFilterName << " - " << optDomainName;
}
if (opt_desc.length())
oss << ": " << opt_desc;
#ifdef HAVE_MPI
TimerCollectiveMPI timerBcast (mpiWorld.getComm());
+ mpiWorld.BcastString (optPhmName);
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);
mpiWorld.getComm().Bcast (&opt_nx, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&opt_ny, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&opt_nsample, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_phmnum, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_filter, 1, MPI::INT, 0);
- mpiWorld.getComm().Bcast (&opt_filter_domain, 1, MPI::INT, 0);
mpiWorld.getComm().Bcast (&opt_filter_param, 1, MPI::DOUBLE, 0);
mpiWorld.getComm().Bcast (&opt_filter_bw, 1, MPI::DOUBLE, 0);
+ mpiWorld.BcastString (optFilterName);
+ mpiWorld.BcastString (optDomainName);
+
if (opt_verbose)
timerBcast.timerEndAndReport ("Time to broadcast variables");
mpiWorld.setTotalWorkUnits (opt_nx);
+ if (mpiWorld.getRank() > 0 && optPhmName != "")
+ phm.createFromPhantom (optPhmName.c_str());
+
if (mpiWorld.getRank() == 0) {
imGlobal = new ImageFile (opt_outfile, opt_nx, opt_ny);
imGlobal->fileCreate();
imGlobal->fileCreate ();
#endif
- if (opt_phmnum >= 0)
- phm.create (opt_phmnum);
-
-#ifdef HAVE_MPI
- else {
- if (mpiWorld.getRank() == 0)
- cerr << "phmnum < 0" << endl;
- return (1);
- }
-#endif
-
ImageFileArray v;
#ifdef HAVE_MPI
if (mpiWorld.getRank() == 0)
if (mpiWorld.getRank() == 0) {
v[opt_nx/2][opt_ny/2] = 1.;
}
- } else if (opt_filter != -1) {
+ } else if (optFilterName != "") {
if (mpiWorld.getRank() == 0) {
- image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
+ image_filter_response (*imGlobal, optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param, opt_trace);
}
} else {
TimerCollectiveMPI timerRasterize (mpiWorld.getComm());
v = imGlobal->getArray ();
if (phm.getComposition() == P_UNIT_PULSE) {
v[opt_nx/2][opt_ny/2] = 1.;
- } else if (opt_filter != -1) {
- image_filter_response (*imGlobal, opt_filter_domain, opt_filter_bw, opt_filter, opt_filter_param, opt_trace);
+ } else if (optFilterName != "") {
+ image_filter_response (*imGlobal, optDomainName.c_str(), opt_filter_bw, optFilterName.c_str(), opt_filter_param, opt_trace);
} else {
#if HAVE_SGP
if (opt_trace >= TRACE_PHM)
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: phm2pj.cpp,v 1.2 2000/06/19 17:58:13 kevin Exp $
+** $Id: phm2pj.cpp,v 1.3 2000/06/22 10:17:28 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
cout << " nview Number of rotated views" << endl;
cout << " --phantom Phantom to use for projection" << endl;
cout << " herman Herman head phantom" << endl;
+ cout << " bherman Bordered herman head phantom" << endl;
cout << " rowland Rowland head phantom" << endl;
cout << " browland Bordered Rowland head phantom" << endl;
cout << " unitpulse Unit pulse phantom" << endl;
int opt_nview;
int opt_nray = 1;
int opt_trace = 0;
- int optPhmNum = -1;
+ string optPhmName;
int opt_verbose = 0;
int opt_debug = 0;
double opt_rotangle = 1;
switch (c) {
case O_PHANTOM:
- if ((optPhmNum = opt_set_phantom (optarg)) < 0) {
+ optPhmName = optarg;
+ if (! phm.createFromPhantom (optPhmName.c_str())) {
+ cout << "ERROR: Invalid phantom name " << optPhmName << endl << endl;
phm2pj_usage(argv[0]);
return (1);
}
- phm.create (optPhmNum);
break;
case O_PHMFILE:
#ifdef HAVE_MPI
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 (optPhmNum != -1) {
- desc << "Phantom=" << name_of_phantom(optPhmNum);
+ } else if (optPhmName != "") {
+ desc << "Phantom=" << optPhmName;
}
if (opt_desc.length()) {
desc << ": " << opt_desc;
#ifdef HAVE_MPI
TimerCollectiveMPI timerBcast(mpiWorld.getComm());
+ mpiWorld.BcastString (optPhmName);
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 (&optPhmNum, 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);
timerBcast.timerEndAndReport ("Time to broadcast variables");
- if (mpiWorld.getRank() > 0 && optPhmNum >= 0)
- phm.create (optPhmNum);
+ if (mpiWorld.getRank() > 0 && optPhmName != "")
+ phm.createFromPhantom (optPhmName.c_str());
#endif
opt_rotangle *= PI;
projects / ctsim.git / commitdiff