r11223: Oct 15, 2006 Version 4.5.0

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index b7d320d49aab76f31668b96e3efbfd178f00289a..2d0d60cc718263746866659b727ba5ed7b28feeb 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -270,12 +270,11 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType
   detInc = proj.detInc();
   nDet = proj.nDet();
   iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection 
-  rotScale = proj.rotInc();
   
   if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
-    rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
+    rotScale = PI / proj.nView(); // scale by number of PI rotations
   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
-    rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
+    rotScale =  (2 * PI) / proj.nView(); // scale by number of 2PI rotations
   else
     sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
   
@@ -527,37 +526,48 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   double det_dy = yInc * sin (theta) / detInc;
   
   // calculate detPosition for first point in image (ix=0, iy=0) 
-  double detPosColStart = start_r * cos (theta - start_phi) / detInc;
+  double detPosColStart = iDetCenter + start_r * cos (theta - start_phi) / detInc;
   
   CubicPolyInterpolator* pCubicInterp = NULL;
-  if (interpType == Backprojector::INTERP_CUBIC)
+  double* deltaFilteredProj = NULL;  
+  if (interpType == Backprojector::INTERP_LINEAR) {
+    // precalculate scaled difference for linear interpolation
+    deltaFilteredProj = new double [nDet];
+    for (int i = 0; i < nDet - 1; i++)
+      deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
+    deltaFilteredProj[nDet - 1] = 0;  // last detector
+  } else if (interpType == Backprojector::INTERP_CUBIC) {
     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+  }
   
+  int iLastDet = nDet - 1;
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
     double curDetPos = detPosColStart;
     ImageFileColumn pImCol = v[ix];
     
     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
       if (interpType == Backprojector::INTERP_NEAREST) {
-        int iDetPos = iDetCenter + nearest<int> (curDetPos);   // calc index in the filtered raysum vector 
+        int iDetPos = nearest<int> (curDetPos);        // calc index in the filtered raysum vector 
         
         if (iDetPos >= 0 && iDetPos < nDet)
           *pImCol++ += filteredProj[iDetPos];
       } else if (interpType == Backprojector::INTERP_LINEAR) {
         double detPosFloor = floor (curDetPos);
-        int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
+        int iDetPos = static_cast<int>(detPosFloor);
         double frac = curDetPos - detPosFloor; // fraction distance from det 
-        if (iDetPos > 0 && iDetPos < nDet - 1)
-          *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
+        if (iDetPos >= 0 && iDetPos <= iLastDet)
+           *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
       } else if (interpType == Backprojector::INTERP_CUBIC) {
-        double p = iDetCenter + curDetPos;     // position along detector
+        double p = curDetPos;  // position along detector
         if (p >= 0 && p < nDet)
           *pImCol++  += pCubicInterp->interpolate (p);
       }
     }  // end for y
   }    // end for x
 
-  if (interpType == Backprojector::INTERP_CUBIC)
+  if (interpType == Backprojector::INTERP_LINEAR)
+    delete deltaFilteredProj;
+  else if (interpType == Backprojector::INTERP_CUBIC)
     delete pCubicInterp;
 }
 
@@ -577,7 +587,7 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
 #elif SIZEOF_LONG == 8
   static const int scaleShift = 32;
 #endif
-  static const long scale = (1 << scaleShift);
+  static const long scale = (1L << scaleShift);
   static const long scaleBitmask = scale - 1;
   static const long halfScale = scale / 2;
   static const double dInvScale = 1. / scale;
@@ -602,7 +612,7 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
   
   int iLastDet = nDet - 1;
   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
-    kint32 curDetPos = detPosColStart;
+    long curDetPos = detPosColStart;
     ImageFileColumn pImCol = v[ix];
     
     if (interpType == Backprojector::INTERP_NEAREST) {
@@ -610,19 +620,26 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
         const int iDetPos = (curDetPos + halfScale) >> scaleShift;
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos];
+        else
+          pImCol++;
+
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
         const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos];
+        else
+          pImCol++;
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_LINEAR) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
         const long iDetPos = curDetPos >> scaleShift;
-        const long detRemainder = curDetPos & scaleBitmask;
-        if (iDetPos >= 0 && iDetPos <= iLastDet)
-          *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+        if (iDetPos >= 0 && iDetPos <= iLastDet) {
+         const long detRemainder = curDetPos & scaleBitmask;
+         *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+       } else
+          pImCol++;
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_CUBIC) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {