r117: *** empty log message ***
authorKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 22 Jun 2000 10:17:28 +0000 (10:17 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 22 Jun 2000 10:17:28 +0000 (10:17 +0000)
35 files changed:
ChangeLog
TODO
configure
configure.in
include/backprojectors.h
include/ct.h
include/ctsupport.h
include/filter.h
include/mpiworld.h
include/phantom.h
include/projections.h
include/scanner.h
include/trace.h [new file with mode: 0644]
libctsim/Makefile.am
libctsim/backprojectors.cpp
libctsim/filter.cpp
libctsim/imagefile.cpp
libctsim/options.cpp [deleted file]
libctsim/phantom.cpp
libctsim/projections.cpp
libctsim/reconstr.cpp [deleted file]
libctsupport/Makefile.am
libctsupport/audio.cpp [deleted file]
libctsupport/clip.cpp
libctsupport/consoleio.cpp
libctsupport/filefuncs.cpp
libctsupport/mathfuncs.cpp [new file with mode: 0644]
libctsupport/normangle.cpp [deleted file]
libctsupport/simpson.cpp [deleted file]
src/Makefile.am
src/ctrec.cpp
src/if2img.cpp
src/mpiworld.cpp
src/phm2if.cpp
src/phm2pj.cpp

index dfaae59..58f4e14 100644 (file)
--- 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 (file)
--- 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
 =======================
index f99dc3f..3301d0c 100755 (executable)
--- 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; }
index b48e0f3..4a38772 100644 (file)
@@ -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.
index c28a11b..08b72d4 100644 (file)
@@ -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
 **  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);
index 4dccb37..ead95ad 100644 (file)
@@ -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 <stdint.h>             /* Standard ints on Linux */
 #endif
 
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <iterator>
+#include <algorithm>
+#include <exception>
+#include <stdexcept>
+#include <memory>
+
+using namespace std;
+
 #ifdef HAVE_MPI
 #include "mpi++.h"
 #include "mpiworld.h"
@@ -122,82 +134,19 @@ extern "C" {
   #include "sgp.h"
 #endif
 
-#include <fstream>
-#include <iostream>
-#include <sstream>
-#include <string>
-#include <iterator>
-#include <algorithm>
-#include <exception>
-#include <stdexcept>
-
-using namespace std;
-
 #include "array2d.h"
 #include "fnetorderstream.h"
 #include "imagefile.h"
 #include "phantom.h"
-#include "projections.h"
 #include "scanner.h"
+#include "backprojectors.h"
+#include "filter.h"
+#include "projections.h"
 
 //----------------------------------------------------------------------//
 //                             USER SYMBOLS                            //
 //----------------------------------------------------------------------//
 
-// Filter types
-static const char O_FILTER_ABS_BANDLIMIT_STR[]= "abs_bandlimit";
-static const char O_FILTER_ABS_SINC_STR[]=      "abs_sinc";
-static const char O_FILTER_ABS_COS_STR[]=       "abs_cos";
-static const char O_FILTER_ABS_HAMMING_STR[]=   "abs_hamming";
-static const char O_FILTER_SHEPP_STR[]=         "shepp";
-static const char O_FILTER_BANDLIMIT_STR[]=     "bandlimit";
-static const char O_FILTER_SINC_STR[]=          "sinc";
-static const char O_FILTER_COS_STR[]=           "cos";
-static const char O_FILTER_HAMMING_STR[]=       "hamming";
-static const char O_FILTER_TRIANGLE_STR[]=      "triangle";
-
-typedef enum {            
-  FILTER_BANDLIMIT, 
-  FILTER_SINC,
-  FILTER_G_HAMMING,
-  FILTER_COSINE,
-  FILTER_TRIANGLE,
-  FILTER_ABS_BANDLIMIT,        // filter times |x| 
-  FILTER_ABS_SINC, 
-  FILTER_ABS_G_HAMMING,
-  FILTER_ABS_COSINE,
-  FILTER_SHEPP
-} FilterType;
-
-
-// Function domains 
-static const char D_FREQ_STR[]=    "freq";
-static const char D_SPATIAL_STR[]= "spatial";
-typedef enum {
-  D_FREQ = 1,
-  D_SPATIAL 
-} DomainType;
-
-
-/* interpolation methods */
-static const char O_INTERP_NEAREST_STR[]=  "nearest";
-static const char O_INTERP_LINEAR_STR[]=   "linear";
-static const char O_INTERP_BSPLINE_STR[]=  "bspline";
-
-#undef HAVE_BSPLINE_INTERP
-typedef enum {          // Interpolation methods
-    I_NEAREST = 1,      // Nearest neighbor
-#if HAVE_BSPLINE_INTERP
-    I_BSPLINE,
-    I_1BSPLINE,      // 1st order B-Spline 
-    I_2BSPLINE,
-    I_3BSPLINE,
-#endif
-    I_LINEAR        // Linear interpolation 
-} InterpolationType;
-
-
 // Trace levels
 static const char O_TRACE_NONE_STR[]=     "none";
 static const char O_TRACE_TEXT_STR[]=     "text";
@@ -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
 
index 3cc00f6..f180091 100644 (file)
@@ -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);
 
 
index 274d792..b544a1f 100644 (file)
@@ -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
index 2889821..5818ed3 100644 (file)
@@ -27,7 +27,7 @@
 
 #include <mpi++.h>
 #include <vector.h>
-
+#include <string>
 
 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;
index 32eedf6..b3765c6 100644 (file)
@@ -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<PhantomElement*> m_listPElem;      // pelem lists
     string m_name;
+    PhantomID m_id;
+    bool m_fail;
+
+    // Standard phantomsa
+    static const char PHM_HERMAN_STR[]=    "herman";
+    static const char PHM_BHERMAN_STR[]=   "bherman";
+    static const char PHM_ROWLAND_STR[]=   "rowland";
+    static const char PHM_BROWLAND_STR[]=  "browland";
+    static const char PHM_UNITPULSE_STR[]= "unitpulse";
+    static PhantomID convertNameToPhantomID (const char* const phmName);
+    static const char* convertPhantomIDToName (const PhantomID phmID);
+
+    void init(void);
 };
 
 typedef slist<PhantomElement*>::iterator PElemIterator;
index 104476d..a9970b3 100644 (file)
@@ -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;}
index 005c802..d383498 100644 (file)
@@ -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 (file)
index 0000000..d302d3c
--- /dev/null
@@ -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
+
index 70b8fa5..b932178 100644 (file)
@@ -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
index ef0de98..9e35092 100644 (file)
@@ -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
 
 #include "ct.h"
 
+Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+{
+  m_fail = false;
+  m_pBackprojectImplem = NULL;
+
+  initBackprojector (proj, im, backprojName, interpName);
+}
+
+void 
+Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
+{
+  if (m_pBackprojectImplem)
+    m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
+}
+
+Backprojector::~Backprojector (void)
+{
+  delete m_pBackprojectImplem;
+}
 
 // FUNCTION IDENTIFICATION
 //     Backproject* projector = selectBackprojector (...)
 //     Selects a backprojector based on BackprojType 
 //     and initializes the backprojector
 
-Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
+bool
+Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
 {
-    Backproject* bj = NULL;
-
-    if (bjType == O_BPROJ_TRIG)
-       bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
-    else if (bjType == O_BPROJ_TABLE)
-       bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
-    else if (bjType == O_BPROJ_DIFF)
-       bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
-    else if (bjType == O_BPROJ_DIFF2)
-       bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
-    else if (bjType == O_BPROJ_IDIFF2)
-       bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
-    else 
-      sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
-
-    return (bj);
+  m_nameBackproject = backprojName;
+  m_nameInterpolation = interpName;
+  m_idBackproject = convertBackprojectNameToID (backprojName);
+  m_idInterpolation = convertInterpolationNameToID (interpName);
+  m_pBackprojectImplem = NULL;
+
+  if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
+    m_fail = true;
+    return false;
+  }
+
+  if (m_idBackproject == BPROJ_TRIG)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_TABLE)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_DIFF)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_DIFF2)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_IDIFF2)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+  else {
+    m_fail = true;
+    return false;
+  }
+
+  return true;
+}
+
+
+const Backprojector::BackprojectID
+Backprojector::convertBackprojectNameToID (const char* const backprojName)
+{
+  BackprojectID backprojID = BPROJ_INVALID;
+
+  if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
+    backprojID = BPROJ_TRIG;
+  else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
+    backprojID = BPROJ_TABLE;
+  else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
+    backprojID = BPROJ_DIFF;
+  else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
+    backprojID = BPROJ_DIFF2;
+  else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
+    backprojID = BPROJ_IDIFF2;
+
+  return (backprojID);
+}
+
+const char*
+Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
+{
+  const char *bprojName = "";
+
+  if (bprojID == BPROJ_TRIG)
+    bprojName = BPROJ_TRIG_STR;
+  else if (bprojID == BPROJ_TABLE)
+    bprojName = BPROJ_TABLE_STR;
+  else if (bprojID == BPROJ_DIFF)
+    bprojName = BPROJ_DIFF_STR;
+  else if (bprojID == BPROJ_DIFF2)
+    bprojName = BPROJ_DIFF2_STR;
+  else if (bprojID == BPROJ_IDIFF2)
+    bprojName = BPROJ_IDIFF2_STR;
+
+  return (bprojName);
+}
+
+
+
+const Backprojector::InterpolationID
+Backprojector::convertInterpolationNameToID (const char* const interpName)
+{
+  InterpolationID interpID = INTERP_INVALID;
+
+  if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
+    interpID = INTERP_NEAREST;
+  else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
+    interpID = INTERP_LINEAR;
+#if HAVE_BSPLINE_INTERP
+  else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
+    interpID = INTERP_BSPLINE;
+#endif
+    
+  return (interpID);
+}
+
+
+/* NAME
+ *     name_of_interp                  Return name of interpolation method
+ *
+ * SYNOPSIS
+ *     name = name_of_interp (interp_type)
+ *     char *name                      Name of interpolation method
+ *     int interp_type                 Method of interpolation
+ *
+ * NOTES
+ *     Returns NULL if interp_type is invalid
+ */
+
+const char*
+Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
+{
+  if (interpID == INTERP_NEAREST)
+    return (INTERP_NEAREST_STR);
+  else if (interpID == INTERP_LINEAR)
+    return (INTERP_LINEAR_STR);
+#if HAVE_BSPLINE_INTERP
+  else if (interpID == INTERP_BSPLINE)
+    return (INTERP_BSPLINE_STR);
+#endif
+  else
+    return ("");
 }
 
 
@@ -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<int> (L / detInc); // calc'd index in the filter raysum array
 
        if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
            errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
        else
          v[ix][iy] += rotInc * filteredProj[iDetPos];
-      } else if (interpType == I_LINEAR) {
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
          double p = L / detInc;        // position along detector
          double pFloor = floor (p);
          int iDetPos = iDetCenter + static_cast<int>(pFloor);
@@ -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<int>(L / detInc);    // calc index in the filtered raysum vector 
 
        if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
          errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
        else
          pImCol[iy] += filteredProj[iDetPos];
-      } else if (interpType == I_LINEAR) {
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
        double dPos = L / detInc;               // position along detector 
        double dPosFloor = floor (dPos);
        int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
@@ -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<int>(curDetPos / detInc);    // calc index in the filtered raysum vector 
 
        if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
            errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
        else
          pImCol[iy] += filteredProj[iDetPos];
-      } else if (interpType == I_LINEAR) {
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
        double detPos = curDetPos / detInc;             // position along detector 
        double detPosFloor = floor (detPos);
        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
@@ -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<int>(L))]);
 #endif
-      if (interpType == I_NEAREST) {
+      if (interpType == Backprojector::INTERP_NEAREST) {
        int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector 
        
        if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
            errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
        else
          *pImCol++ += filteredProj[iDetPos];
-      } else if (interpType == I_LINEAR) {
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
        double detPosFloor = floor (curDetPos);
        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
        double frac = curDetPos - detPosFloor;  // fraction distance from det 
@@ -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) {
index b2d3062..c62eb24 100644 (file)
@@ -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
 
 
 /* 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);
 }
index 4cae5d2..bcb07e2 100644 (file)
@@ -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<int>((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 (file)
index 9bbb585..0000000
+++ /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<InterpolationType>(-1);
-  }
-    
-  return (opt);
-}
-
-FilterType 
-opt_set_filter (const char *optarg)
-{
-  FilterType opt;
-
-  if (strcmp(optarg, O_FILTER_BANDLIMIT_STR) == 0)
-    opt = FILTER_BANDLIMIT;
-  else if (strcmp(optarg, O_FILTER_HAMMING_STR) == 0)
-    opt = FILTER_G_HAMMING;
-  else if (strcmp(optarg, O_FILTER_SINC_STR) == 0)
-    opt = FILTER_SINC;
-  else if (strcmp(optarg, O_FILTER_COS_STR) == 0)
-    opt = FILTER_COSINE;
-  else if (strcmp(optarg, O_FILTER_TRIANGLE_STR) == 0)
-    opt = FILTER_TRIANGLE;
-  else if (strcmp(optarg, O_FILTER_ABS_BANDLIMIT_STR) == 0)
-    opt = FILTER_ABS_BANDLIMIT;
-  else if (strcmp(optarg, O_FILTER_ABS_HAMMING_STR) == 0)
-    opt = FILTER_ABS_G_HAMMING;
-  else if (strcmp(optarg, O_FILTER_ABS_SINC_STR) == 0)
-    opt = FILTER_ABS_SINC;
-  else if (strcmp(optarg, O_FILTER_ABS_COS_STR) == 0)
-    opt = FILTER_ABS_COSINE;
-  else if (strcmp(optarg, O_FILTER_SHEPP_STR) == 0)
-    opt = FILTER_SHEPP;
-  else {
-    sys_error(ERR_WARNING, "Invalid filter type %s\n", optarg);
-    opt = static_cast<FilterType>(-1);
-  }
-
-  return (opt);
-}
-
-const char *
-name_of_filter (const int filter)
-{
-  const char *name = "Unknown filter";
-
-  if (filter == FILTER_SHEPP)
-    name = O_FILTER_SHEPP_STR;
-  else if (filter == FILTER_ABS_COSINE)
-    name = O_FILTER_ABS_COS_STR;
-  else if (filter == FILTER_ABS_SINC)
-    name = O_FILTER_ABS_SINC_STR;
-  else if (filter == FILTER_ABS_G_HAMMING)
-    name = O_FILTER_ABS_HAMMING_STR;
-  else if (filter == FILTER_ABS_BANDLIMIT)
-    name = O_FILTER_ABS_BANDLIMIT_STR;
-  else if (filter == FILTER_COSINE)
-    name = O_FILTER_COS_STR;
-  else if (filter == FILTER_SINC)
-    name = O_FILTER_SINC_STR;
-  else if (filter == FILTER_G_HAMMING)
-    name = O_FILTER_HAMMING_STR;
-  else if (filter == FILTER_BANDLIMIT)
-    name = O_FILTER_BANDLIMIT_STR;
-  else if (filter == FILTER_TRIANGLE)
-    name = O_FILTER_TRIANGLE_STR;
-           
-  return (name);
-}
-      
-DomainType
-opt_set_filter_domain (const char *optarg)
-{
-  DomainType opt;
-
-  if (strcmp(optarg, D_SPATIAL_STR) == 0)
-    opt = D_SPATIAL;
-  else if (strcmp(optarg, D_FREQ_STR) == 0)
-    opt = D_FREQ;
-  else {
-    sys_error(ERR_WARNING, "Invalid filter domain %s\n", optarg);
-    opt = static_cast<DomainType>(-1);
-  }
-
-  return (opt);
-}
-
-const char *
-name_of_filter_domain (const DomainType domain)
-{
-  const char *name = "Unknown domain";
-
-  if (domain == D_SPATIAL)
-    return(D_SPATIAL_STR);
-  else if (domain == D_FREQ)
-    return(D_FREQ_STR);
-
-  return (name);
-}
-
-
-BackprojType
-opt_set_backproj (const char *optarg)
-{
-  BackprojType opt;
-
-  if (strcmp(optarg, O_BPROJ_TRIG_STR) == 0)
-    opt = O_BPROJ_TRIG;
-  else if (strcmp(optarg, O_BPROJ_TABLE_STR) == 0)
-    opt = O_BPROJ_TABLE;
-  else if (strcmp(optarg, O_BPROJ_DIFF_STR) == 0)
-    opt = O_BPROJ_DIFF;
-  else if (strcmp(optarg, O_BPROJ_DIFF2_STR) == 0)
-    opt = O_BPROJ_DIFF2;
-  else if (strcmp(optarg, O_BPROJ_IDIFF2_STR) == 0)
-    opt = O_BPROJ_IDIFF2;
-  else {
-    sys_error(ERR_WARNING, "Invalid backprojection method %s\n", optarg);
-    opt = static_cast<BackprojType>(-1);
-  }
-
-  return (opt);
-}
-
-const char *
-name_of_backproj(const BackprojType bproj)
-{
-  const char *name = "Unknown backprojection method";
-
-  if (bproj == O_BPROJ_TRIG)
-    name = O_BPROJ_TRIG_STR;
-  else if (bproj == O_BPROJ_TABLE)
-    name = O_BPROJ_TABLE_STR;
-  else if (bproj == O_BPROJ_DIFF)
-    name = O_BPROJ_DIFF_STR;
-  else if (bproj == O_BPROJ_DIFF2)
-    name = O_BPROJ_DIFF2_STR;
-  else if (bproj == O_BPROJ_IDIFF2)
-    name = O_BPROJ_IDIFF2_STR;
-
-  return (name);
-}
-
-
-
-/* NAME
- *     name_of_interp                  Return name of interpolation method
- *
- * SYNOPSIS
- *     name = name_of_interp (interp_type)
- *     char *name                      Name of interpolation method
- *     int interp_type                 Method of interpolation
- *
- * NOTES
- *     Returns NULL if interp_type is invalid
- */
-
-const char *
-name_of_interpolation (int interp_type)
-{
-  if (interp_type == I_NEAREST)
-    return (O_INTERP_NEAREST_STR);
-  else if (interp_type == I_LINEAR)
-    return (O_INTERP_LINEAR_STR);
-#if HAVE_BSPLINE_INTERP
-  else if (interp_type == I_BSPLINE)
-    return (O_INTERP_BSPLINE_STR);
-#endif
-  else
-    return ("Unknown interpolation method");
-}
-
-
index ec5d0dc..c080493 100644 (file)
@@ -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
 //
 
 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<int> (POINTS_PER_CIRCLE * theta / TWOPI + 1.5);
 }
index 54f8f7d..5ce0a8e 100644 (file)
@@ -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="<<filterName<< ", interp="<<interpName<<", backproject="<<backprojectName<<endl;
+
+#if HAVE_SGP
+  SGP_ID gid;
+  double plot_xaxis [n_vec_filter];                    // array for plotting 
+
+  if (trace > TRACE_TEXT)  {
+    int i;
+    double f;
+    double filterInc = (filterMax - filterMin) / (n_vec_filter - 1);
+    for (i = 0, f = filterMin; i < n_vec_filter; i++, f += filterInc)
+      plot_xaxis[i] = f;
+      
+    gid = ezplot (plot_xaxis, filter.getFilter(), n_vec_filter);
+    cio_put_str ("Press any key to continue");
+    cio_kb_getc ();
+    sgp2_close (gid);
+  }
+  if (trace >= TRACE_TEXT) {
+    printf ("nview=%d, ndet=%d, det_start=%.4f, det_inc=%.4f\n", m_nView, m_nDet, m_detStart, m_detInc);
+  }
+#endif  //HAVE_SGP
+
+  Backprojector bj (*this, im, backprojectName, interpName);
+  if (bj.fail())
+    return false;
+
+  for (int iview = 0; iview < m_nView; iview++)  {
+    if (trace >= TRACE_TEXT) 
+      printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
+      
+    DetectorArray& darray = getDetectorArray (iview);
+    DetectorValue* detval = darray.detValues();
+
+    for (int j = 0; j < m_nDet; j++)
+      filtered_proj[j] = filter.convolve (detval, det_inc, j, m_nDet);
+
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT)  {
+      ezset  ("clear.");
+      ezset  ("xticks major 5.");
+      ezset  ("xlabel ");
+      ezset  ("ylabel ");
+      ezset  ("xlength .5.");
+      ezset  ("box.");
+      ezset  ("grid.");
+      ezset  ("ufinish yes.");
+      ezplot (detval, plot_xaxis, m_nDet);
+      ezset  ("clear.");
+      ezset  ("xticks major 5.");
+      ezset  ("xlabel ");
+      ezset  ("ylabel ");
+      ezset  ("ustart yes.");
+      ezset  ("xporigin .5.");
+      ezset  ("xlength .5.");
+      ezset ("box");
+      ezset ("grid");
+      gid = ezplot (filtered_proj, plot_xaxis, m_nDet);
+    }
+#endif  //HAVE_SGP
+
+#ifdef HAVE_BSPLINE_INTERP
+    if (interp_type == I_BSPLINE) 
+       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+    
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT && interp_type == I_BSPLINE) {
+       bspline (m_nDet, zoom_factor, spline_order, filtered_proj, filtered_proj);
+      ezplot_1d (filtered_proj, n_filtered_proj);
+    }
+#endif
+#endif
+
+    bj.BackprojectView (filtered_proj, darray.viewAngle());
+
+#ifdef HAVE_SGP
+    if (trace >= TRACE_PLOT) {
+      char str[256];
+      printf ("Do you want to exit with current pic (y/n) -- ");
+      fgets(str, sizeof(str), stdin);
+      sgp2_close (sgp2_get_active_win());
+      if (tolower(str[0]) == 'y') {
+       break;
+      }
+    } 
+#endif  //HAVE_SGP
+  }
+
+  return true;
+}
+
diff --git a/libctsim/reconstr.cpp b/libctsim/reconstr.cpp
deleted file mode 100644 (file)
index 9c88780..0000000
+++ /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);
-}
-
index df80a70..0082f1f 100644 (file)
@@ -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 (file)
index 2def9d2..0000000
+++ /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
-}
index 9c2abbc..ebf94c2 100644 (file)
@@ -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);
index bc6df45..fd895e8 100644 (file)
@@ -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
+}
index 511251d..ab26634 100644 (file)
@@ -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: <R> - Retry | <N> - New name | <A> - Abort program --> ";
-      c = cio_kb_waitc ("RrNnAa", TRUE);
-      c = tolower(c);
-      cerr << c << endl;
-
-      if (c == 'r')            // Retry -- Nothing to do here 
-       ;                       
-      else if (c == 'a')       // Abort -- Exit to OS
-       exit (1);
-      else if (c == 'n') {     // New name - get from console
-       cerr << "Enter new file name -- ";
-       fgets (fname, sizeof(fname), stdin);
-       str_wrm_tail (fname);
-       cerr << endl;
-      } 
-    }
-  } while (fp == NULL);
-
-  return (fp);
-}
diff --git a/libctsupport/mathfuncs.cpp b/libctsupport/mathfuncs.cpp
new file mode 100644 (file)
index 0000000..bc3b383
--- /dev/null
@@ -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 (file)
index 1263801..0000000
+++ /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 (file)
index e7a5270..0000000
+++ /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);
-}
index 42c90f2..383818f 100644 (file)
@@ -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
index 2d94be9..39688cf 100644 (file)
@@ -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
index c3a4bfe..a230974 100644 (file)
@@ -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);
index f5631f1..3b4b993 100644 (file)
@@ -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;
+}
index 60e504b..75119c9 100644 (file)
@@ -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<FilterType>(-1);
-  DomainType opt_filter_domain = D_SPATIAL;
+  string optPhmName;
+  string optFilterName;
+  string optDomainName = "spatial";
   char *opt_outfile = NULL;
   int opt_debug = 0;
   string opt_desc;
@@ -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)
index fa8d0a9..a4510c4 100644 (file)
@@ -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;