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.25 2001/02/09 01:54: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_CUBIC = 2;
63 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
64 #if HAVE_BSPLINE_INTERP
65 const int Backprojector::INTERP_BSPLINE = 4;
66 const int Backprojector::INTERP_1BSPLINE = 5;
67 const int Backprojector::INTERP_2BSPLINE = 6;
68 const int Backprojector::INTERP_3BSPLINE = 7;
71 const char* Backprojector::s_aszInterpName[] =
76 {"freq_preinterpolationj"},
77 #if HAVE_BSPLINE_INTERP
85 const char* Backprojector::s_aszInterpTitle[] =
90 {"Frequency Preinterpolation"},
91 #if HAVE_BSPLINE_INTERP
93 {"B-Spline 1st Order"},
94 {"B-Spline 2nd Order"},
95 {"B-Spline 3rd Order"},
99 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
103 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
106 m_pBackprojectImplem = NULL;
108 initBackprojector (proj, im, backprojName, interpName, interpFactor);
112 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
114 if (m_pBackprojectImplem != NULL)
115 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
118 Backprojector::~Backprojector ()
120 delete m_pBackprojectImplem;
123 // FUNCTION IDENTIFICATION
124 // Backproject* projector = selectBackprojector (...)
127 // Selects a backprojector based on BackprojType
128 // and initializes the backprojector
131 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
133 m_nameBackproject = backprojName;
134 m_nameInterpolation = interpName;
135 m_pBackprojectImplem = NULL;
136 m_idBackproject = convertBackprojectNameToID (backprojName);
137 if (m_idBackproject == BPROJ_INVALID) {
139 m_failMessage = "Invalid backprojection name ";
140 m_failMessage += backprojName;
142 m_idInterpolation = convertInterpNameToID (interpName);
143 if (m_idInterpolation == INTERP_INVALID) {
145 m_failMessage = "Invalid interpolation name ";
146 m_failMessage += interpName;
149 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
154 if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
155 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
156 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
157 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
158 else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
159 if (m_idBackproject == BPROJ_TRIG)
160 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
161 else if (m_idBackproject == BPROJ_TABLE)
162 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
163 else if (m_idBackproject == BPROJ_DIFF)
164 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
165 else if (m_idBackproject == BPROJ_DIFF2)
166 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
167 else if (m_idBackproject == BPROJ_IDIFF2)
168 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
169 else if (m_idBackproject == BPROJ_IDIFF3)
170 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
173 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
182 Backprojector::convertBackprojectNameToID (const char* const backprojName)
184 int backprojID = BPROJ_INVALID;
186 for (int i = 0; i < s_iBackprojectCount; i++)
187 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
196 Backprojector::convertBackprojectIDToName (int bprojID)
198 static const char *bprojName = "";
200 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
201 return (s_aszBackprojectName[bprojID]);
207 Backprojector::convertBackprojectIDToTitle (const int bprojID)
209 static const char *bprojTitle = "";
211 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
212 return (s_aszBackprojectTitle[bprojID]);
219 Backprojector::convertInterpNameToID (const char* const interpName)
221 int interpID = INTERP_INVALID;
223 for (int i = 0; i < s_iInterpCount; i++)
224 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
233 Backprojector::convertInterpIDToName (const int interpID)
235 static const char *interpName = "";
237 if (interpID >= 0 && interpID < s_iInterpCount)
238 return (s_aszInterpName[interpID]);
244 Backprojector::convertInterpIDToTitle (const int interpID)
246 static const char *interpTitle = "";
248 if (interpID >= 0 && interpID < s_iInterpCount)
249 return (s_aszInterpTitle[interpID]);
251 return (interpTitle);
256 // CLASS IDENTICATION
260 // Pure virtual base class for all backprojectors.
262 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
263 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
265 detInc = proj.detInc();
267 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
268 rotScale = proj.rotInc();
270 if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
271 rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
272 else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
273 rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
275 sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
282 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
283 xMax = xMin + proj.phmLen();
284 yMin = -proj.phmLen() / 2;
285 yMax = yMin + proj.phmLen();
287 xInc = (xMax - xMin) / nx; // size of cells
288 yInc = (yMax - yMin) / ny;
290 m_dFocalLength = proj.focalLength();
293 Backproject::~Backproject ()
297 Backproject::ScaleImageByRotIncrement ()
299 for (int ix = 0; ix < nx; ix++)
300 for (int iy = 0; iy < ny; iy++)
301 v[ix][iy] *= rotScale;
304 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
306 sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
307 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
310 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
313 std::ostringstream os;
314 os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
315 os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
316 os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
317 os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
318 os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
320 sys_error (ERR_WARNING, os.str().c_str());
325 // CLASS IDENTICATION
329 // Uses trigometric functions at each point in image for backprojection.
332 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
334 double theta = view_angle;
336 CubicInterpolator* pCubicInterp = NULL;
337 if (interpType == Backprojector::INTERP_CUBIC)
338 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
340 double x = xMin + xInc / 2; // Rectang coords of center of pixel
341 for (int ix = 0; ix < nx; x += xInc, ix++) {
342 double y = yMin + yInc / 2;
343 for (int iy = 0; iy < ny; y += yInc, iy++) {
344 double r = sqrt (x * x + y * y); // distance of cell from center
345 double phi = atan2 (y, x); // angle of cell from center
346 double L = r * cos (theta - phi); // position on detector
348 if (interpType == Backprojector::INTERP_NEAREST) {
349 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
351 if (iDetPos >= 0 && iDetPos < nDet)
352 v[ix][iy] += rotScale * filteredProj[iDetPos];
353 } else if (interpType == Backprojector::INTERP_LINEAR) {
354 double p = L / detInc; // position along detector
355 double pFloor = floor (p);
356 int iDetPos = iDetCenter + static_cast<int>(pFloor);
357 double frac = p - pFloor; // fraction distance from det
358 if (iDetPos >= 0 && iDetPos < nDet - 1)
359 v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
360 } else if (interpType = Backprojector::INTERP_CUBIC) {
361 double p = iDetCenter + (L / detInc); // position along detector
362 if (p >= 0 && p < nDet)
363 v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
368 if (interpType == Backprojector::INTERP_CUBIC)
373 // CLASS IDENTICATION
377 // Precalculates trigometric function value for each point in image for backprojection.
379 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
380 : Backproject (proj, im, interpType, interpFactor)
382 arrayR.initSetSize (im.nx(), im.ny());
383 arrayPhi.initSetSize (im.nx(), im.ny());
384 r = arrayR.getArray();
385 phi = arrayPhi.getArray();
387 double x, y; // Rectang coords of center of pixel
389 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
390 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
391 r[ix][iy] = sqrt (x * x + y * y);
392 phi[ix][iy] = atan2 (y, x);
396 BackprojectTable::~BackprojectTable ()
398 ScaleImageByRotIncrement();
402 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
404 double theta = view_angle;
406 CubicInterpolator* pCubicInterp = NULL;
407 if (interpType == Backprojector::INTERP_CUBIC)
408 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
410 for (int ix = 0; ix < nx; ix++) {
411 ImageFileColumn pImCol = v[ix];
413 for (int iy = 0; iy < ny; iy++) {
414 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
416 if (interpType == Backprojector::INTERP_NEAREST) {
417 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
419 if (iDetPos >= 0 && iDetPos < nDet)
420 pImCol[iy] += filteredProj[iDetPos];
421 } else if (interpType == Backprojector::INTERP_LINEAR) {
422 double dPos = L / detInc; // position along detector
423 double dPosFloor = floor (dPos);
424 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
425 double frac = dPos - dPosFloor; // fraction distance from det
426 if (iDetPos >= 0 && iDetPos < nDet - 1)
427 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
428 } else if (interpType = Backprojector::INTERP_CUBIC) {
429 double p = iDetCenter + (L / detInc); // position along detector
430 if (p >= 0 && p < nDet)
431 pImCol[iy] += pCubicInterp->interpolate (p);
436 if (interpType == Backprojector::INTERP_CUBIC)
441 // CLASS IDENTICATION
445 // Backprojects by precalculating the change in L position for each x & y step in the image.
446 // Iterates in x & y direction by adding difference in L position
448 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
449 : Backproject (proj, im, interpType, interpFactor)
451 // calculate center of first pixel v[0][0]
452 double x = xMin + xInc / 2;
453 double y = yMin + yInc / 2;
454 start_r = sqrt (x * x + y * y);
455 start_phi = atan2 (y, x);
460 BackprojectDiff::~BackprojectDiff()
462 ScaleImageByRotIncrement();
466 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
468 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
469 double det_dx = xInc * cos (theta);
470 double det_dy = yInc * sin (theta);
471 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
473 CubicInterpolator* pCubicInterp = NULL;
474 if (interpType == Backprojector::INTERP_CUBIC)
475 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
477 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
478 double curDetPos = lColStart;
479 ImageFileColumn pImCol = v[ix];
481 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
483 printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
485 if (interpType == Backprojector::INTERP_NEAREST) {
486 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
488 if (iDetPos >= 0 && iDetPos < nDet)
489 pImCol[iy] += filteredProj[iDetPos];
490 } else if (interpType == Backprojector::INTERP_LINEAR) {
491 double detPos = curDetPos / detInc; // position along detector
492 double detPosFloor = floor (detPos);
493 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
494 double frac = detPos - detPosFloor; // fraction distance from det
495 if (iDetPos >= 0 && iDetPos < nDet - 1)
496 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
497 } else if (interpType = Backprojector::INTERP_CUBIC) {
498 double p = iDetCenter + (curDetPos / detInc); // position along detector
499 if (p >= 0 && p < nDet)
500 pImCol[iy] += pCubicInterp->interpolate (p);
505 if (interpType == Backprojector::INTERP_CUBIC)
510 // CLASS IDENTICATION
514 // Optimized version of BackprojectDiff
517 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
519 double theta = view_angle;
521 CubicInterpolator* pCubicInterp = NULL;
522 if (interpType == Backprojector::INTERP_CUBIC)
523 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
525 // Distance betw. 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 = start_r * cos (theta - start_phi) / detInc;
533 printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
535 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
536 double curDetPos = detPosColStart;
537 ImageFileColumn pImCol = v[ix];
539 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
541 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
543 if (interpType == Backprojector::INTERP_NEAREST) {
544 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
546 if (iDetPos >= 0 && iDetPos < nDet)
547 *pImCol++ += filteredProj[iDetPos];
548 } else if (interpType == Backprojector::INTERP_LINEAR) {
549 double detPosFloor = floor (curDetPos);
550 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
551 double frac = curDetPos - detPosFloor; // fraction distance from det
552 if (iDetPos > 0 && iDetPos < nDet - 1)
553 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
554 } else if (interpType = Backprojector::INTERP_CUBIC) {
555 double p = iDetCenter + curDetPos; // position along detector
556 if (p >= 0 && p < nDet)
557 *pImCol++ += pCubicInterp->interpolate (p);
562 if (interpType == Backprojector::INTERP_CUBIC)
566 // CLASS IDENTICATION
567 // BackprojectIntDiff2
570 // Integer version of BackprojectDiff2
573 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
575 double theta = view_angle;
577 CubicInterpolator* pCubicInterp = NULL;
578 if (interpType == Backprojector::INTERP_CUBIC)
579 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
581 static const kint32 scale = 1 << 16;
582 static const double dScale = scale;
583 static const kint32 halfScale = scale / 2;
585 const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
586 const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
588 // calculate L for first point in image (0, 0)
589 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
591 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
592 kint32 curDetPos = detPosColStart;
593 ImageFileColumn pImCol = v[ix];
595 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
596 if (interpType == Backprojector::INTERP_NEAREST) {
597 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
598 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
600 if (iDetPos < 0 || iDetPos >= nDet) // check for index outside of raysum pos
601 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
603 *pImCol++ += filteredProj[iDetPos];
604 } else if (interpType == Backprojector::INTERP_LINEAR) {
605 kint32 detPosFloor = curDetPos / scale;
606 kint32 detPosRemainder = curDetPos % scale;
607 if (detPosRemainder < 0) {
609 detPosRemainder += scale;
611 int iDetPos = iDetCenter + detPosFloor;
612 double frac = detPosRemainder / dScale;
613 if (iDetPos < 0 || iDetPos >= nDet - 1)
614 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
616 *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
617 } else if (interpType = Backprojector::INTERP_CUBIC) {
618 double p = iDetCenter + (static_cast<double>(curDetPos) / scale); // position along detector
619 if (p >= 0 && p < nDet)
620 *pImCol++ += pCubicInterp->interpolate (p);
625 if (interpType == Backprojector::INTERP_CUBIC)
629 // CLASS IDENTICATION
630 // BackprojectIntDiff3
633 // Highly optimized version of BackprojectIntDiff2
636 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
638 double theta = view_angle; // add half PI to view angle to get perpendicular theta angle
639 static const int scaleShift = 16;
640 static const kint32 scale = (1 << scaleShift);
641 static const kint32 scaleBitmask = scale - 1;
642 static const kint32 halfScale = scale / 2;
643 static const double dInvScale = 1. / scale;
645 const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
646 const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
648 // calculate L for first point in image (0, 0)
649 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
651 double* deltaFilteredProj = NULL;
652 CubicInterpolator* pCubicInterp = NULL;
653 if (interpType == Backprojector::INTERP_LINEAR) {
654 // precalculate scaled difference for linear interpolation
655 deltaFilteredProj = new double [nDet];
656 for (int i = 0; i < nDet - 1; i++)
657 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
658 deltaFilteredProj[nDet - 1] = 0; // last detector
659 } else if (interpType == Backprojector::INTERP_CUBIC) {
660 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
663 int iLastDet = nDet - 1;
664 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
665 kint32 curDetPos = detPosColStart;
666 ImageFileColumn pImCol = v[ix];
668 if (interpType == Backprojector::INTERP_NEAREST) {
669 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
670 const int iDetPos = (curDetPos + halfScale) >> 16;
671 if (iDetPos >= 0 && iDetPos <= iLastDet)
672 *pImCol++ += filteredProj[iDetPos];
674 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
675 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
676 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
677 if (iDetPos >= 0 && iDetPos <= iLastDet)
678 *pImCol++ += filteredProj[iDetPos];
680 } else if (interpType == Backprojector::INTERP_LINEAR) {
681 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
682 const kint32 iDetPos = curDetPos >> scaleShift;
683 const kint32 detRemainder = curDetPos & scaleBitmask;
684 if (iDetPos >= 0 && iDetPos <= iLastDet)
685 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
687 } else if (interpType = Backprojector::INTERP_CUBIC) {
688 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
689 *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
694 if (interpType == Backprojector::INTERP_LINEAR)
695 delete deltaFilteredProj;
696 else if (interpType == Backprojector::INTERP_CUBIC)
702 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
704 double beta = view_angle;
706 CubicInterpolator* pCubicInterp = NULL;
707 if (interpType == Backprojector::INTERP_CUBIC)
708 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
710 for (int ix = 0; ix < nx; ix++) {
711 ImageFileColumn pImCol = v[ix];
713 for (int iy = 0; iy < ny; iy++) {
714 double dAngleDiff = beta - phi[ix][iy];
715 double rcos_t = r[ix][iy] * cos (dAngleDiff);
716 double rsin_t = r[ix][iy] * sin (dAngleDiff);
717 double dFLPlusSin = m_dFocalLength + rsin_t;
718 double gamma = atan (rcos_t / dFLPlusSin);
719 double dPos = gamma / detInc; // position along detector
720 double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
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] / dL2;
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])) / dL2;
732 } else if (interpType == Backprojector::INTERP_CUBIC) {
733 double d = iDetCenter + dPos; // position along detector
734 if (d >= 0 && d < nDet)
735 pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
740 if (interpType == Backprojector::INTERP_CUBIC)
745 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
747 double beta = view_angle;
749 CubicInterpolator* pCubicInterp = NULL;
750 if (interpType == Backprojector::INTERP_CUBIC)
751 pCubicInterp = new CubicInterpolator (filteredProj, nDet);
753 for (int ix = 0; ix < nx; ix++) {
754 ImageFileColumn pImCol = v[ix];
756 for (int iy = 0; iy < ny; iy++) {
757 double dAngleDiff = beta - phi[ix][iy];
758 double rcos_t = r[ix][iy] * cos (dAngleDiff);
759 double rsin_t = r[ix][iy] * sin (dAngleDiff);
761 double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
762 double dDetPos = rcos_t / dU;
763 // double to scale for imaginary detector that passes through origin
764 // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also
765 // located focal-length away from the origin.
767 double dPos = dDetPos / detInc; // position along detector array
769 if (interpType == Backprojector::INTERP_NEAREST) {
770 int iDetPos = iDetCenter + nearest<int>(dPos); // calc index in the filtered raysum vector
771 if (iDetPos >= 0 && iDetPos < nDet)
772 pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
773 } else if (interpType == Backprojector::INTERP_LINEAR) {
774 double dPosFloor = floor (dPos);
775 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
776 double frac = dPos - dPosFloor; // fraction distance from det
777 if (iDetPos >= 0 && iDetPos < nDet - 1)
778 pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
780 } else if (interpType == Backprojector::INTERP_CUBIC) {
781 double d = iDetCenter + dPos; // position along detector
782 if (d >= 0 && d < nDet)
783 pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
788 if (interpType == Backprojector::INTERP_CUBIC)