** 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.8 2000/07/13 07:03:21 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
#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);
}
// 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;
}
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]";
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;
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);
// 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();
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)
// 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);
// 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;
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
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;
// 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