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