** 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.37 2001/01/04 21:28:41 kevin Exp $
+** $Id: projections.cpp,v 1.46 2001/01/28 19:10:18 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
{
{"nearest"},
{"bilinear"},
- {"bicubic"},
+// {"bicubic"},
};
const char* Projections::s_aszInterpTitle[] =
{
{"Nearest"},
{"Bilinear"},
- {"Bicubic"},
+// {"Bicubic"},
};
const int Projections::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(char*);
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;
}
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)
calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
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);
return true;
}
+
+bool
+Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
+{
+ unsigned int nx = rIF.nx();
+ unsigned int ny = rIF.ny();
+ ImageFileArray v = rIF.getArray();
+ if (! rIF.isComplex())
+ rIF.convertRealToComplex();
+ ImageFileArray vImag = rIF.getImaginaryArray();
+
+ 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();
+ double** ppdDet = adDet.getArray();
+
+ std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
+ unsigned int iView;
+ double* pdDet = new double [m_nDet];
+ fftw_complex* pcIn = new fftw_complex [m_nDet];
+ fftw_plan plan = fftw_create_plan (m_nDet, FFTW_FORWARD, FFTW_IN_PLACE);
+
+ for (iView = 0; iView < m_nView; iView++) {
+ 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);
+ }
+
+ fftw_destroy_plan (plan);
+ delete [] pcIn;
+
+ calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
+
+ 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;
+#endif
+}
+
+
void
Projections::calcArrayPolarCoordinates (unsigned int nx, unsigned int ny, double** ppdView, double** ppdDet)
{
if (m_geometry != Scanner::GEOMETRY_PARALLEL) {
sys_error (ERR_WARNING, "convertPolar supports Parallel only");
+ return;
}
+ // 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;
if (iInterpolationID == POLAR_INTERP_NEAREST) {
int iView = nearest<int> (ppdView[ix][iy]);
int iDet = nearest<int> (ppdDet[ix][iy]);
- if (iView == nView)
- iView = nView - 1;
+ if (iView == nView) {
+ iView = 0;
+ // iDet = m_nDet - iDet;
+ }
if (iDet >= 0 && iDet < nDet && iView >= 0 && iView < nView) {
v[ix][iy] = ppcDetValue[iView][iDet].real();
if (vImag)
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) {
if (iFloorView < nView - 1)
v2 = ppcDetValue[iFloorView + 1][iFloorDet];
else
- v2 = v1;
+ v2 = ppcDetValue[0][iFloorDet];
if (iFloorDet < nDet - 1)
v4 = ppcDetValue[iFloorView][iFloorDet+1];
else
if (iFloorView < nView - 1 && iFloorDet < nDet - 1)
v3 = ppcDetValue [iFloorView+1][iFloorDet+1];
else if (iFloorView < nView - 1)
- v4 = v2;
+ v3 = v2;
else
- v4 = v1;
- std::complex<double> vInterp = v1 + dFracView * (v2 - v1) + dFracDet * (v4 - v1);
+ v3 = ppcDetValue[0][iFloorDet+1];
+ std::complex<double> vInterp = (1 - dFracView) * (1 - dFracDet) * v1 +
+ dFracView * (1 - dFracDet) * v2 + dFracView * dFracDet * v3 +
+ dFracDet * (1 - dFracView) * v4;
v[ix][iy] = vInterp.real();
if (vImag)
vImag[ix][iy] = vInterp.imag();
}
}
-bool
-Projections::convertFFTPolar (ImageFile& rIF, int iInterpolationID, int iZeropad)
-{
- unsigned int nx = rIF.nx();
- unsigned int ny = rIF.ny();
- ImageFileArray v = rIF.getArray();
- if (! rIF.isComplex())
- rIF.convertRealToComplex();
- ImageFileArray vImag = rIF.getImaginaryArray();
-
- if (! v || nx == 0 || ny == 0)
- return false;
-
- Array2d<double> adView (nx, ny);
- Array2d<double> adDet (nx, ny);
- double** ppdView = adView.getArray();
- double** ppdDet = adDet.getArray();
-
- std::complex<double>** ppcDetValue = new std::complex<double>* [m_nView];
- unsigned int iView;
- double* pdDet = new double [m_nDet];
- fftw_complex* pcIn = new fftw_complex [m_nDet];
- 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);
- 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);
- }
-
- fftw_destroy_plan (plan);
-
- delete [] pcIn;
-
- calcArrayPolarCoordinates (nx, ny, ppdView, ppdDet);
-
- 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;
-}
-
-
projects / ctsim.git / blobdiff