1 /*****************************************************************************
4 ** Name: backprojectors.cpp Classes for backprojection
5 ** Programmer: Kevin Rosenberg
6 ** Date Started: June 2000
8 ** This is part of the CTSim program
9 ** Copyright (c) 1983-2009 Kevin Rosenberg
11 ** This program is free software; you can redistribute it and/or modify
12 ** it under the terms of the GNU General Public License (version 2) as
13 ** published by the Free Software Foundation.
15 ** This program is distributed in the hope that it will be useful,
16 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
17 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 ** GNU General Public License for more details.
20 ** You should have received a copy of the GNU General Public License
21 ** along with this program; if not, write to the Free Software
22 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
23 ******************************************************************************/
26 #include "interpolator.h"
28 const int Backprojector::BPROJ_INVALID = -1;
29 const int Backprojector::BPROJ_TRIG = 0;
30 const int Backprojector::BPROJ_TABLE = 1;
31 const int Backprojector::BPROJ_DIFF = 2;
32 const int Backprojector::BPROJ_IDIFF = 3;
34 const char* const Backprojector::s_aszBackprojectName[] =
42 const char* const Backprojector::s_aszBackprojectTitle[] =
46 "Difference Iteration",
47 "Integer Difference Iteration",
50 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
52 const int Backprojector::INTERP_INVALID = -1;
53 const int Backprojector::INTERP_NEAREST = 0;
54 const int Backprojector::INTERP_LINEAR = 1;
55 const int Backprojector::INTERP_CUBIC = 2;
56 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
57 #if HAVE_BSPLINE_INTERP
58 const int Backprojector::INTERP_BSPLINE = 4;
59 const int Backprojector::INTERP_1BSPLINE = 5;
60 const int Backprojector::INTERP_2BSPLINE = 6;
61 const int Backprojector::INTERP_3BSPLINE = 7;
64 const char* const Backprojector::s_aszInterpName[] =
69 #if HAVE_FREQ_PREINTERP
70 "freq_preinterpolationj",
72 #if HAVE_BSPLINE_INTERP
80 const char* const Backprojector::s_aszInterpTitle[] =
85 #if HAVE_FREQ_PREINTERP
86 "Frequency Preinterpolation",
88 #if HAVE_BSPLINE_INTERP
96 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
100 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
101 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
104 m_pBackprojectImplem = NULL;
106 initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
110 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
112 if (m_pBackprojectImplem != NULL)
113 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
117 Backprojector::PostProcessing()
119 if (m_pBackprojectImplem != NULL)
120 m_pBackprojectImplem->PostProcessing();
123 Backprojector::~Backprojector ()
125 delete m_pBackprojectImplem;
128 // FUNCTION IDENTIFICATION
129 // Backproject* projector = selectBackprojector (...)
132 // Selects a backprojector based on BackprojType
133 // and initializes the backprojector
136 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
137 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
139 m_nameBackproject = backprojName;
140 m_nameInterpolation = interpName;
141 m_pBackprojectImplem = NULL;
142 m_idBackproject = convertBackprojectNameToID (backprojName);
143 if (m_idBackproject == BPROJ_INVALID) {
145 m_failMessage = "Invalid backprojection name ";
146 m_failMessage += backprojName;
148 m_idInterpolation = convertInterpNameToID (interpName);
149 if (m_idInterpolation == INTERP_INVALID) {
151 m_failMessage = "Invalid interpolation name ";
152 m_failMessage += interpName;
155 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
160 if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
161 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI));
162 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
163 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI));
164 else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
165 if (m_idBackproject == BPROJ_TRIG)
166 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI));
167 else if (m_idBackproject == BPROJ_TABLE)
168 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI));
169 else if (m_idBackproject == BPROJ_DIFF)
170 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI));
171 else if (m_idBackproject == BPROJ_IDIFF)
172 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI));
175 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
184 Backprojector::convertBackprojectNameToID (const char* const backprojName)
186 int backprojID = BPROJ_INVALID;
188 for (int i = 0; i < s_iBackprojectCount; i++) {
189 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
198 Backprojector::convertBackprojectIDToName (int bprojID)
200 static const char *bprojName = "";
202 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
203 return (s_aszBackprojectName[bprojID]);
209 Backprojector::convertBackprojectIDToTitle (const int bprojID)
211 static const char *bprojTitle = "";
213 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
214 return (s_aszBackprojectTitle[bprojID]);
221 Backprojector::convertInterpNameToID (const char* const interpName)
223 int interpID = INTERP_INVALID;
225 for (int i = 0; i < s_iInterpCount; i++) {
226 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
235 Backprojector::convertInterpIDToName (const int interpID)
237 static const char *interpName = "";
239 if (interpID >= 0 && interpID < s_iInterpCount)
240 return (s_aszInterpName[interpID]);
246 Backprojector::convertInterpIDToTitle (const int interpID)
248 static const char *interpTitle = "";
250 if (interpID >= 0 && interpID < s_iInterpCount)
251 return (s_aszInterpTitle[interpID]);
253 return (interpTitle);
258 // CLASS IDENTICATION
262 // Pure virtual base class for all backprojectors.
264 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor,
265 const ReconstructionROI* pROI)
266 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false)
268 detInc = proj.detInc();
270 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
272 if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
273 rotScale = PI / proj.nView(); // scale by number of PI rotations
274 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
275 rotScale = (2 * PI) / proj.nView(); // scale by number of 2PI rotations
277 sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
284 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
285 xMax = xMin + proj.phmLen();
286 yMin = -proj.phmLen() / 2;
287 yMax = yMin + proj.phmLen();
290 if (pROI->m_dXMin > xMin)
291 xMin = pROI->m_dXMin;
292 if (pROI->m_dXMax < xMax)
293 xMax = pROI->m_dXMax;
294 if (pROI->m_dYMin > yMin)
295 yMin = pROI->m_dYMin;
296 if (pROI->m_dYMax < yMax)
297 yMax = pROI->m_dYMax;
311 xInc = (xMax - xMin) / nx; // size of cells
312 yInc = (yMax - yMin) / ny;
314 im.setAxisIncrement (xInc, yInc);
315 im.setAxisExtent (xMin, xMax, yMin, yMax);
317 m_dFocalLength = proj.focalLength();
318 m_dSourceDetectorLength = proj.sourceDetectorLength();
321 Backproject::~Backproject ()
325 Backproject::PostProcessing()
327 m_bPostProcessingDone = true;
331 Backproject::ScaleImageByRotIncrement ()
333 for (int ix = 0; ix < nx; ix++)
334 for (int iy = 0; iy < ny; iy++)
335 v[ix][iy] *= rotScale;
338 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
340 sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
341 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
344 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
347 std::ostringstream os;
348 os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
349 os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
350 os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
351 os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
352 os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
354 sys_error (ERR_WARNING, os.str().c_str());
359 // CLASS IDENTICATION
363 // Uses trigometric functions at each point in image for backprojection.
366 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
368 double theta = view_angle;
370 CubicPolyInterpolator* pCubicInterp = NULL;
371 if (interpType == Backprojector::INTERP_CUBIC)
372 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
374 double xstart = xMin + xInc / 2; // Rectang coords of center of pixel
376 #pragma omp parallel for
378 for (int ix = 0; ix < nx; ix++) {
379 double x = xstart + (ix * xInc);
381 double y = yMin + yInc / 2;
382 for (int iy = 0; iy < ny; y += yInc, iy++) {
383 double r = sqrt (x * x + y * y); // distance of cell from center
384 double phi = atan2 (y, x); // angle of cell from center
385 double L = r * cos (theta - phi); // position on detector
387 if (interpType == Backprojector::INTERP_NEAREST) {
388 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
390 if (iDetPos >= 0 && iDetPos < nDet)
391 v[ix][iy] += rotScale * filteredProj[iDetPos];
392 } else if (interpType == Backprojector::INTERP_LINEAR) {
393 double p = L / detInc; // position along detector
394 double pFloor = floor (p);
395 int iDetPos = iDetCenter + static_cast<int>(pFloor);
396 double frac = p - pFloor; // fraction distance from det
397 if (iDetPos >= 0 && iDetPos < nDet - 1)
398 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
399 } else if (interpType == Backprojector::INTERP_CUBIC) {
400 double p = iDetCenter + (L / detInc); // position along detector
401 if (p >= 0 && p < nDet)
402 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
407 if (interpType == Backprojector::INTERP_CUBIC)
412 // CLASS IDENTICATION
416 // Precalculates trigometric function value for each point in image for backprojection.
418 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
419 const int interpFactor, const ReconstructionROI* pROI)
420 : Backproject (proj, im, interpType, interpFactor, pROI)
422 arrayR.initSetSize (im.nx(), im.ny());
423 arrayPhi.initSetSize (im.nx(), im.ny());
424 r = arrayR.getArray();
425 phi = arrayPhi.getArray();
427 double xstart = xMin + xInc / 2;
430 #pragma omp parallel for
432 for (int ix = 0; ix < nx; ix++) {
433 double x = xstart + (ix * xInc);
434 double y = yMin + yInc / 2;
435 for (int iy = 0; iy < ny; iy++, y += yInc) {
436 r[ix][iy] = sqrt (x * x + y * y);
437 phi[ix][iy] = atan2 (y, x);
442 BackprojectTable::~BackprojectTable ()
447 BackprojectTable::PostProcessing()
449 if (! m_bPostProcessingDone) {
450 ScaleImageByRotIncrement();
451 m_bPostProcessingDone = true;
456 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
458 double theta = view_angle;
460 CubicPolyInterpolator* pCubicInterp = NULL;
461 if (interpType == Backprojector::INTERP_CUBIC)
462 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
465 #pragma omp parallel for
467 for (int ix = 0; ix < nx; ix++) {
468 ImageFileColumn pImCol = v[ix];
470 for (int iy = 0; iy < ny; iy++) {
471 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
473 if (interpType == Backprojector::INTERP_NEAREST) {
474 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
476 if (iDetPos >= 0 && iDetPos < nDet) {
477 pImCol[iy] += filteredProj[iDetPos];
479 } else if (interpType == Backprojector::INTERP_LINEAR) {
480 double dPos = L / detInc; // position along detector
481 double dPosFloor = floor (dPos);
482 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
483 double frac = dPos - dPosFloor; // fraction distance from det
484 if (iDetPos >= 0 && iDetPos < nDet - 1) {
485 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
487 } else if (interpType == Backprojector::INTERP_CUBIC) {
488 double p = iDetCenter + (L / detInc); // position along detector
489 if (p >= 0 && p < nDet) {
490 pImCol[iy] += pCubicInterp->interpolate (p);
496 if (interpType == Backprojector::INTERP_CUBIC)
501 // CLASS IDENTICATION
505 // Backprojects by precalculating the change in L position for each x & y step in the image.
506 // Iterates in x & y direction by adding difference in L position
508 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
509 const int interpFactor, const ReconstructionROI* pROI)
510 : Backproject (proj, im, interpType, interpFactor, pROI)
512 // calculate center of first pixel v[0][0]
513 double x = xMin + xInc / 2;
514 double y = yMin + yInc / 2;
515 start_r = sqrt (x * x + y * y);
516 start_phi = atan2 (y, x);
521 BackprojectDiff::~BackprojectDiff ()
526 BackprojectDiff::PostProcessing()
528 if (! m_bPostProcessingDone) {
529 ScaleImageByRotIncrement();
530 m_bPostProcessingDone = true;
535 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
537 double theta = view_angle;
539 // Distance between detectors for an angle given in units of detectors
540 double det_dx = xInc * cos (theta) / detInc;
541 double det_dy = yInc * sin (theta) / detInc;
543 // calculate detPosition for first point in image (ix=0, iy=0)
544 double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc;
546 CubicPolyInterpolator* pCubicInterp = NULL;
547 double* deltaFilteredProj = NULL;
548 if (interpType == Backprojector::INTERP_LINEAR) {
549 // precalculate scaled difference for linear interpolation
550 deltaFilteredProj = new double [nDet];
551 for (int i = 0; i < nDet - 1; i++)
552 deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
553 deltaFilteredProj[nDet - 1] = 0; // last detector
554 } else if (interpType == Backprojector::INTERP_CUBIC) {
555 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
558 int iLastDet = nDet - 1;
561 #pragma omp parallel for
563 for (int ix = 0; ix < nx; ix++) {
564 double detPos = detPosColBase + (ix * det_dx);
565 ImageFileColumn pImCol = v[ix];
567 for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
568 if (interpType == Backprojector::INTERP_NEAREST) {
569 int iDetPos = nearest<int> (detPos); // calc index in the filtered raysum vector
570 if (iDetPos >= 0 && iDetPos < nDet) {
571 *pImCol++ += filteredProj[iDetPos];
573 } else if (interpType == Backprojector::INTERP_LINEAR) {
574 double detPosFloor = floor (detPos);
575 int iDetPos = static_cast<int>(detPosFloor);
576 double frac = detPos - detPosFloor; // fraction distance from det
577 if (iDetPos >= 0 && iDetPos <= iLastDet) {
578 *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
580 } else if (interpType == Backprojector::INTERP_CUBIC) {
581 double p = detPos; // position along detector
582 if (p >= 0 && p < nDet) {
583 *pImCol++ += pCubicInterp->interpolate (p);
589 if (interpType == Backprojector::INTERP_LINEAR)
590 delete deltaFilteredProj;
591 else if (interpType == Backprojector::INTERP_CUBIC)
596 // CLASS IDENTICATION
597 // BackprojectIntDiff
600 // Highly optimized and integer version of BackprojectDiff
603 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
605 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
607 static const int scaleShift = 16;
608 #elif SIZEOF_LONG == 8
609 static const int scaleShift = 32;
611 static const long scale = (1L << scaleShift);
612 static const long scaleBitmask = scale - 1;
613 static const long halfScale = scale / 2;
614 static const double dInvScale = 1. / scale;
616 const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
617 const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
619 // calculate L for first point in image (0, 0)
620 long detPosColBase = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
622 double* deltaFilteredProj = NULL;
623 CubicPolyInterpolator* pCubicInterp = NULL;
624 if (interpType == Backprojector::INTERP_LINEAR) {
625 // precalculate scaled difference for linear interpolation
626 deltaFilteredProj = new double [nDet];
627 for (int i = 0; i < nDet - 1; i++)
628 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
629 deltaFilteredProj[nDet - 1] = 0; // last detector
630 } else if (interpType == Backprojector::INTERP_CUBIC) {
631 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
634 int iLastDet = nDet - 1;
636 #pragma omp parallel for
638 for (int ix = 0; ix < nx; ix++) {
639 long detPos = detPosColBase + (ix * det_dx);
640 ImageFileColumn pImCol = v[ix];
642 for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
643 if (interpType == Backprojector::INTERP_NEAREST) {
644 const int iDetPos = (detPos + halfScale) >> scaleShift;
645 if (iDetPos >= 0 && iDetPos <= iLastDet) {
646 *pImCol++ += filteredProj[iDetPos];
649 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
650 const int iDetPos = ((detPos + halfScale) >> scaleShift) * m_interpFactor;
651 if (iDetPos >= 0 && iDetPos <= iLastDet) {
652 *pImCol++ += filteredProj[iDetPos];
655 } else if (interpType == Backprojector::INTERP_LINEAR) {
656 const long iDetPos = detPos >> scaleShift;
657 if (iDetPos >= 0 && iDetPos <= iLastDet) {
658 const long detRemainder = detPos & scaleBitmask;
659 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
662 } else if (interpType == Backprojector::INTERP_CUBIC) {
663 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(detPos) / scale);
668 if (interpType == Backprojector::INTERP_LINEAR)
669 delete deltaFilteredProj;
670 else if (interpType == Backprojector::INTERP_CUBIC)
676 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
678 double beta = view_angle;
680 CubicPolyInterpolator* pCubicInterp = NULL;
681 if (interpType == Backprojector::INTERP_CUBIC)
682 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
685 #pragma omp parallel for
687 for (int ix = 0; ix < nx; ix++) {
688 ImageFileColumn pImCol = v[ix];
690 for (int iy = 0; iy < ny; iy++) {
691 double dAngleDiff = beta - phi[ix][iy];
692 double rcos_t = r[ix][iy] * cos (dAngleDiff);
693 double rsin_t = r[ix][iy] * sin (dAngleDiff);
694 double dFLPlusSin = m_dFocalLength + rsin_t;
695 double gamma = atan (rcos_t / dFLPlusSin);
696 double dPos = gamma / detInc; // position along detector
697 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
699 if (interpType == Backprojector::INTERP_NEAREST) {
700 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
701 if (iDetPos >= 0 && iDetPos < nDet)
702 pImCol[iy] += filteredProj[iDetPos] / dL2;
703 } else if (interpType == Backprojector::INTERP_LINEAR) {
704 double dPosFloor = floor (dPos);
705 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
706 double frac = dPos - dPosFloor; // fraction distance from det
707 if (iDetPos >= 0 && iDetPos < nDet - 1)
708 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
709 } else if (interpType == Backprojector::INTERP_CUBIC) {
710 double d = iDetCenter + dPos; // position along detector
711 if (d >= 0 && d < nDet)
712 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
717 if (interpType == Backprojector::INTERP_CUBIC)
722 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
724 double beta = view_angle;
726 CubicPolyInterpolator* pCubicInterp = NULL;
727 if (interpType == Backprojector::INTERP_CUBIC)
728 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
731 #pragma omp parallel for
733 for (int ix = 0; ix < nx; ix++) {
734 ImageFileColumn pImCol = v[ix];
736 for (int iy = 0; iy < ny; iy++) {
737 double dAngleDiff = beta - phi[ix][iy];
738 double rcos_t = r[ix][iy] * cos (dAngleDiff);
739 double rsin_t = r[ix][iy] * sin (dAngleDiff);
741 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
742 double dU2 = dU * dU;
744 double dDetPos = rcos_t / dU;
745 // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
746 dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
747 double dPos = dDetPos / detInc; // position along detector array
749 if (interpType == Backprojector::INTERP_NEAREST) {
750 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
751 if (iDetPos >= 0 && iDetPos < nDet)
752 pImCol[iy] += filteredProj[iDetPos] / dU2;
753 } else if (interpType == Backprojector::INTERP_LINEAR) {
754 double dPosFloor = floor (dPos);
755 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
756 double frac = dPos - dPosFloor; // fraction distance from det
757 if (iDetPos >= 0 && iDetPos < nDet - 1)
758 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2;
759 } else if (interpType == Backprojector::INTERP_CUBIC) {
760 double d = iDetCenter + dPos; // position along detector
761 if (d >= 0 && d < nDet)
762 pImCol[iy] += pCubicInterp->interpolate (d) / dU2;
767 if (interpType == Backprojector::INTERP_CUBIC)