From: Kevin M. Rosenberg Date: Thu, 22 Jun 2000 10:17:28 +0000 (+0000) Subject: r117: *** empty log message *** X-Git-Tag: debian-4.5.3-3~900 X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=2d39e823ba389fc68e5317c422b55be006094252 r117: *** empty log message *** --- diff --git a/ChangeLog b/ChangeLog index dfaae59..58f4e14 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,14 @@ +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 diff --git a/TODO b/TODO index b39cb5e..8bb0008 100644 --- a/TODO +++ b/TODO @@ -4,8 +4,6 @@ MISCELLANEOUS 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 ======================= diff --git a/configure b/configure index f99dc3f..3301d0c 100755 --- a/configure +++ b/configure @@ -711,7 +711,7 @@ fi 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; } diff --git a/configure.in b/configure.in index b48e0f3..4a38772 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/ctrec.cpp) -AM_INIT_AUTOMAKE(ctsim,1.9.4) +AM_INIT_AUTOMAKE(ctsim,1.9.5) AM_CONFIG_HEADER(config.h) dnl Checks for programs. diff --git a/include/backprojectors.h b/include/backprojectors.h index c28a11b..08b72d4 100644 --- a/include/backprojectors.h +++ b/include/backprojectors.h @@ -9,17 +9,7 @@ ** 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 @@ -35,21 +25,85 @@ ** 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; @@ -71,7 +125,7 @@ class Backproject 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) {} @@ -82,7 +136,7 @@ class BackprojectTrig : public Backproject 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); @@ -98,7 +152,7 @@ class BackprojectTable : public Backproject 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); @@ -112,7 +166,7 @@ class BackprojectDiff : public Backproject 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) {} @@ -122,7 +176,7 @@ class BackprojectDiff2 : public BackprojectDiff 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) {} @@ -130,4 +184,3 @@ class BackprojectIntDiff2 : public BackprojectDiff }; -Backproject* selectBackprojector (BackprojType type, const Projections& proj, ImageFile& im, InterpolationType interpType); diff --git a/include/ct.h b/include/ct.h index 4dccb37..ead95ad 100644 --- a/include/ct.h +++ b/include/ct.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ct.h,v 1.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 @@ -109,6 +109,18 @@ extern "C" { #include /* Standard ints on Linux */ #endif +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; + #ifdef HAVE_MPI #include "mpi++.h" #include "mpiworld.h" @@ -122,82 +134,19 @@ extern "C" { #include "sgp.h" #endif -#include -#include -#include -#include -#include -#include -#include -#include - -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"; @@ -215,45 +164,6 @@ enum { 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 @@ -267,27 +177,12 @@ int filter_select (double *filter_param); // 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 diff --git a/include/ctsupport.h b/include/ctsupport.h index 3cc00f6..f180091 100644 --- a/include/ctsupport.h +++ b/include/ctsupport.h @@ -9,7 +9,7 @@ ** 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 @@ -150,9 +150,8 @@ typedef unsigned char kuint8; /* 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); @@ -251,8 +250,6 @@ int clip_sector (double& x1, double& y1, double& x2, double& y2, const double u, 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); @@ -265,7 +262,8 @@ void rotate2d (double x[], double y[], int pts, double angle); 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); diff --git a/include/filter.h b/include/filter.h index 274d792..b544a1f 100644 --- a/include/filter.h +++ b/include/filter.h @@ -1,29 +1,92 @@ +#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; @@ -36,3 +99,5 @@ class SignalFilter { }; + +#endif diff --git a/include/mpiworld.h b/include/mpiworld.h index 2889821..5818ed3 100644 --- a/include/mpiworld.h +++ b/include/mpiworld.h @@ -27,7 +27,7 @@ #include #include - +#include class MPIWorld { @@ -62,7 +62,9 @@ class MPIWorld MPI::Intracomm& getComm() { return m_comm; } - + + void BcastString (string& str); + private: int m_myRank; int m_nProcessors; diff --git a/include/phantom.h b/include/phantom.h index 32eedf6..b3765c6 100644 --- a/include/phantom.h +++ b/include/phantom.h @@ -9,7 +9,7 @@ ** 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 @@ -95,6 +95,9 @@ class PhantomElement 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); @@ -104,7 +107,7 @@ class PhantomElement 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); }; @@ -120,7 +123,19 @@ typedef enum { 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) @@ -129,7 +144,9 @@ class Phantom 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); @@ -141,12 +158,22 @@ class Phantom 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); @@ -170,6 +197,19 @@ class Phantom double m_diameter; // diameter of object mutable slist 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::iterator PElemIterator; diff --git a/include/projections.h b/include/projections.h index 104476d..a9970b3 100644 --- a/include/projections.h +++ b/include/projections.h @@ -9,7 +9,7 @@ ** 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 @@ -51,6 +51,8 @@ class Projections 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;} diff --git a/include/scanner.h b/include/scanner.h index 005c802..d383498 100644 --- a/include/scanner.h +++ b/include/scanner.h @@ -9,7 +9,7 @@ ** 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 @@ -29,7 +29,8 @@ #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 diff --git a/include/trace.h b/include/trace.h new file mode 100644 index 0000000..d302d3c --- /dev/null +++ b/include/trace.h @@ -0,0 +1,48 @@ +#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 + diff --git a/libctsim/Makefile.am b/libctsim/Makefile.am index 70b8fa5..b932178 100644 --- a/libctsim/Makefile.am +++ b/libctsim/Makefile.am @@ -1,5 +1,5 @@ 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 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index ef0de98..9e35092 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** 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 @@ -26,6 +26,25 @@ #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 (...) @@ -34,24 +53,122 @@ // 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(new BackprojectTrig (proj, im, interpType)); - else if (bjType == O_BPROJ_TABLE) - bj = static_cast(new BackprojectTable (proj, im, interpType)); - else if (bjType == O_BPROJ_DIFF) - bj = static_cast(new BackprojectDiff (proj, im, interpType)); - else if (bjType == O_BPROJ_DIFF2) - bj = static_cast(new BackprojectDiff2 (proj, im, interpType)); - else if (bjType == O_BPROJ_IDIFF2) - bj = static_cast(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(new BackprojectTrig (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_TABLE) + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_DIFF) + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_DIFF2) + m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation)); + else if (m_idBackproject == BPROJ_IDIFF2) + m_pBackprojectImplem = static_cast(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 (""); } @@ -61,7 +178,7 @@ Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, // 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(); @@ -82,7 +199,7 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Interpol 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); } @@ -133,14 +250,14 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double 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 (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(pFloor); @@ -160,7 +277,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // 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); @@ -193,14 +310,14 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl 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(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(dPosFloor); @@ -222,7 +339,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Backprojects by precalculating the change in L position for each x & y step in the image. // Iterates in x & y direction by adding difference in L position -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] @@ -255,14 +372,14 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double #ifdef DEBUG printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos); #endif - if (interpType == I_NEAREST) { + if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest(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(detPosFloor); @@ -306,14 +423,14 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl #ifdef DEBUG printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(L))]); #endif - if (interpType == I_NEAREST) { + if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (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(detPosFloor); double frac = curDetPos - detPosFloor; // fraction distance from det @@ -352,7 +469,7 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do 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 @@ -360,7 +477,7 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do 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) { diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp index b2d3062..c62eb24 100644 --- a/libctsim/filter.cpp +++ b/libctsim/filter.cpp @@ -9,7 +9,7 @@ ** 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 @@ -29,63 +29,192 @@ /* 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 @@ -105,15 +234,15 @@ SignalFilter::~SignalFilter (void) 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 { @@ -125,7 +254,7 @@ SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double 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); @@ -148,17 +277,17 @@ SignalFilter::spatialResponseCalc (FilterType fType, double bw, double x, double 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.; @@ -209,7 +338,7 @@ SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double p 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); @@ -234,11 +363,11 @@ SignalFilter::frequencyResponse (FilterType fType, double bw, double u, double p 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; @@ -246,7 +375,7 @@ SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, do double b = PI / bw; double b2 = TWOPI / bw; - switch (fType) { + switch (filterID) { case FILTER_BANDLIMIT: q = bw * sinc(u * w, 1.0); break; @@ -284,7 +413,7 @@ SignalFilter::spatialResponseAnalytic (FilterType fType, double bw, double x, do 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; } @@ -356,68 +485,36 @@ SignalFilter::integral_abscos (double u, double w) */ 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); } diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 4cae5d2..bcb07e2 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -9,7 +9,7 @@ ** 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 @@ -89,17 +89,11 @@ Array2dFileLabel::getDateString (void) const /* 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; @@ -109,10 +103,7 @@ image_filter_response (ImageFile& im, const DomainType domain, double bw, const 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]); } @@ -137,7 +128,7 @@ int image_display_scale (const ImageFile& im, const int scale, const double pmin 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); @@ -149,7 +140,7 @@ int image_display_scale (const ImageFile& im, const int scale, const double pmin 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((v[ix][iy] - pmin) * view_scale); if (cval < 0) cval = 0; else if (cval > 255) @@ -160,7 +151,6 @@ int image_display_scale (const ImageFile& im, const int scale, const double pmin g2_image (id_X11, 0., 0., nx, ny, pens); - delete [] pens; return (id_X11); #endif diff --git a/libctsim/options.cpp b/libctsim/options.cpp deleted file mode 100644 index 9bbb585..0000000 --- a/libctsim/options.cpp +++ /dev/null @@ -1,271 +0,0 @@ -/***************************************************************************** -** 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(-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(-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(-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(-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"); -} - - diff --git a/libctsim/phantom.cpp b/libctsim/phantom.cpp index ec5d0dc..c080493 100644 --- a/libctsim/phantom.cpp +++ b/libctsim/phantom.cpp @@ -9,7 +9,7 @@ ** 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 @@ -33,6 +33,19 @@ // Phantom::Phantom (void) +{ + init (); +} + + +Phantom::Phantom (const char* const phmName) +{ + init (); + createFromPhantom (phmName); +} + +void +Phantom::init (void) { m_nPElem = 0; m_xmin = 1E30; @@ -41,9 +54,10 @@ Phantom::Phantom (void) 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++) { @@ -52,29 +66,89 @@ Phantom::~Phantom (void) } -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; } @@ -238,14 +312,14 @@ Phantom::show (void) const * 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 @@ -308,6 +382,12 @@ Phantom::addStdHerman (void) 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 @@ -408,20 +488,7 @@ PhantomElement::PhantomElement (const char *type, const double cx, const double { 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 @@ -446,6 +513,27 @@ PhantomElement::~PhantomElement (void) 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) { @@ -486,8 +574,6 @@ 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) { @@ -578,7 +664,6 @@ PhantomElement::makeVectorOutline (void) } - /* NAME * calc_arc Calculate outline of a arc of a circle * @@ -642,10 +727,9 @@ PhantomElement::calcEllipsePoints (double x[], double y[], const int pts, const */ 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 (POINTS_PER_CIRCLE * theta / TWOPI + 1.5); } diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 54f8f7d..5ce0a8e 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** 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 @@ -427,3 +427,145 @@ Projections::printScanInfo (void) const } + +/* 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="<= 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; +} + diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp deleted file mode 100644 index 9c88780..0000000 --- a/libctsim/reconstr.cpp +++ /dev/null @@ -1,170 +0,0 @@ -/***************************************************************************** -** 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); -} - diff --git a/libctsupport/Makefile.am b/libctsupport/Makefile.am index df80a70..0082f1f 100644 --- a/libctsupport/Makefile.am +++ b/libctsupport/Makefile.am @@ -1,6 +1,6 @@ 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 diff --git a/libctsupport/audio.cpp b/libctsupport/audio.cpp deleted file mode 100644 index 2def9d2..0000000 --- a/libctsupport/audio.cpp +++ /dev/null @@ -1,55 +0,0 @@ -/***************************************************************************** -** 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 -} diff --git a/libctsupport/clip.cpp b/libctsupport/clip.cpp index 9c2abbc..ebf94c2 100644 --- a/libctsupport/clip.cpp +++ b/libctsupport/clip.cpp @@ -14,7 +14,7 @@ ** 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 @@ -151,8 +151,8 @@ clip_circle (double& x1, double& y1, double& x2, double& y2, const double cx, co 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); diff --git a/libctsupport/consoleio.cpp b/libctsupport/consoleio.cpp index bc6df45..fd895e8 100644 --- a/libctsupport/consoleio.cpp +++ b/libctsupport/consoleio.cpp @@ -2,7 +2,7 @@ ** 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 @@ -144,3 +144,36 @@ cio_kb_waitc (const char *astr, int beep_on_error) } +/* 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 +} diff --git a/libctsupport/filefuncs.cpp b/libctsupport/filefuncs.cpp index 511251d..ab26634 100644 --- a/libctsupport/filefuncs.cpp +++ b/libctsupport/filefuncs.cpp @@ -2,7 +2,7 @@ ** 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 @@ -38,15 +38,15 @@ fileBasename (const char* const filename) /* 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; @@ -63,67 +63,3 @@ file_exists (const char *fname) 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: - Retry | - New name | - 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); -} diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp new file mode 100644 index 0000000..bc3b383 --- /dev/null +++ b/libctsupport/mathfuncs.cpp @@ -0,0 +1,83 @@ +/***************************************************************************** +** 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); +} diff --git a/libctsupport/normangle.cpp b/libctsupport/normangle.cpp deleted file mode 100644 index 1263801..0000000 --- a/libctsupport/normangle.cpp +++ /dev/null @@ -1,42 +0,0 @@ -/***************************************************************************** -** 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); -} diff --git a/libctsupport/simpson.cpp b/libctsupport/simpson.cpp deleted file mode 100644 index e7a5270..0000000 --- a/libctsupport/simpson.cpp +++ /dev/null @@ -1,62 +0,0 @@ -/***************************************************************************** -** 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); -} diff --git a/src/Makefile.am b/src/Makefile.am index 42c90f2..383818f 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -40,10 +40,10 @@ LAM_EXTRA_SRC = mpiworld.cpp 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 diff --git a/src/ctrec.cpp b/src/ctrec.cpp index 2d94be9..39688cf 100644 --- a/src/ctrec.cpp +++ b/src/ctrec.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: ctrec.cpp,v 1.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 @@ -112,10 +112,10 @@ ctrec_main (int argc, char * argv[]) 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; @@ -139,22 +139,13 @@ ctrec_main (int argc, char * argv[]) 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); @@ -203,14 +194,14 @@ ctrec_main (int argc, char * argv[]) 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) @@ -233,14 +224,14 @@ ctrec_main (int argc, char * argv[]) } 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); @@ -280,7 +271,7 @@ ctrec_main (int argc, char * argv[]) #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"); @@ -289,7 +280,7 @@ ctrec_main (int argc, char * argv[]) 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 diff --git a/src/if2img.cpp b/src/if2img.cpp index c3a4bfe..a230974 100644 --- a/src/if2img.cpp +++ b/src/if2img.cpp @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: if2img.cpp,v 1.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 @@ -451,7 +451,7 @@ sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densm 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; @@ -474,7 +474,6 @@ sdf2d_to_pgm (ImageFile& im, char *outfile, int nxcell, int nycell, double densm fprintf(fp, "%c ", rowp[ic]); } } - delete rowp; fclose(fp); } @@ -487,7 +486,7 @@ sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double de 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; @@ -513,7 +512,6 @@ sdf2d_to_pgmasc (ImageFile& im, char *outfile, int nxcell, int nycell, double de fprintf(fp, "\n"); } } - delete rowp; fclose(fp); } @@ -531,7 +529,7 @@ sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell 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; @@ -580,7 +578,6 @@ sdf2d_to_png (ImageFile& im, char *outfile, int bitdepth, int nxcell, int nycell 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); @@ -605,7 +602,7 @@ sdf2d_to_gif (ImageFile& im, char *outfile, int nxcell, int nycell, double densm 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; @@ -628,7 +625,6 @@ sdf2d_to_gif (ImageFile& im, char *outfile, int nxcell, int nycell, double densm } } } - delete rowp; if ((out = fopen(outfile,"w")) == NULL) { sys_error(ERR_FATAL, "Error opening output file %s for writing", outfile); diff --git a/src/mpiworld.cpp b/src/mpiworld.cpp index f5631f1..3b4b993 100644 --- a/src/mpiworld.cpp +++ b/src/mpiworld.cpp @@ -9,7 +9,7 @@ ** 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 @@ -64,3 +64,21 @@ MPIWorld::setTotalWorkUnits(int totalWorkUnits) } +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; +} diff --git a/src/phm2if.cpp b/src/phm2if.cpp index 60e504b..75119c9 100644 --- a/src/phm2if.cpp +++ b/src/phm2if.cpp @@ -9,7 +9,7 @@ ** 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 @@ -61,6 +61,7 @@ phm2if_usage (const char *program) 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; @@ -107,9 +108,9 @@ phm2if_main (int argc, char* argv[]) Phantom phm; int opt_nx = 0, opt_ny = 0; int opt_nsample = 1; - int opt_phmnum = -1; - FilterType opt_filter = static_cast(-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; @@ -139,7 +140,9 @@ phm2if_main (int argc, char* argv[]) 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); } @@ -166,16 +169,10 @@ phm2if_main (int argc, char* argv[]) } 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; @@ -223,7 +220,7 @@ phm2if_main (int argc, char* argv[]) } } - 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); @@ -253,11 +250,10 @@ phm2if_main (int argc, char* argv[]) 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; @@ -271,23 +267,27 @@ phm2if_main (int argc, char* argv[]) #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(); @@ -298,17 +298,6 @@ phm2if_main (int argc, char* argv[]) 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) @@ -318,9 +307,9 @@ phm2if_main (int argc, char* argv[]) 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()); @@ -337,8 +326,8 @@ phm2if_main (int argc, char* argv[]) 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) diff --git a/src/phm2pj.cpp b/src/phm2pj.cpp index fa8d0a9..a4510c4 100644 --- a/src/phm2pj.cpp +++ b/src/phm2pj.cpp @@ -9,7 +9,7 @@ ** 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 @@ -59,6 +59,7 @@ phm2pj_usage (const char *program) 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; @@ -95,7 +96,7 @@ phm2pj_main (int argc, char* argv[]) 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; @@ -119,11 +120,12 @@ phm2pj_main (int argc, char* argv[]) 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 @@ -215,8 +217,8 @@ phm2pj_main (int argc, char* argv[]) 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; @@ -228,18 +230,18 @@ phm2pj_main (int argc, char* argv[]) #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;