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.3 2000/06/22 10:17:28 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)
32 m_pBackprojectImplem = NULL;
34 initBackprojector (proj, im, backprojName, interpName);
38 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
40 if (m_pBackprojectImplem)
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)
59 m_nameBackproject = backprojName;
60 m_nameInterpolation = interpName;
61 m_idBackproject = convertBackprojectNameToID (backprojName);
62 m_idInterpolation = convertInterpolationNameToID (interpName);
63 m_pBackprojectImplem = NULL;
65 if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
70 if (m_idBackproject == BPROJ_TRIG)
71 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
72 else if (m_idBackproject == BPROJ_TABLE)
73 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
74 else if (m_idBackproject == BPROJ_DIFF)
75 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
76 else if (m_idBackproject == BPROJ_DIFF2)
77 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
78 else if (m_idBackproject == BPROJ_IDIFF2)
79 m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
89 const Backprojector::BackprojectID
90 Backprojector::convertBackprojectNameToID (const char* const backprojName)
92 BackprojectID backprojID = BPROJ_INVALID;
94 if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
95 backprojID = BPROJ_TRIG;
96 else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
97 backprojID = BPROJ_TABLE;
98 else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
99 backprojID = BPROJ_DIFF;
100 else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
101 backprojID = BPROJ_DIFF2;
102 else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
103 backprojID = BPROJ_IDIFF2;
109 Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
111 const char *bprojName = "";
113 if (bprojID == BPROJ_TRIG)
114 bprojName = BPROJ_TRIG_STR;
115 else if (bprojID == BPROJ_TABLE)
116 bprojName = BPROJ_TABLE_STR;
117 else if (bprojID == BPROJ_DIFF)
118 bprojName = BPROJ_DIFF_STR;
119 else if (bprojID == BPROJ_DIFF2)
120 bprojName = BPROJ_DIFF2_STR;
121 else if (bprojID == BPROJ_IDIFF2)
122 bprojName = BPROJ_IDIFF2_STR;
129 const Backprojector::InterpolationID
130 Backprojector::convertInterpolationNameToID (const char* const interpName)
132 InterpolationID interpID = INTERP_INVALID;
134 if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
135 interpID = INTERP_NEAREST;
136 else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
137 interpID = INTERP_LINEAR;
138 #if HAVE_BSPLINE_INTERP
139 else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
140 interpID = INTERP_BSPLINE;
148 * name_of_interp Return name of interpolation method
151 * name = name_of_interp (interp_type)
152 * char *name Name of interpolation method
153 * int interp_type Method of interpolation
156 * Returns NULL if interp_type is invalid
160 Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
162 if (interpID == INTERP_NEAREST)
163 return (INTERP_NEAREST_STR);
164 else if (interpID == INTERP_LINEAR)
165 return (INTERP_LINEAR_STR);
166 #if HAVE_BSPLINE_INTERP
167 else if (interpID == INTERP_BSPLINE)
168 return (INTERP_BSPLINE_STR);
175 // CLASS IDENTICATION
179 // Pure virtual base class for all backprojectors.
181 Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType)
182 : proj(proj), im(im), interpType(interpType)
184 detInc = proj.detInc();
186 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
187 rotInc = proj.rotInc();
194 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
195 xMax = xMin + proj.phmLen();
196 yMin = -proj.phmLen() / 2;
197 yMax = yMin + proj.phmLen();
199 xInc = (xMax - xMin) / nx; // size of cells
200 yInc = (yMax - yMin) / ny;
202 if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR)
203 sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
206 Backproject::~Backproject (void)
210 Backproject::ScaleImageByRotIncrement (void)
212 for (int ix = 0; ix < nx; ix++)
213 for (int iy = 0; iy < ny; iy++)
217 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
219 printf ("r=%f, phi=%f\n", r, phi);
220 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
223 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
225 printf ("ix=%d, iy=%d\n", ix, iy);
226 printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
227 printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
228 printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
229 printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
230 sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
234 // CLASS IDENTICATION
238 // Uses trigometric functions at each point in image for backprojection.
241 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
243 double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
245 double x, y; // Rectang coords of center of pixel
247 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
248 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
249 double r = sqrt (x * x + y * y); // distance of cell from center
250 double phi = atan2 (y, x); // angle of cell from center
251 double L = r * cos (theta - phi); // position on detector
253 if (interpType == Backprojector::INTERP_NEAREST) {
254 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
256 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
257 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
259 v[ix][iy] += rotInc * filteredProj[iDetPos];
260 } else if (interpType == Backprojector::INTERP_LINEAR) {
261 double p = L / detInc; // position along detector
262 double pFloor = floor (p);
263 int iDetPos = iDetCenter + static_cast<int>(pFloor);
264 double frac = p - pFloor; // fraction distance from det
265 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
266 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
268 v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
274 // CLASS IDENTICATION
278 // Precalculates trigometric function value for each point in image for backprojection.
280 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
281 : Backproject::Backproject (proj, im, interpType)
283 arrayR.initSetSize (nx, ny);
284 arrayPhi.initSetSize (nx, ny);
285 r = arrayR.getArray();
286 phi = arrayPhi.getArray();
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 r[ix][iy] = sqrt (x * x + y * y);
293 phi[ix][iy] = atan2 (y, x);
297 BackprojectTable::~BackprojectTable (void)
299 ScaleImageByRotIncrement();
303 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
305 double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
307 for (int ix = 0; ix < nx; ix++) {
308 ImageFileColumn pImCol = v[ix];
310 for (int iy = 0; iy < ny; iy++) {
311 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
313 if (interpType == Backprojector::INTERP_NEAREST) {
314 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
316 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
317 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
319 pImCol[iy] += filteredProj[iDetPos];
320 } else if (interpType == Backprojector::INTERP_LINEAR) {
321 double dPos = L / detInc; // position along detector
322 double dPosFloor = floor (dPos);
323 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
324 double frac = dPos - dPosFloor; // fraction distance from det
325 if (iDetPos < 0 || iDetPos >= nDet - 1)
326 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
328 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
335 // CLASS IDENTICATION
339 // Backprojects by precalculating the change in L position for each x & y step in the image.
340 // Iterates in x & y direction by adding difference in L position
342 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
343 : Backproject::Backproject (proj, im, interpType)
345 // calculate center of first pixel v[0][0]
346 double x = xMin + xInc / 2;
347 double y = yMin + yInc / 2;
348 start_r = sqrt (x * x + y * y);
349 start_phi = atan2 (y, x);
354 BackprojectDiff::~BackprojectDiff()
356 ScaleImageByRotIncrement();
360 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
362 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
363 double det_dx = xInc * sin (theta);
364 double det_dy = yInc * cos (theta);
365 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
367 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
368 double curDetPos = lColStart;
369 ImageFileColumn pImCol = v[ix];
371 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
373 printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
375 if (interpType == Backprojector::INTERP_NEAREST) {
376 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
378 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
379 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
381 pImCol[iy] += filteredProj[iDetPos];
382 } else if (interpType == Backprojector::INTERP_LINEAR) {
383 double detPos = curDetPos / detInc; // position along detector
384 double detPosFloor = floor (detPos);
385 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
386 double frac = detPos - detPosFloor; // fraction distance from det
387 if (iDetPos < 0 || iDetPos >= nDet - 1)
388 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
390 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
397 // CLASS IDENTICATION
401 // Optimized version of BackprojectDiff
404 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
406 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
408 // Distance betw. detectors for an angle given in units of detectors
409 double det_dx = xInc * sin (theta) / detInc;
410 double det_dy = yInc * cos (theta) / detInc;
412 // calculate detPosition for first point in image (ix=0, iy=0)
413 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
416 printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
418 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
419 double curDetPos = detPosColStart;
420 ImageFileColumn pImCol = v[ix];
422 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
424 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
426 if (interpType == Backprojector::INTERP_NEAREST) {
427 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
429 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
430 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
432 *pImCol++ += filteredProj[iDetPos];
433 } else if (interpType == Backprojector::INTERP_LINEAR) {
434 double detPosFloor = floor (curDetPos);
435 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
436 double frac = curDetPos - detPosFloor; // fraction distance from det
437 if (iDetPos < 0 || iDetPos >= nDet - 1)
438 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
440 *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
446 // CLASS IDENTICATION
447 // BackprojectIntDiff2
450 // Integer version of BackprojectDiff2
453 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
455 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
457 kint32 scale = 1 << 16;
458 double dScale = scale;
459 kint32 halfScale = scale / 2;
461 kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
462 kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
464 // calculate L for first point in image (0, 0)
465 kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
467 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
468 kint32 curDetPos = detPosColStart;
469 ImageFileColumn pImCol = v[ix];
471 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
472 if (interpType == Backprojector::INTERP_NEAREST) {
473 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
474 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
476 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
477 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
479 *pImCol++ += filteredProj[iDetPos];
480 } else if (interpType == Backprojector::INTERP_LINEAR) {
481 kint32 detPosFloor = curDetPos / scale;
482 kint32 detPosRemainder = curDetPos % scale;
483 if (detPosRemainder < 0) {
485 detPosRemainder += scale;
487 int iDetPos = iDetCenter + detPosFloor;
488 double frac = detPosRemainder / dScale;
489 if (iDetPos < 0 || iDetPos >= nDet - 1)
490 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
492 *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);