r186: *** empty log message ***
authorKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 31 Aug 2000 08:38:58 +0000 (08:38 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 31 Aug 2000 08:38:58 +0000 (08:38 +0000)
31 files changed:
ChangeLog
acconfig.h
config.h.in
configure
configure.in
html/simulate.html.in
include/backprojectors.h
include/ezplot.h
include/procsignal.h
include/projections.h
include/scanner.h
include/sgp.h
libctgraphics/ezplot.cpp
libctgraphics/sgp.cpp
libctsim/backprojectors.cpp
libctsim/procsignal.cpp
libctsim/projections.cpp
libctsim/scanner.cpp
libctsupport/mathfuncs.cpp
libctsupport/syserror.cpp
src/ctsim.cpp
src/dialogs.cpp
src/dialogs.h
src/dlgprojections.cpp
src/dlgprojections.h
src/views.cpp
src/views.h
tools/Makefile.am
tools/phm2pj.cpp
tools/pj2if.cpp
tools/pjrec.cpp

index bd59494..e77b594 100644 (file)
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,10 @@
+2.0.0-b11 - 8/31/00
+   Added reconstruction for Equilinear and Equiangular geometries
+   Changed theta to be -PI/2 to make compliant with Kak-Slaney formulas
+   ctsim: Added projection graph to animation of projection collection 
+   ctsim: Added single stepping to projection collection animation
+   ctsim: improved File/Properties display for projection files
+       
 2.0.0-b10 - 8/25/00
    ctsim: Added animation of projection collection processs
    ctsim: Added Auto Scaling for image windows
index 9c06cd5..c7dffb5 100644 (file)
@@ -74,3 +74,5 @@
 
 /* Define for HDF5 library */
 #undef HAVE_HDF5
+
+#undef HAVE_USLEEP
index 31cc9fd..97f445e 100644 (file)
@@ -90,6 +90,9 @@
 /* Define if you have the strtol function.  */
 #undef HAVE_STRTOL
 
+/* Define if you have the usleep function.  */
+#undef HAVE_USLEEP
+
 /* Define if you have the <X11/extensions/Xdbe.h> header file.  */
 #undef HAVE_X11_EXTENSIONS_XDBE_H
 
index b46197c..299bf43 100755 (executable)
--- a/configure
+++ b/configure
@@ -712,7 +712,7 @@ fi
 
 PACKAGE=ctsim
 
-VERSION=2.0.0b10
+VERSION=2.0.0b11
 
 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; }
@@ -2614,7 +2614,7 @@ fi
 
 fi
 
-for ac_func in strtod strtol snprintf htonl
+for ac_func in strtod strtol snprintf htonl usleep
 do
 echo $ac_n "checking for $ac_func""... $ac_c" 1>&6
 echo "configure:2621: checking for $ac_func" >&5
index 09d8ec9..0b5e270 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/ctsim.cpp)
-AM_INIT_AUTOMAKE(ctsim,2.0.0b10)
+AM_INIT_AUTOMAKE(ctsim,2.0.0b11)
 AM_CONFIG_HEADER(config.h)
 
 dnl Checks for programs.
@@ -99,7 +99,7 @@ AC_STRUCT_TM
 
 dnl Checks for library functions.
 AC_FUNC_VPRINTF
-AC_CHECK_FUNCS(strtod strtol snprintf htonl)
+AC_CHECK_FUNCS(strtod strtol snprintf htonl usleep)
 AC_CHECK_FUNC(basename)
 AC_CHECK_FUNC(setjmp)
 if test "${OSTYPE}" = "cygwin" ; then
index 5039901..bd0745a 100644 (file)
@@ -28,7 +28,12 @@ MPI Supercomputing:<br>
 <INPUT type="radio" name="MPI" value="no" checked>No (Single CPU)<br>
 </td><td valign="top">
 <h3>Simulate X-Ray acquistion</h3>
-Number of detectors: <input type="text" name="PJ_NDet" size="4" value="367"><p>
+Geometry:<br>
+<INPUT type="radio" name="geometry" value="parallel" checked>Parallel<br>
+<INPUT type="radio" name="geometry" value="equiangular">Equiangular<br>
+<INPUT type="radio" name="geometry" value="collinear">Collinear<br>
+<p>
+Number of Detectors: <input type="text" name="PJ_NDet" size="4" value="367"><p>
 Number of Rotations: <input type="text" name="PJ_NRot" size="4" value="320"><p>
 Number of Rays<br>(samples) per detector: <input type="text" name="PJ_NRay" size="2" value="3"><p>
 Rotation Angle<br>as a multiple of PI: <input type="text" name="PJ_RotAngle" size="3" value="1.0"><br>
index 1781fe1..b36dfb8 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.h,v 1.13 2000/08/25 15:59:13 kevin Exp $
+**  $Id: backprojectors.h,v 1.14 2000/08/31 08:38:58 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
@@ -132,6 +132,7 @@ class Backproject
     double xMin, xMax, yMin, yMax;     // Retangular coords of phantom
     double xInc, yInc; // size of cells
     int m_interpFactor;
+    double m_dFocalLength;
 
  private:
     Backproject (const Backproject& rhs);
@@ -158,7 +159,7 @@ class BackprojectTable : public Backproject
 
   void BackprojectView (const double* const t, double view_angle);
 
- private:
+ protected:
   Array2d<kfloat64> arrayR;
   Array2d<kfloat64> arrayPhi;
   kfloat64** r;
@@ -211,5 +212,25 @@ class BackprojectIntDiff3 : public BackprojectDiff
   void BackprojectView (const double* const t, double view_angle);
 };
 
+class BackprojectEquilinear : public BackprojectTable
+{
+ public:
+  BackprojectEquilinear (const Projections& proj, ImageFile& im, int interpID, const int interpFactor)
+      : BackprojectTable::BackprojectTable (proj, im, interpID, interpFactor)
+      {}
+
+  void BackprojectView (const double* const t, double view_angle);
+};
+
+class BackprojectEquiangular : public BackprojectTable
+{
+ public:
+  BackprojectEquiangular (const Projections& proj, ImageFile& im, int interpID, const int interpFactor)
+      : BackprojectTable::BackprojectTable (proj, im, interpID, interpFactor)
+      {}
+
+  void BackprojectView (const double* const t, double view_angle);
+};
+
 
 #endif
index 22c2360..716fd3e 100644 (file)
@@ -7,7 +7,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ezplot.h,v 1.12 2000/07/29 19:50:08 kevin Exp $
+**  $Id: ezplot.h,v 1.13 2000/08/31 08:38:58 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
@@ -221,6 +221,7 @@ class EZPlot {
     void addCurve (const double* x, const float* y, int num);
     void addCurve (const double* x, const double* y, int num);
     void addCurve (const double* y, int n);
+    void addCurve (const float* y, int n);
 
     void plot ();
 };     
index e8e9679..ad28581 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.h,v 1.3 2000/08/27 20:32:54 kevin Exp $
+**  $Id: procsignal.h,v 1.4 2000/08/31 08:38:58 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
@@ -58,14 +58,11 @@ class ProcessSignal {
     static const int FILTER_GENERATION_DIRECT;
     static const int FILTER_GENERATION_INVERSE_FOURIER;
 
-    ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1, const int iTraceLevel = Trace::TRACE_NONE, int iGeometry = Scanner::GEOMETRY_PARALLEL);
-
-    ProcessSignal (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad = 0, const int iPreinterpolationFactor = 1, const int iTraceLevel = Trace::TRACE_NONE);
+    ProcessSignal (const char* szFilterName, const char* szFilterMethodName,double bw, double signalIncrement, int n, double param, const char* szDomainName, const char* szFilterGenerationName, const int zeropad = 0, const int preinterpolationFactor = 1, const int iTraceLevel = Trace::TRACE_NONE, int iGeometry = Scanner::GEOMETRY_PARALLEL, double dFocalLength = 1.);
 
     ~ProcessSignal();
 
-    void filterSignal (const double input[], double output[], int view =0) const;
-    void filterSignal (const float input[], double output[], int view = 0) const;
+    void filterSignal (const float input[], double output[]) const;
 
     bool fail(void) const      {return m_fail;}
     const string& failMessage(void) const {return m_failMessage;}
@@ -133,6 +130,7 @@ class ProcessSignal {
     int m_nOutputPoints;
     int m_iPreinterpolationFactor;
     int m_idGeometry;
+    double m_dFocalLength;
 
     bool m_fail;
     string m_failMessage;
@@ -151,7 +149,7 @@ class ProcessSignal {
     fftw_plan m_complexPlanForward, m_complexPlanBackward;
 #endif
 
-    void init (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, const int iTraceLevel, const int iGeometry);
+    void init (const int idFilter, int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, const int iTraceLevel, const int iGeometry, double dFocalLength);
 
     // transforms that use precalculated trig tables, therefore don't 
     // require number of data points (n) as an argument
@@ -159,6 +157,7 @@ class ProcessSignal {
     void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int direction) const;
     void finiteFourierTransform (const complex<double> input[], double output[], const int direction) const;
 
+    double convolve (const double func[], const double filter[], const double dx, const int n, const int np) 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;
 
index 38828a1..8a91616 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.h,v 1.12 2000/08/27 20:32:55 kevin Exp $
+**  $Id: projections.h,v 1.13 2000/08/31 08:38:58 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -32,7 +32,7 @@
 class Scanner;
 class DetectorArray;
 class Array2dFileLabel;
-
+class ostringstream;
 
 // Projections
 class Projections
@@ -40,16 +40,18 @@ class Projections
  public:
   Projections (const int nView, const int nDet);
   Projections (const Scanner& scanner);
-  Projections (void);
-  ~Projections (void);
+  Projections ();
+  ~Projections ();
 
   void initFromScanner (const Scanner& scanner);
 
-  void printProjectionData (void);
-  void printScanInfo (void) const;
+  void printProjectionData ();
+  void printScanInfo (ostringstream& os) const;
 
+  bool read (const string& fname);
   bool read (const char* fname);
   bool write (const char* fname);
+  bool write (const string& fname);
   bool detarrayRead (fnetorderstream& fs, DetectorArray& darray, const int view_num);
   bool detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, const int view_num);
 
@@ -61,20 +63,24 @@ class Projections
   void setPhmLen (double phmLen) { m_phmLen = phmLen;}
   void setCalcTime (double calcTime) {m_calcTime = calcTime;}
   void setRemark (const char* remark) {m_remark = remark;}
-  void setRemark (const string remark) {m_remark = remark;}
-
-  const double detStart(void) const {return m_detStart;}
-  const double rotStart(void) const {return m_rotStart;}
-  const double calcTime(void) const {return m_calcTime;}
-  const double phmLen(void) const {return m_phmLen;}
-  const char*  remark(void) const {return m_remark.c_str();}
-  const double detInc(void) const {return m_detInc;}
-  const double rotInc(void) const {return m_rotInc;}
-  const int nDet(void) const {return m_nDet;}
-  const int nView(void) const {return m_nView;}
-  const string& getFilename(void) const {return m_filename;}
-  Array2dFileLabel& getLabel(void) {return m_label;}
-  const Array2dFileLabel& getLabel(void) const {return m_label;}
+  void setRemark (const string& remark) {m_remark = remark;}
+
+  double detStart() const {return m_detStart;}
+  double rotStart() const {return m_rotStart;}
+  double calcTime() const {return m_calcTime;}
+  double phmLen() const {return m_phmLen;}
+  const char*  remark() const {return m_remark.c_str();}
+  double detInc() const {return m_detInc;}
+  double rotInc() const {return m_rotInc;}
+  int nDet() const {return m_nDet;}
+  int nView() const {return m_nView;}
+  int geometry() const {return m_geometry;}
+  double focalLength() const {return m_focalLength;}
+  double fieldOfView() const {return m_fieldOfView;}
+
+  const string& getFilename() const {return m_filename;}
+  Array2dFileLabel& getLabel() {return m_label;}
+  const Array2dFileLabel& getLabel() const {return m_label;}
 
   DetectorArray& getDetectorArray (const int iview)
       { return (*m_projData[iview]); }
@@ -108,12 +114,12 @@ class Projections
 
   const static kuint16 m_signature = ('P'*256 + 'J');
 
-  bool headerRead (void);
-  bool headerWrite (void);
+  bool headerRead ();
+  bool headerWrite ();
   bool headerRead (fnetorderstream& fs);
   bool headerWrite (fnetorderstream& fs);
-  void newProjData (void);
-  void deleteProjData (void);
+  void newProjData ();
+  void deleteProjData ();
 
   void init (const int nView, const int nDet);
 
index fa6ea1d..7f8f95e 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.h,v 1.13 2000/08/27 20:32:55 kevin Exp $
+**  $Id: scanner.h,v 1.14 2000/08/31 08:38:58 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -93,6 +93,7 @@ class Scanner
   double detLen() const {return m_detLen;}
   double focalLength() const {return m_dFocalLength;}
   double fieldOfView() const {return m_dFieldOfView;}
+  int geometry() const {return m_idGeometry;}
 
   static int getGeometryCount() {return s_iGeometryCount;}
   static const char** getGeometryNameArray() {return s_aszGeometryName;}
@@ -119,12 +120,15 @@ class Scanner
   double m_phmLen;              // Maximum Length of phantom or area of interest
   double m_dXCenter;            // Center of Phantom
   double m_dYCenter;            
+  double m_dAngularDetIncrement;
+  double m_dAngularDetLen;
 
   int m_trace;
   struct {
     double xd1,yd1,xd2,yd2;    /* Coordinates of detector endpoints */
     double xs1,ys1,xs2,ys2;    /* Coordinates of source endpoints */
     double angle;              /* Starting angle */
+    double dAngularDet;
   } m_initPos;
 
   GRFMTX_2D m_rotmtxIncrement;
index 0a5fe03..b026cdc 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: sgp.h,v 1.15 2000/08/27 20:32:55 kevin Exp $
+**  $Id: sgp.h,v 1.16 2000/08/31 08:38:58 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
@@ -72,7 +72,7 @@ public:
 #endif
 
   SGPDriver (const char* szWinTitle = "", int xsize = 640, int ysize = 480);
-
+  
   ~SGPDriver ();
 
   int getPhysicalXSize () const
@@ -96,6 +96,9 @@ public:
 #ifdef HAVE_WXWINDOWS
   wxDC* idWX () const
     { return m_pDC; }
+
+  void setDC (wxDC* dc)
+      { m_pDC = dc; }
 #endif
 };
 
@@ -105,7 +108,7 @@ class SGP {
 private:
   int m_iPhysicalXSize;   // Physical Window size 
   int m_iPhysicalYSize;
-  const SGPDriver m_driver;
+  SGPDriver m_driver;
 
   double xw_min;    // Window extents 
   double yw_min;
@@ -199,6 +202,7 @@ public:
   void setMarker (int idMarker, int color);
   void setRasterOp (int ro);
 
+  void getViewport (double &xmin, double &ymin, double& xmax, double& ymax);
   void getTextExtent (const char *szText, double* x, double* y);
   double getCharHeight ();
 
@@ -219,6 +223,10 @@ public:
   void stylusNDC (double x, double y, bool beam);
   void pointNDC (double x, double y);
   void markerNDC (double x, double y);
+
+#if HAVE_WXWINDOWS
+  void setDC (wxDC* pDC);
+#endif
 };
 
 
index 7038132..c295757 100644 (file)
@@ -6,7 +6,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ezplot.cpp,v 1.11 2000/08/22 07:02:48 kevin Exp $
+**  $Id: ezplot.cpp,v 1.12 2000/08/31 08:38:58 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
@@ -71,6 +71,18 @@ EZPlot::addCurve (const double *y, int n)
 }
 
 
+void 
+EZPlot::addCurve (const float *y, int n)
+{
+  double yDouble [n];
+
+  for (int i = 0; i < n; i++)
+    yDouble[i] = y[i];
+
+  addCurve (yDouble, n);
+}
+
+
 void
 EZPlot::addCurve (const float x[], const double y[], int num)
 {
@@ -300,7 +312,8 @@ EZPlot::plot ()
   xtl_wid = x_fldwid * charwidth;              /* calc size of tick labels */
   ytl_wid = y_fldwid * charwidth;
   tl_height = charheight;
-  
+
+  //  rSGP.getViewport (xp_min, yp_min, xp_max, yp_max);
   /* calculate the extent of the plot frame */
   xp_min = o_xporigin;
   yp_min = o_yporigin;
index 52a2430..c1004d6 100644 (file)
@@ -7,7 +7,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: sgp.cpp,v 1.11 2000/08/27 20:32:55 kevin Exp $
+**  $Id: sgp.cpp,v 1.12 2000/08/31 08:38:58 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
@@ -53,9 +53,8 @@ int SGP::s_iRGBColorCount = sizeof(s_aRGBColor) / sizeof(class RGBColor);
 
 #ifdef HAVE_WXWINDOWS
 SGPDriver::SGPDriver (wxDC* pDC, int xsize = 640, int ysize = 480)
-  : m_iPhysicalXSize(xsize), m_iPhysicalYSize(ysize), m_idDriver(0)
+  : m_iPhysicalXSize(xsize), m_iPhysicalYSize(ysize), m_idDriver(0), m_pDC(pDC)
 {
-  m_pDC = pDC;
   m_idDriver |= SGPDRIVER_WXWINDOWS;
 }
 #endif
@@ -103,7 +102,6 @@ SGP::SGP (const SGPDriver& driver)
   setTextAngle (0.);
   setTextSize (1. / 25.);
   setColor (C_BLACK);
-
 }
 
 
@@ -199,6 +197,15 @@ SGP::setViewport (double xmin, double ymin, double xmax, double ymax)
   viewNDC[3] = ymax;
 }
 
+void
+SGP::getViewport (double& xmin, double& ymin, double& xmax, double& ymax)
+{
+    xmin = xv_min;
+    ymin = yv_min;
+    xmax = xv_max;
+    ymax = yv_max;
+}
+
 
 // NAME
 //     frameViewport           draw box around viewport
@@ -566,47 +573,46 @@ SGP::drawRect (double xmin, double ymin, double xmax, double ymax)
 void 
 SGP::drawCircle (const double r)
 {
-       drawArc (r, 0.0, 7.0);
+       drawArc (r, 0.0, TWOPI);
 }
 
-// =============================================================
-// draw arc around current center.  pass angles and radius     
+//==============================================================
+// draw arc around current center.  angles in radius   
 //==============================================================
 
 void 
 SGP::drawArc (const double r, double start, double stop)
 {
-  if ((stop-start) > 2 * PI)
-    stop = start + 2 * PI;
-  if ((start-stop) > 2 * PI)
-    stop = start + 2 * PI;
-  while (start >= stop)
-    stop += 2*PI;
+  if (start > stop) {
+    double temp = start;
+    start = stop;
+    stop = temp;
+  }
 
   double x = r * cos ((double) start);
   double y = r * sin ((double) start);
   moveRel (x, y);          // move from center to start of arc 
 
-  double theta = 5 * PI / 180;
-  double c = cos(theta);
-  double s = sin(theta);
+  const double thetaIncrement = (5 * (TWOPI / 360)); 
+  double cosTheta = cos(thetaIncrement);
+  double sinTheta = sin(thetaIncrement);
 
   double angle, xp, yp;
-  for (angle = start; angle < stop - theta; angle += theta) {
-    xp = c * x - s * y;
-    yp = s * x + c * y;
-    lineRel (xp - x, yp - y);
+  for (angle = start; angle < stop - thetaIncrement; angle += thetaIncrement) {
+    xp = cosTheta * x - sinTheta * y; // translate point by thetaIncrement
+    yp = sinTheta * x + cosTheta * y;
+    lineAbs (xp, yp);
     x = xp; y = yp;
   }
 
-  c = cos (stop - angle);
-  s = sin (stop - angle);
+  double c = cos (stop - angle);
+  double s = sin (stop - angle);
   xp = c * x - s * y;
   yp = s * x + c * y;
-  lineRel (xp - x, yp - y);
+  lineAbs (xp, yp);
 
-  x = r * cos ((double) stop);
-  y = r * sin ((double) stop);
+  x = r * cos (stop);
+  y = r * sin (stop);
   moveRel (-x, -y);            // move back to center of circle 
 }
 
@@ -783,3 +789,13 @@ const unsigned char SGP::MARKER_BITMAP[MARK_COUNT][5] =
     {'\076', '\042', '\042', '\042', '\076'},    // big open square 
     {'\010', '\024', '\042', '\024', '\010'},    // big open diamond 
 };
+
+
+#if HAVE_WXWINDOWS
+void
+SGP::setDC (wxDC* pDC)
+{
+  if (m_driver.isWX())
+    m_driver.setDC(pDC);
+}
+#endif
index 71aa8d4..0bcd407 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.12 2000/08/25 15:59:13 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.13 2000/08/31 08:38:58 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
@@ -148,22 +148,27 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
     return false;
   }
 
-  if (m_idBackproject == BPROJ_TRIG)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
-  else if (m_idBackproject == BPROJ_TABLE)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
-  else if (m_idBackproject == BPROJ_DIFF)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
-  else if (m_idBackproject == BPROJ_DIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
-  else if (m_idBackproject == BPROJ_IDIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
-  else if (m_idBackproject == BPROJ_IDIFF3)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
-  else {
-    m_fail = true;
-    m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
-    return false;
+  if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
+  else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) 
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
+  else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
+      if (m_idBackproject == BPROJ_TRIG)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
+      else if (m_idBackproject == BPROJ_TABLE)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
+      else if (m_idBackproject == BPROJ_DIFF)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
+      else if (m_idBackproject == BPROJ_DIFF2)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
+      else if (m_idBackproject == BPROJ_IDIFF2)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
+      else if (m_idBackproject == BPROJ_IDIFF3)
+         m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
+  } else {
+      m_fail = true;
+      m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
+      return false;
   }
 
   return true;
@@ -256,7 +261,7 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const int inte
 {
   detInc = proj.detInc();
   nDet = proj.nDet();
-  iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection 
+  iDetCenter = nDet / 2;       // index refering to L=0 projection 
   rotInc = proj.rotInc();
 
   v = im.getArray();
@@ -271,6 +276,8 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const int inte
 
   xInc = (xMax - xMin) / nx;   // size of cells
   yInc = (yMax - yMin) / ny;
+
+  m_dFocalLength = proj.focalLength();
 }
 
 Backproject::~Backproject ()
@@ -312,12 +319,12 @@ void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, doubl
 void
 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = HALFPI + view_angle;   // Add PI/2 to get perpendicular angle to detector     
-  int ix, iy;
-  double x, y;                 // Rectang coords of center of pixel 
+  double theta = view_angle;
 
-  for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
-    for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+  double x = xMin + xInc / 2;  // Rectang coords of center of pixel 
+  for (int ix = 0; ix < nx; x += xInc, ix++) {
+    double y = yMin + yInc / 2;
+    for (int iy = 0; iy < ny; y += yInc, iy++) {
       double r = sqrt (x * x + y * y);   // distance of cell from center
       double phi = atan2 (y, x);         // angle of cell from center
       double L = r * cos (theta - phi);  // position on detector
@@ -340,6 +347,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
            v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }
+  }
 }  
 
 
@@ -374,7 +382,7 @@ BackprojectTable::~BackprojectTable ()
 void
 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = HALFPI + view_angle;  // add half PI to view angle to get perpendicular theta angle
+  double theta = view_angle;
 
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
@@ -431,9 +439,9 @@ BackprojectDiff::~BackprojectDiff()
 void
 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
-  double det_dx = xInc * sin (theta);
-  double det_dy = yInc * cos (theta);
+  double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
+  double det_dx = xInc * cos (theta);
+  double det_dy = yInc * sin (theta);
   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
        
   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
@@ -475,11 +483,11 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
 void
 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
+  double theta = view_angle;
 
   // Distance betw. detectors for an angle given in units of detectors 
-  double det_dx = xInc * sin (theta) / detInc;
-  double det_dy = yInc * cos (theta) / detInc;
+  double det_dx = xInc * cos (theta) / detInc;
+  double det_dy = yInc * sin (theta) / detInc;
 
   // calculate detPosition for first point in image (ix=0, iy=0) 
   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
@@ -524,14 +532,14 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl
 void
 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
+  double theta = view_angle;
 
   static const kint32 scale = 1 << 16;
   static const double dScale = scale;
   static const kint32 halfScale = scale / 2;
 
-  const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
-  const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+  const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+  const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
 
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
@@ -576,15 +584,15 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
 void
 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
 {
-  double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
+  double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
   static const int scaleShift = 16;
   static const kint32 scale = (1 << scaleShift);
   static const kint32 scaleBitmask = scale - 1;
   static const kint32 halfScale = scale / 2;
   static const double dInvScale = 1. / scale;
 
-  const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
-  const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+  const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+  const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
 
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
@@ -624,3 +632,82 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
     } //end linear
   } // end for ix
 }
+
+
+void
+BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+  double beta = view_angle;
+
+  for (int ix = 0; ix < nx; ix++) {
+    ImageFileColumn pImCol = v[ix];
+
+    for (int iy = 0; iy < ny; iy++) { 
+      double dAngleDiff = beta - phi[ix][iy];
+      double rcos_t = r[ix][iy] * cos (dAngleDiff);
+      double rsin_t = r[ix][iy] * sin (dAngleDiff);
+      double dFLPlusSin = m_dFocalLength + rsin_t;
+      double gamma =  atan (rcos_t / dFLPlusSin);
+      double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
+
+      if (interpType == Backprojector::INTERP_NEAREST) {
+       int iDetPos =iDetCenter + nearest<int>(gamma / detInc); // calc index in the filtered raysum vector 
+
+       if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos  
+         ; //    errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+       } else
+         pImCol[iy] += filteredProj[iDetPos] / dL2;
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
+       double dPos = gamma / detInc;           // position along detector 
+       double dPosFloor = floor (dPos);
+       int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+       double frac = dPos - dPosFloor; // fraction distance from det 
+       if (iDetPos < 0 || iDetPos >= nDet - 1) {
+         ; //      errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+       } else
+         pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / dL2;
+      }
+    }  // end for y 
+  }    // end for x 
+}
+
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+  double beta = view_angle;
+
+  for (int ix = 0; ix < nx; ix++) {
+    ImageFileColumn pImCol = v[ix];
+
+    for (int iy = 0; iy < ny; iy++) {
+      double dAngleDiff = beta - phi[ix][iy];
+      double rcos_t = r[ix][iy] * cos (dAngleDiff);
+      double rsin_t = r[ix][iy] * sin (dAngleDiff);
+
+      double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
+      double dDetPos =  rcos_t / dU;
+      // double to scale for imaginary detector that passes through origin
+      // of phantom, see Kak-Slaney Figure 3.22
+      dDetPos *= 2; 
+
+      if (interpType == Backprojector::INTERP_NEAREST) {
+       int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc);      // calc index in the filtered raysum vector 
+
+       if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos 
+         ; ///   errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+       else
+         pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
+       double dPos = dDetPos / detInc;         // position along detector 
+       double dPosFloor = floor (dPos);
+       int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+       double frac = dPos - dPosFloor; // fraction distance from det 
+       if (iDetPos < 0 || iDetPos >= nDet - 1)
+         ; //      errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+       else
+         pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / (dU * dU);
+      }
+    }  // end for y 
+  }    // end for x 
+}
+
index e38fae2..5c10c9e 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: procsignal.cpp,v 1.4 2000/08/27 20:32:55 kevin Exp $
+**  $Id: procsignal.cpp,v 1.5 2000/08/31 08:38:58 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
@@ -77,7 +77,7 @@ const int ProcessSignal::s_iFilterGenerationCount = sizeof(s_aszFilterGeneration
 // CLASS IDENTIFICATION
 //   ProcessSignal
 //
-ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry)
+ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMethodName, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const char* szDomainName, const char* szFilterGenerationName, int iZeropad, int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength)
     : m_adFourierCosTable(NULL), m_adFourierSinTable(NULL), m_adFilter(NULL), m_fail(false)
 {
   m_idFilterMethod = convertFilterMethodNameToID (szFilterMethodName);
@@ -109,18 +109,19 @@ ProcessSignal::ProcessSignal (const char* szFilterName, const char* szFilterMeth
     return;
   }
 
-  init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry);
+  init (m_idFilter, m_idFilterMethod, dBandwidth, dSignalIncrement, nSignalPoints, dFilterParam, m_idDomain, m_idFilterGeneration, iZeropad, iPreinterpolationFactor, iTraceLevel, iGeometry, dFocalLength);
 }
 
 
 void
-ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry)
+ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandwidth, double dSignalIncrement, int nSignalPoints, double dFilterParam, const int idDomain, const int idFilterGeneration, const int iZeropad, const int iPreinterpolationFactor, int iTraceLevel, int iGeometry, double dFocalLength)
 {
   m_idFilter = idFilter;
   m_idDomain = idDomain;
   m_idFilterMethod = idFilterMethod;
   m_idFilterGeneration = idFilterGeneration;
   m_idGeometry = iGeometry;
+  m_dFocalLength = dFocalLength;
 
   if (m_idFilter == SignalFilter::FILTER_INVALID || m_idDomain == SignalFilter::DOMAIN_INVALID || m_idFilterMethod == FILTER_METHOD_INVALID || m_idFilterGeneration == FILTER_GENERATION_INVALID) {
     m_fail = true;
@@ -136,6 +137,13 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
   m_iZeropad = iZeropad;
   m_iPreinterpolationFactor = iPreinterpolationFactor;
 
+  // scale signalInc/BW to signalInc/2 to adjust for imaginary detector
+  // through origin of phantom, see Kak-Slaney Fig 3.22, for Collinear
+  if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+    m_dSignalInc /= 2;
+    m_dBandwidth *= 2;
+  }
+
   if (m_idFilterMethod == FILTER_METHOD_FFT) {
 #if HAVE_FFTW
     m_idFilterMethod = FILTER_METHOD_RFFTW;
@@ -222,10 +230,24 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
            m_adFilter[i] /= m_dSignalInc;
        }
     }
-  } 
+    if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+       for (int i = 0; i < m_nFilterPoints; i++)
+           m_adFilter[i] *= 0.5;
+    } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+       for (int i = 0; i < m_nFilterPoints; i++) {
+         int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
+         double sinScale = sin (iDetFromZero * m_dSignalInc);
+         if (fabs(sinScale) < 1E-7)
+             sinScale = 1;
+         else
+             sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
+         double dScale = 0.5 * sinScale * sinScale;
+         m_adFilter[i] *= dScale;
+       }
+    } // if (geometry)
+  } // if (spatial filtering)
 
-  // Frequency-based filtering
-  else if (m_bFrequencyFiltering) {
+  else if (m_bFrequencyFiltering) {  // Frequency-based filtering
 
     if (m_idFilterGeneration == FILTER_GENERATION_DIRECT) {
       // calculate number of filter points with zeropadding
@@ -256,6 +278,22 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
       SignalFilter filter (m_idFilter, m_dFilterMin, m_dFilterMax, m_nFilterPoints, m_dBandwidth, m_dFilterParam, SignalFilter::DOMAIN_FREQUENCY);
       m_adFilter = new double [m_nFilterPoints];
       filter.copyFilterData (m_adFilter, 0, m_nFilterPoints);
+
+      if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+       for (int i = 0; i < m_nFilterPoints; i++)
+         m_adFilter[i] *= 0.5;
+      } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+       for (int i = 0; i < m_nFilterPoints; i++) {
+         int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
+         double sinScale = sin (iDetFromZero * m_dSignalInc);
+         if (fabs(sinScale) < 1E-7)
+           sinScale = 1;
+         else
+           sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
+         double dScale = 0.5 * sinScale * sinScale;
+         //      m_adFilter[i] *= dScale;
+       }
+      }
       if (m_traceLevel >= Trace::TRACE_PLOT) {
        SGPDriver sgpDriver ("Frequency Filter: Natural Order");
        SGP sgp (sgpDriver);
@@ -311,6 +349,21 @@ ProcessSignal::init (const int idFilter, const int idFilterMethod, double dBandw
        cio_put_str ("Press any key to continue");
        cio_kb_getc ();
       }
+      if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+       for (int i = 0; i < m_nFilterPoints; i++)
+         adSpatialFilter[i] *= 0.5;
+      } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+       for (int i = 0; i < m_nFilterPoints; i++) {
+         int iDetFromZero = i - ((m_nFilterPoints - 1) / 2);
+         double sinScale = sin (iDetFromZero * m_dSignalInc);
+         if (fabs(sinScale) < 1E-7)
+           sinScale = 1;
+         else
+           sinScale = (iDetFromZero * m_dSignalInc) / sinScale;
+         double dScale = 0.5 * sinScale * sinScale;
+         adSpatialFilter[i] *= dScale;
+       }
+      }
       for (int i = nSpatialPoints; i < m_nFilterPoints; i++)
        adSpatialFilter[i] = 0;
 
@@ -470,11 +523,26 @@ ProcessSignal::convertFilterGenerationIDToTitle (const int fgID)
 }
 
 void
-ProcessSignal::filterSignal (const float input[], double output[], int iView) const
+ProcessSignal::filterSignal (const float constInput[], double output[]) const
 {
+  double input [m_nSignalPoints];
+  for (int i = 0; i < m_nSignalPoints; i++)
+    input[i] = constInput[i];
+
+  if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+    for (int i = 0; i < m_nSignalPoints; i++) {
+      int iDetFromCenter = i - (m_nSignalPoints / 2);
+      input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
+    }
+  } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+    for (int i = 0; i < m_nSignalPoints; i++) {
+      int iDetFromCenter = i - (m_nSignalPoints / 2);
+      input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
+    }
+  }
   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
-    for (int i = 0; i < m_nSignalPoints; i++)
-      output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
+      for (int i = 0; i < m_nSignalPoints; i++)
+       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
     double inputSignal[m_nFilterPoints];
     for (int i = 0; i < m_nSignalPoints; i++)
index da298f0..be8e829 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.23 2000/08/27 20:32:55 kevin Exp $
+**  $Id: projections.cpp,v 1.24 2000/08/31 08:38:58 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
@@ -80,10 +80,17 @@ Projections::initFromScanner (const Scanner& scanner)
   m_phmLen = scanner.phmLen();
   m_rotInc = scanner.rotInc();
   m_detInc = scanner.detInc();
+  m_geometry = scanner.geometry();
   m_focalLength = scanner.focalLength();
   m_fieldOfView = scanner.fieldOfView();
   m_rotStart = 0;
-  m_detStart =  -(scanner.detLen() / 2) + (scanner.detInc() / 2);
+  m_detStart =  -(scanner.detLen() / 2);
+#if 0
+  if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) {
+    m_detInc /= 2;
+    cout << "Kludge: detInc /= 2 in Projections::initFromScanner" << endl;
+  }
+#endif
 }
 
 void
@@ -286,6 +293,13 @@ Projections::headerRead (fnetorderstream& fs)
   return true;
 }
 
+bool
+Projections::read (const string& filename)
+{
+  return read (filename.c_str());
+}
+
+
 bool
 Projections::read (const char* filename)
 {
@@ -311,6 +325,12 @@ Projections::read (const char* filename)
 }
 
 
+bool
+Projections::write (const string& filename)
+{
+  return write (filename.c_str());
+}
+
 bool
 Projections::write (const char* filename)
 {
@@ -442,14 +462,16 @@ Projections::detarrayWrite (fnetorderstream& fs, const DetectorArray& darray, co
 void
 Projections::printProjectionData (void)
 {
-  printf("Projections Print\n\n");
+  printf("Projections Data\n\n");
   printf("Description: %s\n", m_remark.c_str());
-  printf("nView = %8d          nDet = %8d\n", m_nView, m_nDet);
-  printf("rotStart = %8.4f     rotInc = %8.4f\n", m_rotStart, m_rotInc);
-  printf("detStart = %8.4f     detInc = %8.4f\n", m_detStart, m_detInc);
+  printf("Geometry: %s\n", Scanner::convertGeometryIDToName (m_geometry));
+  printf("nView       = %8d           nDet = %8d\n", m_nView, m_nDet);
   printf("focalLength = %8.4f  fieldOfView = %8.4f\n", m_focalLength, m_fieldOfView);
+  printf("rotStart    = %8.4f       rotInc = %8.4f\n", m_rotStart, m_rotInc);
+  printf("detStart    = %8.4f       detInc = %8.4f\n", m_detStart, m_detInc);
   if (m_projData != NULL) {
       for (int ir = 0; ir < m_nView; ir++) {
+         printf("View %d: angle %f\n", ir, m_projData[ir]->viewAngle());
          DetectorValue* detval = m_projData[ir]->detValues();
          for (int id = 0; id < m_projData[ir]->nDet(); id++)
              printf("%8.4f  ", detval[id]);
@@ -459,16 +481,19 @@ Projections::printProjectionData (void)
 }
 
 void 
-Projections::printScanInfo (void) const
+Projections::printScanInfo (ostringstream& os) const
 {
-  printf ("Number of detectors: %d\n", m_nDet);
-  printf ("    Number of views: %d\n", m_nView);
-  printf ("             Remark: %s\n", m_remark.c_str());
-  printf ("             phmLen: %f\n", m_phmLen);
-  printf ("           detStart: %f\n", m_detStart);
-  printf ("             detInc: %f\n", m_detInc);
-  printf ("           rotStart: %f\n", m_rotStart);
-  printf ("             rotInc: %f\n", m_rotInc);
+  os << "Number of detectors: " << m_nDet << "\n";
+  os << "    Number of views: " << m_nView<< "\n";
+  os << "             Remark: " << m_remark.c_str()<< "\n";
+  os << "           Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
+  os << "       Focal Length: " << m_focalLength<< "\n";
+  os << "      Field Of View: " << m_fieldOfView<< "\n";
+  os << "             phmLen: " << m_phmLen<< "\n";
+  os << "           detStart: " << m_detStart<< "\n";
+  os << "             detInc: " << m_detInc<< "\n";
+  os << "           rotStart: " << m_rotStart<< "\n";
+  os << "             rotInc: " << m_rotInc<< "\n";
 }
 
 
@@ -512,7 +537,7 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 #endif
 
   double filterBW = 1. / detInc;
-  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace, m_geometry);
+  ProcessSignal processSignal (filterName, filterMethodName, filterBW, m_detInc, m_nDet, filt_param, "spatial", filterGenerationName, zeropad, interpFactor, trace, m_geometry, m_focalLength);
 
   if (processSignal.fail()) {
       sys_error (ERR_SEVERE, "%s [Projections::reconstruct]", processSignal.failMessage().c_str());
@@ -558,12 +583,12 @@ Projections::reconstruct (ImageFile& im, const char* const filterName, double fi
 
   for (int iview = 0; iview < m_nView; iview++)  {
     if (trace >= Trace::TRACE_CONSOLE) 
-      printf ("Reconstructing view %d (last = %d)\n",  iview, m_nView - 1);
+      cout <<"Reconstructing view " << iview << "(last = " << m_nView - 1 << ")\n";
       
     const DetectorArray& darray = getDetectorArray (iview);
     const DetectorValue* detval = darray.detValues();
 
-    processSignal.filterSignal (detval, filteredProj, iview);
+    processSignal.filterSignal (detval, filteredProj);
 
 #ifdef HAVE_BSPLINE_INTERP
     if (interp_type == I_BSPLINE) 
index c45d4fd..ef72c5b 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.12 2000/08/27 20:32:55 kevin Exp $
+**  $Id: scanner.cpp,v 1.13 2000/08/31 08:38:58 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
 
 const int Scanner::GEOMETRY_INVALID = -1;
 const int Scanner::GEOMETRY_PARALLEL = 0;
-const int Scanner::GEOMETRY_EQUILINEAR = 1;
-const int Scanner::GEOMETRY_EQUIANGULAR = 2;
+const int Scanner::GEOMETRY_EQUIANGULAR = 1;
+const int Scanner::GEOMETRY_EQUILINEAR = 2;
 
 const char* Scanner::s_aszGeometryName[] = 
 {
   {"parallel"},
-  {"equilinear"},
   {"equiangular"},
+  {"equilinear"},
 };
 
 const char* Scanner::s_aszGeometryTitle[] = 
 {
   {"Parallel"},
-  {"Equilinear"},
   {"Equiangular"},
+  {"Equilinear"},
 };
 
 const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
@@ -101,8 +101,6 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_nSample = 1;
   if (nDet < 1)
     nDet = 1;
-  //  if ((nDet % 2) == 0)
-  //    ++nDet;                // ensure odd number of detectors
 
   m_nDet     = nDet;
   m_nView    = nView;
@@ -114,23 +112,26 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
 
   m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
   m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
+  m_rotLen  = rot_anglen;
+  m_rotInc  = m_rotLen / m_nView;
   if (m_idGeometry == GEOMETRY_PARALLEL) {
     m_detLen   = m_dFieldOfView;
     m_detInc  = m_detLen / m_nDet;
-    m_rotLen  = rot_anglen;
-    m_rotInc  = m_rotLen / m_nView;
     
     double dHalfDetLen = m_detLen / 2;
     m_initPos.xs1 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter - dHalfDetLen;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
+    m_initPos.ys1 = m_dYCenter + dHalfDetLen;
+    m_initPos.xs2 = m_dXCenter + m_dFocalLength;
     m_initPos.ys2 = m_dYCenter + dHalfDetLen;
-    m_initPos.xd1 = m_dXCenter + m_dFocalLength;
+    m_initPos.xd1 = m_dXCenter - m_dFocalLength;
     m_initPos.yd1 = m_dYCenter - dHalfDetLen;
     m_initPos.xd2 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd2 = m_dYCenter + dHalfDetLen;
+    m_initPos.yd2 = m_dYCenter - dHalfDetLen;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
+#if 0
+    double dAngle = (m_dFieldOfView / 2) / cos (asin (m_dFieldOfView / 2 / m_dFocalLength));
+#else
     double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
     if (dFocalPastPhm <= 0.) {
@@ -139,23 +140,26 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
       return;
     }
     double dAngle = atan( dHalfSquare / dFocalPastPhm );
+#endif
     double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
 
     m_detLen = dHalfDetLen * 2;
     m_detInc  = m_detLen / m_nDet;
-    m_rotLen  = rot_anglen;
-    m_rotInc  = m_rotLen / m_nView;
     
-    m_initPos.xs1 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys2 = m_dYCenter;
-    m_initPos.xd1 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd1 = m_dYCenter - dHalfDetLen;
-    m_initPos.xd2 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd2 = m_dYCenter + dHalfDetLen;
+    m_initPos.angle = 0.0;
+    m_initPos.xs1 = m_dXCenter;
+    m_initPos.ys1 = m_dYCenter + m_dFocalLength;
+    m_initPos.xs2 = m_dXCenter;
+    m_initPos.ys2 = m_dYCenter + m_dFocalLength;
+    m_initPos.xd1 = m_dXCenter - dHalfDetLen;
+    m_initPos.yd1 = m_dYCenter - m_dFocalLength;
+    m_initPos.xd2 = m_dXCenter + dHalfDetLen;
+    m_initPos.yd2 = m_dYCenter - m_dFocalLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+#if 0
+    double dAngle = atan ((m_dFieldOfView / 2) / m_dFocalLength);
+#else
     double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
     if (dFocalPastPhm <= 0.) {
@@ -163,17 +167,20 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
       m_failMessage = "Focal Point inside of phantom";
       return;
     }
-    double dAngle = 2 * atan( dHalfSquare / dFocalPastPhm );
+    double dAngle =  atan ( dHalfSquare / dFocalPastPhm ); 
+#endif
     m_detLen = 2 * dAngle;
     m_detInc = m_detLen / m_nDet;
-    m_rotLen = rot_anglen;
-    m_rotInc = m_rotLen / m_nView;
 
-    m_initPos.xs1 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys2 = m_dYCenter;
-    m_initPos.angle = -dAngle;
+    m_dAngularDetIncrement = m_detInc * 2; // Angular Position 2x gamma angle
+    m_dAngularDetLen = m_detLen * 2;
+    m_initPos.dAngularDet = -m_dAngularDetLen / 2;
+
+    m_initPos.angle = 0;
+    m_initPos.xs1 = m_dXCenter;
+    m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
+    m_initPos.xs2 = m_dXCenter;
+    m_initPos.ys2 = m_dYCenter + m_dFocalLength;
   }
 
   // Calculate incrementatal rotation matrix 
@@ -260,10 +267,8 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 
-  double xd1=0, yd1=0, xd2=0, yd2=0, dDetAngle=0;
-  if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-      dDetAngle = m_initPos.angle + start_angle;
-  else {
+  double xd1=0, yd1=0, xd2=0, yd2=0;
+  if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
       xd1 = m_initPos.xd1;
       yd1 = m_initPos.yd1;
       xd2 = m_initPos.xd2;
@@ -303,12 +308,12 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
     m_pSGP->eraseWindow ();
     m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
     m_pSGP->setRasterOp (RO_COPY);
-    m_pSGP->setColor (C_BLUE);
+    m_pSGP->setColor (C_RED);
     m_pSGP->moveAbs (0., 0.);
     m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen);
     m_pSGP->moveAbs (0., 0.);
     m_pSGP->drawCircle (m_dFocalLength);
-    m_pSGP->setColor (C_RED);
+    m_pSGP->setColor (C_BLUE);
     phm.draw (*m_pSGP);
     m_dTextHeight = m_pSGP->getCharHeight ();
 
@@ -348,15 +353,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
        m_pSGP->lineAbs (xs2, ys2);
        m_pSGP->setPenWidth (2);
        m_pSGP->moveAbs (0., 0.);
-       m_pSGP->drawArc (m_dFocalLength, start_angle + m_initPos.angle, start_angle - m_initPos.angle);
+       m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
       }
       m_pSGP->setPenWidth (1);
     }
     if (m_trace >= Trace::TRACE_CONSOLE)
-      traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / m_nView * 100.);
+      traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast<double>(m_nView) * 100.);
 #endif
            
-    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, dDetAngle);
+    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
     detArray.setViewAngle (viewAngle);
       
 #ifdef HAVE_SGP
@@ -366,9 +371,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
 #endif
     xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
     xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
-    if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-       dDetAngle += m_detInc;
-    else {
+    if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
        xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints 
        xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
     }
@@ -411,10 +414,10 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
   double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0;
   double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0;
   if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-      dAngleInc = m_detInc;
+      dAngleInc = m_dAngularDetIncrement;
       dAngleSampleInc = dAngleInc / m_nSample;
       dAngleSampleOffset = dAngleSampleInc / 2;
-      dAngleMajor = dDetAngle + dAngleSampleOffset;
+      dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
   } else {
       ddx = (xd2 - xd1) / detArray.nDet();  // change in coords 
       ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
@@ -456,10 +459,9 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
 #ifdef HAVE_SGP
        if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
          m_pSGP->setColor (C_YELLOW);
-         m_pSGP->setRasterOp (RO_OR_REVERSE);
+         m_pSGP->setRasterOp (RO_AND);
          m_pSGP->moveAbs (xs, ys);
          m_pSGP->lineAbs (xd, yd);
-         m_pSGP->setRasterOp (RO_SET);
        }
 #endif
 
@@ -470,14 +472,6 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
          traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
          traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
        }
-
-       //      if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
-       //        m_pSGP->setColor (C_YELLOW);
-       //        m_pSGP->setRasterOp (RO_XOR);
-       //        m_pSGP->moveAbs (xs, ys);
-       //        m_pSGP->lineAbs (xd, yd);
-       //        m_pSGP->setRasterOp (RO_SET);
-       //      }
 #endif
        if (m_idGeometry == GEOMETRY_EQUIANGULAR)
            dAngle += dAngleSampleInc;
index bc3b383..1773c8b 100644 (file)
@@ -2,7 +2,7 @@
 **  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 $
+**  $Id: mathfuncs.cpp,v 1.2 2000/08/31 08:38:58 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
@@ -63,10 +63,10 @@ integrateSimpson (const double xmin, const double xmax, const double *y, const i
 
 
 /* NAME
- *    norm_angle       Normalize angle to 0 to 2 * PI range
+ *    normalizeAngle       Normalize angle to 0 to 2 * PI range
  *
  * SYNOPSIS
- *    t = norm_angle (theta)
+ *    t = normalizeAngle (theta)
  *    double t        Normalized angle
  *    double theta     Input angle
  */
index e59449f..939f024 100644 (file)
@@ -2,7 +2,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: syserror.cpp,v 1.4 2000/08/25 15:59:13 kevin Exp $
+**  $Id: syserror.cpp,v 1.5 2000/08/31 08:38:58 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
@@ -66,7 +66,7 @@ void sys_verror (int severity, const char *msg, va_list arg)
       return;
     else if (nErrorCount == MAX_ERROR_COUNT) {
       cout << "*****************************************************************" << endl;
-      cout << "***     M A X I M U M  E R R O R  C O U N T  R E A C H E D    ***" << endl;
+      cout << "***   M A X I M U M   E R R O R   C O U N T   R E A C H E D   ***" << endl;
       cout << "***                                                           ***" << endl;
       cout << "***           No further errors will be reported              ***" << endl;
       cout << "*****************************************************************" << endl;
index bb3b841..385fbee 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsim.cpp,v 1.9 2000/07/23 01:49:03 kevin Exp $
+**  $Id: ctsim.cpp,v 1.10 2000/08/31 08:38:58 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
@@ -53,7 +53,7 @@
 #include <getopt.h>
 #endif
 
-static const char* rcsindent = "$Id";
+static const char* rcsindent = "$Id: ctsim.cpp,v 1.10 2000/08/31 08:38:58 kevin Exp $";
 
 class CTSimApp* theApp = NULL;
 
index a3c27f7..57c95f8 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.cpp,v 1.10 2000/08/27 20:32:55 kevin Exp $
+**  $Id: dialogs.cpp,v 1.11 2000/08/31 08:38:58 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
@@ -630,7 +630,7 @@ DialogGetReconstructionParameters::getFilterGenerationName (void)
 }
 
 
-DialogAutoScaleParameters::DialogAutoScaleParameters (wxFrame *pParent, const ImageFile& rIF)
+DialogAutoScaleParameters::DialogAutoScaleParameters (wxFrame *pParent, const ImageFile& rIF, double dDefaultScaleFactor)
   : wxDialog (pParent, -1, "Auto Scale Parameters", wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION), m_rImageFile(rIF)
 {
   wxBoxSizer* pTopSizer = new wxBoxSizer (wxVERTICAL);
@@ -650,7 +650,9 @@ DialogAutoScaleParameters::DialogAutoScaleParameters (wxFrame *pParent, const Im
 
   wxGridSizer *pGridSizer = new wxGridSizer (2);
   pGridSizer->Add (new wxStaticText (this, -1, "Standard Deviation Factor"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
-  m_pTextCtrlStdDevFactor = new wxTextCtrl (this, -1, "1.0", wxDefaultPosition, wxSize(100, 25), 0);
+  ostringstream osDefaultFactor;
+  osDefaultFactor << dDefaultScaleFactor;
+  m_pTextCtrlStdDevFactor = new wxTextCtrl (this, -1, osDefaultFactor.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
   pGridSizer->Add (m_pTextCtrlStdDevFactor, 0, wxALIGN_CENTER_VERTICAL);
   pTopSizer->Add (pGridSizer, 1, wxALL, 10);
 
@@ -694,3 +696,15 @@ DialogAutoScaleParameters::getMinMax (double* pMin, double* pMax)
   *pMax = dCenter + dHalfWidth;
   *theApp->getLog() << "Setting minimum to " << *pMin << " and maximum to " << *pMax << "\n";
 }
+
+double
+DialogAutoScaleParameters::getAutoScaleFactor ()
+{
+  wxString sStddevFactor = m_pTextCtrlStdDevFactor->GetValue();
+  double dValue = 1.;
+  if (! sStddevFactor.ToDouble (&dValue)) {
+    *theApp->getLog() << "Error: Non-numeric Standard Deviation Factor of " << sStddevFactor << "\n";
+  }
+
+  return dValue;
+}
index 77b7cbf..7271339 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dialogs.h,v 1.11 2000/08/27 20:32:55 kevin Exp $
+**  $Id: dialogs.h,v 1.12 2000/08/31 08:38:58 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
@@ -185,10 +185,11 @@ class DialogGetReconstructionParameters : public wxDialog
 class DialogAutoScaleParameters : public wxDialog
 {
  public:
-    DialogAutoScaleParameters (wxFrame* pParent, const ImageFile& rImageFile);
+    DialogAutoScaleParameters (wxFrame* pParent, const ImageFile& rImageFile, double dDefaultScaleFactor = 1.);
     virtual ~DialogAutoScaleParameters() {}
 
     void getMinMax (double* pMin, double* pMax);
+    double getAutoScaleFactor ();
 
  private:
     const ImageFile& m_rImageFile;
index 8f66dc3..4f79334 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dlgprojections.cpp,v 1.1 2000/08/27 20:32:55 kevin Exp $
+**  $Id: dlgprojections.cpp,v 1.2 2000/08/31 08:38:58 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
@@ -58,17 +58,20 @@ static const int LAYOUT_Y_MARGIN = 4;
 
 BEGIN_EVENT_TABLE(ProjectionsDialog, wxDialog)
    EVT_BUTTON(wxID_CANCEL, ProjectionsDialog::OnCancel)
+   EVT_BUTTON(ID_BTN_PAUSE, ProjectionsDialog::OnPause)
+   EVT_BUTTON(ID_BTN_STEP, ProjectionsDialog::OnStep)
    EVT_CLOSE(ProjectionsDialog::OnClose)
+   EVT_PAINT(ProjectionsDialog::OnPaint)
 END_EVENT_TABLE()
 
 IMPLEMENT_CLASS(ProjectionsDialog, wxDialog)
 
 
 ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, const Phantom& rPhantom, const int iTrace, wxWindow *parent)
-  : wxDialog(parent, -1, "Collect Projections"), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL)
+    : wxDialog(parent, -1, "Collect Projections"), m_rScanner(rScanner), m_rProjections(rProj), m_rPhantom(rPhantom), m_pSGPDriver(NULL), m_pSGP(NULL), m_iTrace(iTrace), m_pDC(NULL), m_btnAbort(0), m_btnPause(0), m_btnStep(0)
 {
     m_state = Continue;
-
+    m_iLastView = -1;
     m_parentTop = parent;
     while ( m_parentTop && m_parentTop->GetParent() )
         m_parentTop = m_parentTop->GetParent();
@@ -84,6 +87,22 @@ ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, con
 
     m_btnAbort->SetConstraints(c);
 
+    m_btnPause = new wxButton (this, ID_BTN_PAUSE, wxString("Pause"));
+    wxLayoutConstraints* cPause = new wxLayoutConstraints;
+    cPause->right.SameAs(this, wxRight, 3*LAYOUT_X_MARGIN + sizeBtn.x);
+    cPause->bottom.SameAs(this, wxBottom, 2*LAYOUT_Y_MARGIN);
+    cPause->width.Absolute(sizeBtn.x);
+    cPause->height.Absolute(sizeBtn.y);
+    m_btnPause->SetConstraints(cPause);
+
+    m_btnStep = new wxButton (this, ID_BTN_STEP, wxString("Step"));
+    wxLayoutConstraints* cStep = new wxLayoutConstraints;
+    cStep->right.SameAs(this, wxRight, 5*LAYOUT_X_MARGIN + sizeBtn.x * 2);
+    cStep->bottom.SameAs(this, wxBottom, 2*LAYOUT_Y_MARGIN);
+    cStep->width.Absolute(sizeBtn.x);
+    cStep->height.Absolute(sizeBtn.y);
+    m_btnStep->SetConstraints(cStep);
+
     SetAutoLayout(TRUE);
     Layout();
 
@@ -92,7 +111,13 @@ ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, con
       sizeDlg.x = max(sizeDlg.x,sizeDlg.y);
       sizeDlg.y = max(sizeDlg.x,sizeDlg.y);
     }
+    if (m_iTrace >= Trace::TRACE_PLOT)
+      sizeDlg.x += 200;
+
+    m_iClientX = sizeDlg.x;
+    m_iClientY = sizeDlg.y;
     SetClientSize(sizeDlg);
+    m_bitmap.Create (m_iClientX, m_iClientY); // save a copy of screen
 
     Centre(wxCENTER_FRAME | wxBOTH);
 
@@ -110,34 +135,64 @@ ProjectionsDialog::ProjectionsDialog (Scanner& rScanner, Projections& rProj, con
 
     wxYield();     // Update the display
 
+    m_pSGPDriver->idWX()->SetFont(*wxSWISS_FONT);
 #ifdef __WXMAC__
     MacUpdateImmediately();
 #endif
 }
 
+void
+ProjectionsDialog::showView (int iViewNumber)
+{
+    if ( iViewNumber < m_rProjections.nView() ) {
+       wxYield();        // update the display
+       m_iLastView = iViewNumber;
+       if (m_iTrace >= Trace::TRACE_PLOT)
+           m_pSGP->setViewport (0, 0, 0.75, 1);
+       m_rScanner.collectProjections (m_rProjections, m_rPhantom, iViewNumber, 1, true, m_iTrace, m_pSGP);
+       if (m_iTrace >= Trace::TRACE_PLOT) {
+           const DetectorArray& detArray = m_rProjections.getDetectorArray (iViewNumber);
+           const DetectorValue* detValues = detArray.detValues();
+           double detPos [detArray.nDet()];
+           for (int i = 0; i < detArray.nDet(); i++)
+               detPos[i] = i;
+           EZPlot ezplot (*m_pSGP);
+           ezplot.ezset("xporigin 0.75");
+           ezplot.ezset("yporigin 0.10");
+           ezplot.ezset("xlength  0.25");
+           ezplot.ezset("ylength  0.90");
+           ezplot.ezset("grid");
+           ezplot.ezset("box");
+           ezplot.addCurve (detValues, detPos, detArray.nDet());
+           ezplot.plot();
+       }
+    }
+}
+
 bool
 ProjectionsDialog::projectView (int iViewNumber)
 {
-    if ( iViewNumber < m_rProjections.nView() ) {
-      m_rScanner.collectProjections (m_rProjections, m_rPhantom, iViewNumber, 1, true, m_iTrace, m_pSGP);
-      wxYield();        // update the display
-    } else {
-      m_state = Finished;    // so that we return TRUE below and 
-                             // that [Cancel] handler knew what to do
+    showView (iViewNumber);
+    wxYield();        // update the display
+    if (m_iTrace >= Trace::TRACE_PLOT)
+       sleep(1);
+    else {
+       m_state = Finished;    // so that we return TRUE below and 
+       // that [Cancel] handler knew what to do
 #if 0
-      if ( m_btnAbort ) 
-       m_btnAbort->SetLabel(_("Close")); // tell the user what he should do...
-      wxYield();
-      
-      (void)ShowModal();
+       if ( m_btnAbort ) 
+           m_btnAbort->SetLabel(_("Close")); // tell the user what he should do...
+       wxYield();
+       
+       (void)ShowModal();
 #endif
     }
-      
+    
 #ifdef __WXMAC__
     MacUpdateImmediately();
 #endif
 
-    return m_state != Canceled;
+    return m_state != Cancelled;
 }
 
 
@@ -152,7 +207,7 @@ void ProjectionsDialog::OnCancel (wxCommandEvent& event)
   } else {
     // request to cancel was received, the next time Update() is called we
     // will handle it
-    m_state = Canceled;
+    m_state = Cancelled;
 
     // update the button state immediately so that the user knows that the
     // request has been noticed
@@ -160,14 +215,71 @@ void ProjectionsDialog::OnCancel (wxCommandEvent& event)
   }
 }
 
+
+void 
+ProjectionsDialog::OnPause (wxCommandEvent& event)
+{
+  if ( m_state == Finished ) {
+    // this means that the count down is already finished and we're being
+    // shown as a modal dialog - so just let the default handler do the job
+    event.Skip();
+  } else {
+      if (m_state == Continue) {
+         m_memoryDC.SelectObject (m_bitmap);       // in memoryDC
+         m_pSGP->setDC (&m_memoryDC);
+         m_memoryDC.SetFont (*wxSWISS_FONT);
+         showView (m_iLastView);
+         m_pSGP->setDC (m_pDC);
+         m_memoryDC.SelectObject(wxNullBitmap);
+         m_state = Paused;
+         m_btnPause->SetLabel (wxString("Resume"));
+      } else if (m_state == Paused) {
+         m_state = Continue;
+         m_btnPause->SetLabel (wxString("Pause"));
+      }
+  }
+}
+
+void 
+ProjectionsDialog::OnStep (wxCommandEvent& event)
+{
+  if ( m_state == Finished ) {
+    // this means that the count down is already finished and we're being
+    // shown as a modal dialog - so just let the default handler do the job
+    event.Skip();
+  } else {
+      if (m_state == Continue) {
+         m_memoryDC.SelectObject (m_bitmap);       // in memoryDC
+         m_pSGP->setDC (&m_memoryDC);
+         m_memoryDC.SetFont (*wxSWISS_FONT);
+         m_rScanner.collectProjections (m_rProjections, m_rPhantom, m_iLastView, 1, true, m_iTrace, m_pSGP);
+         m_pSGP->setDC (m_pDC);
+         m_memoryDC.SelectObject(wxNullBitmap);
+         m_state = Paused;
+         m_btnPause->SetLabel (wxString("Resume"));
+      } else if (m_state == Paused) {
+         projectView (m_iLastView + 1);
+      }
+  }
+}
+
 void ProjectionsDialog::OnClose(wxCloseEvent& event)
 {
-    if ( m_state == Uncancelable )
+    if ( m_state == Uncancellable )
       event.Veto(TRUE);    // can't close this dialog
     else if ( m_state == Finished )
       event.Skip(); // let the default handler close the window as we already terminated
     else
-      m_state = Canceled;          // next Update() will notice it
+      m_state = Cancelled;          // next Update() will notice it
+}
+
+void
+ProjectionsDialog::OnPaint (wxPaintEvent& event)
+{
+  wxPaintDC paintDC (this);
+  if (m_state == Paused) {
+    paintDC.DrawBitmap(m_bitmap, 0, 0, false);
+  }
 }
 
 
index 440d2ff..b8798bb 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: dlgprojections.h,v 1.1 2000/08/27 20:32:55 kevin Exp $
+**  $Id: dlgprojections.h,v 1.2 2000/08/31 08:38:58 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
@@ -34,6 +34,7 @@
 
 #include "wx/setup.h"
 #include "wx/dialog.h"
+#include "wx/dcmemory.h"
 
 class wxButton;
 class wxStaticText;
@@ -69,30 +70,51 @@ public:
    void OnCancel(wxCommandEvent& event);
        // callback to disable "hard" window closing
    void OnClose(wxCloseEvent& event);
+   void OnPaint(wxPaintEvent& event);
+
+   void OnPause(wxCommandEvent& event);
+   void OnStep(wxCommandEvent& event);
+
+   bool isPaused() const {return m_state == Paused;}
+
+   bool isCancelled() const {return m_state == Cancelled;}
 
 private:
    // parent top level window (may be NULL)
    wxWindow *m_parentTop;
+   int m_iLastView;
+   int m_iClientX;   // size of client window
+   int m_iClientY;
+
+   Scanner& m_rScanner;
+   Projections& m_rProjections;
+   const Phantom& m_rPhantom;
+   SGPDriver* m_pSGPDriver;
+   SGP* m_pSGP;
+   const int m_iTrace;
+   wxDC* m_pDC;
+
+   wxButton *m_btnAbort;    // the abort button (or NULL if none)
+   wxButton *m_btnPause; 
+   wxButton *m_btnStep;
+
+   wxMemoryDC m_memoryDC;  // for restoring image on OnPaint
+   wxBitmap m_bitmap;
 
    // continue processing or not (return value for Update())
    enum
    {
-      Uncancelable = -1,   // dialog can't be canceled
-      Canceled,            // can be cancelled and, in fact, was
+      Uncancellable = -1,   // dialog can't be canceled
+      Paused,
+      Cancelled,            // can be cancelled and, in fact, was
       Continue,            // can be cancelled but wasn't
       Finished             // finished, waiting to be removed from screen
    } m_state;
 
-   // the abort button (or NULL if none)
-   wxButton *m_btnAbort;
+   const static int ID_BTN_PAUSE = 19998;
+   const static int ID_BTN_STEP = 19999;
 
-   Scanner& m_rScanner;
-   Projections& m_rProjections;
-   const Phantom& m_rPhantom;
-   SGPDriver* m_pSGPDriver;
-   SGP* m_pSGP;
-   const int m_iTrace;
-   wxDC* m_pDC;
+   void showView (int iViewNumber);
 
    DECLARE_EVENT_TABLE()
 };
index edc35e6..3df51fa 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.17 2000/08/27 20:32:55 kevin Exp $
+**  $Id: views.cpp,v 1.18 2000/08/31 08:38:58 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
@@ -152,7 +152,7 @@ void
 ImageFileView::OnScaleAuto (wxCommandEvent& event)
 {
     const ImageFile& rIF = GetDocument()->getImageFile();
-    DialogAutoScaleParameters dialogAutoScale (m_frame, rIF);
+    DialogAutoScaleParameters dialogAutoScale (m_frame, rIF, m_dAutoScaleFactor);
     int iRetVal = dialogAutoScale.ShowModal();
     if (iRetVal == wxID_OK) {
       m_bMinSpecified = true;
@@ -161,6 +161,7 @@ ImageFileView::OnScaleAuto (wxCommandEvent& event)
       dialogAutoScale.getMinMax (&dMin, &dMax);
       m_dMinPixel = dMin;
       m_dMaxPixel = dMax;
+      m_dAutoScaleFactor = dialogAutoScale.getAutoScaleFactor();
       OnUpdate (this, NULL);
     }
 }
@@ -256,6 +257,7 @@ ImageFileView::OnCreate(wxDocument *doc, long WXUNUSED(flags) )
     
     m_bMinSpecified = false;
     m_bMaxSpecified = false;
+    m_dAutoScaleFactor = 1.;
 
     int width, height;
     m_frame->GetClientSize(&width, &height);
@@ -430,14 +432,21 @@ PhantomView::OnProjections (wxCommandEvent& event)
       if (iTrace > Trace::TRACE_CONSOLE) {
        ProjectionsDialog dialogProjections (theScanner, rProj, rPhantom, iTrace, dynamic_cast<wxWindow*>(m_frame));
        for (int iView = 0; iView < rProj.nView(); iView++) {
-         if (! dialogProjections.projectView (iView)) {
+         ::wxYield();
+         ::wxYield();
+         if (dialogProjections.isCancelled() || ! dialogProjections.projectView (iView)) {
            pProjectionDoc->DeleteAllViews();
            return;
          }
          ::wxYield();
+         ::wxYield();
+         while (dialogProjections.isPaused()) {
+             ::wxYield();
+             ::wxUsleep(50);
+         }
        }
       } else
-       theScanner.collectProjections (rProj, rPhantom);
+       theScanner.collectProjections (rProj, rPhantom, iTrace);
 
       pProjectionDoc->Modify(true);
       pProjectionDoc->UpdateAllViews(this);
@@ -447,6 +456,7 @@ PhantomView::OnProjections (wxCommandEvent& event)
        }
       ostringstream os;
       os << "Projections for " << rPhantom.name() << ": nDet=" << nDet << ", nView=" << nView << ", nSamples=" << nSamples << ", RotAngle=" << dRotAngle << ", FocalLengthRatio=" << dFocalLengthRatio << ", FieldOfViewRatio=" << dFieldOfViewRatio << ", Geometry=" << sGeometry.c_str() << "\n";
+      rProj.setRemark (os.str());
       *theApp->getLog() << os.str().c_str();
     }
   }
@@ -642,7 +652,7 @@ ProjectionFileView::OnProperties (wxCommandEvent& event)
 {
   const Projections& rProj = GetDocument()->getProjections();
   ostringstream os;
-  os << "ProjectionFile " << rProj.getFilename() << ":  Number of Detectors = " << rProj.nDet() << ", Number of Views = " << rProj.nView() << "\n";
+  rProj.printScanInfo(os);
   *theApp->getLog() << os.str().c_str();
 }
 
index a11cfaf..71efce2 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: views.h,v 1.7 2000/08/25 15:59:13 kevin Exp $
+**  $Id: views.h,v 1.8 2000/08/31 08:38:58 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
@@ -54,6 +54,7 @@ private:
     bool m_bMaxSpecified;
     double m_dMinPixel;
     double m_dMaxPixel;
+    double m_dAutoScaleFactor;
 
 public:
     ImageFileView(void);
index a79c7e2..ff51444 100644 (file)
@@ -1,4 +1,4 @@
-bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@  ifinfo phm2if if-1 if-2 if2img 
+bin_PROGRAMS = pjrec phm2pj pj2if @lamprograms@  ifinfo phm2if if-1 if-2 if2img pjinfo
 bin_SCRIPTS = sample-ctsim.sh
 EXTRA_PROGRAMS = pjrec-lam phm2pj-lam phm2if-lam
 INCLUDES=@my_includes@
@@ -25,6 +25,9 @@ if2img_DEPENDENCIES=$(SOURCE_DEPEND)
 pj2if_SOURCES = pj2if.cpp
 pj2if_LDADD=@ctlibs@
 pj2if_DEPENDENCIES=$(SOURCE_DEPEND)
+pjinfo_SOURCES = pjinfo.cpp
+pjinfo_LDADD=@ctlibs@
+pjinfo_DEPENDENCIES=$(SOURCE_DEPEND)
 if_1_SOURCES=if-1.cpp 
 if_1_LDADD=@ctlibs@
 if_1_DEPENDENCIES=$(SOURCE_DEPEND)
index 0459f53..c42f0af 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phm2pj.cpp,v 1.11 2000/08/27 20:32:55 kevin Exp $
+**  $Id: phm2pj.cpp,v 1.12 2000/08/31 08:38:58 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
@@ -50,7 +50,7 @@ static struct option phm2pj_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: phm2pj.cpp,v 1.11 2000/08/27 20:32:55 kevin Exp $";
+static const char* g_szIdStr = "$Id: phm2pj.cpp,v 1.12 2000/08/31 08:38:58 kevin Exp $";
 
 
 void 
@@ -235,7 +235,7 @@ phm2pj_main (int argc, char* argv[])
     }
 
     ostringstream desc;
-    desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", ";
+    desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", Geometry=" << optGeometryName << ", ";
     if (optPhmFileName.length()) {
       desc << "PhantomFile=" << optPhmFileName;
     } else if (optPhmName != "") {
@@ -296,8 +296,11 @@ phm2pj_main (int argc, char* argv[])
   if (mpiWorld.getRank() == 0)
     pjGlobal.initFromScanner (scanner);
   
-  if (opt_verbose)
-    pjGlobal.printScanInfo();
+  if (opt_verbose) {
+    ostringstream os;
+    pjGlobal.printScanInfo(os);
+    cout << os;
+  }
 
   Projections pjLocal (scanner);
   pjLocal.setNView (mpiWorld.getMyLocalWorkUnits());
@@ -348,8 +351,9 @@ phm2pj_main (int argc, char* argv[])
       if (opt_verbose) {
        phm.print();
        cout << endl;
-       pjGlobal.printScanInfo();
-       cout << endl;
+       ostringstream os;
+       pjGlobal.printScanInfo (os);
+       cout << os.str() << endl;
        cout << "  Remark: " << pjGlobal.remark() << endl;
        cout << "Run time: " << pjGlobal.calcTime() << " seconds\n";
       }
index f65d55d..b9bd2e9 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pj2if.cpp,v 1.3 2000/08/03 09:57:29 kevin Exp $
+**  $Id: pj2if.cpp,v 1.4 2000/08/31 08:38:58 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 "timer.h"
 
 
-enum { O_VERBOSE, O_HELP, O_VERSION };
+enum { O_VERBOSE, O_DUMP, O_HELP, O_VERSION };
 
 static struct option my_options[] =
 {
   {"verbose", 0, 0, O_VERBOSE},
+  {"dump", 0, 0, O_DUMP},
   {"help", 0, 0, O_HELP},
   {"version", 0, 0, O_VERSION},
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: pj2if.cpp,v 1.3 2000/08/03 09:57:29 kevin Exp $";
+static const char* g_szIdStr = "$Id: pj2if.cpp,v 1.4 2000/08/31 08:38:58 kevin Exp $";
 
 void 
 pj2if_usage (const char *program)
@@ -55,6 +56,7 @@ pj2if_usage (const char *program)
   cout << "Converts a projection file to a IF file" << endl;
   cout << endl;
   cout << "   --verbose   Verbose mode" << endl;
+  cout << "   --dump      Dump all scan data" << endl;
   cout << "   --version   Print version" << endl;
   cout << "   --help      Print this help message" << endl;
 }
@@ -66,6 +68,7 @@ pj2if_main (const int argc, char *const argv[])
 {
   char *pj_name, *im_name;
   bool optVerbose = false;
+  bool optDump = false;
   extern int optind;
   Timer timerProgram;
 
@@ -80,6 +83,9 @@ pj2if_main (const int argc, char *const argv[])
        case O_VERBOSE:
          optVerbose = true;
          break;
+       case O_DUMP:
+         optDump = true;
+         break;
         case O_VERSION:
 #ifdef VERSION
          cout << "Version " << VERSION << endl << g_szIdStr << endl;
@@ -111,8 +117,13 @@ pj2if_main (const int argc, char *const argv[])
     return (1);
   }
 
-  if (optVerbose)
-      pj.printScanInfo();
+  if (optDump)
+    pj.printProjectionData();
+  else if (optVerbose) {
+    ostringstream os;
+    pj.printScanInfo (os);
+    cout << os.str();
+  }
   
   ImageFile im (pj.nDet(), pj.nView());
   
index 157d658..9fcb890 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: pjrec.cpp,v 1.15 2000/08/27 20:32:55 kevin Exp $
+**  $Id: pjrec.cpp,v 1.16 2000/08/31 08:38:58 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
@@ -49,7 +49,7 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.15 2000/08/27 20:32:55 kevin Exp $";
+static const char* g_szIdStr = "$Id: pjrec.cpp,v 1.16 2000/08/31 08:38:58 kevin Exp $";
 
 void 
 pjrec_usage (const char *program)
@@ -260,8 +260,11 @@ pjrec_main (int argc, char * argv[])
 #ifdef HAVE_MPI
   if (mpiWorld.getRank() == 0) {
     projGlobal.read (pszFilenameProj);
-    if (bOptVerbose)
-      projGlobal.printScanInfo();
+    if (bOptVerbose) {
+      ostringstream os;
+      projGlobal.printScanInfo (os);
+      cout << os.str();
+    }
 
     mpi_ndet = projGlobal.nDet();
     mpi_nview = projGlobal.nView();
@@ -310,8 +313,11 @@ pjrec_main (int argc, char * argv[])
   imLocal = new ImageFile (nx, ny);
 #else
   projGlobal.read (pszFilenameProj);
-  if (bOptVerbose)
-    projGlobal.printScanInfo();
+  if (bOptVerbose) {
+    ostringstream os;
+    projGlobal.printScanInfo(os);
+    cout << os.str();
+  }
 
   imGlobal = new ImageFile (nx, ny);
 #endif