** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: projections.cpp,v 1.75 2001/09/24 11:20:08 kevin Exp $
+** $Id: projections.cpp,v 1.80 2002/06/27 03:19:23 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 "ct.h"
+#include <ctime>\r
const kuint16 Projections::m_signature = ('P'*256 + 'J');
// is there enough data in the data set? Should have 2(Pi+fanAngle)
// coverage minimum
- if ( m_nView < static_cast<int>((2*( M_PI + fanAngle ) ) / dbeta) -1 ){
+ if ( m_nView < static_cast<int>((2*( PI + fanAngle ) ) / dbeta) -1 ){
std::cerr << "Data set does not include 360 +2*FanBeamAngle views"
<< std::endl;
return (1);
}
- if (interpView < 0) // use default position at M_PI+fanAngle
+ if (interpView < 0) // use default position at PI+fanAngle
{
- interpView = static_cast<int> ((M_PI+fanAngle)/dbeta);
+ interpView = static_cast<int> ((PI+fanAngle)/dbeta);
}
else
{
- // check if there is M_PI+fanAngle data on either side of the
+ // check if there is PI+fanAngle data on either side of the
// of the specified image plane
- if ( interpView*dbeta < M_PI+fanAngle ||
- interpView*dbeta + M_PI + fanAngle > m_nView*dbeta)
+ if ( interpView*dbeta < PI+fanAngle ||
+ interpView*dbeta + PI + fanAngle > m_nView*dbeta)
{
std::cerr << "There isn't PI+fanAngle of data on either side of the requested interpolation view" << std::endl;
return(1);
}
- offsetView = interpView - static_cast<int>((M_PI+fanAngle)/dbeta);
+ offsetView = interpView - static_cast<int>((PI+fanAngle)/dbeta);
}
- int last_interp_view = static_cast<int> ((M_PI+fanAngle)/dbeta);
+ int last_interp_view = static_cast<int> ((PI+fanAngle)/dbeta);
// make a new array for data...
for ( int iDet = 0; iDet < m_nDet; iDet++) {
double gamma = (iDet -(m_nDet-1)/2)* dgamma ;
int newiView, newiDet;
- if (beta < M_PI+fanAngle) { //if (M_PI +fanAngle - beta > dbeta )
+ if (beta < PI+fanAngle) { //if (PI +fanAngle - beta > dbeta )
//newbeta = beta;
//newgamma = gamma;
newiDet = iDet;
newiView = iView;
}
- else // (beta > M_PI+fanAngle)
+ else // (beta > PI+fanAngle)
{
//newbeta = beta +2*gamma - 180;
//newgamma = -gamma;
newiDet = -iDet + (m_nDet -1);
- // newiView = nearest<int>((beta + 2*gamma - M_PI)/dbeta);
- //newiView = static_cast<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
- newiView = nearest<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
+ // newiView = nearest<int>((beta + 2*gamma - PI)/dbeta);
+ //newiView = static_cast<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta);
+ newiView = nearest<int>(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta);
}
#ifdef DEBUG
#endif
if ( ( beta > fanAngle - 2*gamma)
- && ( beta < 2*M_PI + fanAngle -2*gamma) )
+ && ( beta < 2*PI + fanAngle -2*gamma) )
{ // not in region 1 or 8
DetectorValue* detval = (m_projData[iView+offsetView])->detValues();
DetectorValue* newdetval = (newdetarray[newiView])->detValues();
if ( beta > fanAngle - 2*gamma
&& beta <= 2*fanAngle ) { // in region 2
newdetval[newiDet] +=
- (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+ (beta +2*gamma - fanAngle)/(PI+2*gamma)
* detval[iDet];
} else if ( beta > 2*fanAngle
- && beta <= M_PI - 2*gamma) { // in region 3
+ && beta <= PI - 2*gamma) { // in region 3
newdetval[newiDet] +=
- (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+ (beta +2*gamma - fanAngle)/(PI+2*gamma)
* detval[iDet];
}
- else if ( beta > M_PI -2*gamma
- && beta <= M_PI + fanAngle ) { // in region 4
+ else if ( beta > PI -2*gamma
+ && beta <= PI + fanAngle ) { // in region 4
newdetval[newiDet] +=
- (beta +2*gamma - fanAngle)/(M_PI+2*gamma)
+ (beta +2*gamma - fanAngle)/(PI+2*gamma)
* detval[iDet];
}
- else if ( beta > M_PI + fanAngle
- && beta <= M_PI +2*fanAngle -2*gamma) { // in region 5
+ else if ( beta > PI + fanAngle
+ && beta <= PI +2*fanAngle -2*gamma) { // in region 5
newdetval[newiDet] +=
- (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma)
+ (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma)
*detval[iDet];
}
- else if ( beta > M_PI +2*fanAngle -2*gamma
- && beta <= 2*M_PI) { // in region 6
+ else if ( beta > PI +2*fanAngle -2*gamma
+ && beta <= 2*PI) { // in region 6
newdetval[newiDet] +=
- (2*M_PI - beta - 2*gamma + fanAngle)/(M_PI-2*gamma)
+ (2*PI - beta - 2*gamma + fanAngle)/(PI-2*gamma)
*detval[iDet];
}
- else if ( beta > 2*M_PI
- && beta <= 2*M_PI + fanAngle -2*gamma){ // in region 7
+ else if ( beta > 2*PI
+ && beta <= 2*PI + fanAngle -2*gamma){ // in region 7
newdetval[newiDet] +=
- (2*M_PI - beta -2*gamma + fanAngle)/(M_PI-2*gamma)
+ (2*PI - beta -2*gamma + fanAngle)/(PI-2*gamma)
*detval[iDet];
}
else
double fanAngle = m_dFanBeamAngle;
// is there enough data?
- if ( m_nView != static_cast<int>(( M_PI+fanAngle ) / dbeta) +1 ){
+ if ( m_nView != static_cast<int>(( PI+fanAngle ) / dbeta) +1 ){
std::cerr << "Data set does seem have enough data to be a halfscan data set" << std::endl;
return (1);
}
double beta2 = iView2 * dbeta;
for ( int iDet2 = 0; iDet2 < m_nDet; iDet2++) {
double gamma2 = (iDet2 -(m_nDet-1)/2)* dgamma ;
- if ( ( beta2 >= M_PI - 2*gamma2) ) { // in redundant data region
+ if ( ( beta2 >= PI - 2*gamma2) ) { // in redundant data region
int iView1, iDet1;
iDet1 = (m_nDet -1) - iDet2;
- //iView1 = nearest<int>((beta2 + 2*gamma2 - M_PI)/dbeta);
+ //iView1 = nearest<int>((beta2 + 2*gamma2 - PI)/dbeta);
iView1 = nearest<int>(( (iView2*dbeta)
- + 2*(iDet2-(m_nDet-1)/2)*dgamma - M_PI)/dbeta);
+ + 2*(iDet2-(m_nDet-1)/2)*dgamma - PI)/dbeta);
DetectorValue* detval2 = (m_projData[iView2])->detValues();
gamma1 = -gamma2;
if ( beta1 <= (fanAngle - 2*gamma1) )
x = beta1 / ( fanAngle - 2*gamma1);
- else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= M_PI - 2*gamma1)
+ else if ( (fanAngle - 2*gamma1 <= beta1 ) && beta1 <= PI - 2*gamma1)
x = 1;
- else if ( (M_PI - 2*gamma1 <= beta1 ) && ( beta1 <=M_PI + fanAngle) )
- x = (M_PI +fanAngle - beta1)/(fanAngle + 2*gamma1);
+ else if ( (PI - 2*gamma1 <= beta1 ) && ( beta1 <=PI + fanAngle) )
+ x = (PI +fanAngle - beta1)/(fanAngle + 2*gamma1);
else {
std::cerr << "Shouldn't be here!"<< std::endl;
return(4);
}
}
// heuristic scaling, why this factor?
- double scalefactor = m_nView * m_rotInc / M_PI;
+ double scalefactor = m_nView * m_rotInc / PI;
for ( int iView = 0 ; iView < m_nView; iView++) {
DetectorValue* detval = (m_projData[iView])->detValues();
for ( int iDet = 0; iDet < m_nDet; iDet++) {
#ifdef MSVC
frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary);
#else
- frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary | std::ios::nocreate);
+ frnetorderstream fileRead (m_filename.c_str(), std::ios::in | std::ios::binary); // | std::ios::nocreate);
#endif
if (fileRead.fail())
double** ppdDet = adDet.getArray();
std::complex<double>** ppcDetValue = new std::complex<double>* [pProj->m_nView];
- unsigned int iView;
+ int iView;
for (iView = 0; iView < pProj->m_nView; iView++) {
ppcDetValue[iView] = new std::complex<double> [pProj->m_nDet];
DetectorValue* detval = pProj->getDetectorArray (iView).detValues();
- for (unsigned int iDet = 0; iDet < pProj->m_nDet; iDet++)
+ for (int iDet = 0; iDet < pProj->m_nDet; iDet++)
ppcDetValue[iView][iDet] = std::complex<double>(detval[iDet], 0);
}
int iZerosAdded = iNumInterpDetWithZeros - iInterpDet;
// For each view, interpolate to nx length, shift to center at origin, and FFt transform
- for (unsigned int iView = 0; iView < m_nView; iView++) {
+ for (int iView = 0; iView < m_nView; iView++) {
DetectorValue* detval = pProj->getDetectorArray(iView).detValues();
LinearInterpolator<DetectorValue> projInterp (detval, pProj->m_nDet);
- for (unsigned int iDet = 0; iDet < iInterpDet; iDet++) {
+ for (int iDet = 0; iDet < iInterpDet; iDet++) {
double dInterpPos = (m_nDet / 2.) + (iDet - dMidPoint) * dInterpScale;
pcIn[iDet].re = projInterp.interpolate (dInterpPos) * dInterpScale;
pcIn[iDet].im = 0;
Fourier::shuffleFourierToNaturalOrder (pcIn, iInterpDet);
if (iZerosAdded > 0) {
- for (unsigned int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++)
+ for (int iDet1 = iMidPoint; iDet1 < iInterpDet; iDet1++)
pcIn[iDet1+iZerosAdded] = pcIn[iDet1];
- for (unsigned int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++)
+ for (int iDet2 = iMidPoint; iDet2 < iMidPoint + iZerosAdded; iDet2++)
pcIn[iDet2].re = pcIn[iDet2].im = 0;
}
fftw_one (plan, pcIn, NULL);
ppcDetValue[iView] = new std::complex<double> [iNumInterpDetWithZeros];
- for (unsigned int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
+ for (int iD = 0; iD < iNumInterpDetWithZeros; iD++) {
ppcDetValue[iView][iD] = std::complex<double> (pcIn[iD].re * dFFTScale, pcIn[iD].im * dFFTScale);
}
{
typedef std::complex<double> complexValue;
- BilinearInterpolator<complexValue>* pBilinear;
+ BilinearInterpolator<complexValue>* pBilinear = NULL;
if (iInterpolationID == POLAR_INTERP_BILINEAR)
pBilinear = new BilinearInterpolator<complexValue> (ppcDetValue, nView, nDetWithZeros);
for (int iv = 0; iv < iNViews; iv++) {
unsigned char* pArgBase = pData + lDataPos;
unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p);
- long lProjNumber = *reinterpret_cast<long*>(p);
+ // long lProjNumber = *reinterpret_cast<long*>(p);
p = pArgBase+20; SwapBytes4IfLittleEndian (p);
long lEscale = *reinterpret_cast<long*>(p);
p = pArgBase+28; SwapBytes4IfLittleEndian (p);
- long lTime = *reinterpret_cast<long*>(p);
+ // long lTime = *reinterpret_cast<long*>(p);
p = pArgBase + 4; SwapBytes4IfLittleEndian (p);
double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
p = pArgBase+12; SwapBytes4IfLittleEndian (p);
- double dAlign = *reinterpret_cast<float*>(p);
+ // double dAlign = *reinterpret_cast<float*>(p);
p = pArgBase + 16; SwapBytes4IfLittleEndian (p);
- double dMaxValue = *reinterpret_cast<float*>(p);
+ // double dMaxValue = *reinterpret_cast<float*>(p);
DetectorArray& detArray = getDetectorArray (iv);
detArray.setViewAngle (dAlpha);
double dViewAngle = m_rotStart;
int iLastFloor = -1;
+ LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
for (int iV = 0; iV < pProjNew->nView(); iV++, dViewAngle += pProjNew->m_rotInc) {
DetectorValue* detValues = pProjNew->getDetectorArray (iV).detValues();
- LinearInterpolator<double> interp (pdThetaValuesForT, pdRaysumsForT, pProjNew->nView(), false);
detValues[iD] = interp.interpolate (dViewAngle, &iLastFloor);
}
}
///////////////////////////////////////////////////////////////////////////////
ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange)
-: m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
- m_iThetaRange (iThetaRange), m_pCoordinates(NULL)
+: m_pCoordinates(NULL), m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
+ m_iThetaRange (iThetaRange)
{
int iGeometry = pProjections->geometry();
double dDetInc = pProjections->detInc();