r7929: fix backprojection

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index b7d320d49aab76f31668b96e3efbfd178f00289a..d38b7189f6b6c878b09ec0549dca5daf242ce24e 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -527,37 +527,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;
 }
 
@@ -602,7 +613,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) {
@@ -622,7 +633,8 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
         const long iDetPos = curDetPos >> scaleShift;
         const long detRemainder = curDetPos & scaleBitmask;
         if (iDetPos >= 0 && iDetPos <= iLastDet)
-          *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+           
       }        // end for iy
     } else if (interpType == Backprojector::INTERP_CUBIC) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {