r636: Optimized Rebinning, Added Reconstruct with Rebinning option

[ctsim.git] / libctsim / projections.cpp

diff --git a/libctsim/projections.cpp b/libctsim/projections.cpp
index 60f9467e06cbb2426ed2db355b9c163248ffdbd7..3eab99c4fe9a8a8725bc2108492eb86302e50a64 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.62 2001/03/13 04:44:25 kevin Exp $
+**  $Id: projections.cpp,v 1.63 2001/03/13 08:24:41 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
@@ -895,59 +895,53 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa
   m_dFocalLength = 510;
   m_dSourceDetectorLength = 890;
   m_detInc = convertDegreesToRadians (3.06976 / 60);
+  m_dFanBeamAngle = (iNDets + 1) * m_detInc;
   m_detStart = -(m_dFanBeamAngle / 2);
   m_rotInc = TWOPI / static_cast<double>(iNViews);
   m_rotStart = HALFPI;
-  m_dFanBeamAngle = (iNDets + 1) * m_detInc;
   m_dViewDiameter = sin (m_dFanBeamAngle / 2) * m_dFocalLength * 2;
 
   if (! ((iNViews == 750 && lDataLength == 1560000L) || (iNViews == 950 && lDataLength == 1976000L) 
                 || (iNViews == 1500 && lDataLength == 3120000)))
     return false;
 
-  int iCenter = (iNDets - 1) / 2; // change from (Nm+1)/2 because of 0 vs. 1 indexing
+  double dCenter = (iNDets - 1.) / 2.; // change from (Nm+1)/2 because of 0 vs. 1 indexing
   double* pdCosScale = new double [iNDets];
   for (int i = 0; i < iNDets; i++)
-    pdCosScale[i] = cos ((i - iCenter) * m_detInc);
+    pdCosScale[i] = 1. / cos ((i - dCenter) * m_detInc);
 
   long lDataPos = 0;
   for (int iv = 0; iv < iNViews; iv++) {
     unsigned char* pArgBase = pData + lDataPos;
-    unsigned char* p = pArgBase+0;
-    SwapBytes4IfLittleEndian (p);
+    unsigned char* p = pArgBase+0; SwapBytes4IfLittleEndian (p);
     long lProjNumber = *reinterpret_cast<long*>(p);
 
-    p = pArgBase+20;
-    SwapBytes4IfLittleEndian (p);
+    p = pArgBase+20;  SwapBytes4IfLittleEndian (p);
     long lEscale = *reinterpret_cast<long*>(p);
 
-    p = pArgBase+28;
-    SwapBytes4IfLittleEndian (p);
+    p = pArgBase+28;  SwapBytes4IfLittleEndian (p);
     long lTime = *reinterpret_cast<long*>(p);
 
-    p = pArgBase + 4;
-    SwapBytes4IfLittleEndian (p);
+    p = pArgBase + 4; SwapBytes4IfLittleEndian (p);
     double dAlpha = *reinterpret_cast<float*>(p) + HALFPI;
 
-    p = pArgBase+12;
-    SwapBytes4IfLittleEndian (p);
+    p = pArgBase+12; SwapBytes4IfLittleEndian (p);
     double dAlign = *reinterpret_cast<float*>(p);
 
-    p = pArgBase + 16;
-    SwapBytes4IfLittleEndian (p);
+    p = pArgBase + 16; SwapBytes4IfLittleEndian (p);
     double dMaxValue = *reinterpret_cast<float*>(p);
 
     DetectorArray& detArray = getDetectorArray (iv);
     detArray.setViewAngle (dAlpha);
     DetectorValue* detval = detArray.detValues();
 
-    double dViewScale = 2294.4871 * ::pow (2.0, -lEscale);
+    double dViewScale = 1. / (2294.4871 * ::pow (2.0, -lEscale));
     lDataPos += 32;
     for (int id = 0; id < iNDets; id++) {
-      int iV = pData[lDataPos+1] + 256 * pData[lDataPos];
+      int iV = pData[lDataPos+1] + (pData[lDataPos] << 8);
       if (iV > 32767)   // two's complement signed conversion
         iV = iV - 65536;
-      detval[id] = iV / (dViewScale * pdCosScale[id]);
+      detval[id] = iV * dViewScale * pdCosScale[id];
       lDataPos += 2;
     }
   }
@@ -957,10 +951,10 @@ Projections::initFromSomatomAR_STAR (int iNViews, int iNDets, unsigned char* pDa
 }
 
 Projections*
-Projections::interpolateToParallel ()
+Projections::interpolateToParallel () const
 {
   if (m_geometry == Scanner::GEOMETRY_PARALLEL)
-    return this;
+    return const_cast<Projections*>(this);
 
   int nDet = m_nDet;
   int nView = m_nView;
@@ -1045,9 +1039,9 @@ Projections::interpolateToParallel ()
 //
 ///////////////////////////////////////////////////////////////////////////////
 
-ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange)
+ParallelRaysums::ParallelRaysums (const Projections* pProjections, int iThetaRange)
 : m_iNumCoordinates(0), m_iNumView(pProjections->nView()), m_iNumDet(pProjections->nDet()),
-  m_iThetaRange (iThetaRange)
+  m_iThetaRange (iThetaRange), m_pCoordinates(NULL)
 {
   int iGeometry = pProjections->geometry();
   double dDetInc = pProjections->detInc();
@@ -1055,9 +1049,10 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange)
   double dFocalLength = pProjections->focalLength();
 
   m_iNumCoordinates =  m_iNumView * m_iNumDet;
+  m_pCoordinates = new ParallelRaysumCoordinate [m_iNumCoordinates];
   m_vecpCoordinates.reserve (m_iNumCoordinates);
   for (int i = 0; i < m_iNumCoordinates; i++)
-    m_vecpCoordinates[i] = new ParallelRaysumCoordinate;
+    m_vecpCoordinates[i] = m_pCoordinates + i;
 
   int iCoordinate = 0;
   for (int iV = 0; iV < m_iNumView; iV++) {
@@ -1096,8 +1091,7 @@ ParallelRaysums::ParallelRaysums (Projections* pProjections, int iThetaRange)
 
 ParallelRaysums::~ParallelRaysums()
 {
-  for (int i = 0; i < m_iNumCoordinates; i++)
-    delete m_vecpCoordinates[i];
+  delete m_pCoordinates;
 }
 
 ParallelRaysums::CoordinateContainer&