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