** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.2 2000/06/19 19:07:33 kevin Exp $
+** $Id: backprojectors.cpp,v 1.11 2000/07/23 01:49:03 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"
+const int Backprojector::BPROJ_INVALID = -1;
+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 char* Backprojector::s_aszBackprojectName[] =
+{
+ {"trig"},
+ {"table"},
+ {"diff"},
+ {"diff2"},
+ {"idiff2"},
+ {"idiff3"},
+};
+
+const char* Backprojector::s_aszBackprojectTitle[] =
+{
+ {"Direct Trigometric"},
+ {"Trigometric Table"},
+ {"Difference Iteration"},
+ {"Difference Iteration Optimized"},
+ {"Integer Difference Iteration Optimized"},
+ {"Integer Difference Iteration Highly-Optimized"},
+};
+
+const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
+
+const int Backprojector::INTERP_INVALID = -1;
+const int Backprojector::INTERP_NEAREST = 0;
+const int Backprojector::INTERP_LINEAR = 1;
+const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 2;
+#if HAVE_BSPLINE_INTERP
+const int Backprojector::INTERP_BSPLINE = 3;
+const int Backprojector::INTERP_1BSPLINE = 4;
+const int Backprojector::INTERP_2BSPLINE = 5;
+const int Backprojector::INTERP_3BSPLINE = 6;
+#endif
+
+const char* Backprojector::s_aszInterpName[] =
+{
+ {"nearest"},
+ {"linear"},
+ {"freq_preinterpolationj"},
+#if HAVE_BSPLINE_INTERP
+ {"bspline"},
+ {"1bspline"},
+ {"2bspline"},
+ {"3bspline"},
+#endif
+};
+
+const char* Backprojector::s_aszInterpTitle[] =
+{
+ {"Nearest"},
+ {"Linear"},
+ {"Frequency Preinterpolationj"},
+#if HAVE_BSPLINE_INTERP
+ {"B-Spline"},
+ {"B-Spline 1st Order"},
+ {"B-Spline 2nd Order"},
+ {"B-Spline 3rd Order"},
+#endif
+};
+
+const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
+
+
+
+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, interpFactor);
+}
+
+void
+Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
+{
+ if (m_pBackprojectImplem != NULL)
+ m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
+}
+
+Backprojector::~Backprojector ()
+{
+ delete m_pBackprojectImplem;
+}
// FUNCTION IDENTIFICATION
// Backproject* projector = selectBackprojector (...)
// Selects a backprojector based on BackprojType
// and initializes the backprojector
-Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
+bool
+Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
+{
+ m_nameBackproject = backprojName;
+ m_nameInterpolation = interpName;
+ m_pBackprojectImplem = NULL;
+ m_idBackproject = convertBackprojectNameToID (backprojName);
+ if (m_idBackproject == BPROJ_INVALID) {
+ m_fail = true;
+ m_failMessage = "Invalid backprojection name ";
+ m_failMessage += backprojName;
+ }
+ m_idInterpolation = convertInterpNameToID (interpName);
+ if (m_idInterpolation == INTERP_INVALID) {
+ m_fail = true;
+ m_failMessage = "Invalid interpolation name ";
+ m_failMessage += interpName;
+ }
+
+ if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
+ m_fail = true;
+ return false;
+ }
+
+ if (m_idBackproject == BPROJ_TRIG)
+ 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, 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 {
+ m_fail = true;
+ m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
+ return false;
+ }
+
+ return true;
+}
+
+
+int
+Backprojector::convertBackprojectNameToID (const char* const backprojName)
{
- Backproject* bj = NULL;
-
- if (bjType == O_BPROJ_TRIG)
- bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
- else if (bjType == O_BPROJ_TABLE)
- bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
- else if (bjType == O_BPROJ_DIFF)
- bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
- else if (bjType == O_BPROJ_DIFF2)
- bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
- else if (bjType == O_BPROJ_IDIFF2)
- bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
- else
- sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
-
- return (bj);
+ int backprojID = BPROJ_INVALID;
+
+ for (int i = 0; i < s_iBackprojectCount; i++)
+ if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
+ backprojID = i;
+ break;
+ }
+
+ return (backprojID);
}
+const char*
+Backprojector::convertBackprojectIDToName (int bprojID)
+{
+ static const char *bprojName = "";
+
+ if (bprojID >= 0 && bprojID < s_iBackprojectCount)
+ return (s_aszBackprojectName[bprojID]);
+
+ return (bprojName);
+}
+
+const char*
+Backprojector::convertBackprojectIDToTitle (const int bprojID)
+{
+ static const char *bprojTitle = "";
+
+ if (bprojID >= 0 && bprojID < s_iBackprojectCount)
+ return (s_aszBackprojectTitle[bprojID]);
+
+ return (bprojTitle);
+}
+
+
+int
+Backprojector::convertInterpNameToID (const char* const interpName)
+{
+ int interpID = INTERP_INVALID;
+
+ for (int i = 0; i < s_iInterpCount; i++)
+ if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
+ interpID = i;
+ break;
+ }
+
+ return (interpID);
+}
+
+const char*
+Backprojector::convertInterpIDToName (const int interpID)
+{
+ static const char *interpName = "";
+
+ if (interpID >= 0 && interpID < s_iInterpCount)
+ return (s_aszInterpName[interpID]);
+
+ return (interpName);
+}
+
+const char*
+Backprojector::convertInterpIDToTitle (const int interpID)
+{
+ static const char *interpTitle = "";
+
+ if (interpID >= 0 && interpID < s_iInterpCount)
+ return (s_aszInterpTitle[interpID]);
+
+ return (interpTitle);
+}
+
+
// CLASS IDENTICATION
// Backproject
// PURPOSE
// Pure virtual base class for all backprojectors.
-Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType)
- : proj(proj), im(im), interpType(interpType)
+Backproject::Backproject (const Projections& proj, ImageFile& im, const int 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 != I_NEAREST && interpType != I_LINEAR)
- sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
}
-Backproject::~Backproject (void)
+Backproject::~Backproject ()
{}
void
-Backproject::ScaleImageByRotIncrement (void)
+Backproject::ScaleImageByRotIncrement ()
{
for (int ix = 0; ix < nx; ix++)
for (int iy = 0; iy < ny; iy++)
double phi = atan2 (y, x); // angle of cell from center
double L = r * cos (theta - phi); // position on detector
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
else
v[ix][iy] += rotInc * filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double p = L / detInc; // position along detector
double pFloor = floor (p);
int iDetPos = iDetCenter + static_cast<int>(pFloor);
// PURPOSE
// Precalculates trigometric function value for each point in image for backprojection.
-BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType)
- : Backproject::Backproject (proj, im, interpType)
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
+ : Backproject::Backproject (proj, im, interpType, interpFactor)
{
arrayR.initSetSize (nx, ny);
arrayPhi.initSetSize (nx, ny);
}
}
-BackprojectTable::~BackprojectTable (void)
+BackprojectTable::~BackprojectTable ()
{
ScaleImageByRotIncrement();
}
for (int iy = 0; iy < ny; iy++) {
double L = r[ix][iy] * cos (theta - phi[ix][iy]);
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
else
pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double dPos = L / detInc; // position along detector
double dPosFloor = floor (dPos);
int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
// 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, InterpolationType interpType)
- : Backproject::Backproject (proj, im, interpType)
+BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
+ : Backproject::Backproject (proj, im, interpType, interpFactor)
{
// calculate center of first pixel v[0][0]
double x = xMin + xInc / 2;
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
#ifdef DEBUG
- printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
+ printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
#endif
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double detPos = curDetPos / detInc; // position along detector
double detPosFloor = floor (detPos);
int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
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>(L))]);
+ printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
#endif
- if (interpType == I_NEAREST) {
+ if (interpType == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
double detPosFloor = floor (curDetPos);
int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
double frac = curDetPos - detPosFloor; // fraction distance from det
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
{
double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
- kint32 scale = 1 << 16;
- double dScale = scale;
- kint32 halfScale = scale / 2;
+ static const kint32 scale = 1 << 16;
+ static const double dScale = scale;
+ static const kint32 halfScale = scale / 2;
- kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
// calculate L for first point in image (0, 0)
kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
ImageFileColumn pImCol = v[ix];
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
- if (interpType == I_NEAREST) {
+ 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
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += filteredProj[iDetPos];
- } else if (interpType == I_LINEAR) {
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
kint32 detPosFloor = curDetPos / scale;
kint32 detPosRemainder = curDetPos % scale;
if (detPosRemainder < 0) {
if (iDetPos < 0 || iDetPos >= nDet - 1)
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
- *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
}
} // end for y
} // end for x
}
+
+// CLASS IDENTICATION
+// BackprojectIntDiff3
+//
+// PURPOSE
+// Highly optimized version of BackprojectIntDiff2
+
+void
+BackprojectIntDiff3::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;
+ static const kint32 scale = (1 << scaleShift);
+ static const kint32 scaleBitmask = scale - 1;
+ static const kint32 halfScale = scale / 2;
+ static const double dInvScale = 1. / scale;
+
+ const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
+ const kint32 det_dy = nearest<kint32> (yInc * cos (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);
+
+ // 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) {
+ 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) {
+ const kint32 iDetPos = curDetPos >> scaleShift;
+ const kint32 detRemainder = curDetPos & scaleBitmask;
+ assert(iDetPos >= 0 && iDetPos < nDet - 1);
+ *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+ } // end for iy
+ } //end linear
+ } // end for ix
+}