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.24 2001/01/27 21:02:20 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_DIFF2 = 3;
34 const int Backprojector::BPROJ_IDIFF2 = 4;
35 const int Backprojector::BPROJ_IDIFF3 = 5;
37 const char* Backprojector::s_aszBackprojectName[] =
47 const char* Backprojector::s_aszBackprojectTitle[] =
49 {"Direct Trigometric"},
50 {"Trigometric Table"},
51 {"Difference Iteration"},
52 {"Difference Iteration Optimized"},
53 {"Integer Difference Iteration Optimized"},
54 {"Integer Difference Iteration Highly-Optimized"},
57 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
59 const int Backprojector::INTERP_INVALID = -1;
60 const int Backprojector::INTERP_NEAREST = 0;
61 const int Backprojector::INTERP_LINEAR = 1;
62 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 2;
63 #if HAVE_BSPLINE_INTERP
64 const int Backprojector::INTERP_BSPLINE = 3;
65 const int Backprojector::INTERP_1BSPLINE = 4;
66 const int Backprojector::INTERP_2BSPLINE = 5;
67 const int Backprojector::INTERP_3BSPLINE = 6;
70 const char* Backprojector::s_aszInterpName[] =
74 {"freq_preinterpolationj"},
75 #if HAVE_BSPLINE_INTERP
83 const char* Backprojector::s_aszInterpTitle[] =
87 {"Frequency Preinterpolation"},
88 #if HAVE_BSPLINE_INTERP
90 {"B-Spline 1st Order"},
91 {"B-Spline 2nd Order"},
92 {"B-Spline 3rd Order"},
96 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
100 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
103 m_pBackprojectImplem = NULL;
105 initBackprojector (proj, im, backprojName, interpName, interpFactor);
109 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
111 if (m_pBackprojectImplem != NULL)
112 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
115 Backprojector::~Backprojector ()
117 delete m_pBackprojectImplem;
120 // FUNCTION IDENTIFICATION
121 // Backproject* projector = selectBackprojector (...)
124 // Selects a backprojector based on BackprojType
125 // and initializes the backprojector
128 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
130 m_nameBackproject = backprojName;
131 m_nameInterpolation = interpName;
132 m_pBackprojectImplem = NULL;
133 m_idBackproject = convertBackprojectNameToID (backprojName);
134 if (m_idBackproject == BPROJ_INVALID) {
136 m_failMessage = "Invalid backprojection name ";
137 m_failMessage += backprojName;
139 m_idInterpolation = convertInterpNameToID (interpName);
140 if (m_idInterpolation == INTERP_INVALID) {
142 m_failMessage = "Invalid interpolation name ";
143 m_failMessage += interpName;
146 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
151 if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
152 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
153 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
154 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
155 else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
156 if (m_idBackproject == BPROJ_TRIG)
157 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
158 else if (m_idBackproject == BPROJ_TABLE)
159 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
160 else if (m_idBackproject == BPROJ_DIFF)
161 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
162 else if (m_idBackproject == BPROJ_DIFF2)
163 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
164 else if (m_idBackproject == BPROJ_IDIFF2)
165 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
166 else if (m_idBackproject == BPROJ_IDIFF3)
167 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
170 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
179 Backprojector::convertBackprojectNameToID (const char* const backprojName)
181 int backprojID = BPROJ_INVALID;
183 for (int i = 0; i < s_iBackprojectCount; i++)
184 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
193 Backprojector::convertBackprojectIDToName (int bprojID)
195 static const char *bprojName = "";
197 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
198 return (s_aszBackprojectName[bprojID]);
204 Backprojector::convertBackprojectIDToTitle (const int bprojID)
206 static const char *bprojTitle = "";
208 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
209 return (s_aszBackprojectTitle[bprojID]);
216 Backprojector::convertInterpNameToID (const char* const interpName)
218 int interpID = INTERP_INVALID;
220 for (int i = 0; i < s_iInterpCount; i++)
221 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
230 Backprojector::convertInterpIDToName (const int interpID)
232 static const char *interpName = "";
234 if (interpID >= 0 && interpID < s_iInterpCount)
235 return (s_aszInterpName[interpID]);
241 Backprojector::convertInterpIDToTitle (const int interpID)
243 static const char *interpTitle = "";
245 if (interpID >= 0 && interpID < s_iInterpCount)
246 return (s_aszInterpTitle[interpID]);
248 return (interpTitle);
253 // CLASS IDENTICATION
257 // Pure virtual base class for all backprojectors.
259 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
260 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
262 detInc = proj.detInc();
264 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
265 rotScale = proj.rotInc();
267 if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
268 rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
269 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
270 rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
272 sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
279 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
280 xMax = xMin + proj.phmLen();
281 yMin = -proj.phmLen() / 2;
282 yMax = yMin + proj.phmLen();
284 xInc = (xMax - xMin) / nx; // size of cells
285 yInc = (yMax - yMin) / ny;
287 m_dFocalLength = proj.focalLength();
290 Backproject::~Backproject ()
294 Backproject::ScaleImageByRotIncrement ()
296 for (int ix = 0; ix < nx; ix++)
297 for (int iy = 0; iy < ny; iy++)
298 v[ix][iy] *= rotScale;
301 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
303 sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
304 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
307 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
310 std::ostringstream os;
311 os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
312 os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
313 os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
314 os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
315 os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
317 sys_error (ERR_WARNING, os.str().c_str());
322 // CLASS IDENTICATION
326 // Uses trigometric functions at each point in image for backprojection.
329 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
331 double theta = view_angle;
333 double x = xMin + xInc / 2; // Rectang coords of center of pixel
334 for (int ix = 0; ix < nx; x += xInc, ix++) {
335 double y = yMin + yInc / 2;
336 for (int iy = 0; iy < ny; y += yInc, iy++) {
337 double r = sqrt (x * x + y * y); // distance of cell from center
338 double phi = atan2 (y, x); // angle of cell from center
339 double L = r * cos (theta - phi); // position on detector
341 if (interpType == Backprojector::INTERP_NEAREST) {
342 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
344 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
345 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
347 v[ix][iy] += rotScale * filteredProj[iDetPos];
348 } else if (interpType == Backprojector::INTERP_LINEAR) {
349 double p = L / detInc; // position along detector
350 double pFloor = floor (p);
351 int iDetPos = iDetCenter + static_cast<int>(pFloor);
352 double frac = p - pFloor; // fraction distance from det
353 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
354 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
356 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
363 // CLASS IDENTICATION
367 // Precalculates trigometric function value for each point in image for backprojection.
369 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
370 : Backproject (proj, im, interpType, interpFactor)
372 arrayR.initSetSize (im.nx(), im.ny());
373 arrayPhi.initSetSize (im.nx(), im.ny());
374 r = arrayR.getArray();
375 phi = arrayPhi.getArray();
377 double x, y; // Rectang coords of center of pixel
379 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
380 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
381 r[ix][iy] = sqrt (x * x + y * y);
382 phi[ix][iy] = atan2 (y, x);
386 BackprojectTable::~BackprojectTable ()
388 ScaleImageByRotIncrement();
392 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
394 double theta = view_angle;
396 for (int ix = 0; ix < nx; ix++) {
397 ImageFileColumn pImCol = v[ix];
399 for (int iy = 0; iy < ny; iy++) {
400 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
402 if (interpType == Backprojector::INTERP_NEAREST) {
403 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
405 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
406 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
408 pImCol[iy] += filteredProj[iDetPos];
409 } else if (interpType == Backprojector::INTERP_LINEAR) {
410 double dPos = L / detInc; // position along detector
411 double dPosFloor = floor (dPos);
412 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
413 double frac = dPos - dPosFloor; // fraction distance from det
414 if (iDetPos < 0 || iDetPos >= nDet - 1)
415 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
417 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
424 // CLASS IDENTICATION
428 // Backprojects by precalculating the change in L position for each x & y step in the image.
429 // Iterates in x & y direction by adding difference in L position
431 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
432 : Backproject (proj, im, interpType, interpFactor)
434 // calculate center of first pixel v[0][0]
435 double x = xMin + xInc / 2;
436 double y = yMin + yInc / 2;
437 start_r = sqrt (x * x + y * y);
438 start_phi = atan2 (y, x);
443 BackprojectDiff::~BackprojectDiff()
445 ScaleImageByRotIncrement();
449 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
451 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
452 double det_dx = xInc * cos (theta);
453 double det_dy = yInc * sin (theta);
454 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
456 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
457 double curDetPos = lColStart;
458 ImageFileColumn pImCol = v[ix];
460 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
462 printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
464 if (interpType == Backprojector::INTERP_NEAREST) {
465 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
467 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
468 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
470 pImCol[iy] += filteredProj[iDetPos];
471 } else if (interpType == Backprojector::INTERP_LINEAR) {
472 double detPos = curDetPos / detInc; // position along detector
473 double detPosFloor = floor (detPos);
474 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
475 double frac = detPos - detPosFloor; // fraction distance from det
476 if (iDetPos < 0 || iDetPos >= nDet - 1)
477 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
479 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
486 // CLASS IDENTICATION
490 // Optimized version of BackprojectDiff
493 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
495 double theta = view_angle;
497 // Distance betw. detectors for an angle given in units of detectors
498 double det_dx = xInc * cos (theta) / detInc;
499 double det_dy = yInc * sin (theta) / detInc;
501 // calculate detPosition for first point in image (ix=0, iy=0)
502 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
505 printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
507 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
508 double curDetPos = detPosColStart;
509 ImageFileColumn pImCol = v[ix];
511 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
513 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
515 if (interpType == Backprojector::INTERP_NEAREST) {
516 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
518 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
519 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
521 *pImCol++ += filteredProj[iDetPos];
522 } else if (interpType == Backprojector::INTERP_LINEAR) {
523 double detPosFloor = floor (curDetPos);
524 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
525 double frac = curDetPos - detPosFloor; // fraction distance from det
526 if (iDetPos < 0 || iDetPos >= nDet - 1)
527 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
529 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
535 // CLASS IDENTICATION
536 // BackprojectIntDiff2
539 // Integer version of BackprojectDiff2
542 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
544 double theta = view_angle;
546 static const kint32 scale = 1 << 16;
547 static const double dScale = scale;
548 static const kint32 halfScale = scale / 2;
550 const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
551 const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
553 // calculate L for first point in image (0, 0)
554 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
556 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
557 kint32 curDetPos = detPosColStart;
558 ImageFileColumn pImCol = v[ix];
560 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
561 if (interpType == Backprojector::INTERP_NEAREST) {
562 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
563 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
565 if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos
566 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
568 *pImCol++ += filteredProj[iDetPos];
569 } else if (interpType == Backprojector::INTERP_LINEAR) {
570 kint32 detPosFloor = curDetPos / scale;
571 kint32 detPosRemainder = curDetPos % scale;
572 if (detPosRemainder < 0) {
574 detPosRemainder += scale;
576 int iDetPos = iDetCenter + detPosFloor;
577 double frac = detPosRemainder / dScale;
578 if (iDetPos < 0 || iDetPos >= nDet - 1)
579 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
581 *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
587 // CLASS IDENTICATION
588 // BackprojectIntDiff3
591 // Highly optimized version of BackprojectIntDiff2
594 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
596 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
597 static const int scaleShift = 16;
598 static const kint32 scale = (1 << scaleShift);
599 static const kint32 scaleBitmask = scale - 1;
600 static const kint32 halfScale = scale / 2;
601 static const double dInvScale = 1. / scale;
603 const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
604 const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
606 // calculate L for first point in image (0, 0)
607 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
609 // precalculate scaled difference for linear interpolation
610 double* deltaFilteredProj = new double [nDet];
611 if (interpType == Backprojector::INTERP_LINEAR) {
612 for (int i = 0; i < nDet - 1; i++)
613 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
615 deltaFilteredProj[nDet - 1] = 0; // last detector
617 int iLastDet = nDet - 1;
618 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
619 kint32 curDetPos = detPosColStart;
620 ImageFileColumn pImCol = v[ix];
622 if (interpType == Backprojector::INTERP_NEAREST) {
623 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
624 const int iDetPos = (curDetPos + halfScale) >> 16;
625 if (iDetPos >= 0 && iDetPos <= iLastDet)
626 *pImCol++ += filteredProj[iDetPos];
628 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
629 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
630 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
631 if (iDetPos >= 0 && iDetPos <= iLastDet)
632 *pImCol++ += filteredProj[iDetPos];
634 } else if (interpType == Backprojector::INTERP_LINEAR) {
635 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
636 const kint32 iDetPos = curDetPos >> scaleShift;
637 const kint32 detRemainder = curDetPos & scaleBitmask;
638 if (iDetPos >= 0 && iDetPos <= iLastDet)
639 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
644 delete deltaFilteredProj;
649 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
651 double beta = view_angle;
653 for (int ix = 0; ix < nx; ix++) {
654 ImageFileColumn pImCol = v[ix];
656 for (int iy = 0; iy < ny; iy++) {
657 double dAngleDiff = beta - phi[ix][iy];
658 double rcos_t = r[ix][iy] * cos (dAngleDiff);
659 double rsin_t = r[ix][iy] * sin (dAngleDiff);
660 double dFLPlusSin = m_dFocalLength + rsin_t;
661 double gamma = atan (rcos_t / dFLPlusSin);
662 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
664 if (interpType == Backprojector::INTERP_NEAREST) {
665 int iDetPos =iDetCenter + nearest<int>(gamma / detInc); // calc index in the filtered raysum vector
667 if (iDetPos < 0 || iDetPos >= nDet) { // check for impossible: index outside of raysum pos
668 ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
670 pImCol[iy] += filteredProj[iDetPos] / dL2;
671 } else if (interpType == Backprojector::INTERP_LINEAR) {
672 double dPos = gamma / detInc; // position along detector
673 double dPosFloor = floor (dPos);
674 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
675 double frac = dPos - dPosFloor; // fraction distance from det
676 if (iDetPos < 0 || iDetPos >= nDet - 1) {
677 ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], gamma, iDetPos);
679 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
686 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
688 double beta = view_angle;
690 for (int ix = 0; ix < nx; ix++) {
691 ImageFileColumn pImCol = v[ix];
693 for (int iy = 0; iy < ny; iy++) {
694 double dAngleDiff = beta - phi[ix][iy];
695 double rcos_t = r[ix][iy] * cos (dAngleDiff);
696 double rsin_t = r[ix][iy] * sin (dAngleDiff);
698 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
699 double dDetPos = rcos_t / dU;
700 // double to scale for imaginary detector that passes through origin
701 // of phantom, see Kak-Slaney Figure 3.22
704 if (interpType == Backprojector::INTERP_NEAREST) {
705 int iDetPos = iDetCenter + nearest<int>(dDetPos / detInc); // calc index in the filtered raysum vector
707 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
708 ; /// errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
710 pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
711 } else if (interpType == Backprojector::INTERP_LINEAR) {
712 double dPos = dDetPos / detInc; // position along detector
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 ; // errorIndexOutsideDetector (ix, iy, beta, r[ix][iy], phi[ix][iy], dDetPos, iDetPos);
719 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / (dU * dU);