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-2000 Kevin Rosenberg
11 ** $Id: backprojectors.cpp,v 1.11 2000/07/23 01:49:03 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 Preinterpolationj"},
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 (m_idBackproject == BPROJ_TRIG)
152 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
153 else if (m_idBackproject == BPROJ_TABLE)
154 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
155 else if (m_idBackproject == BPROJ_DIFF)
156 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
157 else if (m_idBackproject == BPROJ_DIFF2)
158 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
159 else if (m_idBackproject == BPROJ_IDIFF2)
160 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
161 else if (m_idBackproject == BPROJ_IDIFF3)
162 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
165 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
174 Backprojector::convertBackprojectNameToID (const char* const backprojName)
176 int backprojID = BPROJ_INVALID;
178 for (int i = 0; i < s_iBackprojectCount; i++)
179 if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
188 Backprojector::convertBackprojectIDToName (int bprojID)
190 static const char *bprojName = "";
192 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
193 return (s_aszBackprojectName[bprojID]);
199 Backprojector::convertBackprojectIDToTitle (const int bprojID)
201 static const char *bprojTitle = "";
203 if (bprojID >= 0 && bprojID < s_iBackprojectCount)
204 return (s_aszBackprojectTitle[bprojID]);
211 Backprojector::convertInterpNameToID (const char* const interpName)
213 int interpID = INTERP_INVALID;
215 for (int i = 0; i < s_iInterpCount; i++)
216 if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
225 Backprojector::convertInterpIDToName (const int interpID)
227 static const char *interpName = "";
229 if (interpID >= 0 && interpID < s_iInterpCount)
230 return (s_aszInterpName[interpID]);
236 Backprojector::convertInterpIDToTitle (const int interpID)
238 static const char *interpTitle = "";
240 if (interpID >= 0 && interpID < s_iInterpCount)
241 return (s_aszInterpTitle[interpID]);
243 return (interpTitle);
248 // CLASS IDENTICATION
252 // Pure virtual base class for all backprojectors.
254 Backproject::Backproject (const Projections& proj, ImageFile& im, const int interpType, const int interpFactor)
255 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
257 detInc = proj.detInc();
259 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
260 rotInc = proj.rotInc();
267 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
268 xMax = xMin + proj.phmLen();
269 yMin = -proj.phmLen() / 2;
270 yMax = yMin + proj.phmLen();
272 xInc = (xMax - xMin) / nx; // size of cells
273 yInc = (yMax - yMin) / ny;
276 Backproject::~Backproject ()
280 Backproject::ScaleImageByRotIncrement ()
282 for (int ix = 0; ix < nx; ix++)
283 for (int iy = 0; iy < ny; iy++)
287 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
289 printf ("r=%f, phi=%f\n", r, phi);
290 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
293 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
295 printf ("ix=%d, iy=%d\n", ix, iy);
296 printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
297 printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
298 printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
299 printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
300 sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
304 // CLASS IDENTICATION
308 // Uses trigometric functions at each point in image for backprojection.
311 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
313 double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
315 double x, y; // Rectang coords of center of pixel
317 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
318 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
319 double r = sqrt (x * x + y * y); // distance of cell from center
320 double phi = atan2 (y, x); // angle of cell from center
321 double L = r * cos (theta - phi); // position on detector
323 if (interpType == Backprojector::INTERP_NEAREST) {
324 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
326 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
327 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
329 v[ix][iy] += rotInc * filteredProj[iDetPos];
330 } else if (interpType == Backprojector::INTERP_LINEAR) {
331 double p = L / detInc; // position along detector
332 double pFloor = floor (p);
333 int iDetPos = iDetCenter + static_cast<int>(pFloor);
334 double frac = p - pFloor; // fraction distance from det
335 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
336 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
338 v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
344 // CLASS IDENTICATION
348 // Precalculates trigometric function value for each point in image for backprojection.
350 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
351 : Backproject::Backproject (proj, im, interpType, interpFactor)
353 arrayR.initSetSize (nx, ny);
354 arrayPhi.initSetSize (nx, ny);
355 r = arrayR.getArray();
356 phi = arrayPhi.getArray();
358 double x, y; // Rectang coords of center of pixel
360 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
361 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
362 r[ix][iy] = sqrt (x * x + y * y);
363 phi[ix][iy] = atan2 (y, x);
367 BackprojectTable::~BackprojectTable ()
369 ScaleImageByRotIncrement();
373 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
375 double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
377 for (int ix = 0; ix < nx; ix++) {
378 ImageFileColumn pImCol = v[ix];
380 for (int iy = 0; iy < ny; iy++) {
381 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
383 if (interpType == Backprojector::INTERP_NEAREST) {
384 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
386 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
387 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
389 pImCol[iy] += filteredProj[iDetPos];
390 } else if (interpType == Backprojector::INTERP_LINEAR) {
391 double dPos = L / detInc; // position along detector
392 double dPosFloor = floor (dPos);
393 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
394 double frac = dPos - dPosFloor; // fraction distance from det
395 if (iDetPos < 0 || iDetPos >= nDet - 1)
396 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
398 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
405 // CLASS IDENTICATION
409 // Backprojects by precalculating the change in L position for each x & y step in the image.
410 // Iterates in x & y direction by adding difference in L position
412 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
413 : Backproject::Backproject (proj, im, interpType, interpFactor)
415 // calculate center of first pixel v[0][0]
416 double x = xMin + xInc / 2;
417 double y = yMin + yInc / 2;
418 start_r = sqrt (x * x + y * y);
419 start_phi = atan2 (y, x);
424 BackprojectDiff::~BackprojectDiff()
426 ScaleImageByRotIncrement();
430 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
432 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
433 double det_dx = xInc * sin (theta);
434 double det_dy = yInc * cos (theta);
435 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
437 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
438 double curDetPos = lColStart;
439 ImageFileColumn pImCol = v[ix];
441 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
443 printf ("[%2d,%2d]: %8.5f ", ix, iy, curDetPos);
445 if (interpType == Backprojector::INTERP_NEAREST) {
446 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
448 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
449 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
451 pImCol[iy] += filteredProj[iDetPos];
452 } else if (interpType == Backprojector::INTERP_LINEAR) {
453 double detPos = curDetPos / detInc; // position along detector
454 double detPosFloor = floor (detPos);
455 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
456 double frac = detPos - detPosFloor; // fraction distance from det
457 if (iDetPos < 0 || iDetPos >= nDet - 1)
458 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
460 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
467 // CLASS IDENTICATION
471 // Optimized version of BackprojectDiff
474 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
476 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
478 // Distance betw. detectors for an angle given in units of detectors
479 double det_dx = xInc * sin (theta) / detInc;
480 double det_dy = yInc * cos (theta) / detInc;
482 // calculate detPosition for first point in image (ix=0, iy=0)
483 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
486 printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
488 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
489 double curDetPos = detPosColStart;
490 ImageFileColumn pImCol = v[ix];
492 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
494 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
496 if (interpType == Backprojector::INTERP_NEAREST) {
497 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
499 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
500 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
502 *pImCol++ += filteredProj[iDetPos];
503 } else if (interpType == Backprojector::INTERP_LINEAR) {
504 double detPosFloor = floor (curDetPos);
505 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
506 double frac = curDetPos - detPosFloor; // fraction distance from det
507 if (iDetPos < 0 || iDetPos >= nDet - 1)
508 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
510 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
516 // CLASS IDENTICATION
517 // BackprojectIntDiff2
520 // Integer version of BackprojectDiff2
523 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
525 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
527 static const kint32 scale = 1 << 16;
528 static const double dScale = scale;
529 static const kint32 halfScale = scale / 2;
531 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
532 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
534 // calculate L for first point in image (0, 0)
535 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
537 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
538 kint32 curDetPos = detPosColStart;
539 ImageFileColumn pImCol = v[ix];
541 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
542 if (interpType == Backprojector::INTERP_NEAREST) {
543 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
544 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
546 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
547 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
549 *pImCol++ += filteredProj[iDetPos];
550 } else if (interpType == Backprojector::INTERP_LINEAR) {
551 kint32 detPosFloor = curDetPos / scale;
552 kint32 detPosRemainder = curDetPos % scale;
553 if (detPosRemainder < 0) {
555 detPosRemainder += scale;
557 int iDetPos = iDetCenter + detPosFloor;
558 double frac = detPosRemainder / dScale;
559 if (iDetPos < 0 || iDetPos >= nDet - 1)
560 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
562 *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
568 // CLASS IDENTICATION
569 // BackprojectIntDiff3
572 // Highly optimized version of BackprojectIntDiff2
575 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
577 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
578 static const int scaleShift = 16;
579 static const kint32 scale = (1 << scaleShift);
580 static const kint32 scaleBitmask = scale - 1;
581 static const kint32 halfScale = scale / 2;
582 static const double dInvScale = 1. / scale;
584 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
585 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
587 // calculate L for first point in image (0, 0)
588 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
590 // precalculate scaled difference for linear interpolation
591 double deltaFilteredProj [nDet - 1];
592 if (interpType == Backprojector::INTERP_LINEAR) {
593 for (int i = 0; i < nDet - 1; i++)
594 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
597 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
598 kint32 curDetPos = detPosColStart;
599 ImageFileColumn pImCol = v[ix];
601 if (interpType == Backprojector::INTERP_NEAREST) {
602 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
603 const int iDetPos = (curDetPos + halfScale) >> 16;
604 assert(iDetPos >= 0 && iDetPos < nDet);
605 *pImCol++ += filteredProj[iDetPos];
607 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
608 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
609 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
610 assert(iDetPos >= 0 && iDetPos < nDet);
611 *pImCol++ += filteredProj[iDetPos];
613 } else if (interpType == Backprojector::INTERP_LINEAR) {
614 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
615 const kint32 iDetPos = curDetPos >> scaleShift;
616 const kint32 detRemainder = curDetPos & scaleBitmask;
617 assert(iDetPos >= 0 && iDetPos < nDet - 1);
618 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);