r142: *** empty log message ***

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index d3cb2035dc125af3415006a2d204cc61b5cc7907..ed388addbfb292718db205f1e5e13b2561a25c5f 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.6 2000/07/09 08:16:17 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.7 2000/07/11 10:32:44 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
@@ -26,18 +26,18 @@
 
 #include "ct.h"
 
-Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 {
   m_fail = false;
   m_pBackprojectImplem = NULL;
 
-  initBackprojector (proj, im, backprojName, interpName);
+  initBackprojector (proj, im, backprojName, interpName, interpFactor);
 }
 
 void 
 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
 {
-  if (m_pBackprojectImplem)
+  if (m_pBackprojectImplem != NULL)
     m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
 }
 
@@ -54,7 +54,7 @@ Backprojector::~Backprojector (void)
 //     and initializes the backprojector
 
 bool
-Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 {
   m_nameBackproject = backprojName;
   m_nameInterpolation = interpName;
@@ -78,17 +78,17 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
   }
 
   if (m_idBackproject == BPROJ_TRIG)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_TABLE)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_DIFF)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_DIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_IDIFF2)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
   else if (m_idBackproject == BPROJ_IDIFF3)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
   else {
     m_fail = true;
     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
@@ -152,8 +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;
+  else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0)
+    interpID = INTERP_FREQ_PREINTERPOLATION;
 #if HAVE_BSPLINE_INTERP
   else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
     interpID = INTERP_BSPLINE;
@@ -182,8 +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);
+  else if (interpID == INTERP_FREQ_PREINTERPOLATION)
+    return (INTERP_FREQ_PREINTERPOLATION_STR);
 #if HAVE_BSPLINE_INTERP
   else if (interpID == INTERP_BSPLINE)
     return (INTERP_BSPLINE_STR);
@@ -199,8 +199,8 @@ Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
 // PURPOSE
 //   Pure virtual base class for all backprojectors.
 
-Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType)
-  : proj(proj), im(im), interpType(interpType)
+Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor)
+    : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
 {
   detInc = proj.detInc();
   nDet = proj.nDet();
@@ -219,9 +219,6 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, const Backproj
 
   xInc = (xMax - xMin) / nx;   // size of cells
   yInc = (yMax - yMin) / ny;
-
-  if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR)
-    sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
 }
 
 Backproject::~Backproject (void)
@@ -298,8 +295,8 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
 // PURPOSE
 //   Precalculates trigometric function value for each point in image for backprojection.
 
-BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
-  : Backproject::Backproject (proj, im, interpType)
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
+  : Backproject::Backproject (proj, im, interpType, interpFactor)
 {
   arrayR.initSetSize (nx, ny);
   arrayPhi.initSetSize (nx, ny);
@@ -360,8 +357,8 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 //   Iterates in x & y direction by adding difference in L position
 
-BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
-  :  Backproject::Backproject (proj, im, interpType)
+BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
+  :  Backproject::Backproject (proj, im, interpType, interpFactor)
 {
   // calculate center of first pixel v[0][0] 
   double x = xMin + xInc / 2;
@@ -458,7 +455,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl
        if (iDetPos < 0 || iDetPos >= nDet - 1)
            errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
        else
-         *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+         *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
       }
     }  // end for y
   }    // end for x
@@ -539,29 +536,35 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do
   // calculate L for first point in image (0, 0) 
   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
        
+  // precalculate scaled difference for linear interpolation
+  double deltaFilteredProj [nDet - 1];
+  if (interpType == Backprojector::INTERP_LINEAR) {
+    for (int i = 0; i < nDet - 1; i++)
+      deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
+  }
+
   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;
+       const int iDetPos = (curDetPos + halfScale) >> 16;
+       assert(iDetPos >= 0 && iDetPos < nDet);
+       *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;
        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.);
+       const kint32 iDetPos = curDetPos >> scaleShift;
+       const kint32 detRemainder = curDetPos & scaleBitmask;
        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
+       *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
       }        // end for iy
     } //end linear
   } // end for ix