r4228: Interpolation fixes

[ctsim.git] / libctsim / backprojectors.cpp

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index bab1dc7b18fc75e349cc961c43165f884545ed87..7551ff9a1d1e99a1e65336b4e726af20c9f26a12 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -8,7 +8,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: backprojectors.cpp,v 1.29 2001/03/01 07:30:49 kevin Exp $
+**  $Id: backprojectors.cpp,v 1.33 2003/03/23 18:37:42 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
@@ -25,6 +25,7 @@
 ******************************************************************************/
 
 #include "ct.h"
+#include "interpolator.h"
 
 const int Backprojector::BPROJ_INVALID = -1;
 const int Backprojector::BPROJ_TRIG = 0;
@@ -98,12 +99,13 @@ const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const
 
 
 
-Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
+Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, 
+                              const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
 {
   m_fail = false;
   m_pBackprojectImplem = NULL;
   
-  initBackprojector (proj, im, backprojName, interpName, interpFactor);
+  initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
 }
 
 void 
@@ -133,7 +135,8 @@ Backprojector::~Backprojector ()
 //     and initializes the backprojector
 
 bool
-Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
+Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, 
+                                  const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
 {
   m_nameBackproject = backprojName;
   m_nameInterpolation = interpName;
@@ -157,18 +160,18 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const
   }
   
   if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI));
   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) 
-    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
+    m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI));
   else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
     if (m_idBackproject == BPROJ_TRIG)
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI));
     else if (m_idBackproject == BPROJ_TABLE)
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI));
     else if (m_idBackproject == BPROJ_DIFF)
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI));
     else if (m_idBackproject == BPROJ_IDIFF)
-      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor));
+      m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI));
   } else {
     m_fail = true;
     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
@@ -260,7 +263,8 @@ Backprojector::convertInterpIDToTitle (const int interpID)
 // PURPOSE
 //   Pure virtual base class for all backprojectors.
 
-Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
+Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, 
+                          const ReconstructionROI* pROI)
 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false)
 {
   detInc = proj.detInc();
@@ -285,9 +289,34 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType
   yMin = -proj.phmLen() / 2;
   yMax = yMin + proj.phmLen();
   
+  if (pROI) {
+    if (pROI->m_dXMin > xMin)
+      xMin = pROI->m_dXMin;
+    if (pROI->m_dXMax < xMax)
+      xMax = pROI->m_dXMax;
+    if (pROI->m_dYMin > yMin)
+      yMin = pROI->m_dYMin;
+    if (pROI->m_dYMax < yMax)
+      yMax = pROI->m_dYMax;
+
+    if (xMin > xMax) {
+      double temp = xMin;
+      xMin = xMax;
+      xMax = temp;
+    }
+    if (yMin > yMax) {
+      double temp = yMin;
+      yMin = yMax;
+      yMax = temp;
+    }
+  }
+
   xInc = (xMax - xMin) / nx;   // size of cells
   yInc = (yMax - yMin) / ny;
   
+  im.setAxisIncrement (xInc, yInc);
+  im.setAxisExtent (xMin, xMax, yMin, yMax);
+
   m_dFocalLength = proj.focalLength();
   m_dSourceDetectorLength = proj.sourceDetectorLength();
 }
@@ -365,7 +394,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
         double frac = p - pFloor;      // fraction distance from det
         if (iDetPos >= 0 && iDetPos < nDet - 1)
           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
-      } else if (interpType = Backprojector::INTERP_CUBIC) {
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
         double p = iDetCenter + (L / detInc);  // position along detector
         if (p >= 0 && p < nDet)
           v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
@@ -384,8 +413,9 @@ 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, int interpType, const int interpFactor)
-: Backproject (proj, im, interpType, interpFactor)
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, 
+                                    const int interpFactor, const ReconstructionROI* pROI)
+: Backproject (proj, im, interpType, interpFactor, pROI)
 {
   arrayR.initSetSize (im.nx(), im.ny());
   arrayPhi.initSetSize (im.nx(), im.ny());
@@ -441,7 +471,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
         double frac = dPos - dPosFloor;        // fraction distance from det 
         if (iDetPos >= 0 && iDetPos < nDet - 1)
           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
-      } else if (interpType = Backprojector::INTERP_CUBIC) {
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
         double p = iDetCenter + (L / detInc);  // position along detector
         if (p >= 0 && p < nDet)
           pImCol[iy] += pCubicInterp->interpolate (p);
@@ -461,8 +491,9 @@ 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, int interpType, const int interpFactor)
-:  Backproject (proj, im, interpType, interpFactor)
+BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, 
+                                  const int interpFactor, const ReconstructionROI* pROI)
+:  Backproject (proj, im, interpType, interpFactor, pROI)
 {
   // calculate center of first pixel v[0][0] 
   double x = xMin + xInc / 2;
@@ -518,7 +549,7 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
         double frac = curDetPos - detPosFloor; // fraction distance from det 
         if (iDetPos > 0 && iDetPos < nDet - 1)
           *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
-      } else if (interpType = Backprojector::INTERP_CUBIC) {
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
         double p = iDetCenter + curDetPos;     // position along detector
         if (p >= 0 && p < nDet)
           *pImCol++  += pCubicInterp->interpolate (p);
@@ -589,7 +620,7 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
         if (iDetPos >= 0 && iDetPos <= iLastDet)
           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
       }        // end for iy
-    } else if (interpType = Backprojector::INTERP_CUBIC) {
+    } else if (interpType == Backprojector::INTERP_CUBIC) {
       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
       }