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
11 ** $Id: backprojectors.cpp,v 1.31 2001/03/11 15:27:30 kevin Exp $
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 ******************************************************************************/
29 const int Backprojector::BPROJ_INVALID = -1;
30 const int Backprojector::BPROJ_TRIG = 0;
31 const int Backprojector::BPROJ_TABLE = 1;
32 const int Backprojector::BPROJ_DIFF = 2;
33 const int Backprojector::BPROJ_IDIFF = 3;
35 const char* const Backprojector::s_aszBackprojectName[] =
43 const char* const Backprojector::s_aszBackprojectTitle[] =
45 {"Direct Trigometric"},
46 {"Trigometric Table"},
47 {"Difference Iteration"},
48 {"Integer Difference Iteration"},
51 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
53 const int Backprojector::INTERP_INVALID = -1;
54 const int Backprojector::INTERP_NEAREST = 0;
55 const int Backprojector::INTERP_LINEAR = 1;
56 const int Backprojector::INTERP_CUBIC = 2;
57 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
58 #if HAVE_BSPLINE_INTERP
59 const int Backprojector::INTERP_BSPLINE = 4;
60 const int Backprojector::INTERP_1BSPLINE = 5;
61 const int Backprojector::INTERP_2BSPLINE = 6;
62 const int Backprojector::INTERP_3BSPLINE = 7;
65 const char* const Backprojector::s_aszInterpName[] =
70 #if HAVE_FREQ_PREINTERP
71 {"freq_preinterpolationj"},
73 #if HAVE_BSPLINE_INTERP
81 const char* const Backprojector::s_aszInterpTitle[] =
86 #if HAVE_FREQ_PREINTERP
87 {"Frequency Preinterpolation"},
89 #if HAVE_BSPLINE_INTERP
91 {"B-Spline 1st Order"},
92 {"B-Spline 2nd Order"},
93 {"B-Spline 3rd Order"},
97 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
101 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
102 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
105 m_pBackprojectImplem = NULL;
107 initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
111 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
113 if (m_pBackprojectImplem != NULL)
114 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
118 Backprojector::PostProcessing()
120 if (m_pBackprojectImplem != NULL)
121 m_pBackprojectImplem->PostProcessing();
124 Backprojector::~Backprojector ()
126 delete m_pBackprojectImplem;
129 // FUNCTION IDENTIFICATION
130 // Backproject* projector = selectBackprojector (...)
133 // Selects a backprojector based on BackprojType
134 // and initializes the backprojector
137 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
138 const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
140 m_nameBackproject = backprojName;
141 m_nameInterpolation = interpName;
142 m_pBackprojectImplem = NULL;
143 m_idBackproject = convertBackprojectNameToID (backprojName);
144 if (m_idBackproject == BPROJ_INVALID) {
146 m_failMessage = "Invalid backprojection name ";
147 m_failMessage += backprojName;
149 m_idInterpolation = convertInterpNameToID (interpName);
150 if (m_idInterpolation == INTERP_INVALID) {
152 m_failMessage = "Invalid interpolation name ";
153 m_failMessage += interpName;
156 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
161 if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
162 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI));
163 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
164 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI));
165 else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
166 if (m_idBackproject == BPROJ_TRIG)
167 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI));
168 else if (m_idBackproject == BPROJ_TABLE)
169 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI));
170 else if (m_idBackproject == BPROJ_DIFF)
171 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI));
172 else if (m_idBackproject == BPROJ_IDIFF)
173 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI));
176 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
185 Backprojector::convertBackprojectNameToID (const char* const backprojName)
187 int backprojID = BPROJ_INVALID;
189 for (int i = 0; i < s_iBackprojectCount; i++)
190 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
199 Backprojector::convertBackprojectIDToName (int bprojID)
201 static const char *bprojName = "";
203 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
204 return (s_aszBackprojectName[bprojID]);
210 Backprojector::convertBackprojectIDToTitle (const int bprojID)
212 static const char *bprojTitle = "";
214 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
215 return (s_aszBackprojectTitle[bprojID]);
222 Backprojector::convertInterpNameToID (const char* const interpName)
224 int interpID = INTERP_INVALID;
226 for (int i = 0; i < s_iInterpCount; i++)
227 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
236 Backprojector::convertInterpIDToName (const int interpID)
238 static const char *interpName = "";
240 if (interpID >= 0 && interpID < s_iInterpCount)
241 return (s_aszInterpName[interpID]);
247 Backprojector::convertInterpIDToTitle (const int interpID)
249 static const char *interpTitle = "";
251 if (interpID >= 0 && interpID < s_iInterpCount)
252 return (s_aszInterpTitle[interpID]);
254 return (interpTitle);
259 // CLASS IDENTICATION
263 // Pure virtual base class for all backprojectors.
265 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor,
266 const ReconstructionROI* pROI)
267 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false)
269 detInc = proj.detInc();
271 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
272 rotScale = proj.rotInc();
274 if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
275 rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
276 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
277 rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
279 sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
286 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
287 xMax = xMin + proj.phmLen();
288 yMin = -proj.phmLen() / 2;
289 yMax = yMin + proj.phmLen();
292 if (pROI->m_dXMin > xMin)
293 xMin = pROI->m_dXMin;
294 if (pROI->m_dXMax < xMax)
295 xMax = pROI->m_dXMax;
296 if (pROI->m_dYMin > yMin)
297 yMin = pROI->m_dYMin;
298 if (pROI->m_dYMax < yMax)
299 yMax = pROI->m_dYMax;
313 xInc = (xMax - xMin) / nx; // size of cells
314 yInc = (yMax - yMin) / ny;
316 m_dFocalLength = proj.focalLength();
317 m_dSourceDetectorLength = proj.sourceDetectorLength();
320 Backproject::~Backproject ()
324 Backproject::PostProcessing()
326 m_bPostProcessingDone = true;
330 Backproject::ScaleImageByRotIncrement ()
332 for (int ix = 0; ix < nx; ix++)
333 for (int iy = 0; iy < ny; iy++)
334 v[ix][iy] *= rotScale;
337 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
339 sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
340 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
343 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
346 std::ostringstream os;
347 os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
348 os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
349 os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
350 os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
351 os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
353 sys_error (ERR_WARNING, os.str().c_str());
358 // CLASS IDENTICATION
362 // Uses trigometric functions at each point in image for backprojection.
365 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
367 double theta = view_angle;
369 CubicPolyInterpolator* pCubicInterp = NULL;
370 if (interpType == Backprojector::INTERP_CUBIC)
371 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
373 double x = xMin + xInc / 2; // Rectang coords of center of pixel
374 for (int ix = 0; ix < nx; x += xInc, ix++) {
375 double y = yMin + yInc / 2;
376 for (int iy = 0; iy < ny; y += yInc, iy++) {
377 double r = sqrt (x * x + y * y); // distance of cell from center
378 double phi = atan2 (y, x); // angle of cell from center
379 double L = r * cos (theta - phi); // position on detector
381 if (interpType == Backprojector::INTERP_NEAREST) {
382 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
384 if (iDetPos >= 0 && iDetPos < nDet)
385 v[ix][iy] += rotScale * filteredProj[iDetPos];
386 } else if (interpType == Backprojector::INTERP_LINEAR) {
387 double p = L / detInc; // position along detector
388 double pFloor = floor (p);
389 int iDetPos = iDetCenter + static_cast<int>(pFloor);
390 double frac = p - pFloor; // fraction distance from det
391 if (iDetPos >= 0 && iDetPos < nDet - 1)
392 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
393 } else if (interpType == Backprojector::INTERP_CUBIC) {
394 double p = iDetCenter + (L / detInc); // position along detector
395 if (p >= 0 && p < nDet)
396 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
401 if (interpType == Backprojector::INTERP_CUBIC)
406 // CLASS IDENTICATION
410 // Precalculates trigometric function value for each point in image for backprojection.
412 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
413 const int interpFactor, const ReconstructionROI* pROI)
414 : Backproject (proj, im, interpType, interpFactor, pROI)
416 arrayR.initSetSize (im.nx(), im.ny());
417 arrayPhi.initSetSize (im.nx(), im.ny());
418 r = arrayR.getArray();
419 phi = arrayPhi.getArray();
421 double x, y; // Rectang coords of center of pixel
423 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
424 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
425 r[ix][iy] = sqrt (x * x + y * y);
426 phi[ix][iy] = atan2 (y, x);
430 BackprojectTable::~BackprojectTable ()
435 BackprojectTable::PostProcessing()
437 if (! m_bPostProcessingDone) {
438 ScaleImageByRotIncrement();
439 m_bPostProcessingDone = true;
444 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
446 double theta = view_angle;
448 CubicPolyInterpolator* pCubicInterp = NULL;
449 if (interpType == Backprojector::INTERP_CUBIC)
450 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
452 for (int ix = 0; ix < nx; ix++) {
453 ImageFileColumn pImCol = v[ix];
455 for (int iy = 0; iy < ny; iy++) {
456 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
458 if (interpType == Backprojector::INTERP_NEAREST) {
459 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
461 if (iDetPos >= 0 && iDetPos < nDet)
462 pImCol[iy] += filteredProj[iDetPos];
463 } else if (interpType == Backprojector::INTERP_LINEAR) {
464 double dPos = L / detInc; // position along detector
465 double dPosFloor = floor (dPos);
466 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
467 double frac = dPos - dPosFloor; // fraction distance from det
468 if (iDetPos >= 0 && iDetPos < nDet - 1)
469 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
470 } else if (interpType == Backprojector::INTERP_CUBIC) {
471 double p = iDetCenter + (L / detInc); // position along detector
472 if (p >= 0 && p < nDet)
473 pImCol[iy] += pCubicInterp->interpolate (p);
478 if (interpType == Backprojector::INTERP_CUBIC)
483 // CLASS IDENTICATION
487 // Backprojects by precalculating the change in L position for each x & y step in the image.
488 // Iterates in x & y direction by adding difference in L position
490 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
491 const int interpFactor, const ReconstructionROI* pROI)
492 : Backproject (proj, im, interpType, interpFactor, pROI)
494 // calculate center of first pixel v[0][0]
495 double x = xMin + xInc / 2;
496 double y = yMin + yInc / 2;
497 start_r = sqrt (x * x + y * y);
498 start_phi = atan2 (y, x);
503 BackprojectDiff::~BackprojectDiff ()
508 BackprojectDiff::PostProcessing()
510 if (! m_bPostProcessingDone) {
511 ScaleImageByRotIncrement();
512 m_bPostProcessingDone = true;
517 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
519 double theta = view_angle;
521 // Distance between detectors for an angle given in units of detectors
522 double det_dx = xInc * cos (theta) / detInc;
523 double det_dy = yInc * sin (theta) / detInc;
525 // calculate detPosition for first point in image (ix=0, iy=0)
526 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
528 CubicPolyInterpolator* pCubicInterp = NULL;
529 if (interpType == Backprojector::INTERP_CUBIC)
530 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
532 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
533 double curDetPos = detPosColStart;
534 ImageFileColumn pImCol = v[ix];
536 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
537 if (interpType == Backprojector::INTERP_NEAREST) {
538 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
540 if (iDetPos >= 0 && iDetPos < nDet)
541 *pImCol++ += filteredProj[iDetPos];
542 } else if (interpType == Backprojector::INTERP_LINEAR) {
543 double detPosFloor = floor (curDetPos);
544 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
545 double frac = curDetPos - detPosFloor; // fraction distance from det
546 if (iDetPos > 0 && iDetPos < nDet - 1)
547 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
548 } else if (interpType == Backprojector::INTERP_CUBIC) {
549 double p = iDetCenter + curDetPos; // position along detector
550 if (p >= 0 && p < nDet)
551 *pImCol++ += pCubicInterp->interpolate (p);
556 if (interpType == Backprojector::INTERP_CUBIC)
561 // CLASS IDENTICATION
562 // BackprojectIntDiff
565 // Highly optimized and integer version of BackprojectDiff
568 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
570 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
571 static const int scaleShift = 16;
572 static const kint32 scale = (1 << scaleShift);
573 static const kint32 scaleBitmask = scale - 1;
574 static const kint32 halfScale = scale / 2;
575 static const double dInvScale = 1. / scale;
577 const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
578 const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
580 // calculate L for first point in image (0, 0)
581 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
583 double* deltaFilteredProj = NULL;
584 CubicPolyInterpolator* pCubicInterp = NULL;
585 if (interpType == Backprojector::INTERP_LINEAR) {
586 // precalculate scaled difference for linear interpolation
587 deltaFilteredProj = new double [nDet];
588 for (int i = 0; i < nDet - 1; i++)
589 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
590 deltaFilteredProj[nDet - 1] = 0; // last detector
591 } else if (interpType == Backprojector::INTERP_CUBIC) {
592 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
595 int iLastDet = nDet - 1;
596 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
597 kint32 curDetPos = detPosColStart;
598 ImageFileColumn pImCol = v[ix];
600 if (interpType == Backprojector::INTERP_NEAREST) {
601 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
602 const int iDetPos = (curDetPos + halfScale) >> 16;
603 if (iDetPos >= 0 && iDetPos <= iLastDet)
604 *pImCol++ += filteredProj[iDetPos];
606 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
607 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
608 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
609 if (iDetPos >= 0 && iDetPos <= iLastDet)
610 *pImCol++ += filteredProj[iDetPos];
612 } else if (interpType == Backprojector::INTERP_LINEAR) {
613 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
614 const kint32 iDetPos = curDetPos >> scaleShift;
615 const kint32 detRemainder = curDetPos & scaleBitmask;
616 if (iDetPos >= 0 && iDetPos <= iLastDet)
617 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
619 } else if (interpType == Backprojector::INTERP_CUBIC) {
620 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
621 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
626 if (interpType == Backprojector::INTERP_LINEAR)
627 delete deltaFilteredProj;
628 else if (interpType == Backprojector::INTERP_CUBIC)
634 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
636 double beta = view_angle;
638 CubicPolyInterpolator* pCubicInterp = NULL;
639 if (interpType == Backprojector::INTERP_CUBIC)
640 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
642 for (int ix = 0; ix < nx; ix++) {
643 ImageFileColumn pImCol = v[ix];
645 for (int iy = 0; iy < ny; iy++) {
646 double dAngleDiff = beta - phi[ix][iy];
647 double rcos_t = r[ix][iy] * cos (dAngleDiff);
648 double rsin_t = r[ix][iy] * sin (dAngleDiff);
649 double dFLPlusSin = m_dFocalLength + rsin_t;
650 double gamma = atan (rcos_t / dFLPlusSin);
651 double dPos = gamma / detInc; // position along detector
652 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
654 if (interpType == Backprojector::INTERP_NEAREST) {
655 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
656 if (iDetPos >= 0 && iDetPos < nDet)
657 pImCol[iy] += filteredProj[iDetPos] / dL2;
658 } else if (interpType == Backprojector::INTERP_LINEAR) {
659 double dPosFloor = floor (dPos);
660 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
661 double frac = dPos - dPosFloor; // fraction distance from det
662 if (iDetPos >= 0 && iDetPos < nDet - 1)
663 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
664 } else if (interpType == Backprojector::INTERP_CUBIC) {
665 double d = iDetCenter + dPos; // position along detector
666 if (d >= 0 && d < nDet)
667 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
672 if (interpType == Backprojector::INTERP_CUBIC)
677 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
679 double beta = view_angle;
681 CubicPolyInterpolator* pCubicInterp = NULL;
682 if (interpType == Backprojector::INTERP_CUBIC)
683 pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
685 for (int ix = 0; ix < nx; ix++) {
686 ImageFileColumn pImCol = v[ix];
688 for (int iy = 0; iy < ny; iy++) {
689 double dAngleDiff = beta - phi[ix][iy];
690 double rcos_t = r[ix][iy] * cos (dAngleDiff);
691 double rsin_t = r[ix][iy] * sin (dAngleDiff);
693 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
694 double dDetPos = rcos_t / dU;
695 // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
696 dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
697 double dPos = dDetPos / detInc; // position along detector array
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] / (dU * dU));
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]))
710 } else if (interpType == Backprojector::INTERP_CUBIC) {
711 double d = iDetCenter + dPos; // position along detector
712 if (d >= 0 && d < nDet)
713 pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
718 if (interpType == Backprojector::INTERP_CUBIC)