** Date Started: June 2000
**
** This is part of the CTSim program
-** Copyright (C) 1983-2000 Kevin Rosenberg
+** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.5 2000/07/07 15:30:59 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
******************************************************************************/
#include "ct.h"
+#include "interpolator.h"
-Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
+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_IDIFF = 3;
+
+const char* const Backprojector::s_aszBackprojectName[] =
+{
+ {"trig"},
+ {"table"},
+ {"diff"},
+ {"idiff"},
+};
+
+const char* const Backprojector::s_aszBackprojectTitle[] =
+{
+ {"Direct Trigometric"},
+ {"Trigometric Table"},
+ {"Difference Iteration"},
+ {"Integer Difference Iteration"},
+};
+
+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_CUBIC = 2;
+const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
+#if HAVE_BSPLINE_INTERP
+const int Backprojector::INTERP_BSPLINE = 4;
+const int Backprojector::INTERP_1BSPLINE = 5;
+const int Backprojector::INTERP_2BSPLINE = 6;
+const int Backprojector::INTERP_3BSPLINE = 7;
+#endif
+
+const char* const Backprojector::s_aszInterpName[] =
+{
+ {"nearest"},
+ {"linear"},
+ {"cubic"},
+#if HAVE_FREQ_PREINTERP
+ {"freq_preinterpolationj"},
+#endif
+#if HAVE_BSPLINE_INTERP
+ {"bspline"},
+ {"1bspline"},
+ {"2bspline"},
+ {"3bspline"},
+#endif
+};
+
+const char* const Backprojector::s_aszInterpTitle[] =
+{
+ {"Nearest"},
+ {"Linear"},
+ {"Cubic"},
+#if HAVE_FREQ_PREINTERP
+ {"Frequency Preinterpolation"},
+#endif
+#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, const ReconstructionROI* pROI)
{
m_fail = false;
m_pBackprojectImplem = NULL;
-
- initBackprojector (proj, im, backprojName, interpName);
+
+ initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
}
void
Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
{
- if (m_pBackprojectImplem)
+ if (m_pBackprojectImplem != NULL)
m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
}
-Backprojector::~Backprojector (void)
+void
+Backprojector::PostProcessing()
+{
+ if (m_pBackprojectImplem != NULL)
+ m_pBackprojectImplem->PostProcessing();
+}
+
+Backprojector::~Backprojector ()
{
delete m_pBackprojectImplem;
}
// 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, const ReconstructionROI* pROI)
{
m_nameBackproject = backprojName;
m_nameInterpolation = interpName;
m_failMessage = "Invalid backprojection name ";
m_failMessage += backprojName;
}
- m_idInterpolation = convertInterpolationNameToID (interpName);
+ 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));
- else if (m_idBackproject == BPROJ_TABLE)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
- else if (m_idBackproject == BPROJ_DIFF)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
- else if (m_idBackproject == BPROJ_DIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
- else if (m_idBackproject == BPROJ_IDIFF2)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
- else if (m_idBackproject == BPROJ_IDIFF3)
- m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation));
- else {
+
+ if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
+ 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, 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, pROI));
+ else if (m_idBackproject == BPROJ_TABLE)
+ 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, pROI));
+ else if (m_idBackproject == BPROJ_IDIFF)
+ 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]";
return false;
}
-
+
return true;
}
-const Backprojector::BackprojectID
+int
Backprojector::convertBackprojectNameToID (const char* const backprojName)
{
- BackprojectID backprojID = BPROJ_INVALID;
-
- if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
- backprojID = BPROJ_TRIG;
- else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
- backprojID = BPROJ_TABLE;
- else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
- backprojID = BPROJ_DIFF;
- else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
- backprojID = BPROJ_DIFF2;
- else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
- backprojID = BPROJ_IDIFF2;
- else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0)
- backprojID = BPROJ_IDIFF3;
-
- return (backprojID);
+ 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 (const BackprojectID bprojID)
-{
- const char *bprojName = "";
-
- if (bprojID == BPROJ_TRIG)
- bprojName = BPROJ_TRIG_STR;
- else if (bprojID == BPROJ_TABLE)
- bprojName = BPROJ_TABLE_STR;
- else if (bprojID == BPROJ_DIFF)
- bprojName = BPROJ_DIFF_STR;
- else if (bprojID == BPROJ_DIFF2)
- bprojName = BPROJ_DIFF2_STR;
- else if (bprojID == BPROJ_IDIFF2)
- bprojName = BPROJ_IDIFF2_STR;
- else if (bprojID == BPROJ_IDIFF3)
- bprojName = BPROJ_IDIFF3_STR;
-
+Backprojector::convertBackprojectIDToName (int bprojID)
+{
+ static const char *bprojName = "";
+
+ if (bprojID >= 0 && bprojID < s_iBackprojectCount)
+ return (s_aszBackprojectName[bprojID]);
+
return (bprojName);
}
-
-
-const Backprojector::InterpolationID
-Backprojector::convertInterpolationNameToID (const char* const interpName)
+const char*
+Backprojector::convertBackprojectIDToTitle (const int bprojID)
{
- InterpolationID interpID = INTERP_INVALID;
+ static const char *bprojTitle = "";
+
+ if (bprojID >= 0 && bprojID < s_iBackprojectCount)
+ return (s_aszBackprojectTitle[bprojID]);
+
+ return (bprojTitle);
+}
- if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
- interpID = INTERP_NEAREST;
- else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
- interpID = INTERP_LINEAR;
-#if HAVE_BSPLINE_INTERP
- else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
- interpID = INTERP_BSPLINE;
-#endif
- return (interpID);
+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);
}
-
-/* NAME
- * name_of_interp Return name of interpolation method
- *
- * SYNOPSIS
- * name = name_of_interp (interp_type)
- * char *name Name of interpolation method
- * int interp_type Method of interpolation
- *
- * NOTES
- * Returns NULL if interp_type is invalid
- */
+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::convertInterpolationIDToName (const InterpolationID interpID)
+Backprojector::convertInterpIDToTitle (const int interpID)
{
- if (interpID == INTERP_NEAREST)
- return (INTERP_NEAREST_STR);
- else if (interpID == INTERP_LINEAR)
- return (INTERP_LINEAR_STR);
-#if HAVE_BSPLINE_INTERP
- else if (interpID == INTERP_BSPLINE)
- return (INTERP_BSPLINE_STR);
-#endif
- else
- return ("");
+ 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 Backprojector::InterpolationID interpType)
- : proj(proj), im(im), interpType(interpType)
+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();
nDet = proj.nDet();
iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
- rotInc = proj.rotInc();
-
+ rotScale = proj.rotInc();
+
+ if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
+ rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
+ else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
+ rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
+ else
+ sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
+
v = im.getArray();
nx = im.nx();
ny = im.ny();
im.arrayDataClear();
-
+
xMin = -proj.phmLen() / 2; // Retangular coords of phantom
xMax = xMin + proj.phmLen();
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);
- if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR)
- sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
+ m_dFocalLength = proj.focalLength();
+ m_dSourceDetectorLength = proj.sourceDetectorLength();
}
-Backproject::~Backproject (void)
+Backproject::~Backproject ()
{}
void
-Backproject::ScaleImageByRotIncrement (void)
+Backproject::PostProcessing()
+{
+ m_bPostProcessingDone = true;
+}
+
+void
+Backproject::ScaleImageByRotIncrement ()
{
for (int ix = 0; ix < nx; ix++)
for (int iy = 0; iy < ny; iy++)
- v[ix][iy] *= rotInc;
+ v[ix][iy] *= rotScale;
}
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
{
- printf ("r=%f, phi=%f\n", r, phi);
- errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
+ sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
+ errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
}
void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
{
- printf ("ix=%d, iy=%d\n", ix, iy);
- printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
- printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
- printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
- printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
- sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
+#if 1
+ std::ostringstream os;
+ os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
+ os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
+ os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
+ os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
+ os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
+
+ sys_error (ERR_WARNING, os.str().c_str());
+#endif
}
void
BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
{
- double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
- int ix, iy;
- double x, y; // Rectang coords of center of pixel
-
- for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
- for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+ double theta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ 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 y = yMin + yInc / 2;
+ for (int iy = 0; iy < ny; y += yInc, iy++) {
double r = sqrt (x * x + y * y); // distance of cell from center
double phi = atan2 (y, x); // angle of cell from center
double L = r * cos (theta - phi); // position on detector
-
+
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];
+ int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
+
+ if (iDetPos >= 0 && iDetPos < nDet)
+ v[ix][iy] += rotScale * filteredProj[iDetPos];
} else if (interpType == Backprojector::INTERP_LINEAR) {
- double p = L / detInc; // position along detector
- double pFloor = floor (p);
- int iDetPos = iDetCenter + static_cast<int>(pFloor);
- double frac = p - pFloor; // fraction distance from det
- if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
- errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
- else
- v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ double p = L / detInc; // position along detector
+ double pFloor = floor (p);
+ int iDetPos = iDetCenter + static_cast<int>(pFloor);
+ 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) {
+ double p = iDetCenter + (L / detInc); // position along detector
+ if (p >= 0 && p < nDet)
+ v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
}
}
+ }
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
// 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, int interpType,
+ const int interpFactor, const ReconstructionROI* pROI)
+: Backproject (proj, im, interpType, interpFactor, pROI)
{
- arrayR.initSetSize (nx, ny);
- arrayPhi.initSetSize (nx, ny);
+ arrayR.initSetSize (im.nx(), im.ny());
+ arrayPhi.initSetSize (im.nx(), im.ny());
r = arrayR.getArray();
phi = arrayPhi.getArray();
-
+
double x, y; // Rectang coords of center of pixel
int ix, iy;
for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
}
}
-BackprojectTable::~BackprojectTable (void)
+BackprojectTable::~BackprojectTable ()
{
- ScaleImageByRotIncrement();
}
void
-BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
+BackprojectTable::PostProcessing()
{
- double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
+ if (! m_bPostProcessingDone) {
+ ScaleImageByRotIncrement();
+ m_bPostProcessingDone = true;
+ }
+}
+void
+BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double theta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
-
+
for (int iy = 0; iy < ny; iy++) {
double L = r[ix][iy] * cos (theta - phi[ix][iy]);
-
+
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];
+ int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
+
+ if (iDetPos >= 0 && iDetPos < nDet)
+ pImCol[iy] += filteredProj[iDetPos];
} else if (interpType == Backprojector::INTERP_LINEAR) {
- double dPos = L / detInc; // position along detector
- double dPosFloor = floor (dPos);
- int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
- double frac = dPos - dPosFloor; // fraction distance from det
- if (iDetPos < 0 || iDetPos >= nDet - 1)
- errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
- else
- pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+ double dPos = L / detInc; // position along detector
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ 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) {
+ double p = iDetCenter + (L / 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;
}
// 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, 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 y = yMin + yInc / 2;
start_r = sqrt (x * x + y * y);
start_phi = atan2 (y, x);
-
+
im.arrayDataClear();
}
-BackprojectDiff::~BackprojectDiff()
+BackprojectDiff::~BackprojectDiff ()
{
- ScaleImageByRotIncrement();
}
void
-BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
+BackprojectDiff::PostProcessing()
{
- double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
- double det_dx = xInc * sin (theta);
- double det_dy = yInc * cos (theta);
- double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
-
- 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.5lf ", 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) // check for impossible: index outside of raysum pos
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
- 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)
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
- pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
- }
- } // end for y
- } // end for x
+ if (! m_bPostProcessingDone) {
+ ScaleImageByRotIncrement();
+ m_bPostProcessingDone = true;
+ }
}
-
-// 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; // add half PI to view angle to get perpendicular theta angle
-
- // Distance betw. detectors for an angle given in units of detectors
- double det_dx = xInc * sin (theta) / detInc;
- double det_dy = yInc * cos (theta) / detInc;
-
+ double theta = view_angle;
+
+ // 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, rotInc=%8.5f\n", start_r, start_phi, rotInc);
-#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>(L))]);
-#endif
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];
+ int iDetPos = iDetCenter + 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);
- 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]);
+ 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)
+ *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
+ } else if (interpType == Backprojector::INTERP_CUBIC) {
+ double p = iDetCenter + curDetPos; // 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
-// BackprojectIntDiff2
+// BackprojectIntDiff
//
// PURPOSE
-// Integer version of BackprojectDiff2
+// Highly optimized and integer version of BackprojectDiff
void
-BackprojectIntDiff2::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
-
- kint32 scale = 1 << 16;
- double dScale = scale;
- kint32 halfScale = scale / 2;
-
- kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
-
+ 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 * 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);
-
+ kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
+
+ double* deltaFilteredProj = NULL;
+ CubicPolyInterpolator* pCubicInterp = 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]) * dInvScale;
+ 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) {
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 impossible: 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]);
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+ const int iDetPos = (curDetPos + halfScale) >> 16;
+ 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;
+ 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]);
+ } // 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);
}
- } // end for y
- } // end for x
+ } // end Cubic
+ } // end for ix
+
+ if (interpType == Backprojector::INTERP_LINEAR)
+ delete deltaFilteredProj;
+ else if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
-// CLASS IDENTICATION
-// BackprojectIntDiff3
-//
-// PURPOSE
-// Highly optimized version of BackprojectIntDiff2
void
-BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
+BackprojectEquiangular::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 double dScale = scale;
- static const kint32 halfScale = scale / 2;
+ double beta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+ double dFLPlusSin = m_dFocalLength + rsin_t;
+ double gamma = atan (rcos_t / dFLPlusSin);
+ double dPos = gamma / detInc; // position along detector
+ double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
+
+ if (interpType == Backprojector::INTERP_NEAREST) {
+ int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
+ if (iDetPos >= 0 && iDetPos < nDet)
+ pImCol[iy] += filteredProj[iDetPos] / dL2;
+ } else if (interpType == Backprojector::INTERP_LINEAR) {
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos >= 0 && iDetPos < nDet - 1)
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
+ } else if (interpType == Backprojector::INTERP_CUBIC) {
+ double d = iDetCenter + dPos; // position along detector
+ if (d >= 0 && d < nDet)
+ pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
+ }
+ } // end for y
+ } // end for x
- const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
- const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
+}
- // calculate L for first point in image (0, 0)
- kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
-
- for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
- kint32 curDetPos = detPosColStart;
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
+ for (int ix = 0; ix < nx; ix++) {
ImageFileColumn pImCol = v[ix];
+
+ for (int iy = 0; iy < ny; iy++) {
+ double dAngleDiff = beta - phi[ix][iy];
+ double rcos_t = r[ix][iy] * cos (dAngleDiff);
+ double rsin_t = r[ix][iy] * sin (dAngleDiff);
+
+ double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
+ double dDetPos = rcos_t / dU;
+ // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
+ dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
+ double dPos = dDetPos / detInc; // position along detector array
- for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
if (interpType == Backprojector::INTERP_NEAREST) {
- int iDetPos = (curDetPos + halfScale) >> scaleShift;
-
- if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
- errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
- else
- *pImCol++ += filteredProj[iDetPos];
+ int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
+ if (iDetPos >= 0 && iDetPos < nDet)
+ pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
} 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]);
+ double dPosFloor = floor (dPos);
+ int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
+ double frac = dPos - dPosFloor; // fraction distance from det
+ if (iDetPos >= 0 && iDetPos < nDet - 1)
+ pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
+ / (dU * dU);
+ } else if (interpType == Backprojector::INTERP_CUBIC) {
+ double d = iDetCenter + dPos; // position along detector
+ if (d >= 0 && d < nDet)
+ pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
}
- } // end for y
- } // end for x
+ } // end for y
+ } // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
}
+
projects / ctsim.git / blobdiff