r627: no message

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index b3d6d35b2df647b5595d60c01bebcf30771af957..d6a484335943945d8f4ffa02a6bb77e495412fb7 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.55 2001/03/10 23:14:16 kevin Exp $
+**  $Id: projections.cpp,v 1.56 2001/03/10 23:56:58 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
@@ -970,7 +970,7 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa
 
 
 ParallelRaysums::ParallelRaysums (Projections* pProjections)
-: m_ppCoordinates(NULL), m_iNumCoordinates(0)
+: m_iNumCoordinates(0)
 {
   int nDet = pProjections->nDet();
   int nView = pProjections->nView();
@@ -980,30 +980,33 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections)
   double dFocalLength = pProjections->focalLength();
 
   m_iNumCoordinates =  nDet * nView;
-  m_ppCoordinates = new ParallelRaysumCoordinate* [m_iNumCoordinates];
+  m_vecpCoordinates.reserve (m_iNumCoordinates);
   for (int i = 0; i < m_iNumCoordinates; i++)
-    m_ppCoordinates[i] = new ParallelRaysumCoordinate;
-
-  ParallelRaysumCoordinate** ppCurrentCoordinate = m_ppCoordinates;
+    m_vecpCoordinates[i] = new ParallelRaysumCoordinate;
 
+  int iCoordinate = 0;
   for (int iV = 0; iV < nView; iV++) {
     double dViewAngle = pProjections->getDetectorArray(iV).viewAngle();
+
+    double dDetPos = dDetStart;
     for (int iD = 0; iD < nDet; iD++) {
-      ParallelRaysumCoordinate* pC = *ppCurrentCoordinate;
+      ParallelRaysumCoordinate* pC = m_vecpCoordinates[iCoordinate++];
+
       if (iGeometry == Scanner::GEOMETRY_PARALLEL) {
-        pC->m_dTheta = dViewAngle;
-        pC->m_dT = dDetStart + (iD * dDetInc);
+        pC->m_dTheta = normalizeAngle (dViewAngle);
+        pC->m_dT = dDetPos;
+
       } else if (iGeometry == Scanner::GEOMETRY_EQUILINEAR) {
-        double dDetPos = dDetStart + (iD * dDetInc);
         double dFanAngle = atan (dDetPos / pProjections->sourceDetectorLength());
-        pC->m_dTheta = dViewAngle + dFanAngle;
+        pC->m_dTheta = normalizeAngle (dViewAngle + dFanAngle);
         pC->m_dT = dFocalLength * sin(dFanAngle);        
+
       } else if (iGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
-        double dFanAngle = dDetStart + (iD * dDetInc);
-        pC->m_dTheta = dViewAngle + dFanAngle;
-        pC->m_dT = dFocalLength * sin(dFanAngle);        
+        // fan angle is same as dDetPos
+        pC->m_dTheta = normalizeAngle (dViewAngle + dDetPos);
+        pC->m_dT = dFocalLength * sin (dDetPos);        
       }
-      ++ppCurrentCoordinate;
+      dDetPos += dDetInc;
     }
   }
 }
@@ -1011,23 +1014,49 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections)
 ParallelRaysums::~ParallelRaysums()
 {
   for (int i = 0; i < m_iNumCoordinates; i++)
-    delete m_ppCoordinates[i];
+    delete m_vecpCoordinates[i];
+}
 
-  delete m_ppCoordinates;
+ParallelRaysums::CoordinateContainer&
+ParallelRaysums::getSortedByTheta()
+{
+  if (m_vecpSortedByTheta.size() == 0) {
+    m_vecpSortedByTheta.reserve (m_iNumCoordinates);
+    for (int i = 0; i < m_iNumCoordinates; i++)
+      m_vecpSortedByTheta[i] = m_vecpCoordinates[i];
+    std::sort (m_vecpSortedByTheta.begin(), m_vecpSortedByTheta.end(), ParallelRaysumCoordinate::compareByTheta);
+  }
+
+  return m_vecpSortedByTheta;
+}
+
+ParallelRaysums::CoordinateContainer&
+ParallelRaysums::getSortedByT()
+{
+  if (m_vecpSortedByT.size() == 0) {
+    m_vecpSortedByT.reserve (m_iNumCoordinates);
+    for (int i = 0; i < m_iNumCoordinates; i++)
+      m_vecpSortedByT[i] = m_vecpCoordinates[i];
+    std::sort (m_vecpSortedByT.begin(), m_vecpSortedByT.end(), ParallelRaysumCoordinate::compareByT);
+  }
+
+  return m_vecpSortedByT;
 }
 
+
 void
 ParallelRaysums::getLimits (double* dMinT, double* dMaxT, double* dMinTheta, double* dMaxTheta) const
 {
   if (m_iNumCoordinates <= 0)
     return;
 
-  *dMinT = *dMaxT = m_ppCoordinates[0]->m_dT;
-  *dMinTheta = *dMaxTheta = m_ppCoordinates[0]->m_dTheta;
+  *dMinT = *dMaxT = m_vecpCoordinates[0]->m_dT;
+  *dMinTheta = *dMaxTheta = m_vecpCoordinates[0]->m_dTheta;
 
   for (int i = 0; i < m_iNumCoordinates; i++) {
-    double dT = m_ppCoordinates[i]->m_dT;
-    double dTheta = m_ppCoordinates[i]->m_dTheta;
+    double dT = m_vecpCoordinates[i]->m_dT;
+    double dTheta = m_vecpCoordinates[i]->m_dTheta;
+
     if (dT < *dMinT)
       *dMinT = dT;
     else if (dT > *dMaxT)