** This is part of the CTSim program
** Copyright (c) 1983-2001 Kevin Rosenberg
**
-** $Id: backprojectors.cpp,v 1.34 2003/07/04 21:39:40 kevin Exp $
+** $Id$
**
** 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
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;
+ double detPosColStart = iDetCenter + start_r * cos (theta - start_phi) / detInc;
CubicPolyInterpolator* pCubicInterp = NULL;
- if (interpType == Backprojector::INTERP_CUBIC)
+ double* deltaFilteredProj = 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];
+ 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) {
double curDetPos = detPosColStart;
ImageFileColumn pImCol = v[ix];
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
if (interpType == Backprojector::INTERP_NEAREST) {
- int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
+ int iDetPos = 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);
+ int iDetPos = 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]));
+ if (iDetPos >= 0 && iDetPos <= iLastDet)
+ *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
} else if (interpType == Backprojector::INTERP_CUBIC) {
- double p = iDetCenter + curDetPos; // position along detector
+ double p = curDetPos; // position along detector
if (p >= 0 && p < nDet)
*pImCol++ += pCubicInterp->interpolate (p);
}
} // end for y
} // end for x
- if (interpType == Backprojector::INTERP_CUBIC)
+ if (interpType == Backprojector::INTERP_LINEAR)
+ delete deltaFilteredProj;
+ else if (interpType == Backprojector::INTERP_CUBIC)
delete pCubicInterp;
}
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
+#if SIZEOF_LONG == 4
static const int scaleShift = 16;
- static const kint32 scale = (1 << scaleShift);
- static const kint32 scaleBitmask = scale - 1;
- static const kint32 halfScale = scale / 2;
+#elif SIZEOF_LONG == 8
+ static const int scaleShift = 32;
+#endif
+ static const long scale = (1 << scaleShift);
+ static const long scaleBitmask = scale - 1;
+ static const long 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);
+ const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
+ const long det_dy = nearest<long> (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 + iDetCenter) * scale);
+ long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
double* deltaFilteredProj = NULL;
CubicPolyInterpolator* pCubicInterp = NULL;
int iLastDet = nDet - 1;
for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
- kint32 curDetPos = detPosColStart;
+ long 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;
+ const int iDetPos = (curDetPos + halfScale) >> scaleShift;
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;
+ const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * 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;
+ const long iDetPos = curDetPos >> scaleShift;
+ const long detRemainder = curDetPos & scaleBitmask;
if (iDetPos >= 0 && iDetPos <= iLastDet)
- *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+ *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);
+ *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
}
} // end Cubic
} // end for ix
projects / ctsim.git / blobdiff