** Date Started: Aug 84
**
** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
+** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.38 2001/01/06 15:33:15 kevin Exp $
+** $Id: projections.cpp,v 1.48 2001/02/20 17:44:14 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
deleteProjData();
init (scanner.nView(), scanner.nDet());
- 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_dFocalLength = scanner.focalLength();
+ m_dViewDiameter = scanner.viewDiameter();
m_rotStart = 0;
m_detStart = -(scanner.detLen() / 2);
+ m_dFanBeamAngle = scanner.fanBeamAngle();
}
void
kfloat64 _rotInc = m_rotInc;
kfloat64 _detStart = m_detStart;
kfloat64 _detInc = m_detInc;
- kfloat64 _phmLen = m_phmLen;
- kfloat64 _fieldOfView = m_fieldOfView;
- kfloat64 _focalLength = m_focalLength;
-
+ kfloat64 _viewDiameter = m_dViewDiameter;
+ kfloat64 _focalLength = m_dFocalLength;
+ kfloat64 _fanBeamAngle = m_dFanBeamAngle;
+
fs.seekp(0);
if (! fs)
return false;
fs.writeFloat64 (_rotInc);
fs.writeFloat64 (_detStart);
fs.writeFloat64 (_detInc);
- fs.writeFloat64 (_phmLen);
+ fs.writeFloat64 (_viewDiameter);
fs.writeFloat64 (_focalLength);
- fs.writeFloat64 (_fieldOfView);
+ fs.writeFloat64 (_fanBeamAngle);
fs.writeInt16 (_year);
fs.writeInt16 (_month);
fs.writeInt16 (_day);
{
kuint16 _hsize, _signature, _year, _month, _day, _hour, _minute, _second, _remarksize = 0;
kuint32 _nView, _nDet, _geom;
- kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _phmLen, _focalLength, _fieldOfView;
+ kfloat64 _calcTime, _rotStart, _rotInc, _detStart, _detInc, _focalLength, _viewDiameter, _fanBeamAngle;
fs.seekg(0);
if (! fs)
fs.readFloat64 (_rotInc);
fs.readFloat64 (_detStart);
fs.readFloat64 (_detInc);
- fs.readFloat64 (_phmLen);
+ fs.readFloat64 (_viewDiameter);
fs.readFloat64 (_focalLength);
- fs.readFloat64 (_fieldOfView);
+ fs.readFloat64 (_fanBeamAngle);
fs.readInt16 (_year);
fs.readInt16 (_month);
fs.readInt16 (_day);
m_rotInc = _rotInc;
m_detStart = _detStart;
m_detInc = _detInc;
- m_phmLen = _phmLen;
- m_focalLength = _focalLength;
- m_fieldOfView = _fieldOfView;
+ m_dFocalLength = _focalLength;
+ m_dViewDiameter = _viewDiameter;
+ m_dFanBeamAngle = _fanBeamAngle;
m_year = _year;
m_month = _month;
m_day = _day;
bool
Projections::copyViewData (const char* const filename, std::ostream& os, int startView, int endView)
{
- frnetorderstream is (filename, ios::in | ios::binary);
+ frnetorderstream is (filename, std::ios::in | std::ios::binary);
kuint16 sizeHeader, signature;
kuint32 _nView, _nDet;
+ is.seekg (0);
+ if (is.fail()) {
+ sys_error (ERR_SEVERE, "Unable to read projection file %s", filename);
+ return false;
+ }
+
is.readInt16 (sizeHeader);
is.readInt16 (signature);
is.readInt32 (_nView);
int nDet = _nDet;
if (signature != m_signature) {
- sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
+ sys_error (ERR_SEVERE, "Illegal signature in projection file %s", filename);
return false;
}
delete pViewData;
if (is.fail())
- sys_error (ERR_FATAL, "Error reading projection file");
+ sys_error (ERR_SEVERE, "Error reading projection file");
if (os.fail())
- sys_error (ERR_FATAL, "Error writing projection file");
+ sys_error (ERR_SEVERE, "Error writing projection file");
return (! (is.fail() | os.fail()));
}
is.readInt16 (signature);
is.seekg (0);
if (signature != m_signature) {
- sys_error (ERR_FATAL, "Illegal signature in projection file %s", filename);
+ sys_error (ERR_SEVERE, "Illegal signature in projection file %s", filename);
return false;
}
unsigned char* pHdrData = new unsigned char [sizeHeader];
is.read (reinterpret_cast<char*>(pHdrData), sizeHeader);
if (is.fail()) {
- sys_error (ERR_FATAL, "Error reading header");
+ sys_error (ERR_SEVERE, "Error reading header");
return false;
}
os.write (reinterpret_cast<char*>(pHdrData), sizeHeader);
if (os.fail()) {
- sys_error (ERR_FATAL, "Error writing header");
+ sys_error (ERR_SEVERE, "Error writing header");
return false;
}
printf("Description: %s\n", m_remark.c_str());
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("focalLength = %8.4f ViewDiameter = %8.4f\n", m_dFocalLength, m_dViewDiameter);
+ printf("fanBeamAngle= %8.4f\n", convertRadiansToDegrees(m_dFanBeamAngle));
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) {
if (startView < 0)
startView = 0;
+ if (endView < 0)
+ endView = m_nView - 1;
if (startView > m_nView - 1)
startView = m_nView - 1;
if (endView > m_nView - 1)
Projections::printScanInfo (std::ostringstream& os) const
{
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";
+ os << "Number of views: " << m_nView<< "\n";
+ os << "Description: " << m_remark.c_str()<< "\n";
+ os << "Geometry: " << Scanner::convertGeometryIDToName (m_geometry)<< "\n";
+ os << "Focal Length: " << m_dFocalLength<< "\n";
+ os << "View Diameter: " << m_dViewDiameter<< "\n";
+ os << "Fan Beam Angle: " << convertRadiansToDegrees(m_dFanBeamAngle) << "\n";
+ os << "detStart: " << m_detStart<< "\n";
+ os << "detInc: " << m_detInc<< "\n";
+ os << "rotStart: " << m_rotStart<< "\n";
+ os << "rotInc: " << m_rotInc<< "\n";
}
double** ppdView = adView.getArray();
double** ppdDet = adDet.getArray();
- calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
+ if (! calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet))
+ return false;
std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
- for (unsigned int iView = 0; iView < m_nView; iView++) {
+ unsigned int iView;
+ for (iView = 0; iView < m_nView; iView++) {
ppcDetValue[iView] = new std::complex<double> [m_nDet];
for (unsigned int iDet = 0; iDet < m_nDet; iDet++)
ppcDetValue[iView][iDet] = std::complex<double>(getDetectorArray (iView).detValues()[iDet], 0);
if (! v || nx == 0 || ny == 0)
return false;
+#ifndef HAVE_FFT
+ return false;
+#else
Array2d<double> adView (nx, ny);
Array2d<double> adDet (nx, ny);
double** ppdView = adView.getArray();
fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
for (iView = 0; iView < m_nView; iView++) {
- ppcDetValue[iView] = new std::complex<double> [m_nDet];
unsigned int iDet;
for (iDet = 0; iDet < m_nDet; iDet++) {
pcIn[iDet].re = getDetectorArray(iView).detValues()[iDet];
pcIn[iDet].im = 0;
}
fftw_one (plan, pcIn, NULL);
+ ppcDetValue[iView] = new std::complex<double> [m_nDet];
for (iDet = 0; iDet < m_nDet; iDet++)
- ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im);
- Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
+ ppcDetValue[iView][iDet] = std::complex<double> (pcIn[iDet].re, pcIn[iDet].im);
+ Fourier::shuffleFourierToNaturalOrder (ppcDetValue[iView], m_nDet);
}
fftw_destroy_plan (plan);
delete [] pcIn;
- calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
+ bool bError = calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
- interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
+ if (! bError)
+ interpolatePolar (v, vImag, nx, ny, ppcDetValue, ppdView, ppdDet, m_nView, m_nDet, iInterpolationID);
for (iView = 0; iView < m_nView; iView++)
delete [] ppcDetValue[iView];
delete [] ppcDetValue;
- return true;
+ return bError;
+#endif
}
-void
+
+bool
Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet)
{
- double xMin = -m_phmLen / 2;
- double xMax = xMin + m_phmLen;
- double yMin = -m_phmLen / 2;
- double yMax = yMin + m_phmLen;
+ double xMin = -phmLen() / 2;
+ double xMax = xMin + phmLen();
+ double yMin = -phmLen() / 2;
+ double yMax = yMin + phmLen();
double xInc = (xMax - xMin) / nx; // size of cells
double yInc = (yMax - yMin) / ny;
if (m_geometry != Scanner::GEOMETRY_PARALLEL) {
sys_error (ERR_WARNING, "convertPolar supports Parallel only");
+ return false;
}
+ // Calculates polar coordinates (view#, det#) for each point on phantom grid
double x = xMin + xInc / 2; // Rectang coords of center of pixel
for (unsigned int ix = 0; ix < nx; x += xInc, ix++) {
double y = yMin + yInc / 2;
ppdDet[ix][iy] = (r / m_detInc) + iDetCenter;
}
}
+
+ return true;
}
void
v[ix][iy] = 0;
}
} else if (iInterpolationID == POLAR_INTERP_BILINEAR) {
- int iFloorView = ppdView[ix][iy];
+ int iFloorView = static_cast<int>(ppdView[ix][iy]);
double dFracView = ppdView[ix][iy] - iFloorView;
- int iFloorDet = ppdDet[ix][iy];
+ int iFloorDet = static_cast<int>(ppdDet[ix][iy]);
double dFracDet = ppdDet[ix][iy] - iFloorDet;
if (iFloorDet >= 0 && iFloorView >= 0) {