r1018: *** empty log message ***
authorKevin M. Rosenberg <kevin@rosenberg.net>
Mon, 24 Sep 2001 09:40:42 +0000 (09:40 +0000)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Mon, 24 Sep 2001 09:40:42 +0000 (09:40 +0000)
29 files changed:
Makefile.am
Makefile.in
configure
configure.in
helical/Makefile.am [new file with mode: 0644]
helical/dynphm.c [new file with mode: 0644]
helical/sample-helical.sh.in [new file with mode: 0644]
include/imagefile.h
include/projections.h
include/scanner.h
libctsim/imagefile.cpp
libctsim/projections.cpp
libctsim/scanner.cpp
src/ctsim.cpp
src/dialogs.cpp
src/dialogs.h
src/dlgprojections.cpp
src/threadproj.cpp
src/threadproj.h
src/views.cpp
src/views.h
tools/Makefile.am
tools/Makefile.in
tools/Makefile.nt
tools/ctsimtext.cpp
tools/ifexport.cpp
tools/phm2helix.cpp [new file with mode: 0644]
tools/phm2pj.cpp
tools/pjHinterp.cpp [new file with mode: 0644]

index 15cdf07..09164f7 100644 (file)
@@ -16,7 +16,7 @@ else
 EXTRA_DIRS3=
 endif
 
-SUBDIRS=$(EXTRA_DIRS1) $(EXTRA_DIRS2) libctsupport libctsim man doc html cgi-bin include  $(EXTRA_DIRS3) tools
+SUBDIRS=$(EXTRA_DIRS1) $(EXTRA_DIRS2) libctsupport libctsim man doc html cgi-bin include  $(EXTRA_DIRS3) tools helical
 
 EXTRA_DIST=acsite.m4 make.bat msvc/ctsim.dsw msvc/ctsimtext/ctsimtext.dsp msvc/libctsim/libctsim.dsp msvc/ctsim/ctsim.dsp
 
index 2d6f0ff..aba3a41 100644 (file)
@@ -96,7 +96,7 @@ wxlibs = @wxlibs@
 @HAVE_WXWINDOWS_TRUE@EXTRA_DIRS3 = src
 @HAVE_WXWINDOWS_FALSE@EXTRA_DIRS3 = 
 
-SUBDIRS = $(EXTRA_DIRS1) $(EXTRA_DIRS2) libctsupport libctsim man doc html cgi-bin include  $(EXTRA_DIRS3) tools
+SUBDIRS = $(EXTRA_DIRS1) $(EXTRA_DIRS2) libctsupport libctsim man doc html cgi-bin include  $(EXTRA_DIRS3) tools helical
 
 EXTRA_DIST = acsite.m4 make.bat msvc/ctsim.dsw msvc/ctsimtext/ctsimtext.dsp msvc/libctsim/libctsim.dsp msvc/ctsim/ctsim.dsp
 ACLOCAL_M4 = $(top_srcdir)/aclocal.m4
@@ -114,7 +114,7 @@ DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST)
 TAR = gtar
 GZIP_ENV = --best
 DIST_SUBDIRS =  getopt libctgraphics libctsupport libctsim man doc html \
-cgi-bin include src tools
+cgi-bin include src tools helical
 all: all-redirect
 .SUFFIXES:
 $(srcdir)/Makefile.in: Makefile.am $(top_srcdir)/configure.in $(ACLOCAL_M4) 
index c4d45b0..efc5a00 100755 (executable)
--- a/configure
+++ b/configure
@@ -1925,7 +1925,7 @@ if eval "test \"`echo '$''{'ac_cv_lib_$ac_lib_var'+set}'`\" = set"; then
   echo $ac_n "(cached) $ac_c" 1>&6
 else
   ac_save_LIBS="$LIBS"
-LIBS="-lreadline  $LIBS"
+LIBS="-lreadline -lcurses $LIBS"
 cat > conftest.$ac_ext <<EOF
 #line 1931 "configure"
 #include "confdefs.h"
@@ -4883,7 +4883,7 @@ done
 ac_given_srcdir=$srcdir
 ac_given_INSTALL="$INSTALL"
 
-trap 'rm -fr `echo "Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile config.h" | sed "s/:[^ ]*//g"` conftest*; exit 1' 1 2 15
+trap 'rm -fr `echo "Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile helical/Makefile helical/sample-helical.sh config.h" | sed "s/:[^ ]*//g"` conftest*; exit 1' 1 2 15
 EOF
 cat >> $CONFIG_STATUS <<EOF
 
@@ -5014,7 +5014,7 @@ EOF
 
 cat >> $CONFIG_STATUS <<EOF
 
-CONFIG_FILES=\${CONFIG_FILES-"Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile"}
+CONFIG_FILES=\${CONFIG_FILES-"Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile helical/Makefile helical/sample-helical.sh"}
 EOF
 cat >> $CONFIG_STATUS <<\EOF
 for ac_file in .. $CONFIG_FILES; do if test "x$ac_file" != x..; then
index 7463752..bb61c86 100644 (file)
@@ -39,7 +39,7 @@ AC_CHECK_LIB(m, sin)
 AC_CHECK_LIB(curses, main, [curses=true], [curses=false])
 AC_CHECK_LIB(ncurses, main, [ncurses=true], [ncurses=false])
 AC_CHECK_LIB(g2, main, [g2=true], [g2=false])
-AC_CHECK_LIB(readline, main, [readline=true; AC_DEFINE(HAVE_READLINE)], [readline=false])
+AC_CHECK_LIB(readline, main, [readline=true; AC_DEFINE(HAVE_READLINE)], [readline=false], [-lcurses])
 AC_CHECK_LIB(ncurses, main, [ncurses=true; AC_DEFINE(HAVE_NCURSES)], [ncurses=false])
 AC_CHECK_LIB(curses, main, [curses=true; AC_DEFINE(HAVE_CURSES)], [curses=false])
 wxwin=false
@@ -454,4 +454,4 @@ fi
 
 CXXFLAGS="$CFLAGS"
 
-AC_OUTPUT(Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile)
+AC_OUTPUT(Makefile libctgraphics/Makefile libctsupport/Makefile libctsim/Makefile Makefile man/Makefile doc/Makefile cgi-bin/ctsim.cgi cgi-bin/Makefile html/simulate.html html/Makefile include/Makefile getopt/Makefile tools/sample-ctsim.sh cgi-bin/ctsim.conf tools/Makefile src/Makefile helical/Makefile helical/sample-helical.sh)
diff --git a/helical/Makefile.am b/helical/Makefile.am
new file mode 100644 (file)
index 0000000..ce6143e
--- /dev/null
@@ -0,0 +1,5 @@
+bin_SCRIPTS = sample-helical.sh
+
+realclean:
+       rm -f *.pgm *.if *~ *.pj
+
diff --git a/helical/dynphm.c b/helical/dynphm.c
new file mode 100644 (file)
index 0000000..4395f02
--- /dev/null
@@ -0,0 +1,49 @@
+#include <stdio.h>
+#include <math.h>
+int
+main(int argc, char *argv[])
+{
+       int view, nview; 
+       char filename[128];                                                      
+       
+       float mu1=0. , mu2=6., density, afac, s;
+
+       FILE *fp = (FILE *)NULL;
+       if (argc !=4){ 
+               fprintf(stderr, "Usage: %s iview nview phmfilename\n", argv[0]);        
+               exit (1);
+       }
+
+       view = atoi(argv[1]);
+       nview = atoi(argv[2]);
+       sprintf(filename, "%s", argv[3]);
+
+       s = (float)view/((float)(nview-1));
+  
+   if ( s < 7./16. ) 
+               density = mu1; 
+   else if ( s  > 9./16. )
+       density = mu2;
+   else {
+       afac = ( (s - 7./16.) / 2./16.); 
+       density = log(1/((1-afac)*exp(-mu1) + afac * exp(-mu2)));
+   }                                   
+       
+/* 
+       density = mu1 + (mu2-mu1)*s;
+       if (s <=0.5) 
+                       density = mu1; 
+       else
+                       density = mu2; 
+ */
+       if ( (fp = fopen(filename, "w"))  == (FILE *)NULL){
+               fprintf(stderr,"Error, can not open file \"tmpphmfile\"\n"); 
+               exit(2);
+       }
+       fprintf(fp, "rectangle 0 0 11.5 11.5 0 0\n");
+       fprintf(fp, "ellipse   0 0 11.4 11.4 0 1\n");
+       fprintf(fp, "ellipse   0 0 1.25 1.25 0 %f\n", density);
+
+       fclose(fp);     
+       exit(0);
+}
diff --git a/helical/sample-helical.sh.in b/helical/sample-helical.sh.in
new file mode 100644 (file)
index 0000000..9e12f20
--- /dev/null
@@ -0,0 +1,50 @@
+#!/bin/sh
+
+if test "$1" != "" ; then
+  bin=$1
+else
+  bin="@prefix@/bin/"
+fi
+
+if test "$1" = "clean" ; then
+  rm -f sample-phm.png sample-phm.if sample-pj.pj sample-pj.if sample-pj.png sample-rec.if sample-rec.png 540-pj.pj  540-pj.if 540-pj.png  540-rec.if 540-rec.png dynphm
+  exit
+fi
+
+# Generate phantom image
+
+cc -o dynphm dynphm.c -lm
+dynphm 100 200 tmpphm
+${bin}ctsimtext phm2if sample-phm.if 256 256 --nsample 2 --phmfile tmpphm
+rm -f tmpphm 
+if [ -f sample-phm.if ] ; then
+  ${bin}ctsimtext ifexport sample-phm.if sample-phm.png --format png
+fi
+
+# Simulate helical CT data collection and generate raysum sinugram for display
+${bin}ctsimtext phm2helix  sample-pj.pj 367 1080 dynphm --nray 2  --geometry equiangular --rotangle 3
+if [ -f sample-pj.pj ]; then
+  ${bin}ctsimtext pj2if  sample-pj.pj sample-pj.if
+fi
+if [ -f sample-pj.if ]; then
+  ${bin}ctsimtext ifexport sample-pj.if sample-pj.png --format png
+fi
+if [ -f sample-pj.pj ]; then 
+  ${bin}ctsimtext pjHinterp sample-pj.pj 540-pj.pj --interpview 540
+fi
+if [ -f 540-pj.pj ]; then
+  ${bin}ctsimtext pj2if  540-pj.pj 540-pj.if
+fi
+if [ -f sample-pj.if ]; then
+  ${bin}ctsimtext ifexport 540-pj.if 540-pj.png --format png
+fi
+
+# Reconstruct raysums and generate image for display
+${bin}ctsimtext pjrec   540-pj.pj 540-rec.if 256 256 
+if [ -f sample-rec.if ]; then 
+  ${bin}ctsimtext ifexport 540-rec.if 540-rec.png --format png
+
+  # Display comparison statistics
+  ${bin}ctsimtext if2 sample-phm.if 540-rec.if --comp
+fi
+
index 56c255b..6dd8fe7 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: imagefile.h,v 1.35 2001/03/21 21:45:31 kevin Exp $
+**  $Id: imagefile.h,v 1.36 2001/09/24 09:40:42 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
@@ -149,6 +149,7 @@ public:
   static const int EXPORT_FORMAT_DICOM;
   static const int IMPORT_FORMAT_DICOM;
 #endif
+  static const int EXPORT_FORMAT_RAW;
 
   static const int getExportFormatCount() {return s_iExportFormatCount;}
   static const char** getExportFormatNameArray() {return s_aszExportFormatName;}
@@ -231,6 +232,7 @@ public:
   bool writeImagePGM (const char* const outfile, int nxcell, int nycell, double densmin, double densmax);
   bool writeImagePGMASCII (const char* const outfile, int nxcell, int nycell, double densmin, double densmax);
   bool readImagePPM (const char* const pszFile);
+  bool writeImageRaw(const char* const outfile, int nxcell, int nycell);
   bool writeImageText (const char* const outfile);
 
   static double redGrayscaleFactor() {return s_dRedGrayscaleFactor;}
index d2fe89c..b021174 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: projections.h,v 1.35 2001/03/21 21:45:31 kevin Exp $
+**  $Id: projections.h,v 1.36 2001/09/24 09:40:42 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -122,6 +122,10 @@ class Projections
   void printProjectionData ();
   void printScanInfo (std::ostringstream& os) const;
 
+  int Helical180LI(int interpView);
+  int Helical180LI_Equiangular(int interpView);
+  int HalfScanFeather(void);
+
   bool read (const std::string& fname);
   bool read (const char* fname);
   bool write (const char* fname);
index 663ebd7..09ab796 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.h,v 1.21 2001/03/11 06:34:37 kevin Exp $
+**  $Id: scanner.h,v 1.22 2001/09/24 09:40:42 kevin Exp $
 **
 **
 **  This program is free software; you can redistribute it and/or modify
@@ -75,26 +75,26 @@ class Scanner
 
   
   Scanner (const Phantom& phm, const char* const geometryName, int nDet, 
-    int nView, int nSample, const double rot_anglen, 
+    int nView, int iOffsetView, int nSample, const double rot_anglen, 
     double dFocalLengthRatio, double dCenterDetectorRatio, double dViewRatio, double dScanRatio);
   ~Scanner();
   
   void collectProjections (Projections& proj, const Phantom& phm, const int trace = Trace::TRACE_NONE,
     SGP* pSGP = NULL);
 
-  void collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-    bool bStoreAtViewPosition, const int trace = Trace::TRACE_NONE, SGP* pSGP = NULL);
+  void collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, const int iOffsetView, bool bStoreAtViewPosition, const int trace = Trace::TRACE_NONE, SGP* pSGP = NULL);
 
-  void collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-    int iStorageOffset, const int trace = Trace::TRACE_NONE, SGP* pSGP = NULL);
+  void collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, const int iOffsetView, int iStorageOffset, const int trace = Trace::TRACE_NONE, SGP* pSGP = NULL);
 
   void setNView (int nView);
+  void setOffsetView (int iOffsetView);
 
   bool fail() const {return m_fail;}
   const std::string& failMessage() const {return m_failMessage;}
   unsigned int nDet() const {return m_nDet;}
   unsigned int nView() const {return m_nView;}
+  unsigned int offsetView() const {return m_iOffsetView;}
+  unsigned int startView() const {return m_startView;} 
   double rotInc() const {return m_rotInc;}
   double detInc() const {return m_detInc;}
   double detLen() const {return m_detLen;}
@@ -122,6 +122,8 @@ class Scanner
   int m_idGeometry;
   unsigned int m_nDet;         /* Number of detectors in array */
   unsigned int m_nView;                /* Number of rotated views */
+  unsigned int m_iOffsetView; 
+  unsigned int m_startView; 
   unsigned int m_nSample;      /* Number of rays per detector */
   double m_dFocalLength;        // Focal Length, distance from source to center
   double m_dSourceDetectorLength; // Distance from source to detectors
index 356e489..d10e1bd 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: imagefile.cpp,v 1.44 2001/03/30 22:09:09 kevin Exp $
+**  $Id: imagefile.cpp,v 1.45 2001/09/24 09:40:42 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
@@ -46,6 +46,7 @@ const int ImageFile::EXPORT_FORMAT_PNG16 = 4;
 #ifdef HAVE_CTN_DICOM
 const int ImageFile::EXPORT_FORMAT_DICOM = 5;
 #endif
+const int ImageFile::EXPORT_FORMAT_RAW = 6;
 
 const char* ImageFile::s_aszExportFormatName[] = 
 {
@@ -1667,6 +1668,9 @@ ImageFile::exportImage (const char* const pszFormat, const char* const pszFilena
     return bSuccess;
   }
 #endif
+  else if (iFormatID == EXPORT_FORMAT_RAW)  
+        return writeImageRaw(pszFilename, nxcell, nycell);
+       
   
   sys_error (ERR_SEVERE, "Invalid format %s [ImageFile::exportImage]", pszFormat);
   return false;
@@ -1900,3 +1904,26 @@ ImageFile::writeImageGIF (const char* const outfile, int nxcell, int nycell, dou
 }
 #endif
 
+bool
+ImageFile::writeImageRaw (const char* const outfile, int nxcell, int nycell)
+{
+  FILE *fp;
+  int nx = m_nx;
+  int ny = m_ny;
+  ImageFileArray v = getArray();
+  
+  if ((fp = fopen (outfile, "wb")) == NULL)
+    return false;
+  
+  for (int irow = ny - 1; irow >= 0; irow--) {
+    for (int icol = 0; icol < nx; icol++) {
+      float dens = v[icol][irow];
+     fwrite(&dens, sizeof(float), 1, fp);
+    }
+  }
+  
+  fclose(fp);
+  return true;
+}
+
+
index 8440bd4..86976c7 100644 (file)
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: projections.cpp,v 1.73 2001/03/30 19:17:32 kevin Exp $
+**  $Id: projections.cpp,v 1.74 2001/09/24 09:40:42 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
@@ -156,7 +156,7 @@ Projections::initFromScanner (const Scanner& scanner)
   m_dFocalLength = scanner.focalLength();
   m_dSourceDetectorLength = scanner.sourceDetectorLength();
   m_dViewDiameter = scanner.viewDiameter();
-  m_rotStart = 0;
+  m_rotStart = scanner.offsetView()*scanner.rotInc();
   m_dFanBeamAngle = scanner.fanBeamAngle();
 }
 
@@ -167,6 +167,277 @@ Projections::setNView (int nView)  // used by MPI to reduce # of views
   init (nView, m_nDet);
 }
 
+//  Helical 180 Linear Interpolation.
+//  This member function takes a set of helical scan projections and 
+//  performs a linear interpolation between pairs of complementary rays 
+//  to produce a single projection data set approximating what would be
+//  measured at a single axial plane.
+//  Complementary rays are rays which traverse the same path through the 
+//  phantom in opposite directions.
+//
+//  For parallel beam geometry, a ray with a given gantry angle beta and a
+//  detector iDet will have a complementary ray at beta + pi and nDet-iDet
+//
+//  For equiangular or equilinear beam geometry the complementary ray to
+//  gantry angle beta and fan-beam angle gamma is at 
+//  beta-hat = beta +2*gamma + pi, and gamma-hat =  -gamma.
+//  Note that beta-hat - beta depends on gamma and is not constant.
+//
+//  The algorithm used here is from Crawford and King, Med. Phys. 17(6)
+//  1990 p967; what they called method "C", CSH-HH.  It uses interpolation only
+//  between pairs of complementary rays on either side of an image plane.
+//  Input data must sample gantry angles from zero to  
+//  (2*pi + 2* fan-beam-angle).  The data set produced contains gantry
+//  angles from 0 to Pi+fan-beam-angle.  This is a "halfscan" data set,
+//  which still contains redundant data, and can be used with a half scan 
+//  reconstruction to produce an image.
+//  In this particular implementation a lower triangle from (beta,gamma) =
+//  (0,-fanAngle/2)->(2*fanAngle,-fanAngle/2)->(0,fanAngle/2) contains
+//  zeros, but is actually redundant with data contained in the region
+//  (pi+fanAngle,-fanAngle/2)->(pi+fanAngle, fanAngle/2) ->(pi-fanAngle,
+//  fanAngle/2).  
+//
+int 
+Projections::Helical180LI(int interpolation_view)
+{
+   if (m_geometry == Scanner::GEOMETRY_INVALID) 
+   {
+       std::cerr << "Invalid geometry " << m_geometry << std::endl;
+       return (2);
+   } 
+   else if (m_geometry == Scanner::GEOMETRY_PARALLEL) 
+   {
+       std::cerr << "Helical 180LI not yet implemented for PARALLEL geometry"
+                   << std::endl;
+       return (2);
+   }
+   else if (m_geometry == Scanner::GEOMETRY_EQUILINEAR) 
+   {
+       std::cerr << "Helical 180LI not yet implemented for EQUILINEAR geometry"
+                   << std::endl;
+       return (2);
+   }
+   else if (m_geometry == Scanner::GEOMETRY_EQUIANGULAR) 
+   {
+           return Helical180LI_Equiangular(interpolation_view);
+   }
+   else
+   {
+       std::cerr << "Invalid geometry  in projection data file" << m_geometry
+                   << std::endl;
+       return (2);
+   }
+}
+int
+Projections::Helical180LI_Equiangular(int interpView)
+{
+   double dbeta = m_rotInc; 
+   double dgamma =  m_detInc; 
+   double fanAngle = m_dFanBeamAngle;
+   int offsetView=0;
+   
+   // is there enough data in the data set?  Should have 2(Pi+fanAngle)
+   // coverage minimum
+   if ( m_nView <  static_cast<int>((2*( M_PI + fanAngle ) ) / dbeta) -1 ){
+       std::cerr   << "Data set does not include 360 +2*FanBeamAngle views"
+                   << std::endl;
+       return (1);
+   }
+
+   if (interpView < 0)   // use default position at M_PI+fanAngle
+   {
+       interpView = static_cast<int> ((M_PI+fanAngle)/dbeta);
+   }
+   else
+   {
+       // check if there is M_PI+fanAngle data on either side of the 
+       // of the specified image plane
+       if ( interpView*dbeta < M_PI+fanAngle ||            
+            interpView*dbeta + M_PI + fanAngle > m_nView*dbeta) 
+       {
+           std::cerr << "There isn't PI+fanAngle of data on either side of the requested interpolation view" << std::endl;
+           return(1);
+       }
+       offsetView = interpView - static_cast<int>((M_PI+fanAngle)/dbeta);
+
+   }
+   int last_interp_view = static_cast<int> ((M_PI+fanAngle)/dbeta);
+
+   
+// make a new array for data...
+   class DetectorArray ** newdetarray = new DetectorArray * [last_interp_view+1];
+   for ( int i=0 ; i <= last_interp_view ; i++ ){
+       newdetarray[i] = new DetectorArray (m_nDet);
+       newdetarray[i]->setViewAngle((i+offsetView)*dbeta);
+       DetectorValue* newdetval = (newdetarray[i])->detValues();
+       // and initialize the data to zero
+       for (int j=0; j < m_nDet; j++) 
+           newdetval[j] = 0.;
+   }
+
+   int last_acq_view = 2*last_interp_view;
+   for ( int iView = 0 ; iView <= last_acq_view; iView++) {
+       double beta = iView * dbeta; 
+       
+       for ( int iDet = 0; iDet < m_nDet; iDet++) {
+           double gamma = (iDet -(m_nDet-1)/2)* dgamma ;
+           int newiView, newiDet;
+           if (beta < M_PI+fanAngle) { //if (M_PI +fanAngle - beta > dbeta )  
+               //newbeta = beta; 
+               //newgamma = gamma; 
+               newiDet = iDet; 
+               newiView = iView; 
+           }
+           else // (beta > M_PI+fanAngle)
+           {
+               //newbeta = beta +2*gamma - 180;
+               //newgamma = -gamma;
+               newiDet = -iDet + (m_nDet -1);
+               // newiView = nearest<int>((beta + 2*gamma - M_PI)/dbeta);
+               //newiView = static_cast<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
+               newiView = nearest<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
+           } 
+
+#ifdef DEBUG
+//std::cout << beta << " "<< gamma << " " << newbeta << " " << newgamma <<"    " << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl;
+//std::cout << iView-offsetView << " " << iDet << " " << newiView << " " << newiDet << std::endl;
+#endif
+
+           if (   ( beta > fanAngle - 2*gamma) 
+               && ( beta < 2*M_PI + fanAngle -2*gamma)  )
+          {  // not in region  1 or 8
+               DetectorValue* detval = (m_projData[iView+offsetView])->detValues();
+               DetectorValue* newdetval = (newdetarray[newiView])->detValues();
+               if (   beta > fanAngle - 2*gamma  
+                   && beta <= 2*fanAngle ) {  // in region 2
+                   newdetval[newiDet] += 
+                       (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+                               * detval[iDet];
+               } else if ( beta > 2*fanAngle  
+                          && beta <= M_PI - 2*gamma) {  // in region 3
+                   newdetval[newiDet] += 
+                       (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+                               * detval[iDet];
+               } 
+               else if (   beta > M_PI -2*gamma  
+                        && beta <= M_PI + fanAngle ) {  // in region 4
+                   newdetval[newiDet] += 
+                       (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+                               * detval[iDet];
+               } 
+               else if (   beta > M_PI + fanAngle  
+                        && beta <= M_PI +2*fanAngle -2*gamma) { // in region 5
+                   newdetval[newiDet] += 
+                       (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma)
+                               *detval[iDet];
+               } 
+               else if (   beta > M_PI +2*fanAngle -2*gamma 
+                        && beta <= 2*M_PI) {  // in region 6
+                   newdetval[newiDet] += 
+                       (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma)
+                       *detval[iDet];
+               } 
+               else if (   beta > 2*M_PI 
+                        && beta <= 2*M_PI + fanAngle -2*gamma){ // in region 7
+                   newdetval[newiDet] += 
+                       (2*M_PI - beta -2*gamma + fanAngle)/(M_PI-2*gamma)
+                       *detval[iDet];
+               } 
+               else 
+               { 
+                   ; // outside region of interest
+               }
+           }
+       }
+   }
+   deleteProjData();
+   m_projData = newdetarray;
+   m_nView = last_interp_view+1;
+
+   return (0); 
+}
+// HalfScanFeather:
+// A HalfScan Projection Data Set  for equiangular geometry, 
+// covering gantry angles from 0 to  pi+fanBeamAngle 
+// and fan angle gamma from -fanBeamAngle/2 to fanBeamAngle/2
+// contains redundant information.  If one copy of this data is left as
+// zero, (as in the Helical180LI routine above) overweighting is avoided, 
+// but the discontinuity in the data introduces ringing in the image. 
+// This routine makes a copy of the data and applies a weighting to avoid
+// over-representation, as given in Appendix C of Crawford and King, Med
+// Phys 17 1990, p967.
+int
+Projections::HalfScanFeather(void)
+{
+   double dbeta = m_rotInc; 
+   double dgamma =  m_detInc; 
+   double fanAngle = m_dFanBeamAngle;
+
+// is there enough data?  
+   if ( m_nView !=  static_cast<int>(( M_PI+fanAngle ) / dbeta) +1 ){
+       std::cerr   << "Data set does seem have enough data to be a halfscan data set"  << std::endl;
+       return (1);
+   }
+   if (m_geometry == Scanner::GEOMETRY_INVALID) {
+       std::cerr << "Invalid geometry " << m_geometry << std::endl;
+       return (2);
+   }
+
+   if (m_geometry == Scanner::GEOMETRY_PARALLEL) {
+       std::cerr << "HalfScanFeather not yet implemented for PARALLEL geometry"<< std::endl;
+       return (2);
+   }
+
+   for ( int iView2 = 0 ; iView2 < m_nView; iView2++) {
+       double beta2 = iView2 * dbeta; 
+       for ( int iDet2 = 0; iDet2 < m_nDet; iDet2++) {
+           double gamma2 = (iDet2 -(m_nDet-1)/2)* dgamma ;
+           if ( ( beta2 >= M_PI  - 2*gamma2) ) {  // in redundant data region 
+               int iView1, iDet1;
+               iDet1 =  (m_nDet -1) - iDet2;
+               //iView1 = nearest<int>((beta2 + 2*gamma2 - M_PI)/dbeta);
+               iView1 = nearest<int>(( (iView2*dbeta) 
+                               + 2*(iDet2-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
+
+
+               DetectorValue* detval2 = (m_projData[iView2])->detValues();
+               DetectorValue* detval1 = (m_projData[iView1])->detValues();
+
+               detval1[iDet1] = detval2[iDet2] ;
+
+               double x, w1,w2,beta1, gamma1;
+               beta1= iView1*dbeta; 
+               gamma1 = -gamma2;
+               if ( beta1 <= (fanAngle - 2*gamma1) )
+                   x = beta1 / ( fanAngle - 2*gamma1);
+               else if ( (fanAngle  - 2*gamma1 <= beta1 ) && beta1 <= M_PI - 2*gamma1) 
+                   x = 1; 
+               else if ( (M_PI - 2*gamma1 <= beta1 ) && ( beta1 <=M_PI + fanAngle) )  
+                   x = (M_PI +fanAngle - beta1)/(fanAngle + 2*gamma1);
+               else {
+                   std::cerr << "Shouldn't be here!"<< std::endl;
+                   return(4);
+               }
+               w1 = (3*x - 2*x*x)*x;
+               w2 = 1-w1;
+               detval1[iDet1] *= w1; 
+               detval2[iDet2] *= w2;
+
+           } 
+       }
+   }
+   // heuristic scaling, why this factor?  
+   double scalefactor = m_nView * m_rotInc / M_PI;
+   for ( int iView = 0 ; iView < m_nView; iView++) {
+       DetectorValue* detval = (m_projData[iView])->detValues();
+       for ( int iDet = 0; iDet < m_nDet; iDet++) {
+           detval[iDet] *= scalefactor;
+       }
+   }
+
+   return (0); 
+}
+
 // NAME
 // newProjData
 
index 0319912..a90ca44 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.38 2001/04/02 03:49:52 kevin Exp $
+**  $Id: scanner.cpp,v 1.39 2001/09/24 09:40:42 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
@@ -83,8 +83,10 @@ DetectorArray::~DetectorArray (void)
 */
 
 Scanner::Scanner (const Phantom& phm, const char* const geometryName, 
-                  int nDet, int nView, int nSample, const double rot_anglen, 
-                  const double dFocalLengthRatio, const double dCenterDetectorRatio,
+                  int nDet, int nView, int offsetView, 
+                                 int nSample, const double rot_anglen, 
+                  const double dFocalLengthRatio, 
+                                 const double dCenterDetectorRatio,
                   const double dViewRatio, const double dScanRatio)
 {
   m_fail = false;
@@ -106,6 +108,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
   
   m_nDet     = nDet;
   m_nView    = nView;
+  m_iOffsetView = offsetView;
   m_nSample  = nSample;
   m_dFocalLengthRatio = dFocalLengthRatio;
   m_dCenterDetectorRatio = dCenterDetectorRatio;
@@ -144,7 +147,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
-    m_initPos.angle = 0.0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
     m_detLen += dDetectorArrayEndOffset;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
   if (m_dScanDiameter / 2 >= m_dFocalLength) {
@@ -167,7 +170,6 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     }
   
     m_dFanBeamAngle = dAngle * 2;
-    m_initPos.angle = 0.0;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
     m_initPos.xs2 = m_dXCenter;
@@ -176,7 +178,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
-    m_initPos.angle = 0.0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
     if (m_dScanDiameter / 2 > m_dFocalLength) {
       m_fail = true;
@@ -202,7 +204,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
     
     m_dFanBeamAngle = dAngle * 2;
-    m_initPos.angle = 0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
     m_initPos.xs2 = m_dXCenter;
@@ -277,23 +279,22 @@ Scanner::convertGeometryNameToID (const char* const geomName)
 void
 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP)
 {
-  collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP);
+  collectProjections (proj, phm, m_startView, proj.nView(), m_iOffsetView, true, trace, pSGP);
 }
 
 void
-Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-                             bool bStoreAtViewPosition, const int trace, SGP* pSGP)
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, const int iOffsetView,  bool bStoreAtViewPosition, const int trace, SGP* pSGP)
 {
   int iStorageOffset = (bStoreAtViewPosition ? iStartView : 0);
-  collectProjections (proj, phm, iStartView, iNumViews, iStorageOffset, trace, pSGP);
+  collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP);
 }
 
 void
-Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-                             int iStorageOffset, const int trace, SGP* pSGP)
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, const int iOffsetView, 
+      int iStorageOffset, const int trace, SGP* pSGP)
 {
   m_trace = trace;
-  double start_angle = iStartView * proj.rotInc();
+  double start_angle = (iStartView + iOffsetView) * proj.rotInc();
   
   // Calculate initial rotation matrix 
   GRFMTX_2D rotmtx_initial, temp;
index 9248247..2d94fa4 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: ctsim.cpp,v 1.101 2001/03/30 19:17:32 kevin Exp $
+**  $Id: ctsim.cpp,v 1.102 2001/09/24 09:40:42 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
@@ -70,7 +70,7 @@
 #endif
 #endif
 
-static const char* rcsindent = "$Id: ctsim.cpp,v 1.101 2001/03/30 19:17:32 kevin Exp $";
+static const char* rcsindent = "$Id: ctsim.cpp,v 1.102 2001/09/24 09:40:42 kevin Exp $";
 
 struct option CTSimApp::ctsimOptions[] = 
 {
@@ -160,7 +160,7 @@ CTSimApp::OnInit()
   int xDisplay, yDisplay;
   ::wxDisplaySize (&xDisplay, &yDisplay);
   m_pFrame = new MainFrame(m_docManager, (wxFrame *) NULL, -1, "CTSim", wxPoint(0, 0), 
-    wxSize(nearest<int>(xDisplay * .75), nearest<int>(yDisplay * .755)), wxDEFAULT_FRAME_STYLE);
+    wxSize(nearest<int>(xDisplay * .25), nearest<int>(yDisplay * .255)), wxDEFAULT_FRAME_STYLE);
   
   setIconForFrame (m_pFrame);
   m_pFrame->Centre(wxBOTH);
index 8525c62..75596d4 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: dialogs.cpp,v 1.54 2001/03/30 19:17:32 kevin Exp $
+**  $Id: dialogs.cpp,v 1.55 2001/09/24 09:40:42 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
@@ -639,7 +639,7 @@ DialogGetRasterParameters::getViewRatio ()
 
 
 DialogGetProjectionParameters::DialogGetProjectionParameters 
-   (wxWindow* pParent, int iDefaultNDet, int iDefaultNView, int iDefaultNSamples, 
+   (wxWindow* pParent, int iDefaultNDet, int iDefaultNView, int iDefaultOffsetView, int iDefaultNSamples, 
     double dDefaultRotAngle, double dDefaultFocalLength, double dDefaultCenterDetectorLength,
     double dDefaultViewRatio, double dDefaultScanRatio, int iDefaultGeometry, int iDefaultTrace)
 : wxDialog (pParent, -1, _T("Projection Parameters"), wxDefaultPosition, wxDefaultSize, wxDEFAULT_DIALOG_STYLE | wxCAPTION)
@@ -719,6 +719,13 @@ DialogGetProjectionParameters::DialogGetProjectionParameters
     m_pTextCtrlRotAngle = new wxTextCtrl (this, -1, osRotAngle.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
     pText2Sizer->Add (new wxStaticText (this, -1, "Rotation Angle (Fraction of circle)"), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
     pText2Sizer->Add (m_pTextCtrlRotAngle, 0, wxALIGN_CENTER_VERTICAL);
+                       
+    std::ostringstream osOffsetView;
+    osOffsetView << iDefaultOffsetView;
+    m_pTextCtrlOffsetView = new wxTextCtrl (this, -1, osOffsetView.str().c_str(), wxDefaultPosition, wxSize(100, 25), 0);
+    pText2Sizer->Add (new wxStaticText (this, -1, "Gantry offset in units of 'views' "), 0, wxALIGN_RIGHT | wxALIGN_CENTER_VERTICAL);
+    pText2Sizer->Add (m_pTextCtrlOffsetView, 0, wxALIGN_CENTER_VERTICAL);
+
   }
   pGridSizer->Add (pText2Sizer);
 
@@ -772,6 +779,16 @@ DialogGetProjectionParameters::getNView ()
     return (m_iDefaultNView);
 }
 
+unsigned int
+DialogGetProjectionParameters::getOffsetView ()
+{
+       wxString strCtrl = m_pTextCtrlOffsetView->GetValue();
+       unsigned long lValue;
+       if (strCtrl.ToULong (&lValue))
+         return lValue;
+       else
+         return (m_iDefaultOffsetView);
+}
 
 unsigned int
 DialogGetProjectionParameters::getNSamples ()
index 9faf10f..abb71e5 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: dialogs.h,v 1.36 2001/03/30 19:17:32 kevin Exp $
+**  $Id: dialogs.h,v 1.37 2001/09/24 09:40:42 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
@@ -192,13 +192,14 @@ class DialogGetProjectionParameters : public wxDialog
 {
  public:
     DialogGetProjectionParameters (wxWindow* pParent, int iDefaultNDet = 0, 
-      int iDefaultNView = 0, int iDefaultNSamples = 1, double dDefaultRotAngle = 1., 
+      int iDefaultNView = 0, int iDefaultOffsetVIew = 0, int iDefaultNSamples = 1, double dDefaultRotAngle = 1., 
       double dDefaultFocalLength = 1, double dDefaultCenterDetectorLength = 1, double dDefaultViewRatio = 1., 
       double dDefaultScanRatio = 1., int iDefaultGeometry = Scanner::GEOMETRY_PARALLEL, int iDefaultTrace = Trace::TRACE_NONE);
     ~DialogGetProjectionParameters ();
 
     unsigned int getNDet ();
     unsigned int getNView ();
+       unsigned int getOffsetView ();
     unsigned int getNSamples ();
     int getTrace ();
 
@@ -212,6 +213,7 @@ class DialogGetProjectionParameters : public wxDialog
  private:
     wxTextCtrl* m_pTextCtrlNDet;
     wxTextCtrl* m_pTextCtrlNView;
+    wxTextCtrl* m_pTextCtrlOffsetView;
     wxTextCtrl* m_pTextCtrlNSamples;
     wxTextCtrl* m_pTextCtrlRotAngle;
     wxTextCtrl* m_pTextCtrlFocalLength;
@@ -223,6 +225,7 @@ class DialogGetProjectionParameters : public wxDialog
 
     int m_iDefaultNDet;
     int m_iDefaultNView;
+       int m_iDefaultOffsetView;
     int m_iDefaultNSamples;
     int m_iDefaultTrace;
     int m_iDefaultGeometry;
index 76b0ae4..7624beb 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: dlgprojections.cpp,v 1.23 2001/02/22 11:05:38 kevin Exp $
+**  $Id: dlgprojections.cpp,v 1.24 2001/09/24 09:40:42 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,7 +151,7 @@ ProjectionsDialog::showView (int iViewNumber)
                if (m_iTrace >= Trace::TRACE_PLOT)
                        m_pSGP->setViewport (0, 0, 0.66, 1);
                ::wxYield();        // update the display
-               m_rScanner.collectProjections (m_rProjections, m_rPhantom, iViewNumber, 1, true, m_iTrace, m_pSGP);
+               m_rScanner.collectProjections (m_rProjections, m_rPhantom, iViewNumber, 1, m_rScanner.offsetView(), true, m_iTrace, m_pSGP);
                ::wxYield();        // update the display
                if (m_iTrace >= Trace::TRACE_PLOT) {
                        const DetectorArray& detArray = m_rProjections.getDetectorArray (iViewNumber);
index f967269..d61c422 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2001 Kevin Rosenberg
 **
-**  $Id: threadproj.cpp,v 1.15 2001/03/07 21:18:50 kevin Exp $
+**  $Id: threadproj.cpp,v 1.16 2001/09/24 09:40:42 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
 //
 /////////////////////////////////////////////////////////////////////
 
-ProjectorSupervisorThread::ProjectorSupervisorThread (PhantomFileView* pProjView, int iNDet, int iNView, 
+ProjectorSupervisorThread::ProjectorSupervisorThread (PhantomFileView* pProjView, int iNDet, int iNView, int iOffsetView,
    const char* pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
    double dViewRatio, double dScanRatio, const char* const pszLabel)
-: SupervisorThread(), m_pPhantomView(pProjView), m_iNDet(iNDet), m_iNView(iNView), m_strGeometry(pszGeometry), 
+: SupervisorThread(), m_pPhantomView(pProjView), m_iNDet(iNDet), m_iNView(iNView), m_iOffsetView(iOffsetView), m_strGeometry(pszGeometry), 
   m_iNSample(iNSample), m_dRotation(dRotation), m_dFocalLength(dFocalLength), m_dCenterDetectorLength(dCenterDetectorLength),
   m_dViewRatio(dViewRatio), m_dScanRatio(dScanRatio), m_strLabel(pszLabel)
 {
@@ -66,7 +66,7 @@ ProjectorSupervisorThread::ProjectorSupervisorThread (PhantomFileView* pProjView
 wxThread::ExitCode
 ProjectorSupervisorThread::Entry()
 {
-  ProjectorSupervisor projSupervisor (this, m_pPhantomView, m_iNDet, m_iNView, 
+  ProjectorSupervisor projSupervisor (this, m_pPhantomView, m_iNDet, m_iNView, m_iOffsetView, 
    m_strGeometry.c_str(), m_iNSample, m_dRotation, m_dFocalLength, m_dCenterDetectorLength, m_dViewRatio, m_dScanRatio, m_strLabel.c_str());
 
   projSupervisor.start();
@@ -102,17 +102,17 @@ ProjectorSupervisorThread::OnExit()
 //
 /////////////////////////////////////////////////////////////////////
 
-ProjectorSupervisor::ProjectorSupervisor (SupervisorThread* pThread, PhantomFileView* pPhantomView, int iNDet, int iNView, 
+ProjectorSupervisor::ProjectorSupervisor (SupervisorThread* pThread, PhantomFileView* pPhantomView, int iNDet, int iNView, int iOffsetView, 
    const char* pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
    double dViewRatio, double dScanRatio, const char* const pszLabel)
     : BackgroundSupervisor (pThread, pPhantomView->GetFrame(), pPhantomView->GetDocument(), "Projecting", iNView),
       m_pPhantomView(pPhantomView), m_pPhantomDoc(pPhantomView->GetDocument()), 
-      m_iNDet(iNDet), m_iNView(iNView), m_pszGeometry(pszGeometry), m_iNSample(iNSample), 
+      m_iNDet(iNDet), m_iNView(iNView), m_iOffsetView(iOffsetView), m_pszGeometry(pszGeometry), m_iNSample(iNSample), 
       m_dRotation(dRotation), m_dFocalLength(dFocalLength), m_dCenterDetectorLength(dCenterDetectorLength), 
       m_dViewRatio(dViewRatio), m_dScanRatio(dScanRatio), m_pszLabel(pszLabel)
 {
   m_pScanner = new Scanner (m_pPhantomDoc->getPhantom(), m_pszGeometry, m_iNDet,
-                  m_iNView, m_iNSample, m_dRotation, m_dFocalLength, m_dCenterDetectorLength, m_dViewRatio, m_dScanRatio);
+                  m_iNView, m_iOffsetView, m_iNSample, m_dRotation, m_dFocalLength, m_dCenterDetectorLength, m_dViewRatio, m_dScanRatio);
 
   m_vecpChildProjections.reserve (getNumWorkers());
   for (int iThread = 0; iThread < getNumWorkers(); iThread++) {
@@ -138,7 +138,7 @@ ProjectorSupervisor::createWorker (int iThread, int iStartUnit, int iNumUnits)
 {
    ProjectorWorker* pThread = new ProjectorWorker (this, iThread, iStartUnit, iNumUnits);
    m_vecpChildProjections[iThread]->setNView (iNumUnits);
-   pThread->SetParameters (m_pPhantomView, m_vecpChildProjections[iThread], m_pScanner, m_iNDet, m_iNView, 
+   pThread->SetParameters (m_pPhantomView, m_vecpChildProjections[iThread], m_pScanner, m_iNDet, m_iNView, m_iOffsetView, 
      m_pszGeometry, m_iNSample, m_dRotation, m_dFocalLength, m_dCenterDetectorLength, m_dViewRatio, m_dScanRatio);
 
    return pThread;
@@ -198,7 +198,7 @@ ProjectorSupervisor::getProjections()
 
 void
 ProjectorWorker::SetParameters (PhantomFileView* pPhantomView, Projections* pProjections, Scanner* pScanner,
- int iNDet, int iView,
+ int iNDet, int iView, int iOffsetView, 
  const char* pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
  double dViewRatio, double dScanRatio)
 {
index 49212f4..a4f6993 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2001 Kevin Rosenberg
 **
-**  $Id: threadproj.h,v 1.6 2001/03/04 04:27:55 kevin Exp $
+**  $Id: threadproj.h,v 1.7 2001/09/24 09:40:42 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
@@ -48,6 +48,7 @@ private:
   PhantomFileView* m_pPhantomView;
   const int m_iNDet;
   const int m_iNView;
+  const int m_iOffsetView;
   const std::string m_strGeometry;
   const int m_iNSample;
   const double m_dRotation;
@@ -58,7 +59,7 @@ private:
   const std::string m_strLabel;
 
 public:
-  ProjectorSupervisorThread(PhantomFileView* pProjView, int iNDet, int iNView, 
+  ProjectorSupervisorThread(PhantomFileView* pProjView, int iNDet, int iNView, int iOffsetView, 
    const char* pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
    double dViewRatio, double dScanRatio, const char* const pszLabel);
 
@@ -79,6 +80,7 @@ private:
   
   const int m_iNDet;
   const int m_iNView;
+  const int m_iOffsetView;
   const char* const m_pszGeometry;
   const int m_iNSample;
   const double m_dRotation;
@@ -90,7 +92,7 @@ private:
 
 
 public:
-   ProjectorSupervisor (SupervisorThread* pThread, PhantomFileView* pProjView, int iNDet, int iNView, 
+   ProjectorSupervisor (SupervisorThread* pThread, PhantomFileView* pProjView, int iNDet, int iNView,  int iOffsetView,
    const char* pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
    double dViewRatio, double dScanRatio, const char* const pszLabel);
 
@@ -114,6 +116,7 @@ private:
   Scanner* m_pScanner;
   int m_iNDet;
   int m_iNView;
+  int m_iOffsetView;
   const char* m_pszGeometry;
   int m_iNSample;
   double m_dRotation;
@@ -129,7 +132,7 @@ public:
   {}
   
   void SetParameters (PhantomFileView* pPhantomFile, Projections* pProjections, Scanner* pScanner, 
-   int iNDet, int iView, 
+   int iNDet, int iView, int iOffsetView,
    const char* const pszGeometry, int iNSample, double dRotation, double dFocalLength, double dCenterDetectorLength,
    double dViewRatio, double dScanRatio);
 
index afceae4..f78f9e3 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: views.cpp,v 1.145 2001/03/30 21:01:15 kevin Exp $
+**  $Id: views.cpp,v 1.146 2001/09/24 09:40:42 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
@@ -2003,6 +2003,7 @@ PhantomFileView::PhantomFileView()
   m_iDefaultNView = 320;
   m_iDefaultNSample = 2;
 #endif
+  m_iDefaultOffsetView = 0;
   m_dDefaultRotation = 1;
   m_dDefaultFocalLength = 2;
   m_dDefaultCenterDetectorLength = 2;
@@ -2051,7 +2052,7 @@ void
 PhantomFileView::OnProjections (wxCommandEvent& event)
 {
   DialogGetProjectionParameters dialogProjection (getFrameForChild(), 
-    m_iDefaultNDet, m_iDefaultNView, m_iDefaultNSample, m_dDefaultRotation, 
+    m_iDefaultNDet, m_iDefaultNView, m_iDefaultOffsetView, m_iDefaultNSample, m_dDefaultRotation, 
     m_dDefaultFocalLength, m_dDefaultCenterDetectorLength, m_dDefaultViewRatio, m_dDefaultScanRatio, 
     m_iDefaultGeometry, m_iDefaultTrace);
   int retVal = dialogProjection.ShowModal();
@@ -2060,6 +2061,7 @@ PhantomFileView::OnProjections (wxCommandEvent& event)
   
   m_iDefaultNDet = dialogProjection.getNDet();
   m_iDefaultNView = dialogProjection.getNView();
+  m_iDefaultOffsetView = dialogProjection.getOffsetView();
   m_iDefaultNSample = dialogProjection.getNSamples();
   m_iDefaultTrace = dialogProjection.getTrace();
   m_dDefaultRotation = dialogProjection.getRotAngle();
@@ -2076,7 +2078,7 @@ PhantomFileView::OnProjections (wxCommandEvent& event)
     return;
   
   const Phantom& rPhantom = GetDocument()->getPhantom();
-  Scanner theScanner (rPhantom, sGeometry.c_str(), m_iDefaultNDet, m_iDefaultNView, m_iDefaultNSample, 
+  Scanner theScanner (rPhantom, sGeometry.c_str(), m_iDefaultNDet, m_iDefaultNView, m_iDefaultOffsetView, m_iDefaultNSample, 
     dRotationRadians, m_dDefaultFocalLength, m_dDefaultCenterDetectorLength, m_dDefaultViewRatio, m_dDefaultScanRatio);
   if (theScanner.fail()) {
     wxString msg = "Failed making scanner\n";
@@ -2087,13 +2089,18 @@ PhantomFileView::OnProjections (wxCommandEvent& event)
   }
   
   std::ostringstream os;
-  os << "Projections for " << rPhantom.name() << ": nDet=" << m_iDefaultNDet 
-    << ", nView=" << m_iDefaultNView << ", nSamples=" << m_iDefaultNSample 
-    << ", RotAngle=" << m_dDefaultRotation << ", FocalLengthRatio=" << m_dDefaultFocalLength 
+  os << "Projections for " << rPhantom.name() 
+       << ": nDet=" << m_iDefaultNDet 
+    << ", nView=" << m_iDefaultNView 
+       << ", gantry offset=" << m_iDefaultOffsetView 
+       << ", nSamples=" << m_iDefaultNSample 
+    << ", RotAngle=" << m_dDefaultRotation 
+       << ", FocalLengthRatio=" << m_dDefaultFocalLength 
     << ", CenterDetectorLengthRatio=" << m_dDefaultCenterDetectorLength
-    << ", ViewRatio=" << m_dDefaultViewRatio << ", ScanRatio=" << m_dDefaultScanRatio 
-    << ", Geometry=" << sGeometry.c_str() << ", FanBeamAngle=" << 
-    convertRadiansToDegrees (theScanner.fanBeamAngle());
+    << ", ViewRatio=" << m_dDefaultViewRatio 
+       << ", ScanRatio=" << m_dDefaultScanRatio 
+    << ", Geometry=" << sGeometry.c_str() 
+       << ", FanBeamAngle=" << convertRadiansToDegrees (theScanner.fanBeamAngle());
   
   Timer timer;
   Projections* pProj = NULL;
@@ -2118,7 +2125,7 @@ PhantomFileView::OnProjections (wxCommandEvent& event)
 #if HAVE_WXTHREADS
     if (theApp->getUseBackgroundTasks()) {
       ProjectorSupervisorThread* pProjector = new ProjectorSupervisorThread (this, m_iDefaultNDet,
-        m_iDefaultNView, sGeometry.c_str(), m_iDefaultNSample, dRotationRadians,
+        m_iDefaultNView, m_iDefaultOffsetView, sGeometry.c_str(), m_iDefaultNSample, dRotationRadians,
         m_dDefaultFocalLength, m_dDefaultCenterDetectorLength, m_dDefaultViewRatio, m_dDefaultScanRatio, os.str().c_str());
       if (pProjector->Create() != wxTHREAD_NO_ERROR) {
         sys_error (ERR_SEVERE, "Error creating projector thread");
@@ -2135,7 +2142,8 @@ PhantomFileView::OnProjections (wxCommandEvent& event)
       pProj->initFromScanner (theScanner);
       wxProgressDialog dlgProgress (wxString("Projection"), wxString("Projection Progress"), pProj->nView() + 1, getFrameForChild(), wxPD_CAN_ABORT );
       for (int i = 0; i < pProj->nView(); i++) {
-        theScanner.collectProjections (*pProj, rPhantom, i, 1, true, m_iDefaultTrace);
+        //theScanner.collectProjections (*pProj, rPhantom, i, 1, true, m_iDefaultTrace);
+        theScanner.collectProjections (*pProj, rPhantom, i, 1, theScanner.offsetView(), true, m_iDefaultTrace);
         if (! dlgProgress.Update (i+1)) {
           delete pProj;
           return;
index 01f8945..40cc6fd 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: views.h,v 1.55 2001/03/30 21:01:15 kevin Exp $
+**  $Id: views.h,v 1.56 2001/09/24 09:40:42 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
@@ -307,6 +307,7 @@ private:
   
   int m_iDefaultNDet;
   int m_iDefaultNView;
+  int m_iDefaultOffsetView;
   int m_iDefaultNSample;
   int m_iDefaultGeometry;
   int m_iDefaultTrace;
index e214688..951ec21 100644 (file)
@@ -16,6 +16,7 @@ install-exec-hook:
        rm -f $(bindir)/ifinfo
        rm -f $(bindir)/phm2if
        rm -f $(bindir)/phm2pj
+       rm -f $(bindir)/phm2helix
        rm -f $(bindir)/pj2if
        rm -f $(bindir)/pjinfo
        rm -f $(bindir)/pjrec
@@ -25,18 +26,19 @@ install-exec-hook:
        ln -s $(bindir)/ctsimtext $(bindir)/ifinfo
        ln -s $(bindir)/ctsimtext $(bindir)/phm2if
        ln -s $(bindir)/ctsimtext $(bindir)/phm2pj
+       ln -s $(bindir)/ctsimtext $(bindir)/phm2helix
        ln -s $(bindir)/ctsimtext $(bindir)/pj2if
        ln -s $(bindir)/ctsimtext $(bindir)/pjinfo
        ln -s $(bindir)/ctsimtext $(bindir)/pjrec
 
-ctsimtext_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp
+ctsimtext_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp phm2helix.cpp pjHinterp.cpp
 ctsimtext_LDADD=@ctlibs@
 ctsimtext_DEPENDENCIES=$(SOURCE_DEPEND)
 
 realclean:
        rm -f *.pgm *.if *~ *.pj
 
-ctsimtext_lam_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp mpiworld.cpp 
+ctsimtext_lam_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp mpiworld.cpp phm2helix.cpp pjHinterp.cpp
 ctsimtext_lam_LDADD=@ctlamlibs@
 
 if USE_LAM
index 0743b0d..132410b 100644 (file)
@@ -97,11 +97,11 @@ INCLUDES = @my_includes@
 @HAVE_SGP_TRUE@SOURCE_DEPEND = ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a ../libctgraphics/libctgraphics.a
 @HAVE_SGP_FALSE@SOURCE_DEPEND = ../include/ct.h ../libctsim/libctsim.a ../libctsupport/libctsupport.a
 
-ctsimtext_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp
+ctsimtext_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp phm2helix.cpp pjHinterp.cpp
 ctsimtext_LDADD = @ctlibs@
 ctsimtext_DEPENDENCIES = $(SOURCE_DEPEND)
 
-ctsimtext_lam_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp mpiworld.cpp 
+ctsimtext_lam_SOURCES = ctsimtext.cpp if1.cpp if2.cpp ifinfo.cpp ifexport.cpp phm2if.cpp phm2pj.cpp pj2if.cpp pjinfo.cpp pjrec.cpp nographics.cpp mpiworld.cpp phm2helix.cpp pjHinterp.cpp
 ctsimtext_lam_LDADD = @ctlamlibs@
 
 @USE_LAM_TRUE@CC_LAM = $(lamdir)/bin/hcp
@@ -119,10 +119,12 @@ CPPFLAGS = @CPPFLAGS@
 LDFLAGS = @LDFLAGS@
 LIBS = @LIBS@
 ctsimtext_lam_OBJECTS =  ctsimtext.o if1.o if2.o ifinfo.o ifexport.o \
-phm2if.o phm2pj.o pj2if.o pjinfo.o pjrec.o nographics.o mpiworld.o
+phm2if.o phm2pj.o pj2if.o pjinfo.o pjrec.o nographics.o mpiworld.o \
+phm2helix.o pjHinterp.o
 ctsimtext_lam_LDFLAGS = 
 ctsimtext_OBJECTS =  ctsimtext.o if1.o if2.o ifinfo.o ifexport.o \
-phm2if.o phm2pj.o pj2if.o pjinfo.o pjrec.o nographics.o
+phm2if.o phm2pj.o pj2if.o pjinfo.o pjrec.o nographics.o phm2helix.o \
+pjHinterp.o
 ctsimtext_LDFLAGS = 
 SCRIPTS =  $(bin_SCRIPTS)
 
@@ -138,8 +140,9 @@ DISTFILES = $(DIST_COMMON) $(SOURCES) $(HEADERS) $(TEXINFOS) $(EXTRA_DIST)
 TAR = gtar
 GZIP_ENV = --best
 DEP_FILES =  .deps/ctsimtext.P .deps/if1.P .deps/if2.P .deps/ifexport.P \
-.deps/ifinfo.P .deps/mpiworld.P .deps/nographics.P .deps/phm2if.P \
-.deps/phm2pj.P .deps/pj2if.P .deps/pjinfo.P .deps/pjrec.P
+.deps/ifinfo.P .deps/mpiworld.P .deps/nographics.P .deps/phm2helix.P \
+.deps/phm2if.P .deps/phm2pj.P .deps/pj2if.P .deps/pjHinterp.P \
+.deps/pjinfo.P .deps/pjrec.P
 SOURCES = $(ctsimtext_lam_SOURCES) $(ctsimtext_SOURCES)
 OBJECTS = $(ctsimtext_lam_OBJECTS) $(ctsimtext_OBJECTS)
 
@@ -405,6 +408,7 @@ install-exec-hook:
        rm -f $(bindir)/ifinfo
        rm -f $(bindir)/phm2if
        rm -f $(bindir)/phm2pj
+       rm -f $(bindir)/phm2helix
        rm -f $(bindir)/pj2if
        rm -f $(bindir)/pjinfo
        rm -f $(bindir)/pjrec
@@ -414,6 +418,7 @@ install-exec-hook:
        ln -s $(bindir)/ctsimtext $(bindir)/ifinfo
        ln -s $(bindir)/ctsimtext $(bindir)/phm2if
        ln -s $(bindir)/ctsimtext $(bindir)/phm2pj
+       ln -s $(bindir)/ctsimtext $(bindir)/phm2helix
        ln -s $(bindir)/ctsimtext $(bindir)/pj2if
        ln -s $(bindir)/ctsimtext $(bindir)/pjinfo
        ln -s $(bindir)/ctsimtext $(bindir)/pjrec
index 8237a41..c31ba0e 100644 (file)
@@ -44,6 +44,12 @@ ifinfo.exe: ifinfo.obj
 if2img.exe: if2img.obj
         link if2img.obj $(LDFLAGS)
 
+pjHinterp.exe: pjHinterp.obj
+        link pjHinterp.obj $(LDFLAGS)
+
+phm2helix.exe: phm2helix.obj
+        link phm2helix.obj $(LDFLAGS)
+
 clean:
         echo dummy > a.obj
         echo dummy > a.exe
index 75f941f..54236af 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ctsimtext.cpp,v 1.19 2001/04/02 22:50:25 kevin Exp $
+**  $Id: ctsimtext.cpp,v 1.20 2001/09/24 09:40:42 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
@@ -45,7 +45,7 @@ extern "C" {
 // If linked to ctsimtext, but executed as another name, eg pjrec, then program will use that
 // linked name as name of function.
 
-static const char* const g_szIdStr = "$Id: ctsimtext.cpp,v 1.19 2001/04/02 22:50:25 kevin Exp $";
+static const char* const g_szIdStr = "$Id: ctsimtext.cpp,v 1.20 2001/09/24 09:40:42 kevin Exp $";
 static const char* const s_szProgramName = "ctsimtext";
 static const char* const s_szProgramName2 = "ctsimtext.exe";
 static const char* const s_szProgramName3 = "ctsimtext-lam";
@@ -56,6 +56,8 @@ extern int ifexport_main (int argc, char* const argv[]);
 extern int ifinfo_main (int argc, char* const argv[]);
 extern int phm2if_main (int argc, char* const argv[]);
 extern int phm2pj_main (int argc, char* const argv[]);
+extern int phm2helix_main (int argc, char* const argv[]);
+extern int pjHinterp_main (int argc, char* const argv[]);
 extern int pj2if_main (int argc, char* const argv[]);
 extern int pjinfo_main (int argc, char* const argv[]);
 extern int pjrec_main (int argc, char* const argv[]);
@@ -82,6 +84,9 @@ ctsimtext_usage (const char *program)
   std::cout << "  pjrec         Projection reconstruction\n";
   std::cout << "  phm2if        Convert a geometric phantom into an imagefile\n";
   std::cout << "  phm2pj        Take projections of a phantom object\n";
+  std::cout << "  phm2helix     Take projections of a phantom object\n";
+  std::cout << "  pjHinterp    Interpolate helical projections of a phantom object\n";
+
 }
 
 void 
@@ -95,8 +100,10 @@ interactive_usage ()
   std::cout << "  if2           Dual image file conversions\n";
   std::cout << "  phm2if        Convert a geometric phantom into an imagefile\n";
   std::cout << "  phm2pj        Take projections of a phantom object\n";
+  std::cout << "  phm2helix        Take projections of a phantom object\n";
   std::cout << "  pjinfo        Projection file information\n";
   std::cout << "  pj2if         Convert an projection file into an imagefile\n";
+  std::cout << "  pjHinterp      Interpolate helical projections of a phantom object\n";
   std::cout << "  pjrec         Projection reconstruction\n";
   std::cout << "  quit          Quits shell\n";
   std::cout << "All functions accept --help as parameter for online help\n\n";
@@ -250,6 +257,10 @@ processCommand (int argc, char* const argv[])
       return phm2if_main (argc, argv);
     else if (strcasecmp (pszFunction, "phm2pj") == 0)
       return phm2pj_main (argc, argv);
+       else if (strcasecmp (pszFunction, "phm2helix") == 0)
+         return phm2helix_main (argc, argv);
+       else if (strcasecmp (pszFunction, "pjHinterp") == 0)
+         return pjHinterp_main (argc, argv);
     else if (strcasecmp (pszFunction, "pj2if") == 0)
       return pj2if_main (argc, argv);
     else if (strcasecmp (pszFunction, "pjinfo") == 0)
index 8760332..d629422 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: ifexport.cpp,v 1.3 2001/01/09 22:31:47 kevin Exp $
+**  $Id: ifexport.cpp,v 1.4 2001/09/24 09:40:42 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
@@ -48,7 +48,7 @@ static struct option my_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: ifexport.cpp,v 1.3 2001/01/09 22:31:47 kevin Exp $";
+static const char* g_szIdStr = "$Id: ifexport.cpp,v 1.4 2001/09/24 09:40:42 kevin Exp $";
 
 enum { O_AUTO_FULL, O_AUTO_STD0_1, O_AUTO_STD0_5, O_AUTO_STD1, O_AUTO_STD2, O_AUTO_STD3 };
 static const char O_AUTO_FULL_STR[]="full";
@@ -63,11 +63,12 @@ static const char O_CENTER_MEAN_STR[]="mean";
 static const char O_CENTER_MODE_STR[]="mode";
 static const char O_CENTER_MEDIAN_STR[]="median";
 
-enum { O_FORMAT_PNG, O_FORMAT_PNG16, O_FORMAT_PGM, O_FORMAT_PGMASC };
+enum { O_FORMAT_PNG, O_FORMAT_PNG16, O_FORMAT_PGM, O_FORMAT_PGMASC, O_FORMAT_RAW };
 static const char O_FORMAT_PNG_STR[]="png" ;
 static const char O_FORMAT_PNG16_STR[]="png16" ;
 static const char O_FORMAT_PGM_STR[]="pgm";
 static const char O_FORMAT_PGMASC_STR[]="pgmasc";
+static const char O_FORMAT_RAW_STR[]="raw";
 
 void 
 ifexport_usage (const char *program)
@@ -212,6 +213,8 @@ ifexport_main (int argc, char *const argv[])
          else if (strcmp(optarg, O_FORMAT_PNG16_STR) == 0)
            opt_format = O_FORMAT_PNG16;
 #endif
+         else if (strcmp(optarg, O_FORMAT_RAW_STR) == 0)
+               opt_format = O_FORMAT_RAW;
 #if HAVE_GIF
          else if (strcmp(optarg, O_FORMAT_GIF_STR) == 0)
            opt_format = O_FORMAT_GIF;
@@ -347,6 +350,8 @@ ifexport_main (int argc, char *const argv[])
   else if (opt_format == O_FORMAT_GIF) 
     im.writeImageGIF (out_file, opt_scale, opt_scale, densmin, densmax);
 #endif
+  else if (opt_format == O_FORMAT_RAW)
+       im.writeImageRaw (out_file, opt_scale, opt_scale);
   else
     {
       sys_error (ERR_SEVERE, "Internal Error: Invalid format mode %d", opt_format);
diff --git a/tools/phm2helix.cpp b/tools/phm2helix.cpp
new file mode 100644 (file)
index 0000000..8f7f185
--- /dev/null
@@ -0,0 +1,361 @@
+/*****************************************************************************
+** FILE IDENTIFICATION
+**
+**   Name:          phm2helix.cpp
+**   Purpose:       Take projections of a phantom object
+**   Programmer:    Ian Kay
+**   Date Started:  Aug 2001
+**
+**  This is part of the CTSim program
+**  Copyright (C) 1983-2000 Kevin Rosenberg
+**
+**  $Id: phm2helix.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $
+**
+**  This program is free software; you can redistribute it and/or modify
+**  it under the terms of the GNU General Public License (version 2) as
+**  published by the Free Software Foundation.
+**
+**  This program is distributed in the hope that it will be useful,
+**  but WITHOUT ANY WARRANTY; without even the implied warranty of
+**  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+**  GNU General Public License for more details.
+**
+**  You should have received a copy of the GNU General Public License
+**  along with this program; if not, write to the Free Software
+**  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+******************************************************************************/
+
+#include "ct.h"
+#include "timer.h"
+
+
+enum { O_PHANTOMPROG, O_PHMFILE,  O_DESC, O_NRAY, O_ROTANGLE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
+  O_VIEW_RATIO, O_SCAN_RATIO, O_OFFSETVIEW,  O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
+
+static struct option phm2helix_options[] = 
+{
+  {"phantom", 1, 0, O_PHANTOMPROG},
+  {"desc", 1, 0, O_DESC},
+  {"nray", 1, 0, O_NRAY},
+  {"rotangle", 1, 0, O_ROTANGLE},
+  {"geometry", 1, 0, O_GEOMETRY},
+  {"focal-length", 1, 0, O_FOCAL_LENGTH},
+  {"center-detector-length", 1, 0, O_CENTER_DETECTOR_LENGTH},
+  {"view-ratio", 1, 0, O_VIEW_RATIO},
+  {"scan-ratio", 1, 0, O_SCAN_RATIO},
+  {"offsetview", 1, 0, O_OFFSETVIEW},
+  {"trace", 1, 0, O_TRACE},
+  {"verbose", 0, 0, O_VERBOSE},
+  {"help", 0, 0, O_HELP},
+  {"debug", 0, 0, O_DEBUG},
+  {"version", 0, 0, O_VERSION},
+  {"phmfile", 1, 0, O_PHMFILE},
+  {0, 0, 0, 0}
+};
+
+static const char* g_szIdStr = "$Id: phm2helix.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $";
+
+
+void 
+phm2helix_usage (const char *program)
+{
+         std::cout << "usage: " << fileBasename(program) << " outfile ndet nview phmprog [OPTIONS]\n";
+         std::cout << "Calculate (projections) through time varying phantom object  \n\n";
+         std::cout << "     outfile          Name of output file for projectsions\n";
+         std::cout << "     ndet             Number of detectors\n";
+         std::cout << "     nview            Number of rotated views\n";
+         std::cout << "     phmprog          Name of phm generation executable\n";
+         std::cout << "     --phmfile        Temp phantom file name \n";
+         std::cout << "     --desc           Description of raysum\n";
+         std::cout << "     --nray           Number of rays per detector (default = 1)\n";
+         std::cout << "     --rotangle       Angle to rotate view through (fraction of a circle)\n";
+         std::cout << "                      (default = select appropriate for geometry)\n";
+         std::cout << "     --geometry       Geometry of scanning\n";
+         std::cout << "        parallel      Parallel scan beams (default)\n";
+         std::cout << "        equilinear    Equilinear divergent scan beams\n";
+         std::cout << "        equiangular   Equiangular divergent scan beams\n";
+         std::cout << "     --focal-length   Focal length ratio (ratio to radius of phantom)\n";
+         std::cout << "                      (default = 1)\n";
+         std::cout << "     --view-ratio     Length to view (view diameter to phantom diameter)\n";
+         std::cout << "                      (default = 1)\n";
+         std::cout << "     --scan-ratio     Length to scan (scan diameter to view diameter)\n";
+         std::cout << "                      (default = 1)\n";
+         std::cout << "     --offsetview      Intial gantry offset in  'views' (default = 0)\n"; 
+         std::cout << "     --trace          Trace level to use\n";
+         std::cout << "        none          No tracing (default)\n";
+         std::cout << "        console       Trace text level\n";
+         std::cout << "        phantom       Trace phantom image\n";
+         std::cout << "        proj          Trace projections\n";
+         std::cout << "        plot          Trace plot\n";
+         std::cout << "        clipping      Trace clipping\n";
+         std::cout << "     --verbose        Verbose mode\n";
+         std::cout << "     --debug          Debug mode\n";
+         std::cout << "     --version        Print version\n";
+         std::cout << "     --help           Print this help message\n";
+}
+
+
+int 
+phm2helix_main (int argc, char* const argv[])
+{
+       Phantom phm;
+       std::string optGeometryName = Scanner::convertGeometryIDToName(Scanner::GEOMETRY_PARALLEL);
+       char *opt_outfile = NULL;
+       std::string opt_desc;
+       std::string opt_PhmProg;
+       std::string opt_PhmFileName = "tmpphmfile";
+       int opt_ndet;
+       int opt_nview;
+       int opt_offsetview = 0;
+       int opt_nray = 1;
+       double dOptFocalLength = 2.;
+       double dOptCenterDetectorLength = 2;
+       double dOptViewRatio = 1.;
+       double dOptScanRatio = 1.;
+       int opt_trace = Trace::TRACE_NONE;
+       int opt_verbose = 0;
+       int opt_debug = 0;
+       double opt_rotangle = -1;
+       char* endptr = NULL;
+       char* endstr;
+
+       Timer timerProgram;
+
+       while (1) {
+               int c = getopt_long(argc, argv, "", phm2helix_options, NULL);
+
+               if (c == -1)
+               break;
+               
+       switch (c) {
+       case O_VERBOSE:
+               opt_verbose = 1;
+               break;
+       case O_DEBUG:
+               opt_debug = 1;
+               break;
+       case O_TRACE:
+               if ((opt_trace = Trace::convertTraceNameToID(optarg)) 
+                                                               == Trace::TRACE_INVALID) {
+                       phm2helix_usage(argv[0]);
+                       return (1);
+               }
+               break;
+                         case O_PHMFILE:
+                               opt_PhmFileName = optarg;
+                               break;
+                         case O_DESC:
+                               opt_desc = optarg;
+                               break;
+                         case O_ROTANGLE:
+                               opt_rotangle = strtod(optarg, &endptr);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                       std::cerr << "Error setting --rotangle to " << optarg << std::endl;
+                                       phm2helix_usage(argv[0]);
+                                       return (1);
+                               }
+                               break;
+                         case O_GEOMETRY:
+                               optGeometryName = optarg;
+                               break;
+                         case O_FOCAL_LENGTH:
+                               dOptFocalLength = strtod(optarg, &endptr);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                       std::cerr << "Error setting --focal-length to " << optarg << std::endl;
+                                       phm2helix_usage(argv[0]);
+                                       return (1);
+                               }
+                               break;
+                         case  O_CENTER_DETECTOR_LENGTH:
+                               dOptCenterDetectorLength = strtod(optarg, &endptr);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                       std::cerr << "Error setting --center-detector-length to " << optarg << std::endl; 
+                                       phm2helix_usage(argv[0]);
+                                       return (1);
+                               }
+                         break; 
+                         case O_VIEW_RATIO:
+                               dOptViewRatio = strtod(optarg, &endptr);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                       std::cerr << "Error setting --view-ratio to " << optarg << std::endl;
+                                       phm2helix_usage(argv[0]);
+                                       return (1);
+                               }
+                               break;
+                         case O_SCAN_RATIO:
+                               dOptScanRatio = strtod(optarg, &endptr);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                       std::cerr << "Error setting --scan-ratio to " << optarg << std::endl;
+                                       phm2helix_usage(argv[0]);
+                                       return (1);
+                               }
+                               break;
+                         case O_NRAY:
+                               opt_nray = strtol(optarg, &endptr, 10);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) {
+                                 std::cerr << "Error setting --nray to %s" << optarg << std::endl;
+                                 phm2helix_usage(argv[0]);
+                                 return (1);
+                               }
+                               break;
+                         case O_OFFSETVIEW:
+                               opt_offsetview = strtol(optarg, &endptr, 10);
+                               endstr = optarg + strlen(optarg);
+                               if (endptr != endstr) { 
+                                 std::cerr << "Error setting --offsetview to %s" << optarg << std::endl;
+                                 phm2helix_usage(argv[0]);
+                                 return (1);
+                               }
+                               break;
+                         case O_VERSION:
+#ifdef VERSION
+                                 std::cout << "Version: " << VERSION << std::endl << g_szIdStr << std::endl;
+#else
+                                         std::cout << "Unknown version number\n";
+#endif
+                                 return (0);
+                         case O_HELP:
+                         case '?':
+                               phm2helix_usage(argv[0]);
+                               return (0);
+                         default:
+                               phm2helix_usage(argv[0]);
+                               return (1);
+                       } // end of switch 
+         } // end of while loop
+         
+         if (optind + 4 != argc) {
+               phm2helix_usage(argv[0]);
+               return (1);
+         }
+
+         opt_outfile = argv[optind];
+         opt_ndet = strtol(argv[optind+1], &endptr, 10);
+         endstr = argv[optind+1] + strlen(argv[optind+1]);
+         if (endptr != endstr) {
+               std::cerr << "Error setting --ndet to " << argv[optind+1] << std::endl;
+               phm2helix_usage(argv[0]);
+               return (1);
+         }
+         opt_nview = strtol(argv[optind+2], &endptr, 10);
+         endstr = argv[optind+2] + strlen(argv[optind+2]);
+         if (endptr != endstr) {
+               std::cerr << "Error setting --nview to " << argv[optind+2] << std::endl;
+               phm2helix_usage(argv[0]);
+               return (1);
+         }
+         opt_PhmProg = argv[optind+3];
+
+         if (opt_rotangle < 0) {
+               if (optGeometryName.compare ("parallel") == 0)
+                 opt_rotangle = 0.5;
+               else
+                 opt_rotangle = 1.0;
+         }
+
+         std::ostringstream desc;
+         desc << "phm2helix: NDet=" << opt_ndet 
+                  << ", Nview=" << opt_nview 
+                  << ", NRay=" << opt_nray 
+                  << ", RotAngle=" << opt_rotangle  
+                  << ", OffsetView =" << opt_offsetview 
+                  << ", Geometry=" << optGeometryName 
+                  << ", PhantomProg=" << opt_PhmProg
+                  << ", PhmFileName=" << opt_PhmFileName;
+         if (opt_desc.length()) {
+               desc << ": " << opt_desc;
+         }
+         opt_desc = desc.str();
+
+         opt_rotangle *= TWOPI;
+
+         int stat; 
+         char extcommand[100];
+         if(opt_debug != 0)
+                       std::cout  <<  opt_PhmProg  <<  " " << 0 << " " <<  opt_nview << " " << opt_PhmFileName  << std::endl; 
+          //extcommand <<  opt_PhmProg  <<  " " << 0 << " " <<  opt_nview << " " << opt_PhmFileName ; 
+
+         sprintf(extcommand, "%s %d %d %s",    opt_PhmProg.c_str(), 0, opt_nview, opt_PhmFileName.c_str() );
+
+          stat = system( extcommand ); 
+         if (stat != 0 ) 
+                       std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl; 
+               
+         phm.createFromFile (opt_PhmFileName.c_str());
+         remove(opt_PhmFileName.c_str());
+
+         Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, 
+                               opt_offsetview, opt_nray, opt_rotangle, dOptFocalLength, 
+                               dOptCenterDetectorLength, dOptViewRatio, dOptScanRatio);
+         if (scanner.fail()) {
+                std::cout << "Scanner Creation Error: " << scanner.failMessage() 
+                                                       << std::endl;
+                return (1);
+         }
+
+         Projections pjGlobal(scanner);
+
+         
+         for( int iView = 0; iView < opt_nview; iView++ ){
+               if(opt_debug != 0)
+                       std::cout  <<  opt_PhmProg  <<  " " << iView << " " <<  opt_nview << " " << opt_PhmFileName  << std::endl; 
+          //extcommand <<  opt_PhmProg  <<  " " << iView << " " <<  opt_nview << " " << opt_PhmFileName ; 
+                       sprintf(extcommand, "%s %d %d %s",      opt_PhmProg.c_str(), iView, opt_nview, opt_PhmFileName.c_str() );
+               stat = system( extcommand ); 
+
+               if (stat != 0 ) 
+                       std::cerr << "Error executing external phantom program " << opt_PhmProg << " with command " << extcommand << std::endl; 
+               Phantom phmtmp; 
+               phmtmp.createFromFile (opt_PhmFileName.c_str());
+
+               scanner.collectProjections (pjGlobal, phmtmp, iView, 
+                         1, scanner.offsetView(), true, opt_trace);
+               remove(opt_PhmFileName.c_str());
+         }
+
+         
+         pjGlobal.setCalcTime (timerProgram.timerEnd());
+               pjGlobal.setRemark (opt_desc);
+               pjGlobal.write (opt_outfile);
+               if (opt_verbose) {
+                 phm.print (std::cout);
+                 std::cout << std::endl;
+               std::ostringstream os;
+               pjGlobal.printScanInfo (os);
+               std::cout << os.str() << std::endl;
+               std::cout << "  Remark: " << pjGlobal.remark() << std::endl;
+               std::cout << "Run time: " << pjGlobal.calcTime() << " seconds\n";
+         }
+               
+         return (0);
+}
+
+
+#ifndef NO_MAIN
+int 
+main (int argc, char* argv[])
+{
+  int retval = 1;
+
+  try {
+    retval = phm2helix_main(argc, argv);
+#if HAVE_DMALLOC
+    //    dmalloc_shutdown();
+#endif
+  } catch (exception e) {
+    std::cerr << "Exception: " << e.what() << std::endl;
+  } catch (...) {
+    std::cerr << "Unknown exception\n";
+  }
+
+  return (retval);
+}
+#endif
+
index cfae91f..ffa7ec2 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: phm2pj.cpp,v 1.30 2001/04/02 03:49:52 kevin Exp $
+**  $Id: phm2pj.cpp,v 1.31 2001/09/24 09:40:42 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
@@ -30,7 +30,7 @@
 
 
 enum { O_PHANTOM, O_DESC, O_NRAY, O_ROTANGLE, O_PHMFILE, O_GEOMETRY, O_FOCAL_LENGTH, O_CENTER_DETECTOR_LENGTH,
-O_VIEW_RATIO, O_SCAN_RATIO, O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
+O_VIEW_RATIO, O_SCAN_RATIO, O_OFFSETVIEW, O_TRACE, O_VERBOSE, O_HELP, O_DEBUG, O_VERSION };
 
 static struct option phm2pj_options[] = 
 {
@@ -42,6 +42,7 @@ static struct option phm2pj_options[] =
   {"geometry", 1, 0, O_GEOMETRY},
   {"focal-length", 1, 0, O_FOCAL_LENGTH},
   {"center-detector-length", 1, 0, O_CENTER_DETECTOR_LENGTH},
+  {"offsetview", 1, 0, O_OFFSETVIEW},
   {"view-ratio", 1, 0, O_VIEW_RATIO},
   {"scan-ratio", 1, 0, O_SCAN_RATIO},
   {"trace", 1, 0, O_TRACE},
@@ -52,7 +53,7 @@ static struct option phm2pj_options[] =
   {0, 0, 0, 0}
 };
 
-static const char* g_szIdStr = "$Id: phm2pj.cpp,v 1.30 2001/04/02 03:49:52 kevin Exp $";
+static const char* g_szIdStr = "$Id: phm2pj.cpp,v 1.31 2001/09/24 09:40:42 kevin Exp $";
 
 
 void 
@@ -84,6 +85,7 @@ phm2pj_usage (const char *program)
   std::cout << "                      (default = 1)\n";
   std::cout << "     --scan-ratio     Length to scan (scan diameter to view diameter)\n";
   std::cout << "                      (default = 1)\n";
+  std::cout << "     --offsetview     Initial gantry offset in 'views' (default = 0)\n";
   std::cout << "     --trace          Trace level to use\n";
   std::cout << "        none          No tracing (default)\n";
   std::cout << "        console       Trace text level\n";
@@ -112,6 +114,7 @@ phm2pj_main (int argc, char* const argv[])
   std::string optPhmFileName;
   int opt_ndet;
   int opt_nview;
+  int opt_offsetview = 0;
   int opt_nray = 1;
   double dOptFocalLength = 2.;
   double dOptCenterDetectorLength = 2.;
@@ -219,6 +222,16 @@ phm2pj_main (int argc, char* const argv[])
           return (1);
         }
         break;
+         case O_OFFSETVIEW:
+               opt_offsetview = strtol(optarg, &endptr, 10);
+               endstr = optarg + strlen(optarg);
+               if (endptr != endstr) { 
+                 std::cerr << "Error setting --offsetview to %s" << optarg << std::endl;
+                 phm2pj_usage(argv[0]);
+                 return (1);
+               }
+               break;
+
       case O_VERSION:
 #ifdef VERSION
         std::cout << "Version: " << VERSION << std::endl << g_szIdStr << std::endl;
@@ -270,7 +283,7 @@ phm2pj_main (int argc, char* const argv[])
     }
     
     std::ostringstream desc;
-    desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << ", Geometry=" << optGeometryName << ", ";
+    desc << "phm2pj: NDet=" << opt_ndet << ", Nview=" << opt_nview << ", NRay=" << opt_nray << ", RotAngle=" << opt_rotangle << "OffsetView =" << opt_offsetview << ", Geometry=" << optGeometryName << ", ";
     if (optPhmFileName.length()) {
       desc << "PhantomFile=" << optPhmFileName;
     } else if (optPhmName != "") {
@@ -324,8 +337,8 @@ phm2pj_main (int argc, char* const argv[])
 #endif
   
   opt_rotangle *= TWOPI;
-  Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_nray, opt_rotangle, dOptFocalLength
-    dOptCenterDetectorLength, dOptViewRatio, dOptScanRatio);
+  Scanner scanner (phm, optGeometryName.c_str(), opt_ndet, opt_nview, opt_offsetview, opt_nray
+               opt_rotangle, dOptFocalLength, dOptCenterDetectorLength, dOptViewRatio, dOptScanRatio);
   if (scanner.fail()) {
     std::cout << "Scanner Creation Error: " << scanner.failMessage() << std::endl;
     return (1);
@@ -361,7 +374,7 @@ phm2pj_main (int argc, char* const argv[])
   
 #else
   Projections pjGlobal (scanner);
-  scanner.collectProjections (pjGlobal, phm, opt_trace);
+  scanner.collectProjections (pjGlobal, phm, 0, opt_nview, opt_offsetview, true, opt_trace);
 #endif
   
 #ifdef HAVE_MPI
diff --git a/tools/pjHinterp.cpp b/tools/pjHinterp.cpp
new file mode 100644 (file)
index 0000000..3d1403f
--- /dev/null
@@ -0,0 +1,175 @@
+/*****************************************************************************
+* ** FILE IDENTIFICATION
+* **
+* **   Name:          phm2helix.cpp
+* **   Purpose:       Take projections of a phantom object
+* **   Programmer:    Ian Kay
+* **   Date Started:  Aug 2001
+* **
+* **  This is part of the CTSim program
+* **  Copyright (C) 1983-2000 Kevin Rosenberg
+* **
+* **  $Id: pjHinterp.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $
+* **
+* **  This program is free software; you can redistribute it and/or modify
+* **  it under the terms of the GNU General Public License (version 2) as
+* **  published by the Free Software Foundation.
+* **
+* **  This program is distributed in the hope that it will be useful,
+* **  but WITHOUT ANY WARRANTY; without even the implied warranty of
+* **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+* **  GNU General Public License for more details.
+* **
+* **  You should have received a copy of the GNU General Public License
+* **  along with this program; if not, write to the Free Software
+* **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+* ******************************************************************************/
+
+#include "ct.h"
+#include "timer.h"
+
+
+enum { O_INTERPVIEW, O_VERBOSE, O_TRACE, O_HELP, O_DEBUG, O_VERSION};
+
+
+static struct option my_options[] =
+{
+       {"interpview", 1, 0, O_INTERPVIEW},
+       {"trace", 1, 0, O_TRACE},
+       {"debug", 0, 0, O_DEBUG},
+       {"verbose", 0, 0, O_VERBOSE},
+       {"help", 0, 0, O_HELP},
+       {"version", 0, 0, O_VERSION},
+       {0, 0, 0, 0}
+};
+
+static const char* g_szIdStr = "$Id: pjHinterp.cpp,v 1.1 2001/09/24 09:40:42 kevin Exp $";
+
+void pjHinterp_usage ( const char *program )
+{
+       std::cout << "usage: " << fileBasename(program) << " raysum-file interp-raysum-file [OPTIONS]" << std::endl;
+       std::cout << "Interpolation of helical data in raw data space" << std::endl;
+       std::cout << "  raysum-file    Input raysum file" << std::endl;
+    std::cout << "  interp-file    Output interpolated raysum file " << std::endl;
+       std::cout << "  --trace        Set tracing to level" << std::endl;
+       std::cout << "     none        No tracing (default)" << std::endl;
+       std::cout << "     console     Text level tracing" << std::endl;
+       std::cout << "     phantom     Trace phantom" << std::endl;
+       std::cout << "     proj        Trace allrays" << std::endl;
+       std::cout << "     plot        Trace plotting" << std::endl;
+       std::cout << "     clipping    Trace clipping" << std::endl;
+       std::cout << "  --verbose      Turn on verbose mode" << std::endl;
+       std::cout << "  --debug        Turn on debug mode" << std::endl;
+       std::cout << "  --version      Print version" << std::endl;
+       std::cout << "  --help         Print this help message" << std::endl;
+}
+
+
+int
+pjHinterp_main(int argc, char * const argv[])
+{
+       Projections projGlobal;
+       char* pszProjFilename = NULL;
+       char* pszInterpFilename = NULL;
+       bool bOptVerbose = false;
+       bool bOptDebug = 1;
+       int optTrace = Trace::TRACE_NONE;
+       char *endptr = NULL;
+       char *endstr; 
+       int opt_interpview=-1;
+
+       while (1) {
+               int c = getopt_long(argc, argv, "", my_options, NULL);
+
+               if (c == -1)
+                break;
+
+               switch (c) {
+               case O_INTERPVIEW:
+                       opt_interpview = strtol(optarg, &endptr, 10);
+                       endstr = optarg + strlen(optarg);
+                       if (endptr != endstr) { 
+                               std::cerr << "Error setting --interpview to %s" << optarg << std::endl;
+                               pjHinterp_usage(argv[0]);
+                               return(1);
+                       }
+                       break;
+               case O_VERBOSE:
+                       bOptVerbose = true;
+                       break;
+               case O_DEBUG:
+                       bOptDebug = true;
+                       break;
+               case O_TRACE:
+                       if ((optTrace = Trace::convertTraceNameToID(optarg)) 
+                                                       == Trace::TRACE_INVALID) {
+                               pjHinterp_usage(argv[0]);
+                               return (1);
+                       }
+                       break;
+               case O_VERSION:
+#ifdef VERSION
+                       std::cout <<  "Version " <<  VERSION << std::endl << 
+                                                                                               g_szIdStr << std::endl;
+#else
+                       std::cout << "Unknown version number" << std::endl;
+#endif
+                       return (0);
+
+               case O_HELP:
+               case '?':      
+                       pjHinterp_usage(argv[0]);
+                       return (0);
+               default:
+                       pjHinterp_usage(argv[0]);
+                       return (1);
+               } // end switch
+       } // end while
+
+    if (optind + 2 != argc) {
+      pjHinterp_usage(argv[0]);
+      return (1);
+    }
+
+    pszProjFilename = argv[optind];
+
+    pszInterpFilename = argv[optind + 1];
+
+       Projections projections;
+       if ( projections.read(pszProjFilename) != true ){ 
+                       std::cerr << "Error reading input file " << pszProjFilename << std::endl;       
+                       return (1);
+       }
+       if (bOptVerbose) { 
+                       ostringstream os;
+                       projections.printScanInfo(os);
+                       std::cout << os.str();
+       }
+
+       int status = projections.Helical180LI(opt_interpview);
+       if ( status != 0 )  return (1);
+       status = projections.HalfScanFeather();
+       if ( status != 0 )  return (1);
+       projections.write( pszInterpFilename  );
+       
+       return (0);
+}
+
+
+#ifndef NO_MAIN
+int
+main (int argc, char* argv[])
+{
+  int retval = 1;
+
+  try {
+    retval = pjHinterp_main(argc, argv);
+  } catch (exception e) {
+      std::cerr << "Exception: " << e.what() << std::endl;
+  } catch (...) {
+      std::cerr << "Unknown exception" << std::endl;
+  }
+
+  return (retval);
+}
+#endif