- 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];
- } // 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);
- *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]);
- } // end for iy
- } //end linear
- } // end for 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
+
+ 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] / (dU * dU));
+ } 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]))
+ / (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
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;