X-Git-Url: http://git.kpe.io/?a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=ef1ef12e655738f3b1e1b06bd13fbba148f38b28;hb=87ee1eb8f451d1f1b278c160e5f90ff20301560c;hp=26455facaa9facfd9d2629807a9d66ec82a7362a;hpb=8a7697ce57b56cdc43698cd1241ad98d49f9b5ac;p=ctsim.git

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 26455fa..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;
 }
 
@@ -572,17 +583,21 @@ void
 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 {
   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
+#if SIZEOF_LONG == 4
   static const int scaleShift = 16;
-  static const kint32 scale = (1 << scaleShift);
-  static const kint32 scaleBitmask = scale - 1;
-  static const kint32 halfScale = scale / 2;
+#elif SIZEOF_LONG == 8
+  static const int scaleShift = 32;
+#endif
+  static const long scale = (1L << scaleShift);
+  static const long scaleBitmask = scale - 1;
+  static const long halfScale = scale / 2;
   static const double dInvScale = 1. / scale;
   
-  const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
-  const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
+  const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
+  const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
   
   // calculate L for first point in image (0, 0) 
-  kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
+  long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
   
   double* deltaFilteredProj = NULL;  
   CubicPolyInterpolator* pCubicInterp = NULL;
@@ -598,31 +613,32 @@ 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) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        const int iDetPos = (curDetPos + halfScale) >> 16;
+        const int iDetPos = (curDetPos + halfScale) >> scaleShift;
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos];
       }	// 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) >> 16) * m_interpFactor;
+        const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos];
       }	// end for iy
     } else if (interpType == Backprojector::INTERP_LINEAR) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        const kint32 iDetPos = curDetPos >> scaleShift;
-        const kint32 detRemainder = curDetPos & scaleBitmask;
-        if (iDetPos >= 0 && iDetPos <= iLastDet)
-          *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+        const long iDetPos = curDetPos >> scaleShift;
+        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) {
-        *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
+        *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
       }
     } // end Cubic
   } // end for ix