X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fprojections.cpp;h=2ca1ce38e12af9fbf708595028df1f9dc9c3fbc8;hp=d715b73ce04188de6dab100e33f20f56a38f00dd;hb=c0f892798de8f89715266150f7d8e413f2cf29fe;hpb=7c456cf1278a3dfd6149094292dffea81e795cc6 diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp index d715b73..2ca1ce3 100644 --- a/libctsim/projections.cpp +++ b/libctsim/projections.cpp @@ -8,7 +8,7 @@ ** 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.77 2002/05/28 18:43:16 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 @@ -238,30 +238,30 @@ Projections::Helical180LI_Equiangular(int interpView) // is there enough data in the data set? Should have 2(Pi+fanAngle) // coverage minimum - if ( m_nView < static_cast((2*( M_PI + fanAngle ) ) / dbeta) -1 ){ + if ( m_nView < static_cast((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 ((M_PI+fanAngle)/dbeta); + interpView = static_cast ((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((M_PI+fanAngle)/dbeta); + offsetView = interpView - static_cast((PI+fanAngle)/dbeta); } - int last_interp_view = static_cast ((M_PI+fanAngle)/dbeta); + int last_interp_view = static_cast ((PI+fanAngle)/dbeta); // make a new array for data... @@ -282,20 +282,20 @@ Projections::Helical180LI_Equiangular(int interpView) 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((beta + 2*gamma - M_PI)/dbeta); - //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); - newiView = nearest(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - M_PI)/dbeta); + // newiView = nearest((beta + 2*gamma - PI)/dbeta); + //newiView = static_cast(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); + newiView = nearest(( (iView*dbeta) + 2*(iDet-(m_nDet-1)/2)*dgamma - PI)/dbeta); } #ifdef DEBUG @@ -304,43 +304,43 @@ Projections::Helical180LI_Equiangular(int interpView) #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 @@ -374,7 +374,7 @@ Projections::HalfScanFeather(void) double fanAngle = m_dFanBeamAngle; // is there enough data? - if ( m_nView != static_cast(( M_PI+fanAngle ) / dbeta) +1 ){ + if ( m_nView != static_cast(( PI+fanAngle ) / dbeta) +1 ){ std::cerr << "Data set does seem have enough data to be a halfscan data set" << std::endl; return (1); } @@ -392,12 +392,12 @@ Projections::HalfScanFeather(void) 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((beta2 + 2*gamma2 - M_PI)/dbeta); + //iView1 = nearest((beta2 + 2*gamma2 - PI)/dbeta); iView1 = nearest(( (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(); @@ -410,10 +410,10 @@ Projections::HalfScanFeather(void) 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); @@ -427,7 +427,7 @@ Projections::HalfScanFeather(void) } } // 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++) { @@ -650,7 +650,7 @@ Projections::read (const char* filename) #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())