** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.30 2001/03/07 19:02:44 kevin Exp $
+** $Id$
**
** 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"
+#include "interpolator.h"
const int Backprojector::BPROJ_INVALID = -1;
const int Backprojector::BPROJ_TRIG = 0;
const char* const Backprojector::s_aszBackprojectName[] =
{
- {"trig"},
- {"table"},
- {"diff"},
- {"idiff"},
+ "trig",
+ "table",
+ "diff",
+ "idiff",
};
const char* const Backprojector::s_aszBackprojectTitle[] =
{
- {"Direct Trigometric"},
- {"Trigometric Table"},
- {"Difference Iteration"},
- {"Integer Difference Iteration"},
+ "Direct Trigometric",
+ "Trigometric Table",
+ "Difference Iteration",
+ "Integer Difference Iteration",
};
const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
const char* const Backprojector::s_aszInterpName[] =
{
- {"nearest"},
- {"linear"},
- {"cubic"},
+ "nearest",
+ "linear",
+ "cubic",
#if HAVE_FREQ_PREINTERP
- {"freq_preinterpolationj"},
+ "freq_preinterpolationj",
#endif
#if HAVE_BSPLINE_INTERP
- {"bspline"},
- {"1bspline"},
- {"2bspline"},
- {"3bspline"},
+ "bspline",
+ "1bspline",
+ "2bspline",
+ "3bspline",
#endif
};
const char* const Backprojector::s_aszInterpTitle[] =
{
- {"Nearest"},
- {"Linear"},
- {"Cubic"},
+ "Nearest",
+ "Linear",
+ "Cubic",
#if HAVE_FREQ_PREINTERP
- {"Frequency Preinterpolation"},
+ "Frequency Preinterpolation",
#endif
#if HAVE_BSPLINE_INTERP
- {"B-Spline"},
- {"B-Spline 1st Order"},
- {"B-Spline 2nd Order"},
- {"B-Spline 3rd Order"},
+ "B-Spline",
+ "B-Spline 1st Order",
+ "B-Spline 2nd Order",
+ "B-Spline 3rd Order",
#endif
};
-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
// 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;
}
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]";
// 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();
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();
}
// 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());
// 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;
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;
}
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;
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