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.8 2000/07/13 07:03:21 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 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
32 m_pBackprojectImplem = NULL;
34 initBackprojector (proj, im, backprojName, interpName, interpFactor);
38 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
40 if (m_pBackprojectImplem != NULL)
41 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
44 Backprojector::~Backprojector (void)
46 delete m_pBackprojectImplem;
49 // FUNCTION IDENTIFICATION
50 // Backproject* projector = selectBackprojector (...)
53 // Selects a backprojector based on BackprojType
54 // and initializes the backprojector
57 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
59 m_nameBackproject = backprojName;
60 m_nameInterpolation = interpName;
61 m_pBackprojectImplem = NULL;
62 m_idBackproject = convertBackprojectNameToID (backprojName);
63 if (m_idBackproject == BPROJ_INVALID) {
65 m_failMessage = "Invalid backprojection name ";
66 m_failMessage += backprojName;
68 m_idInterpolation = convertInterpolationNameToID (interpName);
69 if (m_idInterpolation == INTERP_INVALID) {
71 m_failMessage = "Invalid interpolation name ";
72 m_failMessage += interpName;
75 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
80 if (m_idBackproject == BPROJ_TRIG)
81 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
82 else if (m_idBackproject == BPROJ_TABLE)
83 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
84 else if (m_idBackproject == BPROJ_DIFF)
85 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
86 else if (m_idBackproject == BPROJ_DIFF2)
87 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
88 else if (m_idBackproject == BPROJ_IDIFF2)
89 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
90 else if (m_idBackproject == BPROJ_IDIFF3)
91 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
94 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
102 const Backprojector::BackprojectID
103 Backprojector::convertBackprojectNameToID (const char* const backprojName)
105 BackprojectID backprojID = BPROJ_INVALID;
107 if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
108 backprojID = BPROJ_TRIG;
109 else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
110 backprojID = BPROJ_TABLE;
111 else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
112 backprojID = BPROJ_DIFF;
113 else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
114 backprojID = BPROJ_DIFF2;
115 else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
116 backprojID = BPROJ_IDIFF2;
117 else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0)
118 backprojID = BPROJ_IDIFF3;
124 Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
126 const char *bprojName = "";
128 if (bprojID == BPROJ_TRIG)
129 bprojName = BPROJ_TRIG_STR;
130 else if (bprojID == BPROJ_TABLE)
131 bprojName = BPROJ_TABLE_STR;
132 else if (bprojID == BPROJ_DIFF)
133 bprojName = BPROJ_DIFF_STR;
134 else if (bprojID == BPROJ_DIFF2)
135 bprojName = BPROJ_DIFF2_STR;
136 else if (bprojID == BPROJ_IDIFF2)
137 bprojName = BPROJ_IDIFF2_STR;
138 else if (bprojID == BPROJ_IDIFF3)
139 bprojName = BPROJ_IDIFF3_STR;
146 const Backprojector::InterpolationID
147 Backprojector::convertInterpolationNameToID (const char* const interpName)
149 InterpolationID interpID = INTERP_INVALID;
151 if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
152 interpID = INTERP_NEAREST;
153 else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
154 interpID = INTERP_LINEAR;
155 else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0)
156 interpID = INTERP_FREQ_PREINTERPOLATION;
157 #if HAVE_BSPLINE_INTERP
158 else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
159 interpID = INTERP_BSPLINE;
167 * name_of_interp Return name of interpolation method
170 * name = name_of_interp (interp_type)
171 * char *name Name of interpolation method
172 * int interp_type Method of interpolation
175 * Returns NULL if interp_type is invalid
179 Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
181 if (interpID == INTERP_NEAREST)
182 return (INTERP_NEAREST_STR);
183 else if (interpID == INTERP_LINEAR)
184 return (INTERP_LINEAR_STR);
185 else if (interpID == INTERP_FREQ_PREINTERPOLATION)
186 return (INTERP_FREQ_PREINTERPOLATION_STR);
187 #if HAVE_BSPLINE_INTERP
188 else if (interpID == INTERP_BSPLINE)
189 return (INTERP_BSPLINE_STR);
196 // CLASS IDENTICATION
200 // Pure virtual base class for all backprojectors.
202 Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor)
203 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
205 detInc = proj.detInc();
207 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
208 rotInc = proj.rotInc();
215 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
216 xMax = xMin + proj.phmLen();
217 yMin = -proj.phmLen() / 2;
218 yMax = yMin + proj.phmLen();
220 xInc = (xMax - xMin) / nx; // size of cells
221 yInc = (yMax - yMin) / ny;
224 Backproject::~Backproject (void)
228 Backproject::ScaleImageByRotIncrement (void)
230 for (int ix = 0; ix < nx; ix++)
231 for (int iy = 0; iy < ny; iy++)
235 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
237 printf ("r=%f, phi=%f\n", r, phi);
238 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
241 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
243 printf ("ix=%d, iy=%d\n", ix, iy);
244 printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
245 printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
246 printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
247 printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
248 sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
252 // CLASS IDENTICATION
256 // Uses trigometric functions at each point in image for backprojection.
259 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
261 double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
263 double x, y; // Rectang coords of center of pixel
265 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
266 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
267 double r = sqrt (x * x + y * y); // distance of cell from center
268 double phi = atan2 (y, x); // angle of cell from center
269 double L = r * cos (theta - phi); // position on detector
271 if (interpType == Backprojector::INTERP_NEAREST) {
272 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
274 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
275 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
277 v[ix][iy] += rotInc * filteredProj[iDetPos];
278 } else if (interpType == Backprojector::INTERP_LINEAR) {
279 double p = L / detInc; // position along detector
280 double pFloor = floor (p);
281 int iDetPos = iDetCenter + static_cast<int>(pFloor);
282 double frac = p - pFloor; // fraction distance from det
283 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
284 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
286 v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
292 // CLASS IDENTICATION
296 // Precalculates trigometric function value for each point in image for backprojection.
298 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
299 : Backproject::Backproject (proj, im, interpType, interpFactor)
301 arrayR.initSetSize (nx, ny);
302 arrayPhi.initSetSize (nx, ny);
303 r = arrayR.getArray();
304 phi = arrayPhi.getArray();
306 double x, y; // Rectang coords of center of pixel
308 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
309 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
310 r[ix][iy] = sqrt (x * x + y * y);
311 phi[ix][iy] = atan2 (y, x);
315 BackprojectTable::~BackprojectTable (void)
317 ScaleImageByRotIncrement();
321 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
323 double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
325 for (int ix = 0; ix < nx; ix++) {
326 ImageFileColumn pImCol = v[ix];
328 for (int iy = 0; iy < ny; iy++) {
329 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
331 if (interpType == Backprojector::INTERP_NEAREST) {
332 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
334 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
335 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
337 pImCol[iy] += filteredProj[iDetPos];
338 } else if (interpType == Backprojector::INTERP_LINEAR) {
339 double dPos = L / detInc; // position along detector
340 double dPosFloor = floor (dPos);
341 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
342 double frac = dPos - dPosFloor; // fraction distance from det
343 if (iDetPos < 0 || iDetPos >= nDet - 1)
344 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
346 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
353 // CLASS IDENTICATION
357 // Backprojects by precalculating the change in L position for each x & y step in the image.
358 // Iterates in x & y direction by adding difference in L position
360 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
361 : Backproject::Backproject (proj, im, interpType, interpFactor)
363 // calculate center of first pixel v[0][0]
364 double x = xMin + xInc / 2;
365 double y = yMin + yInc / 2;
366 start_r = sqrt (x * x + y * y);
367 start_phi = atan2 (y, x);
372 BackprojectDiff::~BackprojectDiff()
374 ScaleImageByRotIncrement();
378 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
380 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
381 double det_dx = xInc * sin (theta);
382 double det_dy = yInc * cos (theta);
383 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
385 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
386 double curDetPos = lColStart;
387 ImageFileColumn pImCol = v[ix];
389 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
391 printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
393 if (interpType == Backprojector::INTERP_NEAREST) {
394 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
396 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
397 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
399 pImCol[iy] += filteredProj[iDetPos];
400 } else if (interpType == Backprojector::INTERP_LINEAR) {
401 double detPos = curDetPos / detInc; // position along detector
402 double detPosFloor = floor (detPos);
403 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
404 double frac = detPos - detPosFloor; // fraction distance from det
405 if (iDetPos < 0 || iDetPos >= nDet - 1)
406 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
408 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
415 // CLASS IDENTICATION
419 // Optimized version of BackprojectDiff
422 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
424 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
426 // Distance betw. detectors for an angle given in units of detectors
427 double det_dx = xInc * sin (theta) / detInc;
428 double det_dy = yInc * cos (theta) / detInc;
430 // calculate detPosition for first point in image (ix=0, iy=0)
431 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
434 printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
436 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
437 double curDetPos = detPosColStart;
438 ImageFileColumn pImCol = v[ix];
440 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
442 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
444 if (interpType == Backprojector::INTERP_NEAREST) {
445 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
447 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
448 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
450 *pImCol++ += filteredProj[iDetPos];
451 } else if (interpType == Backprojector::INTERP_LINEAR) {
452 double detPosFloor = floor (curDetPos);
453 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
454 double frac = curDetPos - detPosFloor; // fraction distance from det
455 if (iDetPos < 0 || iDetPos >= nDet - 1)
456 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
458 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
464 // CLASS IDENTICATION
465 // BackprojectIntDiff2
468 // Integer version of BackprojectDiff2
471 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
473 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
475 static const kint32 scale = 1 << 16;
476 static const double dScale = scale;
477 static const kint32 halfScale = scale / 2;
479 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
480 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
482 // calculate L for first point in image (0, 0)
483 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
485 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
486 kint32 curDetPos = detPosColStart;
487 ImageFileColumn pImCol = v[ix];
489 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
490 if (interpType == Backprojector::INTERP_NEAREST) {
491 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
492 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
494 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
495 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
497 *pImCol++ += filteredProj[iDetPos];
498 } else if (interpType == Backprojector::INTERP_LINEAR) {
499 kint32 detPosFloor = curDetPos / scale;
500 kint32 detPosRemainder = curDetPos % scale;
501 if (detPosRemainder < 0) {
503 detPosRemainder += scale;
505 int iDetPos = iDetCenter + detPosFloor;
506 double frac = detPosRemainder / dScale;
507 if (iDetPos < 0 || iDetPos >= nDet - 1)
508 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
510 *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
516 // CLASS IDENTICATION
517 // BackprojectIntDiff3
520 // Highly optimized version of BackprojectIntDiff2
523 BackprojectIntDiff3::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
526 static const int scaleShift = 16;
527 static const kint32 scale = (1 << scaleShift);
528 static const kint32 scaleBitmask = scale - 1;
529 static const kint32 halfScale = scale / 2;
530 static const double dInvScale = 1. / scale;
532 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
533 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
535 // calculate L for first point in image (0, 0)
536 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
538 // precalculate scaled difference for linear interpolation
539 double deltaFilteredProj [nDet - 1];
540 if (interpType == Backprojector::INTERP_LINEAR) {
541 for (int i = 0; i < nDet - 1; i++)
542 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
545 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
546 kint32 curDetPos = detPosColStart;
547 ImageFileColumn pImCol = v[ix];
549 if (interpType == Backprojector::INTERP_NEAREST) {
550 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
551 const int iDetPos = (curDetPos + halfScale) >> 16;
552 assert(iDetPos >= 0 && iDetPos < nDet);
553 *pImCol++ += filteredProj[iDetPos];
555 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
556 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
557 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
558 assert(iDetPos >= 0 && iDetPos < nDet);
559 *pImCol++ += filteredProj[iDetPos];
561 } else if (interpType == Backprojector::INTERP_LINEAR) {
562 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
563 const kint32 iDetPos = curDetPos >> scaleShift;
564 const kint32 detRemainder = curDetPos & scaleBitmask;
565 assert(iDetPos >= 0 && iDetPos < nDet - 1);
566 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);