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