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.9 2000/07/20 11:17:31 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 char Backprojector::BPROJ_TRIG_STR[]= "trig";
30 const char Backprojector::BPROJ_TABLE_STR[]= "table";
31 const char Backprojector::BPROJ_DIFF_STR[]= "diff";
32 const char Backprojector::BPROJ_DIFF2_STR[]= "diff2";
33 const char Backprojector::BPROJ_IDIFF2_STR[]= "idiff2";
34 const char Backprojector::BPROJ_IDIFF3_STR[]= "idiff3";
36 const char Backprojector::BPROJ_TRIG_TITLE_STR[]= "Direc Trigometric";
37 const char Backprojector::BPROJ_TABLE_TITLE_STR[]= "Trig Table";
38 const char Backprojector::BPROJ_DIFF_TITLE_STR[]= "Diff";
39 const char Backprojector::BPROJ_DIFF2_TITLE_STR[]= "Diff2";
40 const char Backprojector::BPROJ_IDIFF2_TITLE_STR[]= "Integer Diff2";
41 const char Backprojector::BPROJ_IDIFF3_TITLE_STR[]= "Integer Diff3";
43 const char Backprojector::INTERP_NEAREST_STR[]= "nearest";
44 const char Backprojector::INTERP_LINEAR_STR[]= "linear";
45 const char Backprojector::INTERP_BSPLINE_STR[]= "bspline";
46 const char Backprojector::INTERP_FREQ_PREINTERPOLATION_STR[]= "freq_preinterpolation";
48 const char Backprojector::INTERP_NEAREST_TITLE_STR[]= "Nearest";
49 const char Backprojector::INTERP_LINEAR_TITLE_STR[]= "Linear";
50 const char Backprojector::INTERP_BSPLINE_TITLE_STR[]= "B-Spline";
51 const char Backprojector::INTERP_FREQ_PREINTERPOLATION_TITLE_STR[]= "Frequency Preinterpolation";
54 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
57 m_pBackprojectImplem = NULL;
59 initBackprojector (proj, im, backprojName, interpName, interpFactor);
63 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
65 if (m_pBackprojectImplem != NULL)
66 m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
69 Backprojector::~Backprojector (void)
71 delete m_pBackprojectImplem;
74 // FUNCTION IDENTIFICATION
75 // Backproject* projector = selectBackprojector (...)
78 // Selects a backprojector based on BackprojType
79 // and initializes the backprojector
82 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
84 m_nameBackproject = backprojName;
85 m_nameInterpolation = interpName;
86 m_pBackprojectImplem = NULL;
87 m_idBackproject = convertBackprojectNameToID (backprojName);
88 if (m_idBackproject == BPROJ_INVALID) {
90 m_failMessage = "Invalid backprojection name ";
91 m_failMessage += backprojName;
93 m_idInterpolation = convertInterpolationNameToID (interpName);
94 if (m_idInterpolation == INTERP_INVALID) {
96 m_failMessage = "Invalid interpolation name ";
97 m_failMessage += interpName;
100 if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
105 if (m_idBackproject == BPROJ_TRIG)
106 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
107 else if (m_idBackproject == BPROJ_TABLE)
108 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
109 else if (m_idBackproject == BPROJ_DIFF)
110 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
111 else if (m_idBackproject == BPROJ_DIFF2)
112 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
113 else if (m_idBackproject == BPROJ_IDIFF2)
114 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
115 else if (m_idBackproject == BPROJ_IDIFF3)
116 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
119 m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
127 const Backprojector::BackprojectID
128 Backprojector::convertBackprojectNameToID (const char* const backprojName)
130 BackprojectID backprojID = BPROJ_INVALID;
132 if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
133 backprojID = BPROJ_TRIG;
134 else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
135 backprojID = BPROJ_TABLE;
136 else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
137 backprojID = BPROJ_DIFF;
138 else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
139 backprojID = BPROJ_DIFF2;
140 else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
141 backprojID = BPROJ_IDIFF2;
142 else if (strcasecmp (backprojName, BPROJ_IDIFF3_STR) == 0)
143 backprojID = BPROJ_IDIFF3;
149 Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
151 const char *bprojName = "";
153 if (bprojID == BPROJ_TRIG)
154 bprojName = BPROJ_TRIG_STR;
155 else if (bprojID == BPROJ_TABLE)
156 bprojName = BPROJ_TABLE_STR;
157 else if (bprojID == BPROJ_DIFF)
158 bprojName = BPROJ_DIFF_STR;
159 else if (bprojID == BPROJ_DIFF2)
160 bprojName = BPROJ_DIFF2_STR;
161 else if (bprojID == BPROJ_IDIFF2)
162 bprojName = BPROJ_IDIFF2_STR;
163 else if (bprojID == BPROJ_IDIFF3)
164 bprojName = BPROJ_IDIFF3_STR;
171 const Backprojector::InterpolationID
172 Backprojector::convertInterpolationNameToID (const char* const interpName)
174 InterpolationID interpID = INTERP_INVALID;
176 if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
177 interpID = INTERP_NEAREST;
178 else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
179 interpID = INTERP_LINEAR;
180 else if (strcasecmp (interpName, INTERP_FREQ_PREINTERPOLATION_STR) == 0)
181 interpID = INTERP_FREQ_PREINTERPOLATION;
182 #if HAVE_BSPLINE_INTERP
183 else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
184 interpID = INTERP_BSPLINE;
192 * name_of_interp Return name of interpolation method
195 * name = name_of_interp (interp_type)
196 * char *name Name of interpolation method
197 * int interp_type Method of interpolation
200 * Returns NULL if interp_type is invalid
204 Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
206 if (interpID == INTERP_NEAREST)
207 return (INTERP_NEAREST_STR);
208 else if (interpID == INTERP_LINEAR)
209 return (INTERP_LINEAR_STR);
210 else if (interpID == INTERP_FREQ_PREINTERPOLATION)
211 return (INTERP_FREQ_PREINTERPOLATION_STR);
212 #if HAVE_BSPLINE_INTERP
213 else if (interpID == INTERP_BSPLINE)
214 return (INTERP_BSPLINE_STR);
221 // CLASS IDENTICATION
225 // Pure virtual base class for all backprojectors.
227 Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType, const int interpFactor)
228 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
230 detInc = proj.detInc();
232 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
233 rotInc = proj.rotInc();
240 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
241 xMax = xMin + proj.phmLen();
242 yMin = -proj.phmLen() / 2;
243 yMax = yMin + proj.phmLen();
245 xInc = (xMax - xMin) / nx; // size of cells
246 yInc = (yMax - yMin) / ny;
249 Backproject::~Backproject (void)
253 Backproject::ScaleImageByRotIncrement (void)
255 for (int ix = 0; ix < nx; ix++)
256 for (int iy = 0; iy < ny; iy++)
260 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
262 printf ("r=%f, phi=%f\n", r, phi);
263 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
266 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
268 printf ("ix=%d, iy=%d\n", ix, iy);
269 printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
270 printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
271 printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
272 printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
273 sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
277 // CLASS IDENTICATION
281 // Uses trigometric functions at each point in image for backprojection.
284 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
286 double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
288 double x, y; // Rectang coords of center of pixel
290 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
291 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
292 double r = sqrt (x * x + y * y); // distance of cell from center
293 double phi = atan2 (y, x); // angle of cell from center
294 double L = r * cos (theta - phi); // position on detector
296 if (interpType == Backprojector::INTERP_NEAREST) {
297 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
299 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
300 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
302 v[ix][iy] += rotInc * filteredProj[iDetPos];
303 } else if (interpType == Backprojector::INTERP_LINEAR) {
304 double p = L / detInc; // position along detector
305 double pFloor = floor (p);
306 int iDetPos = iDetCenter + static_cast<int>(pFloor);
307 double frac = p - pFloor; // fraction distance from det
308 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
309 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
311 v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
317 // CLASS IDENTICATION
321 // Precalculates trigometric function value for each point in image for backprojection.
323 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
324 : Backproject::Backproject (proj, im, interpType, interpFactor)
326 arrayR.initSetSize (nx, ny);
327 arrayPhi.initSetSize (nx, ny);
328 r = arrayR.getArray();
329 phi = arrayPhi.getArray();
331 double x, y; // Rectang coords of center of pixel
333 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
334 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
335 r[ix][iy] = sqrt (x * x + y * y);
336 phi[ix][iy] = atan2 (y, x);
340 BackprojectTable::~BackprojectTable (void)
342 ScaleImageByRotIncrement();
346 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
348 double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
350 for (int ix = 0; ix < nx; ix++) {
351 ImageFileColumn pImCol = v[ix];
353 for (int iy = 0; iy < ny; iy++) {
354 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
356 if (interpType == Backprojector::INTERP_NEAREST) {
357 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
359 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
360 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
362 pImCol[iy] += filteredProj[iDetPos];
363 } else if (interpType == Backprojector::INTERP_LINEAR) {
364 double dPos = L / detInc; // position along detector
365 double dPosFloor = floor (dPos);
366 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
367 double frac = dPos - dPosFloor; // fraction distance from det
368 if (iDetPos < 0 || iDetPos >= nDet - 1)
369 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
371 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
378 // CLASS IDENTICATION
382 // Backprojects by precalculating the change in L position for each x & y step in the image.
383 // Iterates in x & y direction by adding difference in L position
385 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType, const int interpFactor)
386 : Backproject::Backproject (proj, im, interpType, interpFactor)
388 // calculate center of first pixel v[0][0]
389 double x = xMin + xInc / 2;
390 double y = yMin + yInc / 2;
391 start_r = sqrt (x * x + y * y);
392 start_phi = atan2 (y, x);
397 BackprojectDiff::~BackprojectDiff()
399 ScaleImageByRotIncrement();
403 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
405 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
406 double det_dx = xInc * sin (theta);
407 double det_dy = yInc * cos (theta);
408 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
410 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
411 double curDetPos = lColStart;
412 ImageFileColumn pImCol = v[ix];
414 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
416 printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
418 if (interpType == Backprojector::INTERP_NEAREST) {
419 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
421 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
422 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
424 pImCol[iy] += filteredProj[iDetPos];
425 } else if (interpType == Backprojector::INTERP_LINEAR) {
426 double detPos = curDetPos / detInc; // position along detector
427 double detPosFloor = floor (detPos);
428 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
429 double frac = detPos - detPosFloor; // fraction distance from det
430 if (iDetPos < 0 || iDetPos >= nDet - 1)
431 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
433 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
440 // CLASS IDENTICATION
444 // Optimized version of BackprojectDiff
447 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
449 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
451 // Distance betw. detectors for an angle given in units of detectors
452 double det_dx = xInc * sin (theta) / detInc;
453 double det_dy = yInc * cos (theta) / detInc;
455 // calculate detPosition for first point in image (ix=0, iy=0)
456 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
459 printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
461 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
462 double curDetPos = detPosColStart;
463 ImageFileColumn pImCol = v[ix];
465 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
467 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
469 if (interpType == Backprojector::INTERP_NEAREST) {
470 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
472 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
473 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
475 *pImCol++ += filteredProj[iDetPos];
476 } else if (interpType == Backprojector::INTERP_LINEAR) {
477 double detPosFloor = floor (curDetPos);
478 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
479 double frac = curDetPos - detPosFloor; // fraction distance from det
480 if (iDetPos < 0 || iDetPos >= nDet - 1)
481 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
483 *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
489 // CLASS IDENTICATION
490 // BackprojectIntDiff2
493 // Integer version of BackprojectDiff2
496 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
498 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
500 static const kint32 scale = 1 << 16;
501 static const double dScale = scale;
502 static const kint32 halfScale = scale / 2;
504 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
505 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
507 // calculate L for first point in image (0, 0)
508 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
510 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
511 kint32 curDetPos = detPosColStart;
512 ImageFileColumn pImCol = v[ix];
514 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
515 if (interpType == Backprojector::INTERP_NEAREST) {
516 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
517 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
519 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
520 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
522 *pImCol++ += filteredProj[iDetPos];
523 } else if (interpType == Backprojector::INTERP_LINEAR) {
524 kint32 detPosFloor = curDetPos / scale;
525 kint32 detPosRemainder = curDetPos % scale;
526 if (detPosRemainder < 0) {
528 detPosRemainder += scale;
530 int iDetPos = iDetCenter + detPosFloor;
531 double frac = detPosRemainder / dScale;
532 if (iDetPos < 0 || iDetPos >= nDet - 1)
533 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
535 *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
541 // CLASS IDENTICATION
542 // BackprojectIntDiff3
545 // Highly optimized version of BackprojectIntDiff2
548 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
550 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
551 static const int scaleShift = 16;
552 static const kint32 scale = (1 << scaleShift);
553 static const kint32 scaleBitmask = scale - 1;
554 static const kint32 halfScale = scale / 2;
555 static const double dInvScale = 1. / scale;
557 const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
558 const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
560 // calculate L for first point in image (0, 0)
561 kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
563 // precalculate scaled difference for linear interpolation
564 double deltaFilteredProj [nDet - 1];
565 if (interpType == Backprojector::INTERP_LINEAR) {
566 for (int i = 0; i < nDet - 1; i++)
567 deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
570 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
571 kint32 curDetPos = detPosColStart;
572 ImageFileColumn pImCol = v[ix];
574 if (interpType == Backprojector::INTERP_NEAREST) {
575 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
576 const int iDetPos = (curDetPos + halfScale) >> 16;
577 assert(iDetPos >= 0 && iDetPos < nDet);
578 *pImCol++ += filteredProj[iDetPos];
580 } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
581 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
582 const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
583 assert(iDetPos >= 0 && iDetPos < nDet);
584 *pImCol++ += filteredProj[iDetPos];
586 } else if (interpType == Backprojector::INTERP_LINEAR) {
587 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
588 const kint32 iDetPos = curDetPos >> scaleShift;
589 const kint32 detRemainder = curDetPos & scaleBitmask;
590 assert(iDetPos >= 0 && iDetPos < nDet - 1);
591 *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);