- double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
- double det_dx = xInc * cos (theta);
- double det_dy = yInc * sin (theta);
- double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
-
- CubicInterpolator* pCubicInterp = NULL;
- if (interpType == Backprojector::INTERP_CUBIC)
- pCubicInterp = new CubicInterpolator (filteredProj, nDet);
-
- for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
- double curDetPos = lColStart;
- ImageFileColumn pImCol = v[ix];
-
- for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-#ifdef DEBUG
- printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
-#endif
- if (interpType == Backprojector::INTERP_NEAREST) {
- int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
-
- if (iDetPos >= 0 && iDetPos < nDet)
- pImCol[iy] += filteredProj[iDetPos];
- } else if (interpType == Backprojector::INTERP_LINEAR) {
- double detPos = curDetPos / detInc; // position along detector
- double detPosFloor = floor (detPos);
- int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
- double frac = detPos - detPosFloor; // fraction distance from det
- if (iDetPos >= 0 && iDetPos < nDet - 1)
- pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
- } else if (interpType = Backprojector::INTERP_CUBIC) {
- double p = iDetCenter + (curDetPos / detInc); // position along detector
- if (p >= 0 && p < nDet)
- pImCol[iy] += pCubicInterp->interpolate (p);
- }
- } // end for y
- } // end for x
-
- if (interpType == Backprojector::INTERP_CUBIC)
- delete pCubicInterp;