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 x = xMin + xInc / 2; // Rectang coords of center of pixel
375 for (int ix = 0; ix < nx; x += xInc, ix++) {
376 double y = yMin + yInc / 2;
377 for (int iy = 0; iy < ny; y += yInc, iy++) {
378 double r = sqrt (x * x + y * y); // distance of cell from center
379 double phi = atan2 (y, x); // angle of cell from center
380 double L = r * cos (theta - phi); // position on detector
382 if (interpType == Backprojector::INTERP_NEAREST) {
383 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
385 if (iDetPos >= 0 && iDetPos < nDet)
386 v[ix][iy] += rotScale * filteredProj[iDetPos];
387 } else if (interpType == Backprojector::INTERP_LINEAR) {
388 double p = L / detInc; // position along detector
389 double pFloor = floor (p);
390 int iDetPos = iDetCenter + static_cast<int>(pFloor);
391 double frac = p - pFloor; // fraction distance from det
392 if (iDetPos >= 0 && iDetPos < nDet - 1)
393 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
394 } else if (interpType == Backprojector::INTERP_CUBIC) {
395 double p = iDetCenter + (L / detInc); // position along detector
396 if (p >= 0 && p < nDet)
397 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
402 if (interpType == Backprojector::INTERP_CUBIC)
407 // CLASS IDENTICATION
411 // Precalculates trigometric function value for each point in image for backprojection.
413 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
414 const int interpFactor, const ReconstructionROI* pROI)
415 : Backproject (proj, im, interpType, interpFactor, pROI)
417 arrayR.initSetSize (im.nx(), im.ny());
418 arrayPhi.initSetSize (im.nx(), im.ny());
419 r = arrayR.getArray();
420 phi = arrayPhi.getArray();
422 double x, y; // Rectang coords of center of pixel
424 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
425 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
426 r[ix][iy] = sqrt (x * x + y * y);
427 phi[ix][iy] = atan2 (y, x);
431 BackprojectTable::~BackprojectTable ()
436 BackprojectTable::PostProcessing()
438 if (! m_bPostProcessingDone) {
439 ScaleImageByRotIncrement();
440 m_bPostProcessingDone = true;
445 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
447 double theta = view_angle;
449 CubicPolyInterpolator* pCubicInterp = NULL;
450 if (interpType == Backprojector::INTERP_CUBIC)
451 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
453 for (int ix = 0; ix < nx; ix++) {
454 ImageFileColumn pImCol = v[ix];
456 for (int iy = 0; iy < ny; iy++) {
457 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
459 if (interpType == Backprojector::INTERP_NEAREST) {
460 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
462 if (iDetPos >= 0 && iDetPos < nDet)
463 pImCol[iy] += filteredProj[iDetPos];
464 } else if (interpType == Backprojector::INTERP_LINEAR) {
465 double dPos = L / detInc; // position along detector
466 double dPosFloor = floor (dPos);
467 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
468 double frac = dPos - dPosFloor; // fraction distance from det
469 if (iDetPos >= 0 && iDetPos < nDet - 1)
470 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
471 } else if (interpType == Backprojector::INTERP_CUBIC) {
472 double p = iDetCenter + (L / detInc); // position along detector
473 if (p >= 0 && p < nDet)
474 pImCol[iy] += pCubicInterp->interpolate (p);
479 if (interpType == Backprojector::INTERP_CUBIC)
484 // CLASS IDENTICATION
488 // Backprojects by precalculating the change in L position for each x & y step in the image.
489 // Iterates in x & y direction by adding difference in L position
491 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
492 const int interpFactor, const ReconstructionROI* pROI)
493 : Backproject (proj, im, interpType, interpFactor, pROI)
495 // calculate center of first pixel v[0][0]
496 double x = xMin + xInc / 2;
497 double y = yMin + yInc / 2;
498 start_r = sqrt (x * x + y * y);
499 start_phi = atan2 (y, x);
504 BackprojectDiff::~BackprojectDiff ()
509 BackprojectDiff::PostProcessing()
511 if (! m_bPostProcessingDone) {
512 ScaleImageByRotIncrement();
513 m_bPostProcessingDone = true;
518 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
520 double theta = view_angle;
522 // Distance between detectors for an angle given in units of detectors
523 double det_dx = xInc * cos (theta) / detInc;
524 double det_dy = yInc * sin (theta) / detInc;
526 // calculate detPosition for first point in image (ix=0, iy=0)
527 double detPosColStart = iDetCenter + start_r * cos (theta - start_phi) / detInc;
529 CubicPolyInterpolator* pCubicInterp = NULL;
530 double* deltaFilteredProj = NULL;
531 if (interpType == Backprojector::INTERP_LINEAR) {
532 // precalculate scaled difference for linear interpolation
533 deltaFilteredProj = new double [nDet];
534 for (int i = 0; i < nDet - 1; i++)
535 deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
536 deltaFilteredProj[nDet - 1] = 0; // last detector
537 } else if (interpType == Backprojector::INTERP_CUBIC) {
538 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
541 int iLastDet = nDet - 1;
542 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
543 double curDetPos = detPosColStart;
544 ImageFileColumn pImCol = v[ix];
546 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
547 if (interpType == Backprojector::INTERP_NEAREST) {
548 int iDetPos = nearest<int> (curDetPos); // calc index in the filtered raysum vector
550 if (iDetPos >= 0 && iDetPos < nDet)
551 *pImCol++ += filteredProj[iDetPos];
552 } else if (interpType == Backprojector::INTERP_LINEAR) {
553 double detPosFloor = floor (curDetPos);
554 int iDetPos = static_cast<int>(detPosFloor);
555 double frac = curDetPos - detPosFloor; // fraction distance from det
556 if (iDetPos >= 0 && iDetPos <= iLastDet)
557 *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
558 } else if (interpType == Backprojector::INTERP_CUBIC) {
559 double p = curDetPos; // position along detector
560 if (p >= 0 && p < nDet)
561 *pImCol++ += pCubicInterp->interpolate (p);
566 if (interpType == Backprojector::INTERP_LINEAR)
567 delete deltaFilteredProj;
568 else if (interpType == Backprojector::INTERP_CUBIC)
573 // CLASS IDENTICATION
574 // BackprojectIntDiff
577 // Highly optimized and integer version of BackprojectDiff
580 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
582 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
584 static const int scaleShift = 16;
585 #elif SIZEOF_LONG == 8
586 static const int scaleShift = 32;
588 static const long scale = (1L << scaleShift);
589 static const long scaleBitmask = scale - 1;
590 static const long halfScale = scale / 2;
591 static const double dInvScale = 1. / scale;
593 const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
594 const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
596 // calculate L for first point in image (0, 0)
597 long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
599 double* deltaFilteredProj = NULL;
600 CubicPolyInterpolator* pCubicInterp = NULL;
601 if (interpType == Backprojector::INTERP_LINEAR) {
602 // precalculate scaled difference for linear interpolation
603 deltaFilteredProj = new double [nDet];
604 for (int i = 0; i < nDet - 1; i++)
605 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
606 deltaFilteredProj[nDet - 1] = 0; // last detector
607 } else if (interpType == Backprojector::INTERP_CUBIC) {
608 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
611 int iLastDet = nDet - 1;
612 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
613 long curDetPos = detPosColStart;
614 ImageFileColumn pImCol = v[ix];
616 if (interpType == Backprojector::INTERP_NEAREST) {
617 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
618 const int iDetPos = (curDetPos + halfScale) >> scaleShift;
619 if (iDetPos >= 0 && iDetPos <= iLastDet)
620 *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];
633 } else if (interpType == Backprojector::INTERP_LINEAR) {
634 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
635 const long iDetPos = curDetPos >> scaleShift;
636 if (iDetPos >= 0 && iDetPos <= iLastDet) {
637 const long detRemainder = curDetPos & scaleBitmask;
638 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
642 } else if (interpType == Backprojector::INTERP_CUBIC) {
643 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
644 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
649 if (interpType == Backprojector::INTERP_LINEAR)
650 delete deltaFilteredProj;
651 else if (interpType == Backprojector::INTERP_CUBIC)
657 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
659 double beta = view_angle;
661 CubicPolyInterpolator* pCubicInterp = NULL;
662 if (interpType == Backprojector::INTERP_CUBIC)
663 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
665 for (int ix = 0; ix < nx; ix++) {
666 ImageFileColumn pImCol = v[ix];
668 for (int iy = 0; iy < ny; iy++) {
669 double dAngleDiff = beta - phi[ix][iy];
670 double rcos_t = r[ix][iy] * cos (dAngleDiff);
671 double rsin_t = r[ix][iy] * sin (dAngleDiff);
672 double dFLPlusSin = m_dFocalLength + rsin_t;
673 double gamma = atan (rcos_t / dFLPlusSin);
674 double dPos = gamma / detInc; // position along detector
675 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
677 if (interpType == Backprojector::INTERP_NEAREST) {
678 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
679 if (iDetPos >= 0 && iDetPos < nDet)
680 pImCol[iy] += filteredProj[iDetPos] / dL2;
681 } else if (interpType == Backprojector::INTERP_LINEAR) {
682 double dPosFloor = floor (dPos);
683 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
684 double frac = dPos - dPosFloor; // fraction distance from det
685 if (iDetPos >= 0 && iDetPos < nDet - 1)
686 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
687 } else if (interpType == Backprojector::INTERP_CUBIC) {
688 double d = iDetCenter + dPos; // position along detector
689 if (d >= 0 && d < nDet)
690 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
695 if (interpType == Backprojector::INTERP_CUBIC)
700 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
702 double beta = view_angle;
704 CubicPolyInterpolator* pCubicInterp = NULL;
705 if (interpType == Backprojector::INTERP_CUBIC)
706 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
708 for (int ix = 0; ix < nx; ix++) {
709 ImageFileColumn pImCol = v[ix];
711 for (int iy = 0; iy < ny; iy++) {
712 double dAngleDiff = beta - phi[ix][iy];
713 double rcos_t = r[ix][iy] * cos (dAngleDiff);
714 double rsin_t = r[ix][iy] * sin (dAngleDiff);
716 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
717 double dDetPos = rcos_t / dU;
718 // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
719 dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
720 double dPos = dDetPos / detInc; // position along detector array
722 if (interpType == Backprojector::INTERP_NEAREST) {
723 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
724 if (iDetPos >= 0 && iDetPos < nDet)
725 pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
726 } else if (interpType == Backprojector::INTERP_LINEAR) {
727 double dPosFloor = floor (dPos);
728 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
729 double frac = dPos - dPosFloor; // fraction distance from det
730 if (iDetPos >= 0 && iDetPos < nDet - 1)
731 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
733 } else if (interpType == Backprojector::INTERP_CUBIC) {
734 double d = iDetCenter + dPos; // position along detector
735 if (d >= 0 && d < nDet)
736 pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
741 if (interpType == Backprojector::INTERP_CUBIC)