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
@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
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)
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"
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
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
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
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)
--- /dev/null
+bin_SCRIPTS = sample-helical.sh
+
+realclean:
+ rm -f *.pgm *.if *~ *.pj
+
--- /dev/null
+#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);
+}
--- /dev/null
+#!/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
+
** 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
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;}
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;}
** 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
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);
** 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
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;}
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
** 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
#ifdef HAVE_CTN_DICOM
const int ImageFile::EXPORT_FORMAT_DICOM = 5;
#endif
+const int ImageFile::EXPORT_FORMAT_RAW = 6;
const char* ImageFile::s_aszExportFormatName[] =
{
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;
}
#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;
+}
+
+
** 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
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();
}
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
** 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
*/
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;
m_nDet = nDet;
m_nView = nView;
+ m_iOffsetView = offsetView;
m_nSample = nSample;
m_dFocalLengthRatio = dFocalLengthRatio;
m_dCenterDetectorRatio = dCenterDetectorRatio;
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) {
}
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;
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;
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;
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;
** 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
#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[] =
{
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);
** 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
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)
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);
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 ()
** 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
{
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 ();
private:
wxTextCtrl* m_pTextCtrlNDet;
wxTextCtrl* m_pTextCtrlNView;
+ wxTextCtrl* m_pTextCtrlOffsetView;
wxTextCtrl* m_pTextCtrlNSamples;
wxTextCtrl* m_pTextCtrlRotAngle;
wxTextCtrl* m_pTextCtrlFocalLength;
int m_iDefaultNDet;
int m_iDefaultNView;
+ int m_iDefaultOffsetView;
int m_iDefaultNSamples;
int m_iDefaultTrace;
int m_iDefaultGeometry;
** 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
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);
** 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)
{
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();
//
/////////////////////////////////////////////////////////////////////
-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++) {
{
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;
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)
{
** 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
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;
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);
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;
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);
Scanner* m_pScanner;
int m_iNDet;
int m_iNView;
+ int m_iOffsetView;
const char* m_pszGeometry;
int m_iNSample;
double m_dRotation;
{}
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);
** 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
m_iDefaultNView = 320;
m_iDefaultNSample = 2;
#endif
+ m_iDefaultOffsetView = 0;
m_dDefaultRotation = 1;
m_dDefaultFocalLength = 2;
m_dDefaultCenterDetectorLength = 2;
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();
m_iDefaultNDet = dialogProjection.getNDet();
m_iDefaultNView = dialogProjection.getNView();
+ m_iDefaultOffsetView = dialogProjection.getOffsetView();
m_iDefaultNSample = dialogProjection.getNSamples();
m_iDefaultTrace = dialogProjection.getTrace();
m_dDefaultRotation = dialogProjection.getRotAngle();
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";
}
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;
#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");
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;
** 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
int m_iDefaultNDet;
int m_iDefaultNView;
+ int m_iDefaultOffsetView;
int m_iDefaultNSample;
int m_iDefaultGeometry;
int m_iDefaultTrace;
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
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
@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
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)
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)
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
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
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
** 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
// 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";
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[]);
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
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";
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)
** 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
{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";
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)
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;
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);
--- /dev/null
+/*****************************************************************************
+** 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
+
** 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
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[] =
{
{"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},
{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
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";
std::string optPhmFileName;
int opt_ndet;
int opt_nview;
+ int opt_offsetview = 0;
int opt_nray = 1;
double dOptFocalLength = 2.;
double dOptCenterDetectorLength = 2.;
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;
}
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 != "") {
#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);
#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
--- /dev/null
+/*****************************************************************************
+* ** 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
projects / ctsim.git / commitdiff