** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.25 2001/02/09 01:54:20 kevin Exp $
+** $Id: backprojectors.cpp,v 1.26 2001/02/11 04:56:37 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
const int Backprojector::BPROJ_TRIG = 0;
const int Backprojector::BPROJ_TABLE = 1;
const int Backprojector::BPROJ_DIFF = 2;
-const int Backprojector::BPROJ_DIFF2 = 3;
-const int Backprojector::BPROJ_IDIFF2 = 4;
-const int Backprojector::BPROJ_IDIFF3 = 5;
+const int Backprojector::BPROJ_IDIFF = 3;
const char* Backprojector::s_aszBackprojectName[] =
{
{"trig"},
{"table"},
{"diff"},
- {"diff2"},
- {"idiff2"},
- {"idiff3"},
+ {"idiff"},
};
const char* Backprojector::s_aszBackprojectTitle[] =
{"Direct Trigometric"},
{"Trigometric Table"},
{"Difference Iteration"},
- {"Difference Iteration Optimized"},
- {"Integer Difference Iteration Optimized"},
- {"Integer Difference Iteration Highly-Optimized"},
+ {"Integer Difference Iteration"},
};
const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
{"nearest"},
{"linear"},
{"cubic"},
+#if HAVE_FREQ_PREINTERP
{"freq_preinterpolationj"},
+#endif
#if HAVE_BSPLINE_INTERP
{"bspline"},
{"1bspline"},
{"Nearest"},
{"Linear"},
{"Cubic"},
+#if HAVE_FREQ_PREINTERP
{"Frequency Preinterpolation"},
+#endif
#if HAVE_BSPLINE_INTERP
{"B-Spline"},
{"B-Spline 1st Order"},
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, interpFactor));
- else if (m_idBackproject == BPROJ_DIFF2)
- 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, interpFactor));
- else if (m_idBackproject == BPROJ_IDIFF3)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
+ else if (m_idBackproject == BPROJ_IDIFF)
+ m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor));
} else {
m_fail = true;
m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
{
double theta = view_angle;
- CubicInterpolator* pCubicInterp = NULL;
+ CubicPolyInterpolator* pCubicInterp = NULL;
if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
double x = xMin + xInc / 2; // Rectang coords of center of pixel
for (int ix = 0; ix < nx; x += xInc, ix++) {
{
double theta = view_angle;
- CubicInterpolator* pCubicInterp = NULL;
+ CubicPolyInterpolator* pCubicInterp = NULL;
if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
ScaleImageByRotIncrement();
}
-void
-BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
-{
- double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
- double det_dx = xInc * cos (theta);
- double det_dy = yInc * sin (theta);
- double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
-
- CubicInterpolator* pCubicInterp = NULL;
- if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
-
- for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
- double curDetPos = lColStart;
- ImageFileColumn pImCol = v[ix];
-
- for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-#ifdef DEBUG
- printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
-#endif
- if (interpType == Backprojector::INTERP_NEAREST) {
- int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
-
- if (iDetPos >= 0 && iDetPos < nDet)
- pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == Backprojector::INTERP_LINEAR) {
- double detPos = curDetPos / detInc; // position along detector
- double detPosFloor = floor (detPos);
- int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
- double frac = detPos - detPosFloor; // 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) {
- double p = iDetCenter + (curDetPos / detInc); // position along detector
- if (p >= 0 && p < nDet)
- pImCol[iy] += pCubicInterp->interpolate (p);
- }
- } // end for y
- } // end for x
-
- if (interpType == Backprojector::INTERP_CUBIC)
- delete pCubicInterp;
-}
-
-
-// CLASS IDENTICATION
-// BackprojectDiff2
-//
-// PURPOSE
-// Optimized version of BackprojectDiff
void
-BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
+BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
{
double theta = view_angle;
- CubicInterpolator* pCubicInterp = NULL;
- if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
-
- // Distance betw. detectors for an angle given in units of detectors
+ // Distance between detectors for an angle given in units of detectors
double det_dx = xInc * cos (theta) / detInc;
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;
-#ifdef DEBUG
- printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
-#endif
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
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) {
-#ifdef DEBUG
- printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
-#endif
if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
delete pCubicInterp;
}
-// CLASS IDENTICATION
-// BackprojectIntDiff2
-//
-// PURPOSE
-// Integer version of BackprojectDiff2
-
-void
-BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
-{
- double theta = view_angle;
-
- CubicInterpolator* pCubicInterp = NULL;
- if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
-
- static const kint32 scale = 1 << 16;
- static const double dScale = scale;
- static const kint32 halfScale = scale / 2;
-
- const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (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 * scale);
-
- for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
- kint32 curDetPos = detPosColStart;
- ImageFileColumn pImCol = v[ix];
-
- for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
- if (interpType == Backprojector::INTERP_NEAREST) {
- int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
- int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
-
- if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
- *pImCol++ += filteredProj[iDetPos];
- } else if (interpType == Backprojector::INTERP_LINEAR) {
- kint32 detPosFloor = curDetPos / scale;
- kint32 detPosRemainder = curDetPos % scale;
- if (detPosRemainder < 0) {
- detPosFloor--;
- detPosRemainder += scale;
- }
- int iDetPos = iDetCenter + detPosFloor;
- double frac = detPosRemainder / dScale;
- if (iDetPos < 0 || iDetPos >= nDet - 1)
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
- *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
- } else if (interpType = Backprojector::INTERP_CUBIC) {
- double p = iDetCenter + (static_cast<double>(curDetPos) / scale); // position along detector
- if (p >= 0 && p < nDet)
- *pImCol++ += pCubicInterp->interpolate (p);
- }
- } // end for y
- } // end for x
-
- if (interpType == Backprojector::INTERP_CUBIC)
- delete pCubicInterp;
-}
// CLASS IDENTICATION
-// BackprojectIntDiff3
+// BackprojectIntDiff
//
// PURPOSE
-// Highly optimized version of BackprojectIntDiff2
+// Highly optimized and integer version of BackprojectDiff
void
-BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
+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
static const int scaleShift = 16;
kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
double* deltaFilteredProj = NULL;
- CubicInterpolator* pCubicInterp = NULL;
+ CubicPolyInterpolator* pCubicInterp = NULL;
if (interpType == Backprojector::INTERP_LINEAR) {
// precalculate scaled difference for linear interpolation
deltaFilteredProj = new double [nDet];
deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
deltaFilteredProj[nDet - 1] = 0; // last detector
} else if (interpType == Backprojector::INTERP_CUBIC) {
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
}
int iLastDet = nDet - 1;
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
*pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
}
- } // end cubic
+ } // end Cubic
} // end for ix
if (interpType == Backprojector::INTERP_LINEAR)
{
double beta = view_angle;
- CubicInterpolator* pCubicInterp = NULL;
+ CubicPolyInterpolator* pCubicInterp = NULL;
if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
{
double beta = view_angle;
- CubicInterpolator* pCubicInterp = NULL;
+ CubicPolyInterpolator* pCubicInterp = NULL;
if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
projects / ctsim.git / blobdiff