+
+
+void
+BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[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 dFLPlusSin = m_dFocalLength + rsin_t;
+ double gamma = atan (rcos_t / dFLPlusSin);
+ double dPos = gamma / detInc; // position along detector
+ double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
+
+ 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] / dL2;
+ } 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])) / dL2;
+ } else if (interpType == Backprojector::INTERP_CUBIC) {
+ double d = iDetCenter + dPos; // position along detector
+ if (d >= 0 && d < nDet)
+ pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
+ }
+ } // end for y
+ } // end for x
+
+ if (interpType == Backprojector::INTERP_CUBIC)
+ delete pCubicInterp;
+}
+
+void
+BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
+{
+ double beta = view_angle;
+
+ CubicPolyInterpolator* pCubicInterp = NULL;
+ if (interpType == Backprojector::INTERP_CUBIC)
+ pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
+
+ for (int ix = 0; ix < nx; ix++) {
+ ImageFileColumn pImCol = v[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;
+}
+