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;
}
#elif SIZEOF_LONG == 8
static const int scaleShift = 32;
#endif
- static const long scale = (1 << scaleShift);
+ static const long scale = (1L << scaleShift);
static const long scaleBitmask = scale - 1;
static const long halfScale = scale / 2;
static const double dInvScale = 1. / scale;
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) {
} else if (interpType == Backprojector::INTERP_LINEAR) {
for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
const long iDetPos = curDetPos >> scaleShift;
- const long detRemainder = curDetPos & scaleBitmask;
- if (iDetPos >= 0 && iDetPos <= iLastDet)
- *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
+ if (iDetPos >= 0 && iDetPos <= iLastDet) {
+ const long detRemainder = curDetPos & scaleBitmask;
+ *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) {