libctsim/backprojectors.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:         backprojectors.cpp         Classes for backprojection
   5 **   Programmer:   Kevin Rosenberg
   6 **   Date Started: June 2000
   7 **
   8 **  This is part of the CTSim program
   9 **  Copyright (C) 1983-2000 Kevin Rosenberg
  10 **
  11 **  $Id: backprojectors.cpp,v 1.2 2000/06/19 19:07:33 kevin Exp $
  12 **
  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.
  16 **
  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.
  21 **
  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 ******************************************************************************/
  26
  27 #include "ct.h"
  28
  29
  30 // FUNCTION IDENTIFICATION
  31 //     Backproject* projector = selectBackprojector (...)
  32 //
  33 // PURPOSE
  34 //     Selects a backprojector based on BackprojType
  35 //     and initializes the backprojector
  36
  37 Backproject* selectBackprojector (BackprojType bjType, const Projections& proj, ImageFile& im, InterpolationType interpType)
  38 {
  39     Backproject* bj = NULL;
  40
  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));
  51     else
  52       sys_error (ERR_WARNING, "Illegal backproject type %d [selectBackprojector]");
  53
  54     return (bj);
  55 }
  56
  57
  58 // CLASS IDENTICATION
  59 //   Backproject
  60 //
  61 // PURPOSE
  62 //   Pure virtual base class for all backprojectors.
  63
  64 Backproject::Backproject (const Projections& proj, ImageFile& im, const InterpolationType interpType)
  65   : proj(proj), im(im), interpType(interpType)
  66 {
  67   detInc = proj.detInc();
  68   nDet = proj.nDet();
  69   iDetCenter = (nDet - 1) / 2;  // index refering to L=0 projection
  70   rotInc = proj.rotInc();
  71
  72   v = im.getArray();
  73   nx = im.nx();
  74   ny = im.ny();
  75   im.arrayDataClear();
  76
  77   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
  78   xMax = xMin + proj.phmLen();
  79   yMin = -proj.phmLen() / 2;
  80   yMax = yMin + proj.phmLen();
  81
  82   xInc = (xMax - xMin) / nx;    // size of cells
  83   yInc = (yMax - yMin) / ny;
  84
  85   if (interpType != I_NEAREST && interpType != I_LINEAR)
  86     sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
  87 }
  88
  89 Backproject::~Backproject (void)
  90 {}
  91
  92 void
  93 Backproject::ScaleImageByRotIncrement (void)
  94 {
  95   for (int ix = 0; ix < nx; ix++)
  96     for (int iy = 0; iy < ny; iy++)
  97       v[ix][iy] *= rotInc;
  98 }
  99
 100 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 101 {
 102     printf ("r=%f, phi=%f\n", r, phi);
 103     errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 104 }
 105
 106 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
 107 {
 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);
 114 }
 115
 116
 117 // CLASS IDENTICATION
 118 //   BackprojectTrig
 119 //
 120 // PURPOSE
 121 //   Uses trigometric functions at each point in image for backprojection.
 122
 123 void
 124 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 125 {
 126   double theta = HALFPI + view_angle;   // Add PI/2 to get perpendicular angle to detector
 127   int ix, iy;
 128   double x, y;                  // Rectang coords of center of pixel
 129
 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
 135
 136       if (interpType == I_NEAREST) {
 137         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 138
 139         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 140             errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 141         else
 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);
 150           else
 151             v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 152       }
 153     }
 154 }
 155
 156
 157 // CLASS IDENTICATION
 158 //   BackprojectTable
 159 //
 160 // PURPOSE
 161 //   Precalculates trigometric function value for each point in image for backprojection.
 162
 163 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, InterpolationType interpType)
 164   : Backproject::Backproject (proj, im, interpType)
 165 {
 166   arrayR.initSetSize (nx, ny);
 167   arrayPhi.initSetSize (nx, ny);
 168   r = arrayR.getArray();
 169   phi = arrayPhi.getArray();
 170
 171   double x, y;                  // Rectang coords of center of pixel
 172   int ix, iy;
 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);
 177     }
 178 }
 179
 180 BackprojectTable::~BackprojectTable (void)
 181 {
 182   ScaleImageByRotIncrement();
 183 }
 184
 185 void
 186 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 187 {
 188   double theta = HALFPI + view_angle;  // add half PI to view angle to get perpendicular theta angle
 189
 190   for (int ix = 0; ix < nx; ix++) {
 191     ImageFileColumn pImCol = v[ix];
 192
 193     for (int iy = 0; iy < ny; iy++) {
 194       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 195
 196       if (interpType == I_NEAREST) {
 197         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 198
 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);
 201         else
 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);
 210         else
 211           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 212       }
 213     }   // end for y
 214   }     // end for x
 215 }
 216
 217
 218 // CLASS IDENTICATION
 219 //   BackprojectDiff
 220 //
 221 // PURPOSE
 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
 224
 225 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, InterpolationType interpType)
 226   :  Backproject::Backproject (proj, im, interpType)
 227 {
 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);
 233
 234   im.arrayDataClear();
 235 }
 236
 237 BackprojectDiff::~BackprojectDiff()
 238 {
 239   ScaleImageByRotIncrement();
 240 }
 241
 242 void
 243 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 244 {
 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
 249
 250   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
 251     double curDetPos = lColStart;
 252     ImageFileColumn pImCol = v[ix];
 253
 254     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 255 #ifdef DEBUG
 256       printf ("[%2d,%2d]:  %8.5lf  ", ix, iy, curDetPos);
 257 #endif
 258       if (interpType == I_NEAREST) {
 259         int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);    // calc index in the filtered raysum vector
 260
 261         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 262             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 263         else
 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);
 272         else
 273           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 274       }
 275     }   // end for y
 276   }     // end for x
 277 }
 278
 279
 280 // CLASS IDENTICATION
 281 //   BackprojectDiff2
 282 //
 283 // PURPOSE
 284 //   Optimized version of BackprojectDiff
 285
 286 void
 287 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 288 {
 289   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 290
 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;
 294
 295   // calculate detPosition for first point in image (ix=0, iy=0)
 296   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 297
 298 #ifdef DEBUG
 299   printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
 300 #endif
 301   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 302     double curDetPos = detPosColStart;
 303     ImageFileColumn pImCol = v[ix];
 304
 305     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 306 #ifdef DEBUG
 307       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
 308 #endif
 309       if (interpType == I_NEAREST) {
 310         int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector
 311
 312         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 313             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 314         else
 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);
 322         else
 323           *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 324       }
 325     }   // end for y
 326   }     // end for x
 327 }
 328
 329 // CLASS IDENTICATION
 330 //   BackprojectIntDiff2
 331 //
 332 // PURPOSE
 333 //   Integer version of BackprojectDiff2
 334
 335 void
 336 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 337 {
 338   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 339
 340   kint32 scale = 1 << 16;
 341   double dScale = scale;
 342   kint32 halfScale = scale / 2;
 343
 344   kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
 345   kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
 346
 347   // calculate L for first point in image (0, 0)
 348   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
 349
 350   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 351     kint32 curDetPos = detPosColStart;
 352     ImageFileColumn pImCol = v[ix];
 353
 354     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 355       if (interpType == I_NEAREST) {
 356         int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
 357         int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector
 358
 359         if (iDetPos < 0 || iDetPos >= nDet)  // check for impossible: index outside of raysum pos
 360             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 361         else
 362           *pImCol++ += filteredProj[iDetPos];
 363       } else if (interpType == I_LINEAR) {
 364         kint32 detPosFloor = curDetPos / scale;
 365         kint32 detPosRemainder = curDetPos % scale;
 366         if (detPosRemainder < 0) {
 367           detPosFloor--;
 368           detPosRemainder += scale;
 369         }
 370         int iDetPos = iDetCenter + detPosFloor;
 371         double frac = detPosRemainder / dScale;
 372         if (iDetPos < 0 || iDetPos >= nDet - 1)
 373             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 374         else
 375           *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 376       }
 377     }   // end for y
 378   }     // end for x
 379 }