** This is part of the CTSim program
** Copyright (C) 1983-2000 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.3 2000/06/22 10:17:28 kevin Exp $
+** $Id: backprojectors.cpp,v 1.6 2000/07/09 08:16:17 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
{
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 = convertInterpolationNameToID (interpName);
- m_pBackprojectImplem = NULL;
+ if (m_idInterpolation == INTERP_INVALID) {
+ m_fail = true;
+ m_failMessage = "Invalid interpolation name ";
+ m_failMessage += interpName;
+ }
- if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
+ if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
m_fail = true;
return false;
}
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 {
m_fail = true;
+ m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
return false;
}
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);
}
bprojName = BPROJ_DIFF2_STR;
else if (bprojID == BPROJ_IDIFF2)
bprojName = BPROJ_IDIFF2_STR;
+ else if (bprojID == BPROJ_IDIFF3)
+ bprojName = BPROJ_IDIFF3_STR;
return (bprojName);
}
interpID = INTERP_NEAREST;
else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
interpID = INTERP_LINEAR;
+ else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATE_STR) == 0)
+ interpID = INTERP_FREQ_PREINTERPOLATE;
#if HAVE_BSPLINE_INTERP
else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
interpID = INTERP_BSPLINE;
#endif
-
+
return (interpID);
}
return (INTERP_NEAREST_STR);
else if (interpID == INTERP_LINEAR)
return (INTERP_LINEAR_STR);
+ else if (interpID == INTERP_FREQ_PREINTERPOLATE)
+ return (INTERP_FREQ_PREINTERPOLATE_STR);
#if HAVE_BSPLINE_INTERP
else if (interpID == INTERP_BSPLINE)
return (INTERP_BSPLINE_STR);
{
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);
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 double dScale = scale;
+ 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);
+
+ 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) {
+ int iDetPos = (curDetPos + halfScale) >> 16;
+ 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) {
+ kint32 iDetPos = curDetPos >> scaleShift;
+ kint32 detRemainder = curDetPos & scaleBitmask;
+#if 0
+ double frac = detRemainder * (1. / 65536.);
+ assert(iDetPos >= 0 && iDetPos < nDet - 1);
+ *pImCol++ += ((1.-frac) * filteredProj[iDetPos]) + (frac * filteredProj[iDetPos+1]);
+#else
+ assert(iDetPos >= 0 && iDetPos < nDet - 1);
+ const double* const detPointer = &filteredProj[iDetPos];
+ *pImCol++ += (((scale-detRemainder) * *detPointer) + (detRemainder * *(detPointer+1))) * dInvScale;
+#endif
+ } // end for iy
+ } //end linear
+ } // end for ix
+}