r140: *** empty log message ***

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 1f6e92c1bece4c526d80abdad76e47d1d70a7f75..d3cb2035dc125af3415006a2d204cc61b5cc7907 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.4 2000/06/25 17:32:24 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.6 2000/07/09 08:16:17 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
@@ -87,6 +87,8 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
   else if (m_idBackproject == BPROJ_IDIFF2)
     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+  else if (m_idBackproject == BPROJ_IDIFF3)
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation));
   else {
     m_fail = true;
     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
@@ -112,6 +114,8 @@ Backprojector::convertBackprojectNameToID (const char* const backprojName)
     backprojID = BPROJ_DIFF2;
   else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
     backprojID = BPROJ_IDIFF2;
+  else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0)
+    backprojID = BPROJ_IDIFF3;
 
   return (backprojID);
 }
@@ -131,6 +135,8 @@ Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
     bprojName = BPROJ_DIFF2_STR;
   else if (bprojID == BPROJ_IDIFF2)
     bprojName = BPROJ_IDIFF2_STR;
+  else if (bprojID == BPROJ_IDIFF3)
+    bprojName = BPROJ_IDIFF3_STR;
 
   return (bprojName);
 }
@@ -146,6 +152,8 @@ Backprojector::convertInterpolationNameToID (const char* const interpName)
     interpID = INTERP_NEAREST;
   else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
     interpID = INTERP_LINEAR;
+  else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATE_STR) == 0)
+    interpID = INTERP_FREQ_PREINTERPOLATE;
 #if HAVE_BSPLINE_INTERP
   else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
     interpID = INTERP_BSPLINE;
@@ -174,6 +182,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
     return (INTERP_NEAREST_STR);
   else if (interpID == INTERP_LINEAR)
     return (INTERP_LINEAR_STR);
+  else if (interpID == INTERP_FREQ_PREINTERPOLATE)
+    return (INTERP_FREQ_PREINTERPOLATE_STR);
 #if HAVE_BSPLINE_INTERP
   else if (interpID == INTERP_BSPLINE)
     return (INTERP_BSPLINE_STR);
@@ -465,12 +475,12 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
 {
   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 
-  kint32 scale = 1 << 16;
-  double dScale = scale;
-  kint32 halfScale = scale / 2;
+  static const kint32 scale = 1 << 16;
+  static const double dScale = scale;
+  static const kint32 halfScale = scale / 2;
 
-  kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
-  kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+  const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+  const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
 
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
@@ -500,8 +510,59 @@ BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const do
        if (iDetPos < 0 || iDetPos >= nDet - 1)
            errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
        else
-         *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+         *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
       }
     }  // end for y
   }    // end for x
 }
+
+// CLASS IDENTICATION
+//   BackprojectIntDiff3
+//
+// PURPOSE
+//   Highly optimized version of BackprojectIntDiff2
+
+void
+BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+  double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
+  static const int scaleShift = 16;
+  static const kint32 scale = (1 << scaleShift);
+  static const kint32 scaleBitmask = scale - 1;
+  static const double dScale = scale;
+  static const kint32 halfScale = scale / 2;
+  static const double dInvScale = 1. / scale;
+
+  const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+  const kint32 det_dy = nearest<kint32> (yInc * cos (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);
+       
+  for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
+    kint32 curDetPos = detPosColStart;
+    ImageFileColumn pImCol = v[ix];
+
+    if (interpType == Backprojector::INTERP_NEAREST) {
+      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+       int iDetPos = (curDetPos + halfScale) >> 16;
+       assert(iDetPos >= 0 && iDetPos < nDet);
+       *pImCol++ += filteredProj[iDetPos];
+      }        // end for iy
+    } else if (interpType == Backprojector::INTERP_LINEAR) {
+      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+       kint32 iDetPos = curDetPos >> scaleShift;
+       kint32 detRemainder = curDetPos & scaleBitmask;
+#if 0
+       double frac = detRemainder * (1. / 65536.);
+       assert(iDetPos >= 0 && iDetPos < nDet - 1);
+       *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]);
+#else
+       assert(iDetPos >= 0 && iDetPos < nDet - 1);
+       const double* const detPointer = &filteredProj[iDetPos];
+       *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale;
+#endif
+      }        // end for iy
+    } //end linear
+  } // end for ix
+}