projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r2088: *** empty log message ***
[ctsim.git] / libctsim / projections.cpp
diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index d715b73ce04188de6dab100e33f20f56a38f00dd..2ca1ce38e12af9fbf708595028df1f9dc9c3fbc8 100644 (file)
--- 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<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...
@@ -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<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
@@ -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<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);
    }
@@ -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<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();
@@ -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())
Computed Tomography simulator and reconstructor