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.1 2000/06/19 02:59:34 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 ******************************************************************************/
30 // FUNCTION IDENTIFICATION
31 // Backproject* projector = selectBackprojector (...)
34 // Selects a backprojector based on BackprojType
35 // and initializes the backprojector
37 Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
39 Backproject* bj = NULL;
41 if (bjType == O_BPROJ_TRIG)
42 bj = static_cast<Backproject*>(new BackprojectTrig (proj, im, interpType));
43 else if (bjType == O_BPROJ_TABLE)
44 bj = static_cast<Backproject*>(new BackprojectTable (proj, im, interpType));
45 else if (bjType == O_BPROJ_DIFF)
46 bj = static_cast<Backproject*>(new BackprojectDiff (proj, im, interpType));
47 else if (bjType == O_BPROJ_DIFF2)
48 bj = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, interpType));
49 else if (bjType == O_BPROJ_IDIFF2)
50 bj = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, interpType));
52 sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
62 // Pure virtual base class for all backprojectors.
64 Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType)
65 : proj(proj), im(im), interpType(interpType)
67 detInc = proj.detInc();
69 iDetCenter = (nDet - 1) / 2; // index refering to L=0 projection
70 rotInc = proj.rotInc();
77 xMin = -proj.phmLen() / 2; // Retangular coords of phantom
78 xMax = xMin + proj.phmLen();
79 yMin = -proj.phmLen() / 2;
80 yMax = yMin + proj.phmLen();
82 xInc = (xMax - xMin) / nx; // size of cells
83 yInc = (yMax - yMin) / ny;
85 if (interpType != I_NEAREST && interpType != I_LINEAR)
86 sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
89 Backproject::~Backproject (void)
93 Backproject::ScaleImageByRotIncrement (void)
95 for (int ix = 0; ix < nx; ix++)
96 for (int iy = 0; iy < ny; iy++)
100 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
102 printf ("r=%f, phi=%f\n", r, phi);
103 errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
106 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
108 printf ("ix=%d, iy=%d\n", ix, iy);
109 printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
110 printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
111 printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
112 printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
113 sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
117 // CLASS IDENTICATION
121 // Uses trigometric functions at each point in image for backprojection.
124 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
126 double theta = HALFPI + view_angle; // Add PI/2 to get perpendicular angle to detector
128 double x, y; // Rectang coords of center of pixel
130 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
131 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
132 double r = sqrt (x * x + y * y); // distance of cell from center
133 double phi = atan2 (y, x); // angle of cell from center
134 double L = r * cos (theta - phi); // position on detector
136 if (interpType == I_NEAREST) {
137 int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
139 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
140 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
142 v[ix][iy] += rotInc * filteredProj[iDetPos];
143 } else if (interpType == I_LINEAR) {
144 double p = L / detInc; // position along detector
145 double pFloor = floor (p);
146 int iDetPos = iDetCenter + static_cast<int>(pFloor);
147 double frac = p - pFloor; // fraction distance from det
148 if (iDetPos < 0 || iDetPos >= nDet - 1) // check for impossible: index outside of raysum pos
149 errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
151 v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
157 // CLASS IDENTICATION
161 // Precalculates trigometric function value for each point in image for backprojection.
163 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType)
164 : Backproject::Backproject (proj, im, interpType)
166 arrayR.initSetSize (nx, ny);
167 arrayPhi.initSetSize (nx, ny);
168 r = arrayR.getArray();
169 phi = arrayPhi.getArray();
171 double x, y; // Rectang coords of center of pixel
173 for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
174 for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
175 r[ix][iy] = sqrt (x * x + y * y);
176 phi[ix][iy] = atan2 (y, x);
180 BackprojectTable::~BackprojectTable (void)
182 ScaleImageByRotIncrement();
186 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
188 double theta = HALFPI + view_angle; // add half PI to view angle to get perpendicular theta angle
190 for (int ix = 0; ix < nx; ix++) {
191 ImageFileColumn pImCol = v[ix];
193 for (int iy = 0; iy < ny; iy++) {
194 double L = r[ix][iy] * cos (theta - phi[ix][iy]);
196 if (interpType == I_NEAREST) {
197 int iDetPos = iDetCenter + nearest<int>(L / detInc); // calc index in the filtered raysum vector
199 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
200 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
202 pImCol[iy] += filteredProj[iDetPos];
203 } else if (interpType == I_LINEAR) {
204 double dPos = L / detInc; // position along detector
205 double dPosFloor = floor (dPos);
206 int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
207 double frac = dPos - dPosFloor; // fraction distance from det
208 if (iDetPos < 0 || iDetPos >= nDet - 1)
209 errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
211 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
218 // CLASS IDENTICATION
222 // Backprojects by precalculating the change in L position for each x & y step in the image.
223 // Iterates in x & y direction by adding difference in L position
225 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType)
226 : Backproject::Backproject (proj, im, interpType)
228 // calculate center of first pixel v[0][0]
229 double x = xMin + xInc / 2;
230 double y = yMin + yInc / 2;
231 start_r = sqrt (x * x + y * y);
232 start_phi = atan2 (y, x);
237 BackprojectDiff::~BackprojectDiff()
239 ScaleImageByRotIncrement();
243 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
245 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
246 double det_dx = xInc * sin (theta);
247 double det_dy = yInc * cos (theta);
248 double lColStart = start_r * cos (theta - start_phi); // calculate L for first point in image
250 for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
251 double curDetPos = lColStart;
252 ImageFileColumn pImCol = v[ix];
254 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
256 printf ("[%2d,%2d]: %8.5lf ", ix, iy, curDetPos);
258 if (interpType == I_NEAREST) {
259 int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc); // calc index in the filtered raysum vector
261 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
262 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
264 pImCol[iy] += filteredProj[iDetPos];
265 } else if (interpType == I_LINEAR) {
266 double detPos = curDetPos / detInc; // position along detector
267 double detPosFloor = floor (detPos);
268 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
269 double frac = detPos - detPosFloor; // fraction distance from det
270 if (iDetPos < 0 || iDetPos >= nDet - 1)
271 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
273 pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
280 // CLASS IDENTICATION
284 // Optimized version of BackprojectDiff
287 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
289 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
291 // Distance betw. detectors for an angle given in units of detectors
292 double det_dx = xInc * sin (theta) / detInc;
293 double det_dy = yInc * cos (theta) / detInc;
295 // calculate detPosition for first point in image (ix=0, iy=0)
296 double detPosColStart = start_r * cos (theta - start_phi) / detInc;
299 printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
301 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
302 double curDetPos = detPosColStart;
303 ImageFileColumn pImCol = v[ix];
305 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
307 printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
309 if (interpType == I_NEAREST) {
310 int iDetPos = iDetCenter + nearest<int> (curDetPos); // calc index in the filtered raysum vector
312 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
313 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
315 *pImCol++ += filteredProj[iDetPos];
316 } else if (interpType == I_LINEAR) {
317 double detPosFloor = floor (curDetPos);
318 int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
319 double frac = curDetPos - detPosFloor; // fraction distance from det
320 if (iDetPos < 0 || iDetPos >= nDet - 1)
321 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
323 *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
329 // CLASS IDENTICATION
330 // BackprojectIntDiff2
333 // Integer version of BackprojectDiff2
336 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
338 double theta = - view_angle; // add half PI to view angle to get perpendicular theta angle
341 long int scale = 1 << 32;
343 long int scale = 1 << 16;
345 double dScale = scale;
346 long int halfScale = scale / 2;
348 long int det_dx = nearest<long int> (xInc * sin (theta) / detInc * scale);
349 long int det_dy = nearest<long int> (yInc * cos (theta) / detInc * scale);
351 // calculate L for first point in image (0, 0)
352 long int detPosColStart = nearest<long int> (start_r * cos (theta - start_phi) / detInc * scale);
354 for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
355 long int curDetPos = detPosColStart;
356 ImageFileColumn pImCol = v[ix];
358 for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
359 if (interpType == I_NEAREST) {
360 int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
361 int iDetPos = iDetCenter + detPosNearest; // calc index in the filtered raysum vector
363 if (iDetPos < 0 || iDetPos >= nDet) // check for impossible: index outside of raysum pos
364 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
366 *pImCol++ += filteredProj[iDetPos];
367 } else if (interpType == I_LINEAR) {
368 long int detPosFloor = curDetPos / scale;
369 long int detPosRemainder = curDetPos % scale;
370 if (detPosRemainder < 0) {
372 detPosRemainder += scale;
374 int iDetPos = iDetCenter + detPosFloor;
375 double frac = detPosRemainder / dScale;
376 if (iDetPos < 0 || iDetPos >= nDet - 1)
377 errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
379 *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);