From bd1d464294e037da19ccc80d8cc60475768eb2ca Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Mon, 24 Sep 2001 09:40:42 +0000 Subject: [PATCH] r1018: *** empty log message *** --- Makefile.am | 2 +- Makefile.in | 4 +- configure | 6 +- configure.in | 4 +- helical/Makefile.am | 5 + helical/dynphm.c | 49 +++++ helical/sample-helical.sh.in | 50 +++++ include/imagefile.h | 4 +- include/projections.h | 6 +- include/scanner.h | 16 +- libctsim/imagefile.cpp | 29 ++- libctsim/projections.cpp | 275 +++++++++++++++++++++++++- libctsim/scanner.cpp | 29 +-- src/ctsim.cpp | 6 +- src/dialogs.cpp | 21 +- src/dialogs.h | 7 +- src/dlgprojections.cpp | 4 +- src/threadproj.cpp | 18 +- src/threadproj.h | 11 +- src/views.cpp | 30 +-- src/views.h | 3 +- tools/Makefile.am | 6 +- tools/Makefile.in | 17 +- tools/Makefile.nt | 6 + tools/ctsimtext.cpp | 15 +- tools/ifexport.cpp | 11 +- tools/phm2helix.cpp | 361 +++++++++++++++++++++++++++++++++++ tools/phm2pj.cpp | 27 ++- tools/pjHinterp.cpp | 175 +++++++++++++++++ 29 files changed, 1109 insertions(+), 88 deletions(-) create mode 100644 helical/Makefile.am create mode 100644 helical/dynphm.c create mode 100644 helical/sample-helical.sh.in create mode 100644 tools/phm2helix.cpp create mode 100644 tools/pjHinterp.cpp diff --git a/Makefile.am b/Makefile.am index 15cdf07..09164f7 100644 --- a/Makefile.am +++ b/Makefile.am @@ -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 diff --git a/Makefile.in b/Makefile.in index 2d6f0ff..aba3a41 100644 --- a/Makefile.in +++ b/Makefile.in @@ -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) diff --git a/configure b/configure index c4d45b0..efc5a00 100755 --- 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 <> $CONFIG_STATUS <> $CONFIG_STATUS <> $CONFIG_STATUS <<\EOF for ac_file in .. $CONFIG_FILES; do if test "x$ac_file" != x..; then diff --git a/configure.in b/configure.in index 7463752..bb61c86 100644 --- a/configure.in +++ b/configure.in @@ -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 index 0000000..ce6143e --- /dev/null +++ b/helical/Makefile.am @@ -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 index 0000000..4395f02 --- /dev/null +++ b/helical/dynphm.c @@ -0,0 +1,49 @@ +#include +#include +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 index 0000000..9e12f20 --- /dev/null +++ b/helical/sample-helical.sh.in @@ -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 + diff --git a/include/imagefile.h b/include/imagefile.h index 56c255b..6dd8fe7 100644 --- a/include/imagefile.h +++ b/include/imagefile.h @@ -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;} diff --git a/include/projections.h b/include/projections.h index d2fe89c..b021174 100644 --- a/include/projections.h +++ b/include/projections.h @@ -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); diff --git a/include/scanner.h b/include/scanner.h index 663ebd7..09ab796 100644 --- a/include/scanner.h +++ b/include/scanner.h @@ -9,7 +9,7 @@ ** This is part of the CTSim program ** Copyright (C) 1983-2000 Kevin Rosenberg ** -** $Id: scanner.h,v 1.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 diff --git a/libctsim/imagefile.cpp b/libctsim/imagefile.cpp index 356e489..d10e1bd 100644 --- a/libctsim/imagefile.cpp +++ b/libctsim/imagefile.cpp @@ -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; +} + + diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index 8440bd4..86976c7 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -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((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 ((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((M_PI+fanAngle)/dbeta); + + } + int last_interp_view = static_cast ((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((beta + 2*gamma - M_PI)/dbeta); + //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); + newiView = nearest(( (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(( 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((beta2 + 2*gamma2 - M_PI)/dbeta); + iView1 = nearest(( (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 diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp index 0319912..a90ca44 100644 --- a/libctsim/scanner.cpp +++ b/libctsim/scanner.cpp @@ -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; diff --git a/src/ctsim.cpp b/src/ctsim.cpp index 9248247..2d94fa4 100644 --- a/src/ctsim.cpp +++ b/src/ctsim.cpp @@ -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(xDisplay * .75), nearest(yDisplay * .755)), wxDEFAULT_FRAME_STYLE); + wxSize(nearest(xDisplay * .25), nearest(yDisplay * .255)), wxDEFAULT_FRAME_STYLE); setIconForFrame (m_pFrame); m_pFrame->Centre(wxBOTH); diff --git a/src/dialogs.cpp b/src/dialogs.cpp index 8525c62..75596d4 100644 --- a/src/dialogs.cpp +++ b/src/dialogs.cpp @@ -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 () diff --git a/src/dialogs.h b/src/dialogs.h index 9faf10f..abb71e5 100644 --- a/src/dialogs.h +++ b/src/dialogs.h @@ -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; diff --git a/src/dlgprojections.cpp b/src/dlgprojections.cpp index 76b0ae4..7624beb 100644 --- a/src/dlgprojections.cpp +++ b/src/dlgprojections.cpp @@ -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); diff --git a/src/threadproj.cpp b/src/threadproj.cpp index f967269..d61c422 100644 --- a/src/threadproj.cpp +++ b/src/threadproj.cpp @@ -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 @@ -54,10 +54,10 @@ // ///////////////////////////////////////////////////////////////////// -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) { diff --git a/src/threadproj.h b/src/threadproj.h index 49212f4..a4f6993 100644 --- a/src/threadproj.h +++ b/src/threadproj.h @@ -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); diff --git a/src/views.cpp b/src/views.cpp index afceae4..f78f9e3 100644 --- a/src/views.cpp +++ b/src/views.cpp @@ -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; diff --git a/src/views.h b/src/views.h index 01f8945..40cc6fd 100644 --- a/src/views.h +++ b/src/views.h @@ -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; diff --git a/tools/Makefile.am b/tools/Makefile.am index e214688..951ec21 100644 --- a/tools/Makefile.am +++ b/tools/Makefile.am @@ -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 diff --git a/tools/Makefile.in b/tools/Makefile.in index 0743b0d..132410b 100644 --- a/tools/Makefile.in +++ b/tools/Makefile.in @@ -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 diff --git a/tools/Makefile.nt b/tools/Makefile.nt index 8237a41..c31ba0e 100644 --- a/tools/Makefile.nt +++ b/tools/Makefile.nt @@ -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 diff --git a/tools/ctsimtext.cpp b/tools/ctsimtext.cpp index 75f941f..54236af 100644 --- a/tools/ctsimtext.cpp +++ b/tools/ctsimtext.cpp @@ -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) diff --git a/tools/ifexport.cpp b/tools/ifexport.cpp index 8760332..d629422 100644 --- a/tools/ifexport.cpp +++ b/tools/ifexport.cpp @@ -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 index 0000000..8f7f185 --- /dev/null +++ b/tools/phm2helix.cpp @@ -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 + diff --git a/tools/phm2pj.cpp b/tools/phm2pj.cpp index cfae91f..ffa7ec2 100644 --- a/tools/phm2pj.cpp +++ b/tools/phm2pj.cpp @@ -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 index 0000000..3d1403f --- /dev/null +++ b/tools/pjHinterp.cpp @@ -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 -- 2.34.1