+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
/* Define for HDF5 library */
#undef HAVE_HDF5
+
+#undef HAVE_USLEEP
/* Define if you have the strtol function. */
#undef HAVE_STRTOL
+/* Define if you have the usleep function. */
+#undef HAVE_USLEEP
+
/* Define if you have the <X11/extensions/Xdbe.h> header file. */
#undef HAVE_X11_EXTENSIONS_XDBE_H
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; }
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
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.
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
<INPUT type="radio" name="MPI" value="no" checked>No (Single CPU)<br>
</td><td valign="top">
<h3>Simulate X-Ray acquistion</h3>
-Number of detectors: <input type="text" name="PJ_NDet" size="4" value="367"><p>
+Geometry:<br>
+<INPUT type="radio" name="geometry" value="parallel" checked>Parallel<br>
+<INPUT type="radio" name="geometry" value="equiangular">Equiangular<br>
+<INPUT type="radio" name="geometry" value="collinear">Collinear<br>
+<p>
+Number of Detectors: <input type="text" name="PJ_NDet" size="4" value="367"><p>
Number of Rotations: <input type="text" name="PJ_NRot" size="4" value="320"><p>
Number of Rays<br>(samples) per detector: <input type="text" name="PJ_NRay" size="2" value="3"><p>
Rotation Angle<br>as a multiple of PI: <input type="text" name="PJ_RotAngle" size="3" value="1.0"><br>
** 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
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);
void BackprojectView (const double* const t, double view_angle);
- private:
+ protected:
Array2d<kfloat64> arrayR;
Array2d<kfloat64> arrayPhi;
kfloat64** r;
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
** 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
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 ();
};
** 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
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;}
int m_nOutputPoints;
int m_iPreinterpolationFactor;
int m_idGeometry;
+ double m_dFocalLength;
bool m_fail;
string m_failMessage;
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
void finiteFourierTransform (const complex<double> input[], complex<double> output[], const int direction) const;
void finiteFourierTransform (const complex<double> input[], double output[], const int direction) const;
+ double convolve (const double func[], const double filter[], const double dx, const int n, const int np) const;
double convolve (const double f[], const double dx, const int n, const int np) const;
double convolve (const float f[], const double dx, const int n, const int np) const;
** 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
class Scanner;
class DetectorArray;
class Array2dFileLabel;
-
+class ostringstream;
// Projections
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);
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]); }
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);
** 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
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;}
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;
** 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
#endif
SGPDriver (const char* szWinTitle = "", int xsize = 640, int ysize = 480);
-
+
~SGPDriver ();
int getPhysicalXSize () const
#ifdef HAVE_WXWINDOWS
wxDC* idWX () const
{ return m_pDC; }
+
+ void setDC (wxDC* dc)
+ { m_pDC = dc; }
#endif
};
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;
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 ();
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
};
** 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
}
+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)
{
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;
** 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
#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
setTextAngle (0.);
setTextSize (1. / 25.);
setColor (C_BLACK);
-
}
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
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
}
{'\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
** 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
return false;
}
- if (m_idBackproject == BPROJ_TRIG)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_TABLE)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_DIFF)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_DIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_IDIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
- else if (m_idBackproject == BPROJ_IDIFF3)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
- else {
- m_fail = true;
- m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
- return false;
+ if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
+ else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
+ else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
+ if (m_idBackproject == BPROJ_TRIG)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_TABLE)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_DIFF)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_DIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_IDIFF2)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_IDIFF3)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
+ } else {
+ m_fail = true;
+ m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
+ return false;
}
return true;
{
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();
xInc = (xMax - xMin) / nx; // size of cells
yInc = (yMax - yMin) / ny;
+
+ m_dFocalLength = proj.focalLength();
}
Backproject::~Backproject ()
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
v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
}
}
+ }
}
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];
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) {
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;
void
BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle;
static const kint32 scale = 1 << 16;
static const double dScale = scale;
static const kint32 halfScale = scale / 2;
- const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
// calculate L for first point in image (0, 0)
kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
void
BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
+ double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
static const int scaleShift = 16;
static const kint32 scale = (1 << scaleShift);
static const kint32 scaleBitmask = scale - 1;
static const kint32 halfScale = scale / 2;
static const double dInvScale = 1. / scale;
- const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
// calculate L for first point in image (0, 0)
kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
} //end linear
} // end for ix
}
+
+
+void
+BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+ double dFLPlusSin = m_dFocalLength + rsin_t;
+ double gamma = atan (rcos_t / dFLPlusSin);
+ double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ int iDetPos =iDetCenter + nearest<int>(gamma / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+ } else
+ pImCol[iy] += filteredProj[iDetPos] / dL2;
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ double dPos = gamma / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1) {
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
+ } else
+ pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / dL2;
+ }
+ } // end for y
+ } // end for x
+}
+
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+
+ double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
+ double dDetPos = rcos_t / dU;
+ // double to scale for imaginary detector that passes through origin
+ // of phantom, see Kak-Slaney Figure 3.22
+ dDetPos *= 2;
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
+ ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+ else
+ pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ double dPos = dDetPos / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos < 0 || iDetPos >= nDet - 1)
+ ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
+ else
+ pImCol[iy] += (((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1])) / (dU * dU);
+ }
+ } // end for y
+ } // end for x
+}
+
** 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
// 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);
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;
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;
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
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);
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;
}
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++)
** 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
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
return true;
}
+bool
+Projections::read (const string& filename)
+{
+ return read (filename.c_str());
+}
+
+
bool
Projections::read (const char* filename)
{
}
+bool
+Projections::write (const string& filename)
+{
+ return write (filename.c_str());
+}
+
bool
Projections::write (const char* filename)
{
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]);
}
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";
}
#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());
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)
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: scanner.cpp,v 1.12 2000/08/27 20:32:55 kevin Exp $
+** $Id: scanner.cpp,v 1.13 2000/08/31 08:38:58 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
const int Scanner::GEOMETRY_INVALID = -1;
const int Scanner::GEOMETRY_PARALLEL = 0;
-const int Scanner::GEOMETRY_EQUILINEAR = 1;
-const int Scanner::GEOMETRY_EQUIANGULAR = 2;
+const int Scanner::GEOMETRY_EQUIANGULAR = 1;
+const int Scanner::GEOMETRY_EQUILINEAR = 2;
const char* Scanner::s_aszGeometryName[] =
{
{"parallel"},
- {"equilinear"},
{"equiangular"},
+ {"equilinear"},
};
const char* Scanner::s_aszGeometryTitle[] =
{
{"Parallel"},
- {"Equilinear"},
{"Equiangular"},
+ {"Equilinear"},
};
const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
m_nSample = 1;
if (nDet < 1)
nDet = 1;
- // if ((nDet % 2) == 0)
- // ++nDet; // ensure odd number of detectors
m_nDet = nDet;
m_nView = nView;
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.) {
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.) {
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
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;
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 ();
m_pSGP->lineAbs (xs2, ys2);
m_pSGP->setPenWidth (2);
m_pSGP->moveAbs (0., 0.);
- m_pSGP->drawArc (m_dFocalLength, start_angle + m_initPos.angle, start_angle - m_initPos.angle);
+ m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
}
m_pSGP->setPenWidth (1);
}
if (m_trace >= Trace::TRACE_CONSOLE)
- traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / m_nView * 100.);
+ traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast<double>(m_nView) * 100.);
#endif
- projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, dDetAngle);
+ projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
detArray.setViewAngle (viewAngle);
#ifdef HAVE_SGP
#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);
}
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
#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
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;
** 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
/* 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
*/
** 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
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;
** 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
#include <getopt.h>
#endif
-static const char* rcsindent = "$Id";
+static const char* rcsindent = "$Id: ctsim.cpp,v 1.10 2000/08/31 08:38:58 kevin Exp $";
class CTSimApp* theApp = NULL;
** 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
}
-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);
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);
*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;
+}
** 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
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;
** 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
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();
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();
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);
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;
}
} 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
}
}
+
+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);
+ }
}
** 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
#include "wx/setup.h"
#include "wx/dialog.h"
+#include "wx/dcmemory.h"
class wxButton;
class wxStaticText;
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()
};
** 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
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;
dialogAutoScale.getMinMax (&dMin, &dMax);
m_dMinPixel = dMin;
m_dMaxPixel = dMax;
+ m_dAutoScaleFactor = dialogAutoScale.getAutoScaleFactor();
OnUpdate (this, NULL);
}
}
m_bMinSpecified = false;
m_bMaxSpecified = false;
+ m_dAutoScaleFactor = 1.;
int width, height;
m_frame->GetClientSize(&width, &height);
if (iTrace > Trace::TRACE_CONSOLE) {
ProjectionsDialog dialogProjections (theScanner, rProj, rPhantom, iTrace, dynamic_cast<wxWindow*>(m_frame));
for (int iView = 0; iView < rProj.nView(); iView++) {
- if (! dialogProjections.projectView (iView)) {
+ ::wxYield();
+ ::wxYield();
+ if (dialogProjections.isCancelled() || ! dialogProjections.projectView (iView)) {
pProjectionDoc->DeleteAllViews();
return;
}
::wxYield();
+ ::wxYield();
+ while (dialogProjections.isPaused()) {
+ ::wxYield();
+ ::wxUsleep(50);
+ }
}
} else
- theScanner.collectProjections (rProj, rPhantom);
+ theScanner.collectProjections (rProj, rPhantom, iTrace);
pProjectionDoc->Modify(true);
pProjectionDoc->UpdateAllViews(this);
}
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();
}
}
{
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();
}
** 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
bool m_bMaxSpecified;
double m_dMinPixel;
double m_dMaxPixel;
+ double m_dAutoScaleFactor;
public:
ImageFileView(void);
-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@
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)
** 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
{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
}
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 != "") {
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());
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";
}
** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: pj2if.cpp,v 1.3 2000/08/03 09:57:29 kevin Exp $
+** $Id: pj2if.cpp,v 1.4 2000/08/31 08:38:58 kevin Exp $
**
** This program is free software; you can redistribute it and/or modify
** it under the terms of the GNU General Public License (version 2) as
#include "timer.h"
-enum { O_VERBOSE, O_HELP, O_VERSION };
+enum { O_VERBOSE, O_DUMP, O_HELP, O_VERSION };
static struct option my_options[] =
{
{"verbose", 0, 0, O_VERBOSE},
+ {"dump", 0, 0, O_DUMP},
{"help", 0, 0, O_HELP},
{"version", 0, 0, O_VERSION},
{0, 0, 0, 0}
};
-static const char* g_szIdStr = "$Id: pj2if.cpp,v 1.3 2000/08/03 09:57:29 kevin Exp $";
+static const char* g_szIdStr = "$Id: pj2if.cpp,v 1.4 2000/08/31 08:38:58 kevin Exp $";
void
pj2if_usage (const char *program)
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;
}
{
char *pj_name, *im_name;
bool optVerbose = false;
+ bool optDump = false;
extern int optind;
Timer timerProgram;
case O_VERBOSE:
optVerbose = true;
break;
+ case O_DUMP:
+ optDump = true;
+ break;
case O_VERSION:
#ifdef VERSION
cout << "Version " << VERSION << endl << g_szIdStr << endl;
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());
** 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
{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)
#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();
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