** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.9 2000/07/20 11:17:31 kevin Exp $
+** $Id: backprojectors.cpp,v 1.12 2000/08/25 15:59:13 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 char Backprojector::BPROJ_TRIG_STR[]= "trig";
-const char Backprojector::BPROJ_TABLE_STR[]= "table";
-const char Backprojector::BPROJ_DIFF_STR[]= "diff";
-const char Backprojector::BPROJ_DIFF2_STR[]= "diff2";
-const char Backprojector::BPROJ_IDIFF2_STR[]= "idiff2";
-const char Backprojector::BPROJ_IDIFF3_STR[]= "idiff3";
-
-const char Backprojector::BPROJ_TRIG_TITLE_STR[]= "Direc Trigometric";
-const char Backprojector::BPROJ_TABLE_TITLE_STR[]= "Trig Table";
-const char Backprojector::BPROJ_DIFF_TITLE_STR[]= "Diff";
-const char Backprojector::BPROJ_DIFF2_TITLE_STR[]= "Diff2";
-const char Backprojector::BPROJ_IDIFF2_TITLE_STR[]= "Integer Diff2";
-const char Backprojector::BPROJ_IDIFF3_TITLE_STR[]= "Integer Diff3";
-
-const char Backprojector::INTERP_NEAREST_STR[]= "nearest";
-const char Backprojector::INTERP_LINEAR_STR[]= "linear";
-const char Backprojector::INTERP_BSPLINE_STR[]= "bspline";
-const char Backprojector::INTERP_FREQ_PREINTERPOLATION_STR[]= "freq_preinterpolation";
-
-const char Backprojector::INTERP_NEAREST_TITLE_STR[]= "Nearest";
-const char Backprojector::INTERP_LINEAR_TITLE_STR[]= "Linear";
-const char Backprojector::INTERP_BSPLINE_TITLE_STR[]= "B-Spline";
-const char Backprojector::INTERP_FREQ_PREINTERPOLATION_TITLE_STR[]= "Frequency Preinterpolation";
+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_pBackprojectImplem->BackprojectView (viewData, viewAngle);
}
-Backprojector::~Backprojector (void)
+Backprojector::~Backprojector ()
{
delete m_pBackprojectImplem;
}
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 ";
}
-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;
+ 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)
+Backprojector::convertBackprojectIDToName (int 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;
+ 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);
+}
-const Backprojector::InterpolationID
-Backprojector::convertInterpolationNameToID (const char* const interpName)
+
+int
+Backprojector::convertInterpNameToID (const char* const interpName)
{
- InterpolationID interpID = INTERP_INVALID;
-
- if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
- interpID = INTERP_NEAREST;
- else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
- interpID = INTERP_LINEAR;
- 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;
-#endif
+ 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 = "";
-/* 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
- */
+ 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);
- else if (interpID == INTERP_FREQ_PREINTERPOLATION)
- return (INTERP_FREQ_PREINTERPOLATION_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, const int interpFactor)
+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();
yInc = (yMax - yMin) / ny;
}
-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++)
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);
+ 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);
+ 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());
}
// PURPOSE
// Precalculates trigometric function value for each point in image for backprojection.
-BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
+BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
: Backproject::Backproject (proj, im, interpType, interpFactor)
{
arrayR.initSetSize (nx, ny);
}
}
-BackprojectTable::~BackprojectTable (void)
+BackprojectTable::~BackprojectTable ()
{
ScaleImageByRotIncrement();
}
// 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, const int interpFactor)
+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]
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 == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
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 == Backprojector::INTERP_NEAREST) {
int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
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
+ if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos
errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
else
*pImCol++ += filteredProj[iDetPos];
kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
// precalculate scaled difference for linear interpolation
- double deltaFilteredProj [nDet - 1];
+ double deltaFilteredProj [nDet];
if (interpType == Backprojector::INTERP_LINEAR) {
for (int i = 0; i < nDet - 1; i++)
deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
}
+ deltaFilteredProj[nDet - 1] = 0; // last detector
+ int iLastDet = nDet - 1;
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];
+ 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;
- assert(iDetPos >= 0 && iDetPos < nDet);
+ 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;
- assert(iDetPos >= 0 && iDetPos < nDet - 1);
- *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+ if (iDetPos >= 0 && iDetPos <= iLastDet)
+ *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
} // end for iy
} //end linear
} // end for ix
projects / ctsim.git / blobdiff