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-2001 Kevin Rosenberg
13 ** This program is free software; you can redistribute it and/or modify
14 ** it under the terms of the GNU General Public License (version 2) as
15 ** published by the Free Software Foundation.
17 ** This program is distributed in the hope that it will be useful,
18 ** but WITHOUT ANY WARRANTY; without even the implied warranty of
19 ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 ** GNU General Public License for more details.
22 ** You should have received a copy of the GNU General Public License
23 ** along with this program; if not, write to the Free Software
24 ** Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 ******************************************************************************/
28 #include "interpolator.h"
30 const int Backprojector::BPROJ_INVALID = -1;
31 const int Backprojector::BPROJ_TRIG = 0;
32 const int Backprojector::BPROJ_TABLE = 1;
33 const int Backprojector::BPROJ_DIFF = 2;
34 const int Backprojector::BPROJ_IDIFF = 3;
36 const char* const Backprojector::s_aszBackprojectName[] =
44 const char* const Backprojector::s_aszBackprojectTitle[] =
48 "Difference Iteration",
49 "Integer Difference Iteration",
52 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
54 const int Backprojector::INTERP_INVALID = -1;
55 const int Backprojector::INTERP_NEAREST = 0;
56 const int Backprojector::INTERP_LINEAR = 1;
57 const int Backprojector::INTERP_CUBIC = 2;
58 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
59 #if HAVE_BSPLINE_INTERP
60 const int Backprojector::INTERP_BSPLINE = 4;
61 const int Backprojector::INTERP_1BSPLINE = 5;
62 const int Backprojector::INTERP_2BSPLINE = 6;
63 const int Backprojector::INTERP_3BSPLINE = 7;
66 const char* const Backprojector::s_aszInterpName[] =
71 #if HAVE_FREQ_PREINTERP
72 "freq_preinterpolationj",
74 #if HAVE_BSPLINE_INTERP
82 const char* const Backprojector::s_aszInterpTitle[] =
87 #if HAVE_FREQ_PREINTERP
88 "Frequency Preinterpolation",
90 #if HAVE_BSPLINE_INTERP
98 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
102 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
103 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
106 m_pBackprojectImplem = NULL;
108 initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
112 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
114 if (m_pBackprojectImplem != NULL)
115 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
119 Backprojector::PostProcessing()
121 if (m_pBackprojectImplem != NULL)
122 m_pBackprojectImplem->PostProcessing();
125 Backprojector::~Backprojector ()
127 delete m_pBackprojectImplem;
130 // FUNCTION IDENTIFICATION
131 // Backproject* projector = selectBackprojector (...)
134 // Selects a backprojector based on BackprojType
135 // and initializes the backprojector
138 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
139 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
141 m_nameBackproject = backprojName;
142 m_nameInterpolation = interpName;
143 m_pBackprojectImplem = NULL;
144 m_idBackproject = convertBackprojectNameToID (backprojName);
145 if (m_idBackproject == BPROJ_INVALID) {
147 m_failMessage = "Invalid backprojection name ";
148 m_failMessage += backprojName;
150 m_idInterpolation = convertInterpNameToID (interpName);
151 if (m_idInterpolation == INTERP_INVALID) {
153 m_failMessage = "Invalid interpolation name ";
154 m_failMessage += interpName;
157 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
162 if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
163 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI));
164 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
165 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI));
166 else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
167 if (m_idBackproject == BPROJ_TRIG)
168 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI));
169 else if (m_idBackproject == BPROJ_TABLE)
170 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI));
171 else if (m_idBackproject == BPROJ_DIFF)
172 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI));
173 else if (m_idBackproject == BPROJ_IDIFF)
174 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI));
177 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
186 Backprojector::convertBackprojectNameToID (const char* const backprojName)
188 int backprojID = BPROJ_INVALID;
190 for (int i = 0; i < s_iBackprojectCount; i++)
191 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
200 Backprojector::convertBackprojectIDToName (int bprojID)
202 static const char *bprojName = "";
204 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
205 return (s_aszBackprojectName[bprojID]);
211 Backprojector::convertBackprojectIDToTitle (const int bprojID)
213 static const char *bprojTitle = "";
215 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
216 return (s_aszBackprojectTitle[bprojID]);
223 Backprojector::convertInterpNameToID (const char* const interpName)
225 int interpID = INTERP_INVALID;
227 for (int i = 0; i < s_iInterpCount; i++)
228 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
237 Backprojector::convertInterpIDToName (const int interpID)
239 static const char *interpName = "";
241 if (interpID >= 0 && interpID < s_iInterpCount)
242 return (s_aszInterpName[interpID]);
248 Backprojector::convertInterpIDToTitle (const int interpID)
250 static const char *interpTitle = "";
252 if (interpID >= 0 && interpID < s_iInterpCount)
253 return (s_aszInterpTitle[interpID]);
255 return (interpTitle);
260 // CLASS IDENTICATION
264 // Pure virtual base class for all backprojectors.
266 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor,
267 const ReconstructionROI* pROI)
268 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false)
270 detInc = proj.detInc();
272 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
273 rotScale = proj.rotInc();
275 if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
276 rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
277 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
278 rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
280 sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
287 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
288 xMax = xMin + proj.phmLen();
289 yMin = -proj.phmLen() / 2;
290 yMax = yMin + proj.phmLen();
293 if (pROI->m_dXMin > xMin)
294 xMin = pROI->m_dXMin;
295 if (pROI->m_dXMax < xMax)
296 xMax = pROI->m_dXMax;
297 if (pROI->m_dYMin > yMin)
298 yMin = pROI->m_dYMin;
299 if (pROI->m_dYMax < yMax)
300 yMax = pROI->m_dYMax;
314 xInc = (xMax - xMin) / nx; // size of cells
315 yInc = (yMax - yMin) / ny;
317 im.setAxisIncrement (xInc, yInc);
318 im.setAxisExtent (xMin, xMax, yMin, yMax);
320 m_dFocalLength = proj.focalLength();
321 m_dSourceDetectorLength = proj.sourceDetectorLength();
324 Backproject::~Backproject ()
328 Backproject::PostProcessing()
330 m_bPostProcessingDone = true;
334 Backproject::ScaleImageByRotIncrement ()
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 x = xMin + xInc / 2; // Rectang coords of center of pixel
378 for (int ix = 0; ix < nx; x += xInc, ix++) {
379 double y = yMin + yInc / 2;
380 for (int iy = 0; iy < ny; y += yInc, iy++) {
381 double r = sqrt (x * x + y * y); // distance of cell from center
382 double phi = atan2 (y, x); // angle of cell from center
383 double L = r * cos (theta - phi); // position on detector
385 if (interpType == Backprojector::INTERP_NEAREST) {
386 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
388 if (iDetPos >= 0 && iDetPos < nDet)
389 v[ix][iy] += rotScale * filteredProj[iDetPos];
390 } else if (interpType == Backprojector::INTERP_LINEAR) {
391 double p = L / detInc; // position along detector
392 double pFloor = floor (p);
393 int iDetPos = iDetCenter + static_cast<int>(pFloor);
394 double frac = p - pFloor; // fraction distance from det
395 if (iDetPos >= 0 && iDetPos < nDet - 1)
396 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
397 } else if (interpType == Backprojector::INTERP_CUBIC) {
398 double p = iDetCenter + (L / detInc); // position along detector
399 if (p >= 0 && p < nDet)
400 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
405 if (interpType == Backprojector::INTERP_CUBIC)
410 // CLASS IDENTICATION
414 // Precalculates trigometric function value for each point in image for backprojection.
416 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
417 const int interpFactor, const ReconstructionROI* pROI)
418 : Backproject (proj, im, interpType, interpFactor, pROI)
420 arrayR.initSetSize (im.nx(), im.ny());
421 arrayPhi.initSetSize (im.nx(), im.ny());
422 r = arrayR.getArray();
423 phi = arrayPhi.getArray();
425 double x, y; // Rectang coords of center of pixel
427 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
428 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
429 r[ix][iy] = sqrt (x * x + y * y);
430 phi[ix][iy] = atan2 (y, x);
434 BackprojectTable::~BackprojectTable ()
439 BackprojectTable::PostProcessing()
441 if (! m_bPostProcessingDone) {
442 ScaleImageByRotIncrement();
443 m_bPostProcessingDone = true;
448 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
450 double theta = view_angle;
452 CubicPolyInterpolator* pCubicInterp = NULL;
453 if (interpType == Backprojector::INTERP_CUBIC)
454 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
456 for (int ix = 0; ix < nx; ix++) {
457 ImageFileColumn pImCol = v[ix];
459 for (int iy = 0; iy < ny; iy++) {
460 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
462 if (interpType == Backprojector::INTERP_NEAREST) {
463 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
465 if (iDetPos >= 0 && iDetPos < nDet)
466 pImCol[iy] += filteredProj[iDetPos];
467 } else if (interpType == Backprojector::INTERP_LINEAR) {
468 double dPos = L / detInc; // position along detector
469 double dPosFloor = floor (dPos);
470 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
471 double frac = dPos - dPosFloor; // fraction distance from det
472 if (iDetPos >= 0 && iDetPos < nDet - 1)
473 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
474 } else if (interpType == Backprojector::INTERP_CUBIC) {
475 double p = iDetCenter + (L / detInc); // position along detector
476 if (p >= 0 && p < nDet)
477 pImCol[iy] += pCubicInterp->interpolate (p);
482 if (interpType == Backprojector::INTERP_CUBIC)
487 // CLASS IDENTICATION
491 // Backprojects by precalculating the change in L position for each x & y step in the image.
492 // Iterates in x & y direction by adding difference in L position
494 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
495 const int interpFactor, const ReconstructionROI* pROI)
496 : Backproject (proj, im, interpType, interpFactor, pROI)
498 // calculate center of first pixel v[0][0]
499 double x = xMin + xInc / 2;
500 double y = yMin + yInc / 2;
501 start_r = sqrt (x * x + y * y);
502 start_phi = atan2 (y, x);
507 BackprojectDiff::~BackprojectDiff ()
512 BackprojectDiff::PostProcessing()
514 if (! m_bPostProcessingDone) {
515 ScaleImageByRotIncrement();
516 m_bPostProcessingDone = true;
521 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
523 double theta = view_angle;
525 // Distance between detectors for an angle given in units of detectors
526 double det_dx = xInc * cos (theta) / detInc;
527 double det_dy = yInc * sin (theta) / detInc;
529 // calculate detPosition for first point in image (ix=0, iy=0)
530 double detPosColStart = iDetCenter + start_r * cos (theta - start_phi) / detInc;
532 CubicPolyInterpolator* pCubicInterp = NULL;
533 double* deltaFilteredProj = NULL;
534 if (interpType == Backprojector::INTERP_LINEAR) {
535 // precalculate scaled difference for linear interpolation
536 deltaFilteredProj = new double [nDet];
537 for (int i = 0; i < nDet - 1; i++)
538 deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
539 deltaFilteredProj[nDet - 1] = 0; // last detector
540 } else if (interpType == Backprojector::INTERP_CUBIC) {
541 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
544 int iLastDet = nDet - 1;
545 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
546 double curDetPos = detPosColStart;
547 ImageFileColumn pImCol = v[ix];
549 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
550 if (interpType == Backprojector::INTERP_NEAREST) {
551 int iDetPos = nearest<int> (curDetPos); // calc index in the filtered raysum vector
553 if (iDetPos >= 0 && iDetPos < nDet)
554 *pImCol++ += filteredProj[iDetPos];
555 } else if (interpType == Backprojector::INTERP_LINEAR) {
556 double detPosFloor = floor (curDetPos);
557 int iDetPos = static_cast<int>(detPosFloor);
558 double frac = curDetPos - detPosFloor; // fraction distance from det
559 if (iDetPos >= 0 && iDetPos <= iLastDet)
560 *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
561 } else if (interpType == Backprojector::INTERP_CUBIC) {
562 double p = curDetPos; // position along detector
563 if (p >= 0 && p < nDet)
564 *pImCol++ += pCubicInterp->interpolate (p);
569 if (interpType == Backprojector::INTERP_LINEAR)
570 delete deltaFilteredProj;
571 else if (interpType == Backprojector::INTERP_CUBIC)
576 // CLASS IDENTICATION
577 // BackprojectIntDiff
580 // Highly optimized and integer version of BackprojectDiff
583 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
585 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
587 static const int scaleShift = 16;
588 #elif SIZEOF_LONG == 8
589 static const int scaleShift = 32;
591 static const long scale = (1L << scaleShift);
592 static const long scaleBitmask = scale - 1;
593 static const long halfScale = scale / 2;
594 static const double dInvScale = 1. / scale;
596 const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
597 const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
599 // calculate L for first point in image (0, 0)
600 long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
602 double* deltaFilteredProj = NULL;
603 CubicPolyInterpolator* pCubicInterp = NULL;
604 if (interpType == Backprojector::INTERP_LINEAR) {
605 // precalculate scaled difference for linear interpolation
606 deltaFilteredProj = new double [nDet];
607 for (int i = 0; i < nDet - 1; i++)
608 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
609 deltaFilteredProj[nDet - 1] = 0; // last detector
610 } else if (interpType == Backprojector::INTERP_CUBIC) {
611 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
614 int iLastDet = nDet - 1;
615 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
616 long curDetPos = detPosColStart;
617 ImageFileColumn pImCol = v[ix];
619 if (interpType == Backprojector::INTERP_NEAREST) {
620 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
621 const int iDetPos = (curDetPos + halfScale) >> scaleShift;
622 if (iDetPos >= 0 && iDetPos <= iLastDet)
623 *pImCol++ += filteredProj[iDetPos];
625 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
626 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
627 const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
628 if (iDetPos >= 0 && iDetPos <= iLastDet)
629 *pImCol++ += filteredProj[iDetPos];
631 } else if (interpType == Backprojector::INTERP_LINEAR) {
632 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
633 const long iDetPos = curDetPos >> scaleShift;
634 if (iDetPos >= 0 && iDetPos <= iLastDet) {
635 const long detRemainder = curDetPos & scaleBitmask;
636 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
639 } else if (interpType == Backprojector::INTERP_CUBIC) {
640 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
641 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
646 if (interpType == Backprojector::INTERP_LINEAR)
647 delete deltaFilteredProj;
648 else if (interpType == Backprojector::INTERP_CUBIC)
654 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
656 double beta = view_angle;
658 CubicPolyInterpolator* pCubicInterp = NULL;
659 if (interpType == Backprojector::INTERP_CUBIC)
660 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
662 for (int ix = 0; ix < nx; ix++) {
663 ImageFileColumn pImCol = v[ix];
665 for (int iy = 0; iy < ny; iy++) {
666 double dAngleDiff = beta - phi[ix][iy];
667 double rcos_t = r[ix][iy] * cos (dAngleDiff);
668 double rsin_t = r[ix][iy] * sin (dAngleDiff);
669 double dFLPlusSin = m_dFocalLength + rsin_t;
670 double gamma = atan (rcos_t / dFLPlusSin);
671 double dPos = gamma / detInc; // position along detector
672 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
674 if (interpType == Backprojector::INTERP_NEAREST) {
675 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
676 if (iDetPos >= 0 && iDetPos < nDet)
677 pImCol[iy] += filteredProj[iDetPos] / dL2;
678 } else if (interpType == Backprojector::INTERP_LINEAR) {
679 double dPosFloor = floor (dPos);
680 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
681 double frac = dPos - dPosFloor; // fraction distance from det
682 if (iDetPos >= 0 && iDetPos < nDet - 1)
683 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
684 } else if (interpType == Backprojector::INTERP_CUBIC) {
685 double d = iDetCenter + dPos; // position along detector
686 if (d >= 0 && d < nDet)
687 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
692 if (interpType == Backprojector::INTERP_CUBIC)
697 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
699 double beta = view_angle;
701 CubicPolyInterpolator* pCubicInterp = NULL;
702 if (interpType == Backprojector::INTERP_CUBIC)
703 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
705 for (int ix = 0; ix < nx; ix++) {
706 ImageFileColumn pImCol = v[ix];
708 for (int iy = 0; iy < ny; iy++) {
709 double dAngleDiff = beta - phi[ix][iy];
710 double rcos_t = r[ix][iy] * cos (dAngleDiff);
711 double rsin_t = r[ix][iy] * sin (dAngleDiff);
713 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
714 double dDetPos = rcos_t / dU;
715 // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
716 dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
717 double dPos = dDetPos / detInc; // position along detector array
719 if (interpType == Backprojector::INTERP_NEAREST) {
720 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
721 if (iDetPos >= 0 && iDetPos < nDet)
722 pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
723 } else if (interpType == Backprojector::INTERP_LINEAR) {
724 double dPosFloor = floor (dPos);
725 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
726 double frac = dPos - dPosFloor; // fraction distance from det
727 if (iDetPos >= 0 && iDetPos < nDet - 1)
728 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
730 } else if (interpType == Backprojector::INTERP_CUBIC) {
731 double d = iDetCenter + dPos; // position along detector
732 if (d >= 0 && d < nDet)
733 pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
738 if (interpType == Backprojector::INTERP_CUBIC)