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.12 2000/08/25 15:59:13 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 const int Backprojector::BPROJ_INVALID = -1;
  30 const int Backprojector::BPROJ_TRIG = 0;
  31 const int Backprojector::BPROJ_TABLE = 1;
  32 const int Backprojector::BPROJ_DIFF = 2;
  33 const int Backprojector::BPROJ_DIFF2 = 3;
  34 const int Backprojector::BPROJ_IDIFF2 = 4;
  35 const int Backprojector::BPROJ_IDIFF3 = 5;
  36
  37 const char* Backprojector::s_aszBackprojectName[] =
  38 {
  39   {"trig"},
  40   {"table"},
  41   {"diff"},
  42   {"diff2"},
  43   {"idiff2"},
  44   {"idiff3"},
  45 };
  46
  47 const char* Backprojector::s_aszBackprojectTitle[] =
  48 {
  49   {"Direct Trigometric"},
  50   {"Trigometric Table"},
  51   {"Difference Iteration"},
  52   {"Difference Iteration Optimized"},
  53   {"Integer Difference Iteration Optimized"},
  54   {"Integer Difference Iteration Highly-Optimized"},
  55 };
  56
  57 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
  58
  59 const int Backprojector::INTERP_INVALID = -1;
  60 const int Backprojector::INTERP_NEAREST = 0;
  61 const int Backprojector::INTERP_LINEAR = 1;
  62 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 2;
  63 #if HAVE_BSPLINE_INTERP
  64 const int Backprojector::INTERP_BSPLINE = 3;
  65 const int Backprojector::INTERP_1BSPLINE = 4;
  66 const int Backprojector::INTERP_2BSPLINE = 5;
  67 const int Backprojector::INTERP_3BSPLINE = 6;
  68 #endif
  69
  70 const char* Backprojector::s_aszInterpName[] =
  71 {
  72   {"nearest"},
  73   {"linear"},
  74   {"freq_preinterpolationj"},
  75 #if HAVE_BSPLINE_INTERP
  76   {"bspline"},
  77   {"1bspline"},
  78   {"2bspline"},
  79   {"3bspline"},
  80 #endif
  81 };
  82
  83 const char* Backprojector::s_aszInterpTitle[] =
  84 {
  85   {"Nearest"},
  86   {"Linear"},
  87   {"Frequency Preinterpolationj"},
  88 #if HAVE_BSPLINE_INTERP
  89   {"B-Spline"},
  90   {"B-Spline 1st Order"},
  91   {"B-Spline 2nd Order"},
  92   {"B-Spline 3rd Order"},
  93 #endif
  94 };
  95
  96 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
  97
  98
  99
 100 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 101 {
 102   m_fail = false;
 103   m_pBackprojectImplem = NULL;
 104
 105   initBackprojector (proj, im, backprojName, interpName, interpFactor);
 106 }
 107
 108 void
 109 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
 110 {
 111   if (m_pBackprojectImplem != NULL)
 112     m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
 113 }
 114
 115 Backprojector::~Backprojector ()
 116 {
 117   delete m_pBackprojectImplem;
 118 }
 119
 120 // FUNCTION IDENTIFICATION
 121 //     Backproject* projector = selectBackprojector (...)
 122 //
 123 // PURPOSE
 124 //     Selects a backprojector based on BackprojType
 125 //     and initializes the backprojector
 126
 127 bool
 128 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 129 {
 130   m_nameBackproject = backprojName;
 131   m_nameInterpolation = interpName;
 132   m_pBackprojectImplem = NULL;
 133   m_idBackproject = convertBackprojectNameToID (backprojName);
 134   if (m_idBackproject == BPROJ_INVALID) {
 135     m_fail = true;
 136     m_failMessage = "Invalid backprojection name ";
 137     m_failMessage += backprojName;
 138   }
 139   m_idInterpolation = convertInterpNameToID (interpName);
 140   if (m_idInterpolation == INTERP_INVALID) {
 141     m_fail = true;
 142     m_failMessage = "Invalid interpolation name ";
 143     m_failMessage += interpName;
 144   }
 145
 146   if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
 147     m_fail = true;
 148     return false;
 149   }
 150
 151   if (m_idBackproject == BPROJ_TRIG)
 152     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
 153   else if (m_idBackproject == BPROJ_TABLE)
 154     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
 155   else if (m_idBackproject == BPROJ_DIFF)
 156     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
 157   else if (m_idBackproject == BPROJ_DIFF2)
 158     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
 159   else if (m_idBackproject == BPROJ_IDIFF2)
 160     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
 161   else if (m_idBackproject == BPROJ_IDIFF3)
 162     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
 163   else {
 164     m_fail = true;
 165     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
 166     return false;
 167   }
 168
 169   return true;
 170 }
 171
 172
 173 int
 174 Backprojector::convertBackprojectNameToID (const char* const backprojName)
 175 {
 176   int backprojID = BPROJ_INVALID;
 177
 178   for (int i = 0; i < s_iBackprojectCount; i++)
 179       if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
 180           backprojID = i;
 181           break;
 182       }
 183
 184   return (backprojID);
 185 }
 186
 187 const char*
 188 Backprojector::convertBackprojectIDToName (int bprojID)
 189 {
 190   static const char *bprojName = "";
 191
 192   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 193       return (s_aszBackprojectName[bprojID]);
 194
 195   return (bprojName);
 196 }
 197
 198 const char*
 199 Backprojector::convertBackprojectIDToTitle (const int bprojID)
 200 {
 201   static const char *bprojTitle = "";
 202
 203   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 204       return (s_aszBackprojectTitle[bprojID]);
 205
 206   return (bprojTitle);
 207 }
 208
 209
 210 int
 211 Backprojector::convertInterpNameToID (const char* const interpName)
 212 {
 213   int interpID = INTERP_INVALID;
 214
 215   for (int i = 0; i < s_iInterpCount; i++)
 216       if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
 217           interpID = i;
 218           break;
 219       }
 220
 221   return (interpID);
 222 }
 223
 224 const char*
 225 Backprojector::convertInterpIDToName (const int interpID)
 226 {
 227   static const char *interpName = "";
 228
 229   if (interpID >= 0 && interpID < s_iInterpCount)
 230       return (s_aszInterpName[interpID]);
 231
 232   return (interpName);
 233 }
 234
 235 const char*
 236 Backprojector::convertInterpIDToTitle (const int interpID)
 237 {
 238   static const char *interpTitle = "";
 239
 240   if (interpID >= 0 && interpID < s_iInterpCount)
 241       return (s_aszInterpTitle[interpID]);
 242
 243   return (interpTitle);
 244 }
 245
 246
 247
 248 // CLASS IDENTICATION
 249 //   Backproject
 250 //
 251 // PURPOSE
 252 //   Pure virtual base class for all backprojectors.
 253
 254 Backproject::Backproject (const Projections& proj, ImageFile& im, const int interpType, const int interpFactor)
 255     : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
 256 {
 257   detInc = proj.detInc();
 258   nDet = proj.nDet();
 259   iDetCenter = (nDet - 1) / 2;  // index refering to L=0 projection
 260   rotInc = proj.rotInc();
 261
 262   v = im.getArray();
 263   nx = im.nx();
 264   ny = im.ny();
 265   im.arrayDataClear();
 266
 267   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
 268   xMax = xMin + proj.phmLen();
 269   yMin = -proj.phmLen() / 2;
 270   yMax = yMin + proj.phmLen();
 271
 272   xInc = (xMax - xMin) / nx;    // size of cells
 273   yInc = (yMax - yMin) / ny;
 274 }
 275
 276 Backproject::~Backproject ()
 277 {}
 278
 279 void
 280 Backproject::ScaleImageByRotIncrement ()
 281 {
 282   for (int ix = 0; ix < nx; ix++)
 283     for (int iy = 0; iy < ny; iy++)
 284       v[ix][iy] *= rotInc;
 285 }
 286
 287 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 288 {
 289     sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
 290     errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 291 }
 292
 293 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
 294 {
 295   ostringstream os;
 296   os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
 297   os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
 298   os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
 299   os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
 300   os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
 301
 302   sys_error (ERR_WARNING, os.str().c_str());
 303 }
 304
 305
 306 // CLASS IDENTICATION
 307 //   BackprojectTrig
 308 //
 309 // PURPOSE
 310 //   Uses trigometric functions at each point in image for backprojection.
 311
 312 void
 313 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 314 {
 315   double theta = HALFPI + view_angle;   // Add PI/2 to get perpendicular angle to detector
 316   int ix, iy;
 317   double x, y;                  // Rectang coords of center of pixel
 318
 319   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 320     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 321       double r = sqrt (x * x + y * y);   // distance of cell from center
 322       double phi = atan2 (y, x);         // angle of cell from center
 323       double L = r * cos (theta - phi);  // position on detector
 324
 325       if (interpType == Backprojector::INTERP_NEAREST) {
 326         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 327
 328         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 329             errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 330         else
 331           v[ix][iy] += rotInc * filteredProj[iDetPos];
 332       } else if (interpType == Backprojector::INTERP_LINEAR) {
 333           double p = L / detInc;        // position along detector
 334           double pFloor = floor (p);
 335           int iDetPos = iDetCenter + static_cast<int>(pFloor);
 336           double frac = p - pFloor;     // fraction distance from det
 337           if (iDetPos < 0 || iDetPos >= nDet - 1)       // check for impossible: index outside of raysum pos
 338             errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 339           else
 340             v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 341       }
 342     }
 343 }
 344
 345
 346 // CLASS IDENTICATION
 347 //   BackprojectTable
 348 //
 349 // PURPOSE
 350 //   Precalculates trigometric function value for each point in image for backprojection.
 351
 352 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
 353   : Backproject::Backproject (proj, im, interpType, interpFactor)
 354 {
 355   arrayR.initSetSize (nx, ny);
 356   arrayPhi.initSetSize (nx, ny);
 357   r = arrayR.getArray();
 358   phi = arrayPhi.getArray();
 359
 360   double x, y;                  // Rectang coords of center of pixel
 361   int ix, iy;
 362   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 363     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 364       r[ix][iy] = sqrt (x * x + y * y);
 365       phi[ix][iy] = atan2 (y, x);
 366     }
 367 }
 368
 369 BackprojectTable::~BackprojectTable ()
 370 {
 371   ScaleImageByRotIncrement();
 372 }
 373
 374 void
 375 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 376 {
 377   double theta = HALFPI + view_angle;  // add half PI to view angle to get perpendicular theta angle
 378
 379   for (int ix = 0; ix < nx; ix++) {
 380     ImageFileColumn pImCol = v[ix];
 381
 382     for (int iy = 0; iy < ny; iy++) {
 383       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 384
 385       if (interpType == Backprojector::INTERP_NEAREST) {
 386         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 387
 388         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 389           errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
 390         else
 391           pImCol[iy] += filteredProj[iDetPos];
 392       } else if (interpType == Backprojector::INTERP_LINEAR) {
 393         double dPos = L / detInc;               // position along detector
 394         double dPosFloor = floor (dPos);
 395         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 396         double frac = dPos - dPosFloor; // fraction distance from det
 397         if (iDetPos < 0 || iDetPos >= nDet - 1)
 398             errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
 399         else
 400           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 401       }
 402     }   // end for y
 403   }     // end for x
 404 }
 405
 406
 407 // CLASS IDENTICATION
 408 //   BackprojectDiff
 409 //
 410 // PURPOSE
 411 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 412 //   Iterates in x & y direction by adding difference in L position
 413
 414 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
 415   :  Backproject::Backproject (proj, im, interpType, interpFactor)
 416 {
 417   // calculate center of first pixel v[0][0]
 418   double x = xMin + xInc / 2;
 419   double y = yMin + yInc / 2;
 420   start_r = sqrt (x * x + y * y);
 421   start_phi = atan2 (y, x);
 422
 423   im.arrayDataClear();
 424 }
 425
 426 BackprojectDiff::~BackprojectDiff()
 427 {
 428   ScaleImageByRotIncrement();
 429 }
 430
 431 void
 432 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 433 {
 434   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 435   double det_dx = xInc * sin (theta);
 436   double det_dy = yInc * cos (theta);
 437   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
 438
 439   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
 440     double curDetPos = lColStart;
 441     ImageFileColumn pImCol = v[ix];
 442
 443     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 444 #ifdef DEBUG
 445       printf ("[%2d,%2d]:  %8.5f  ", ix, iy, curDetPos);
 446 #endif
 447       if (interpType == Backprojector::INTERP_NEAREST) {
 448         int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);    // calc index in the filtered raysum vector
 449
 450         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 451             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 452         else
 453           pImCol[iy] += filteredProj[iDetPos];
 454       } else if (interpType == Backprojector::INTERP_LINEAR) {
 455         double detPos = curDetPos / detInc;             // position along detector
 456         double detPosFloor = floor (detPos);
 457         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 458         double frac = detPos - detPosFloor;     // fraction distance from det
 459         if (iDetPos < 0 || iDetPos >= nDet - 1)
 460             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 461         else
 462           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 463       }
 464     }   // end for y
 465   }     // end for x
 466 }
 467
 468
 469 // CLASS IDENTICATION
 470 //   BackprojectDiff2
 471 //
 472 // PURPOSE
 473 //   Optimized version of BackprojectDiff
 474
 475 void
 476 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 477 {
 478   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 479
 480   // Distance betw. detectors for an angle given in units of detectors
 481   double det_dx = xInc * sin (theta) / detInc;
 482   double det_dy = yInc * cos (theta) / detInc;
 483
 484   // calculate detPosition for first point in image (ix=0, iy=0)
 485   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 486
 487 #ifdef DEBUG
 488   printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
 489 #endif
 490   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 491     double curDetPos = detPosColStart;
 492     ImageFileColumn pImCol = v[ix];
 493
 494     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 495 #ifdef DEBUG
 496       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
 497 #endif
 498       if (interpType == Backprojector::INTERP_NEAREST) {
 499         int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector
 500
 501         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 502             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 503         else
 504           *pImCol++ += filteredProj[iDetPos];
 505       } else if (interpType == Backprojector::INTERP_LINEAR) {
 506         double detPosFloor = floor (curDetPos);
 507         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 508         double frac = curDetPos - detPosFloor;  // fraction distance from det
 509         if (iDetPos < 0 || iDetPos >= nDet - 1)
 510             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 511         else
 512           *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
 513       }
 514     }   // end for y
 515   }     // end for x
 516 }
 517
 518 // CLASS IDENTICATION
 519 //   BackprojectIntDiff2
 520 //
 521 // PURPOSE
 522 //   Integer version of BackprojectDiff2
 523
 524 void
 525 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 526 {
 527   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 528
 529   static const kint32 scale = 1 << 16;
 530   static const double dScale = scale;
 531   static const kint32 halfScale = scale / 2;
 532
 533   const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
 534   const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
 535
 536   // calculate L for first point in image (0, 0)
 537   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
 538
 539   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 540     kint32 curDetPos = detPosColStart;
 541     ImageFileColumn pImCol = v[ix];
 542
 543     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 544       if (interpType == Backprojector::INTERP_NEAREST) {
 545         int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
 546         int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector
 547
 548         if (iDetPos < 0 || iDetPos >= nDet)  // check for index outside of raysum pos
 549             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 550         else
 551           *pImCol++ += filteredProj[iDetPos];
 552       } else if (interpType == Backprojector::INTERP_LINEAR) {
 553         kint32 detPosFloor = curDetPos / scale;
 554         kint32 detPosRemainder = curDetPos % scale;
 555         if (detPosRemainder < 0) {
 556           detPosFloor--;
 557           detPosRemainder += scale;
 558         }
 559         int iDetPos = iDetCenter + detPosFloor;
 560         double frac = detPosRemainder / dScale;
 561         if (iDetPos < 0 || iDetPos >= nDet - 1)
 562             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 563         else
 564           *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 565       }
 566     }   // end for y
 567   }     // end for x
 568 }
 569
 570 // CLASS IDENTICATION
 571 //   BackprojectIntDiff3
 572 //
 573 // PURPOSE
 574 //   Highly optimized version of BackprojectIntDiff2
 575
 576 void
 577 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
 578 {
 579   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 580   static const int scaleShift = 16;
 581   static const kint32 scale = (1 << scaleShift);
 582   static const kint32 scaleBitmask = scale - 1;
 583   static const kint32 halfScale = scale / 2;
 584   static const double dInvScale = 1. / scale;
 585
 586   const kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
 587   const kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
 588
 589   // calculate L for first point in image (0, 0)
 590   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 591
 592   // precalculate scaled difference for linear interpolation
 593   double deltaFilteredProj [nDet];
 594   if (interpType == Backprojector::INTERP_LINEAR) {
 595     for (int i = 0; i < nDet - 1; i++)
 596       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 597   }
 598   deltaFilteredProj[nDet - 1] = 0;  // last detector
 599
 600   int iLastDet = nDet - 1;
 601   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 602     kint32 curDetPos = detPosColStart;
 603     ImageFileColumn pImCol = v[ix];
 604
 605     if (interpType == Backprojector::INTERP_NEAREST) {
 606       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 607         const int iDetPos = (curDetPos + halfScale) >> 16;
 608         if (iDetPos >= 0 && iDetPos <= iLastDet)
 609           *pImCol++ += filteredProj[iDetPos];
 610       } // end for iy
 611     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 612       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 613         const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
 614         if (iDetPos >= 0 && iDetPos <= iLastDet)
 615         *pImCol++ += filteredProj[iDetPos];
 616       } // end for iy
 617     } else if (interpType == Backprojector::INTERP_LINEAR) {
 618       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 619         const kint32 iDetPos = curDetPos >> scaleShift;
 620         const kint32 detRemainder = curDetPos & scaleBitmask;
 621         if (iDetPos >= 0 && iDetPos <= iLastDet)
 622           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 623       } // end for iy
 624     } //end linear
 625   } // end for ix
 626 }