From a05f3cb550877e94aa118cc04b361c0c8fdb3dc3 Mon Sep 17 00:00:00 2001 From: "Kevin M. Rosenberg" Date: Thu, 31 Aug 2000 08:38:58 +0000 Subject: [PATCH] r186: *** empty log message *** --- ChangeLog | 7 ++ acconfig.h | 2 + config.h.in | 3 + configure | 4 +- configure.in | 4 +- html/simulate.html.in | 7 +- include/backprojectors.h | 25 +++++- include/ezplot.h | 3 +- include/procsignal.h | 13 ++- include/projections.h | 54 ++++++------ include/scanner.h | 6 +- include/sgp.h | 14 +++- libctgraphics/ezplot.cpp | 17 +++- libctgraphics/sgp.cpp | 66 +++++++++------ libctsim/backprojectors.cpp | 159 ++++++++++++++++++++++++++++-------- libctsim/procsignal.cpp | 88 +++++++++++++++++--- libctsim/projections.cpp | 61 ++++++++++---- libctsim/scanner.cpp | 104 +++++++++++------------ libctsupport/mathfuncs.cpp | 6 +- libctsupport/syserror.cpp | 4 +- src/ctsim.cpp | 4 +- src/dialogs.cpp | 20 ++++- src/dialogs.h | 5 +- src/dlgprojections.cpp | 150 +++++++++++++++++++++++++++++----- src/dlgprojections.h | 46 ++++++++--- src/views.cpp | 20 +++-- src/views.h | 3 +- tools/Makefile.am | 5 +- tools/phm2pj.cpp | 18 ++-- tools/pj2if.cpp | 21 +++-- tools/pjrec.cpp | 18 ++-- 31 files changed, 700 insertions(+), 257 deletions(-) diff --git a/ChangeLog b/ChangeLog index bd59494..e77b594 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,10 @@ +2.0.0-b11 - 8/31/00 + Added reconstruction for Equilinear and Equiangular geometries + Changed theta to be -PI/2 to make compliant with Kak-Slaney formulas + ctsim: Added projection graph to animation of projection collection + ctsim: Added single stepping to projection collection animation + ctsim: improved File/Properties display for projection files + 2.0.0-b10 - 8/25/00 ctsim: Added animation of projection collection processs ctsim: Added Auto Scaling for image windows diff --git a/acconfig.h b/acconfig.h index 9c06cd5..c7dffb5 100644 --- a/acconfig.h +++ b/acconfig.h @@ -74,3 +74,5 @@ /* Define for HDF5 library */ #undef HAVE_HDF5 + +#undef HAVE_USLEEP diff --git a/config.h.in b/config.h.in index 31cc9fd..97f445e 100644 --- a/config.h.in +++ b/config.h.in @@ -90,6 +90,9 @@ /* Define if you have the strtol function. */ #undef HAVE_STRTOL +/* Define if you have the usleep function. */ +#undef HAVE_USLEEP + /* Define if you have the header file. */ #undef HAVE_X11_EXTENSIONS_XDBE_H diff --git a/configure b/configure index b46197c..299bf43 100755 --- a/configure +++ b/configure @@ -712,7 +712,7 @@ fi PACKAGE=ctsim -VERSION=2.0.0b10 +VERSION=2.0.0b11 if test "`cd $srcdir && pwd`" != "`pwd`" && test -f $srcdir/config.status; then { echo "configure: error: source directory already configured; run "make distclean" there first" 1>&2; exit 1; } @@ -2614,7 +2614,7 @@ fi fi -for ac_func in strtod strtol snprintf htonl +for ac_func in strtod strtol snprintf htonl usleep do echo $ac_n "checking for $ac_func""... $ac_c" 1>&6 echo "configure:2621: checking for $ac_func" >&5 diff --git a/configure.in b/configure.in index 09d8ec9..0b5e270 100644 --- a/configure.in +++ b/configure.in @@ -4,7 +4,7 @@ dnl Must reset CDPATH so that bash's cd does not print to stdout dnl CDPATH= AC_INIT(src/ctsim.cpp) -AM_INIT_AUTOMAKE(ctsim,2.0.0b10) +AM_INIT_AUTOMAKE(ctsim,2.0.0b11) AM_CONFIG_HEADER(config.h) dnl Checks for programs. @@ -99,7 +99,7 @@ AC_STRUCT_TM dnl Checks for library functions. AC_FUNC_VPRINTF -AC_CHECK_FUNCS(strtod strtol snprintf htonl) +AC_CHECK_FUNCS(strtod strtol snprintf htonl usleep) AC_CHECK_FUNC(basename) AC_CHECK_FUNC(setjmp) if test "${OSTYPE}" = "cygwin" ; then diff --git a/html/simulate.html.in b/html/simulate.html.in index 5039901..bd0745a 100644 --- a/html/simulate.html.in +++ b/html/simulate.html.in @@ -28,7 +28,12 @@ MPI Supercomputing:
No (Single CPU)

Simulate X-Ray acquistion

-Number of detectors:

+Geometry:
+Parallel
+Equiangular
+Collinear
+

+Number of Detectors:

Number of Rotations:

Number of Rays
(samples) per detector:

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