X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=ef1ef12e655738f3b1e1b06bd13fbba148f38b28;hb=9828c17ccbcc9f0f3fddaddc629dd7124c7cda45;hp=b7d320d49aab76f31668b96e3efbfd178f00289a;hpb=860a2e338a79d7a01d551bc5a70def43f60bdb77;p=ctsim.git

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index b7d320d..ef1ef12 100644
--- 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;
 }
 
@@ -577,7 +588,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 +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) {
@@ -620,9 +631,10 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
     } 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]);
+	}
       }	// end for iy
     } else if (interpType == Backprojector::INTERP_CUBIC) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {