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 ()
334 #pragma omp parallel for
336 for (int ix = 0; ix < nx; ix++)
337 for (int iy = 0; iy < ny; iy++)
338 v[ix][iy] *= rotScale;
341 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
343 sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
344 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
347 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
350 std::ostringstream os;
351 os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
352 os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
353 os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
354 os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
355 os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
357 sys_error (ERR_WARNING, os.str().c_str());
362 // CLASS IDENTICATION
366 // Uses trigometric functions at each point in image for backprojection.
369 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
371 double theta = view_angle;
373 CubicPolyInterpolator* pCubicInterp = NULL;
374 if (interpType == Backprojector::INTERP_CUBIC)
375 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
377 double xstart = xMin + xInc / 2; // Rectang coords of center of pixel
379 #pragma omp parallel for
381 for (int ix = 0; ix < nx; ix++) {
382 double x = xstart + (ix * xInc);
384 double y = yMin + yInc / 2;
385 for (int iy = 0; iy < ny; y += yInc, iy++) {
386 double r = sqrt (x * x + y * y); // distance of cell from center
387 double phi = atan2 (y, x); // angle of cell from center
388 double L = r * cos (theta - phi); // position on detector
390 if (interpType == Backprojector::INTERP_NEAREST) {
391 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
393 if (iDetPos >= 0 && iDetPos < nDet)
394 v[ix][iy] += rotScale * filteredProj[iDetPos];
395 } else if (interpType == Backprojector::INTERP_LINEAR) {
396 double p = L / detInc; // position along detector
397 double pFloor = floor (p);
398 int iDetPos = iDetCenter + static_cast<int>(pFloor);
399 double frac = p - pFloor; // fraction distance from det
400 if (iDetPos >= 0 && iDetPos < nDet - 1)
401 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
402 } else if (interpType == Backprojector::INTERP_CUBIC) {
403 double p = iDetCenter + (L / detInc); // position along detector
404 if (p >= 0 && p < nDet)
405 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
410 if (interpType == Backprojector::INTERP_CUBIC)
415 // CLASS IDENTICATION
419 // Precalculates trigometric function value for each point in image for backprojection.
421 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
422 const int interpFactor, const ReconstructionROI* pROI)
423 : Backproject (proj, im, interpType, interpFactor, pROI)
425 arrayR.initSetSize (im.nx(), im.ny());
426 arrayPhi.initSetSize (im.nx(), im.ny());
427 r = arrayR.getArray();
428 phi = arrayPhi.getArray();
430 double xstart = xMin + xInc / 2;
433 #pragma omp parallel for
435 for (int ix = 0; ix < nx; ix++) {
436 double x = xstart + (ix * xInc);
437 double y = yMin + yInc / 2;
438 for (int iy = 0; iy < ny; iy++, y += yInc) {
439 r[ix][iy] = sqrt (x * x + y * y);
440 phi[ix][iy] = atan2 (y, x);
445 BackprojectTable::~BackprojectTable ()
450 BackprojectTable::PostProcessing()
452 if (! m_bPostProcessingDone) {
453 ScaleImageByRotIncrement();
454 m_bPostProcessingDone = true;
459 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
461 double theta = view_angle;
463 CubicPolyInterpolator* pCubicInterp = NULL;
464 if (interpType == Backprojector::INTERP_CUBIC)
465 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
468 #pragma omp parallel for
470 for (int ix = 0; ix < nx; ix++) {
471 ImageFileColumn pImCol = v[ix];
473 for (int iy = 0; iy < ny; iy++) {
474 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
476 if (interpType == Backprojector::INTERP_NEAREST) {
477 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
479 if (iDetPos >= 0 && iDetPos < nDet) {
480 pImCol[iy] += filteredProj[iDetPos];
482 } else if (interpType == Backprojector::INTERP_LINEAR) {
483 double dPos = L / detInc; // position along detector
484 double dPosFloor = floor (dPos);
485 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
486 double frac = dPos - dPosFloor; // fraction distance from det
487 if (iDetPos >= 0 && iDetPos < nDet - 1) {
488 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
490 } else if (interpType == Backprojector::INTERP_CUBIC) {
491 double p = iDetCenter + (L / detInc); // position along detector
492 if (p >= 0 && p < nDet) {
493 pImCol[iy] += pCubicInterp->interpolate (p);
499 if (interpType == Backprojector::INTERP_CUBIC)
504 // CLASS IDENTICATION
508 // Backprojects by precalculating the change in L position for each x & y step in the image.
509 // Iterates in x & y direction by adding difference in L position
511 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
512 const int interpFactor, const ReconstructionROI* pROI)
513 : Backproject (proj, im, interpType, interpFactor, pROI)
515 // calculate center of first pixel v[0][0]
516 double x = xMin + xInc / 2;
517 double y = yMin + yInc / 2;
518 start_r = sqrt (x * x + y * y);
519 start_phi = atan2 (y, x);
524 BackprojectDiff::~BackprojectDiff ()
529 BackprojectDiff::PostProcessing()
531 if (! m_bPostProcessingDone) {
532 ScaleImageByRotIncrement();
533 m_bPostProcessingDone = true;
538 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
540 double theta = view_angle;
542 // Distance between detectors for an angle given in units of detectors
543 double det_dx = xInc * cos (theta) / detInc;
544 double det_dy = yInc * sin (theta) / detInc;
546 // calculate detPosition for first point in image (ix=0, iy=0)
547 double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc;
549 CubicPolyInterpolator* pCubicInterp = NULL;
550 double* deltaFilteredProj = NULL;
551 if (interpType == Backprojector::INTERP_LINEAR) {
552 // precalculate scaled difference for linear interpolation
553 deltaFilteredProj = new double [nDet];
555 #pragma omp parallel for
557 for (int i = 0; i < nDet - 1; i++)
558 deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
559 deltaFilteredProj[nDet - 1] = 0; // last detector
560 } else if (interpType == Backprojector::INTERP_CUBIC) {
561 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
564 int iLastDet = nDet - 1;
567 #pragma omp parallel for
569 for (int ix = 0; ix < nx; ix++) {
570 double detPos = detPosColBase + (ix * det_dx);
571 ImageFileColumn pImCol = v[ix];
573 for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
574 if (interpType == Backprojector::INTERP_NEAREST) {
575 int iDetPos = nearest<int> (detPos); // calc index in the filtered raysum vector
576 if (iDetPos >= 0 && iDetPos < nDet) {
577 *pImCol++ += filteredProj[iDetPos];
579 } else if (interpType == Backprojector::INTERP_LINEAR) {
580 double detPosFloor = floor (detPos);
581 int iDetPos = static_cast<int>(detPosFloor);
582 double frac = detPos - detPosFloor; // fraction distance from det
583 if (iDetPos >= 0 && iDetPos <= iLastDet) {
584 *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
586 } else if (interpType == Backprojector::INTERP_CUBIC) {
587 double p = detPos; // position along detector
588 if (p >= 0 && p < nDet) {
589 *pImCol++ += pCubicInterp->interpolate (p);
595 if (interpType == Backprojector::INTERP_LINEAR)
596 delete deltaFilteredProj;
597 else if (interpType == Backprojector::INTERP_CUBIC)
602 // CLASS IDENTICATION
603 // BackprojectIntDiff
606 // Highly optimized and integer version of BackprojectDiff
609 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
611 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
613 static const int scaleShift = 16;
614 #elif SIZEOF_LONG == 8
615 static const int scaleShift = 32;
617 static const long scale = (1L << scaleShift);
618 static const long scaleBitmask = scale - 1;
619 static const long halfScale = scale / 2;
620 static const double dInvScale = 1. / scale;
622 const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
623 const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
625 // calculate L for first point in image (0, 0)
626 long detPosColBase = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
628 double* deltaFilteredProj = NULL;
629 CubicPolyInterpolator* pCubicInterp = NULL;
630 if (interpType == Backprojector::INTERP_LINEAR) {
631 // precalculate scaled difference for linear interpolation
632 deltaFilteredProj = new double [nDet];
634 #pragma omp parallel for
636 for (int i = 0; i < nDet - 1; i++)
637 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
638 deltaFilteredProj[nDet - 1] = 0; // last detector
639 } else if (interpType == Backprojector::INTERP_CUBIC) {
640 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
643 int iLastDet = nDet - 1;
645 #pragma omp parallel for
647 for (int ix = 0; ix < nx; ix++) {
648 long detPos = detPosColBase + (ix * det_dx);
649 ImageFileColumn pImCol = v[ix];
651 for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
652 if (interpType == Backprojector::INTERP_NEAREST) {
653 const int iDetPos = (detPos + halfScale) >> scaleShift;
654 if (iDetPos >= 0 && iDetPos <= iLastDet) {
655 *pImCol++ += filteredProj[iDetPos];
658 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
659 const int iDetPos = ((detPos + halfScale) >> scaleShift) * m_interpFactor;
660 if (iDetPos >= 0 && iDetPos <= iLastDet) {
661 *pImCol++ += filteredProj[iDetPos];
664 } else if (interpType == Backprojector::INTERP_LINEAR) {
665 const long iDetPos = detPos >> scaleShift;
666 if (iDetPos >= 0 && iDetPos <= iLastDet) {
667 const long detRemainder = detPos & scaleBitmask;
668 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
671 } else if (interpType == Backprojector::INTERP_CUBIC) {
672 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(detPos) / scale);
677 if (interpType == Backprojector::INTERP_LINEAR)
678 delete deltaFilteredProj;
679 else if (interpType == Backprojector::INTERP_CUBIC)
685 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
687 double beta = view_angle;
689 CubicPolyInterpolator* pCubicInterp = NULL;
690 if (interpType == Backprojector::INTERP_CUBIC)
691 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
694 #pragma omp parallel for
696 for (int ix = 0; ix < nx; ix++) {
697 ImageFileColumn pImCol = v[ix];
699 for (int iy = 0; iy < ny; iy++) {
700 double dAngleDiff = beta - phi[ix][iy];
701 double rcos_t = r[ix][iy] * cos (dAngleDiff);
702 double rsin_t = r[ix][iy] * sin (dAngleDiff);
703 double dFLPlusSin = m_dFocalLength + rsin_t;
704 double gamma = atan (rcos_t / dFLPlusSin);
705 double dPos = gamma / detInc; // position along detector
706 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
708 if (interpType == Backprojector::INTERP_NEAREST) {
709 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
710 if (iDetPos >= 0 && iDetPos < nDet)
711 pImCol[iy] += filteredProj[iDetPos] / dL2;
712 } else if (interpType == Backprojector::INTERP_LINEAR) {
713 double dPosFloor = floor (dPos);
714 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
715 double frac = dPos - dPosFloor; // fraction distance from det
716 if (iDetPos >= 0 && iDetPos < nDet - 1)
717 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
718 } else if (interpType == Backprojector::INTERP_CUBIC) {
719 double d = iDetCenter + dPos; // position along detector
720 if (d >= 0 && d < nDet)
721 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
726 if (interpType == Backprojector::INTERP_CUBIC)
731 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
733 double beta = view_angle;
735 CubicPolyInterpolator* pCubicInterp = NULL;
736 if (interpType == Backprojector::INTERP_CUBIC)
737 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
740 #pragma omp parallel for
742 for (int ix = 0; ix < nx; ix++) {
743 ImageFileColumn pImCol = v[ix];
745 for (int iy = 0; iy < ny; iy++) {
746 double dAngleDiff = beta - phi[ix][iy];
747 double rcos_t = r[ix][iy] * cos (dAngleDiff);
748 double rsin_t = r[ix][iy] * sin (dAngleDiff);
750 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
751 double dU2 = dU * dU;
753 double dDetPos = rcos_t / dU;
754 // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
755 dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
756 double dPos = dDetPos / detInc; // position along detector array
758 if (interpType == Backprojector::INTERP_NEAREST) {
759 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
760 if (iDetPos >= 0 && iDetPos < nDet)
761 pImCol[iy] += filteredProj[iDetPos] / dU2;
762 } else if (interpType == Backprojector::INTERP_LINEAR) {
763 double dPosFloor = floor (dPos);
764 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
765 double frac = dPos - dPosFloor; // fraction distance from det
766 if (iDetPos >= 0 && iDetPos < nDet - 1)
767 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2;
768 } else if (interpType == Backprojector::INTERP_CUBIC) {
769 double d = iDetCenter + dPos; // position along detector
770 if (d >= 0 && d < nDet)
771 pImCol[iy] += pCubicInterp->interpolate (d) / dU2;
776 if (interpType == Backprojector::INTERP_CUBIC)