X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=7551ff9a1d1e99a1e65336b4e726af20c9f26a12;hp=8afe1892054fb7969fdc878b290b1e587e1b55ec;hb=4f15a69a3150f180bd93fcbe1c87dbaca92a77ad;hpb=4c687b64a6e11cc735525f9f65f962b53ba6b595 diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index 8afe189..7551ff9 100644 --- a/libctsim/backprojectors.cpp +++ b/libctsim/backprojectors.cpp @@ -8,7 +8,7 @@ ** This is part of the CTSim program ** Copyright (c) 1983-2001 Kevin Rosenberg ** -** $Id: backprojectors.cpp,v 1.25 2001/02/09 01:54:20 kevin Exp $ +** $Id: backprojectors.cpp,v 1.33 2003/03/23 18:37:42 kevin Exp $ ** ** 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 @@ -25,33 +25,28 @@ ******************************************************************************/ #include "ct.h" +#include "interpolator.h" const int Backprojector::BPROJ_INVALID = -1; const int Backprojector::BPROJ_TRIG = 0; const int Backprojector::BPROJ_TABLE = 1; const int Backprojector::BPROJ_DIFF = 2; -const int Backprojector::BPROJ_DIFF2 = 3; -const int Backprojector::BPROJ_IDIFF2 = 4; -const int Backprojector::BPROJ_IDIFF3 = 5; +const int Backprojector::BPROJ_IDIFF = 3; -const char* Backprojector::s_aszBackprojectName[] = +const char* const Backprojector::s_aszBackprojectName[] = { {"trig"}, {"table"}, {"diff"}, - {"diff2"}, - {"idiff2"}, - {"idiff3"}, + {"idiff"}, }; -const char* Backprojector::s_aszBackprojectTitle[] = +const char* const Backprojector::s_aszBackprojectTitle[] = { {"Direct Trigometric"}, {"Trigometric Table"}, {"Difference Iteration"}, - {"Difference Iteration Optimized"}, - {"Integer Difference Iteration Optimized"}, - {"Integer Difference Iteration Highly-Optimized"}, + {"Integer Difference Iteration"}, }; const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*); @@ -68,12 +63,14 @@ const int Backprojector::INTERP_2BSPLINE = 6; const int Backprojector::INTERP_3BSPLINE = 7; #endif -const char* Backprojector::s_aszInterpName[] = +const char* const Backprojector::s_aszInterpName[] = { {"nearest"}, {"linear"}, {"cubic"}, +#if HAVE_FREQ_PREINTERP {"freq_preinterpolationj"}, +#endif #if HAVE_BSPLINE_INTERP {"bspline"}, {"1bspline"}, @@ -82,12 +79,14 @@ const char* Backprojector::s_aszInterpName[] = #endif }; -const char* Backprojector::s_aszInterpTitle[] = +const char* const Backprojector::s_aszInterpTitle[] = { {"Nearest"}, {"Linear"}, {"Cubic"}, +#if HAVE_FREQ_PREINTERP {"Frequency Preinterpolation"}, +#endif #if HAVE_BSPLINE_INTERP {"B-Spline"}, {"B-Spline 1st Order"}, @@ -100,12 +99,13 @@ const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const -Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) +Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, + const char* const interpName, const int interpFactor, const ReconstructionROI* pROI) { m_fail = false; m_pBackprojectImplem = NULL; - initBackprojector (proj, im, backprojName, interpName, interpFactor); + initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI); } void @@ -115,6 +115,13 @@ Backprojector::BackprojectView (const double* const viewData, const double viewA m_pBackprojectImplem->BackprojectView (viewData, viewAngle); } +void +Backprojector::PostProcessing() +{ + if (m_pBackprojectImplem != NULL) + m_pBackprojectImplem->PostProcessing(); +} + Backprojector::~Backprojector () { delete m_pBackprojectImplem; @@ -128,7 +135,8 @@ Backprojector::~Backprojector () // and initializes the backprojector bool -Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor) +Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, + const char* const interpName, const int interpFactor, const ReconstructionROI* pROI) { m_nameBackproject = backprojName; m_nameInterpolation = interpName; @@ -152,22 +160,18 @@ Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const } if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR) - m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR) - m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI)); else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) { if (m_idBackproject == BPROJ_TRIG) - m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI)); else if (m_idBackproject == BPROJ_TABLE) - m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI)); else if (m_idBackproject == BPROJ_DIFF) - m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_DIFF2) - m_pBackprojectImplem = static_cast(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF2) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor)); - else if (m_idBackproject == BPROJ_IDIFF3) - m_pBackprojectImplem = static_cast(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor)); + m_pBackprojectImplem = static_cast(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI)); + else if (m_idBackproject == BPROJ_IDIFF) + m_pBackprojectImplem = static_cast(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI)); } else { m_fail = true; m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]"; @@ -259,8 +263,9 @@ Backprojector::convertInterpIDToTitle (const int interpID) // PURPOSE // Pure virtual base class for all backprojectors. -Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) -: proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor) +Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor, + const ReconstructionROI* pROI) +: proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false) { detInc = proj.detInc(); nDet = proj.nDet(); @@ -284,15 +289,47 @@ Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType yMin = -proj.phmLen() / 2; yMax = yMin + proj.phmLen(); + if (pROI) { + if (pROI->m_dXMin > xMin) + xMin = pROI->m_dXMin; + if (pROI->m_dXMax < xMax) + xMax = pROI->m_dXMax; + if (pROI->m_dYMin > yMin) + yMin = pROI->m_dYMin; + if (pROI->m_dYMax < yMax) + yMax = pROI->m_dYMax; + + if (xMin > xMax) { + double temp = xMin; + xMin = xMax; + xMax = temp; + } + if (yMin > yMax) { + double temp = yMin; + yMin = yMax; + yMax = temp; + } + } + xInc = (xMax - xMin) / nx; // size of cells yInc = (yMax - yMin) / ny; + im.setAxisIncrement (xInc, yInc); + im.setAxisExtent (xMin, xMax, yMin, yMax); + m_dFocalLength = proj.focalLength(); + m_dSourceDetectorLength = proj.sourceDetectorLength(); } Backproject::~Backproject () {} +void +Backproject::PostProcessing() +{ + m_bPostProcessingDone = true; +} + void Backproject::ScaleImageByRotIncrement () { @@ -333,9 +370,9 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double { double theta = view_angle; - CubicInterpolator* pCubicInterp = NULL; + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); double x = xMin + xInc / 2; // Rectang coords of center of pixel for (int ix = 0; ix < nx; x += xInc, ix++) { @@ -357,7 +394,7 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double double frac = p - pFloor; // fraction distance from det if (iDetPos >= 0 && iDetPos < nDet - 1) v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); - } else if (interpType = Backprojector::INTERP_CUBIC) { + } else if (interpType == Backprojector::INTERP_CUBIC) { double p = iDetCenter + (L / detInc); // position along detector if (p >= 0 && p < nDet) v[ix][iy] += rotScale * pCubicInterp->interpolate (p); @@ -376,8 +413,9 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double // PURPOSE // Precalculates trigometric function value for each point in image for backprojection. -BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) -: Backproject (proj, im, interpType, interpFactor) +BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, + const int interpFactor, const ReconstructionROI* pROI) +: Backproject (proj, im, interpType, interpFactor, pROI) { arrayR.initSetSize (im.nx(), im.ny()); arrayPhi.initSetSize (im.nx(), im.ny()); @@ -395,7 +433,15 @@ BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int BackprojectTable::~BackprojectTable () { - ScaleImageByRotIncrement(); +} + +void +BackprojectTable::PostProcessing() +{ + if (! m_bPostProcessingDone) { + ScaleImageByRotIncrement(); + m_bPostProcessingDone = true; + } } void @@ -403,9 +449,9 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl { double theta = view_angle; - CubicInterpolator* pCubicInterp = NULL; + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; @@ -425,7 +471,7 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl double frac = dPos - dPosFloor; // 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) { + } else if (interpType == Backprojector::INTERP_CUBIC) { double p = iDetCenter + (L / detInc); // position along detector if (p >= 0 && p < nDet) pImCol[iy] += pCubicInterp->interpolate (p); @@ -445,8 +491,9 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl // Backprojects by precalculating the change in L position for each x & y step in the image. // Iterates in x & y direction by adding difference in L position -BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor) -: Backproject (proj, im, interpType, interpFactor) +BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, + const int interpFactor, const ReconstructionROI* pROI) +: Backproject (proj, im, interpType, interpFactor, pROI) { // calculate center of first pixel v[0][0] double x = xMin + xInc / 2; @@ -457,89 +504,40 @@ BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int in im.arrayDataClear(); } -BackprojectDiff::~BackprojectDiff() +BackprojectDiff::~BackprojectDiff () { - ScaleImageByRotIncrement(); } void -BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) +BackprojectDiff::PostProcessing() { - 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(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(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; + if (! m_bPostProcessingDone) { + ScaleImageByRotIncrement(); + m_bPostProcessingDone = true; + } } - -// CLASS IDENTICATION -// BackprojectDiff2 -// -// PURPOSE -// Optimized version of BackprojectDiff - void -BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle) +BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle) { double theta = view_angle; - CubicInterpolator* pCubicInterp = NULL; - if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); - - // Distance betw. detectors for an angle given in units of detectors + // Distance between detectors for an angle given in units of detectors double det_dx = xInc * cos (theta) / detInc; 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; -#ifdef DEBUG - printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale); -#endif + CubicPolyInterpolator* pCubicInterp = NULL; + if (interpType == Backprojector::INTERP_CUBIC) + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); + 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) { -#ifdef DEBUG - printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest(curDetPos)]); -#endif if (interpType == Backprojector::INTERP_NEAREST) { int iDetPos = iDetCenter + nearest (curDetPos); // calc index in the filtered raysum vector @@ -551,7 +549,7 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl double frac = curDetPos - detPosFloor; // fraction distance from det if (iDetPos > 0 && iDetPos < nDet - 1) *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])); - } else if (interpType = Backprojector::INTERP_CUBIC) { + } else if (interpType == Backprojector::INTERP_CUBIC) { double p = iDetCenter + curDetPos; // position along detector if (p >= 0 && p < nDet) *pImCol++ += pCubicInterp->interpolate (p); @@ -563,77 +561,15 @@ BackprojectDiff2::BackprojectView (const double* const filteredProj, const doubl delete pCubicInterp; } -// CLASS IDENTICATION -// BackprojectIntDiff2 -// -// PURPOSE -// Integer version of BackprojectDiff2 - -void -BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle) -{ - double theta = view_angle; - - CubicInterpolator* pCubicInterp = NULL; - if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); - - static const kint32 scale = 1 << 16; - static const double dScale = scale; - static const kint32 halfScale = scale / 2; - - const kint32 det_dx = nearest (xInc * cos (theta) / detInc * scale); - const kint32 det_dy = nearest (yInc * sin (theta) / detInc * scale); - - // calculate L for first point in image (0, 0) - kint32 detPosColStart = nearest (start_r * cos (theta - start_phi) / detInc * scale); - - for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) { - kint32 curDetPos = detPosColStart; - ImageFileColumn pImCol = v[ix]; - - for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { - if (interpType == Backprojector::INTERP_NEAREST) { - int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale)); - int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector - - if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += filteredProj[iDetPos]; - } else if (interpType == Backprojector::INTERP_LINEAR) { - kint32 detPosFloor = curDetPos / scale; - kint32 detPosRemainder = curDetPos % scale; - if (detPosRemainder < 0) { - detPosFloor--; - detPosRemainder += scale; - } - int iDetPos = iDetCenter + detPosFloor; - double frac = detPosRemainder / dScale; - if (iDetPos < 0 || iDetPos >= nDet - 1) - errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos); - else - *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]); - } else if (interpType = Backprojector::INTERP_CUBIC) { - double p = iDetCenter + (static_cast(curDetPos) / scale); // position along detector - if (p >= 0 && p < nDet) - *pImCol++ += pCubicInterp->interpolate (p); - } - } // end for y - } // end for x - - if (interpType == Backprojector::INTERP_CUBIC) - delete pCubicInterp; -} // CLASS IDENTICATION -// BackprojectIntDiff3 +// BackprojectIntDiff // // PURPOSE -// Highly optimized version of BackprojectIntDiff2 +// Highly optimized and integer version of BackprojectDiff void -BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle) +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 static const int scaleShift = 16; @@ -649,7 +585,7 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do kint32 detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); double* deltaFilteredProj = NULL; - CubicInterpolator* pCubicInterp = NULL; + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_LINEAR) { // precalculate scaled difference for linear interpolation deltaFilteredProj = new double [nDet]; @@ -657,7 +593,7 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale; deltaFilteredProj[nDet - 1] = 0; // last detector } else if (interpType == Backprojector::INTERP_CUBIC) { - pCubicInterp = new CubicInterpolator (filteredProj, nDet); + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); } int iLastDet = nDet - 1; @@ -684,11 +620,11 @@ BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const do if (iDetPos >= 0 && iDetPos <= iLastDet) *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]); } // end for iy - } else if (interpType = Backprojector::INTERP_CUBIC) { + } else if (interpType == Backprojector::INTERP_CUBIC) { for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) { *pImCol++ += pCubicInterp->interpolate (static_cast(curDetPos) / 65536); } - } // end cubic + } // end Cubic } // end for ix if (interpType == Backprojector::INTERP_LINEAR) @@ -703,9 +639,9 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const { double beta = view_angle; - CubicInterpolator* pCubicInterp = NULL; + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; @@ -746,9 +682,9 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const { double beta = view_angle; - CubicInterpolator* pCubicInterp = NULL; + CubicPolyInterpolator* pCubicInterp = NULL; if (interpType == Backprojector::INTERP_CUBIC) - pCubicInterp = new CubicInterpolator (filteredProj, nDet); + pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet); for (int ix = 0; ix < nx; ix++) { ImageFileColumn pImCol = v[ix]; @@ -760,10 +696,8 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double dU = (m_dFocalLength + rsin_t) / m_dFocalLength; double dDetPos = rcos_t / dU; - // double to scale for imaginary detector that passes through origin - // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also - // located focal-length away from the origin. - dDetPos *= 2; + // 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) {