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.32 2001/03/18 18:08:25 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   im.setAxisIncrement (xInc, yInc);
 317   im.setAxisExtent (xMin, xMax, yMin, yMax);
 318
 319   m_dFocalLength = proj.focalLength();
 320   m_dSourceDetectorLength = proj.sourceDetectorLength();
 321 }
 322
 323 Backproject::~Backproject ()
 324 {}
 325
 326 void
 327 Backproject::PostProcessing()
 328 {
 329   m_bPostProcessingDone = true;
 330 }
 331
 332 void
 333 Backproject::ScaleImageByRotIncrement ()
 334 {
 335   for (int ix = 0; ix < nx; ix++)
 336     for (int iy = 0; iy < ny; iy++)
 337       v[ix][iy] *= rotScale;
 338 }
 339
 340 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 341 {
 342   sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
 343   errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 344 }
 345
 346 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
 347 {
 348 #if 1
 349   std::ostringstream os;
 350   os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
 351   os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
 352   os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
 353   os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
 354   os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
 355
 356   sys_error (ERR_WARNING, os.str().c_str());
 357 #endif
 358 }
 359
 360
 361 // CLASS IDENTICATION
 362 //   BackprojectTrig
 363 //
 364 // PURPOSE
 365 //   Uses trigometric functions at each point in image for backprojection.
 366
 367 void
 368 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 369 {
 370   double theta = view_angle;
 371
 372   CubicPolyInterpolator* pCubicInterp = NULL;
 373   if (interpType == Backprojector::INTERP_CUBIC)
 374     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 375
 376   double x = xMin + xInc / 2;   // Rectang coords of center of pixel
 377   for (int ix = 0; ix < nx; x += xInc, ix++) {
 378     double y = yMin + yInc / 2;
 379     for (int iy = 0; iy < ny; y += yInc, iy++) {
 380       double r = sqrt (x * x + y * y);   // distance of cell from center
 381       double phi = atan2 (y, x);         // angle of cell from center
 382       double L = r * cos (theta - phi);  // position on detector
 383
 384       if (interpType == Backprojector::INTERP_NEAREST) {
 385         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 386
 387         if (iDetPos >= 0 && iDetPos < nDet)
 388           v[ix][iy] += rotScale * filteredProj[iDetPos];
 389       } else if (interpType == Backprojector::INTERP_LINEAR) {
 390         double p = L / detInc;  // position along detector
 391         double pFloor = floor (p);
 392         int iDetPos = iDetCenter + static_cast<int>(pFloor);
 393         double frac = p - pFloor;       // fraction distance from det
 394         if (iDetPos >= 0 && iDetPos < nDet - 1)
 395           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 396       } else if (interpType == Backprojector::INTERP_CUBIC) {
 397         double p = iDetCenter + (L / detInc);   // position along detector
 398         if (p >= 0 && p < nDet)
 399           v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
 400       }
 401     }
 402   }
 403
 404   if (interpType == Backprojector::INTERP_CUBIC)
 405     delete pCubicInterp;
 406 }
 407
 408
 409 // CLASS IDENTICATION
 410 //   BackprojectTable
 411 //
 412 // PURPOSE
 413 //   Precalculates trigometric function value for each point in image for backprojection.
 414
 415 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
 416                                     const int interpFactor, const ReconstructionROI* pROI)
 417 : Backproject (proj, im, interpType, interpFactor, pROI)
 418 {
 419   arrayR.initSetSize (im.nx(), im.ny());
 420   arrayPhi.initSetSize (im.nx(), im.ny());
 421   r = arrayR.getArray();
 422   phi = arrayPhi.getArray();
 423
 424   double x, y;                  // Rectang coords of center of pixel
 425   int ix, iy;
 426   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 427     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 428       r[ix][iy] = sqrt (x * x + y * y);
 429       phi[ix][iy] = atan2 (y, x);
 430     }
 431 }
 432
 433 BackprojectTable::~BackprojectTable ()
 434 {
 435 }
 436
 437 void
 438 BackprojectTable::PostProcessing()
 439 {
 440   if (! m_bPostProcessingDone) {
 441     ScaleImageByRotIncrement();
 442     m_bPostProcessingDone = true;
 443   }
 444 }
 445
 446 void
 447 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 448 {
 449   double theta = view_angle;
 450
 451   CubicPolyInterpolator* pCubicInterp = NULL;
 452   if (interpType == Backprojector::INTERP_CUBIC)
 453     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 454
 455   for (int ix = 0; ix < nx; ix++) {
 456     ImageFileColumn pImCol = v[ix];
 457
 458     for (int iy = 0; iy < ny; iy++) {
 459       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 460
 461       if (interpType == Backprojector::INTERP_NEAREST) {
 462         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 463
 464         if (iDetPos >= 0 && iDetPos < nDet)
 465           pImCol[iy] += filteredProj[iDetPos];
 466       } else if (interpType == Backprojector::INTERP_LINEAR) {
 467         double dPos = L / detInc;               // position along detector
 468         double dPosFloor = floor (dPos);
 469         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 470         double frac = dPos - dPosFloor; // fraction distance from det
 471         if (iDetPos >= 0 && iDetPos < nDet - 1)
 472           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 473       } else if (interpType == Backprojector::INTERP_CUBIC) {
 474         double p = iDetCenter + (L / detInc);   // position along detector
 475         if (p >= 0 && p < nDet)
 476           pImCol[iy] += pCubicInterp->interpolate (p);
 477       }
 478     }   // end for y
 479   }     // end for x
 480
 481   if (interpType == Backprojector::INTERP_CUBIC)
 482     delete pCubicInterp;
 483 }
 484
 485
 486 // CLASS IDENTICATION
 487 //   BackprojectDiff
 488 //
 489 // PURPOSE
 490 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 491 //   Iterates in x & y direction by adding difference in L position
 492
 493 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
 494                                   const int interpFactor, const ReconstructionROI* pROI)
 495 :  Backproject (proj, im, interpType, interpFactor, pROI)
 496 {
 497   // calculate center of first pixel v[0][0]
 498   double x = xMin + xInc / 2;
 499   double y = yMin + yInc / 2;
 500   start_r = sqrt (x * x + y * y);
 501   start_phi = atan2 (y, x);
 502
 503   im.arrayDataClear();
 504 }
 505
 506 BackprojectDiff::~BackprojectDiff ()
 507 {
 508 }
 509
 510 void
 511 BackprojectDiff::PostProcessing()
 512 {
 513   if (! m_bPostProcessingDone) {
 514     ScaleImageByRotIncrement();
 515     m_bPostProcessingDone = true;
 516   }
 517 }
 518
 519 void
 520 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 521 {
 522   double theta = view_angle;
 523
 524   // Distance between detectors for an angle given in units of detectors
 525   double det_dx = xInc * cos (theta) / detInc;
 526   double det_dy = yInc * sin (theta) / detInc;
 527
 528   // calculate detPosition for first point in image (ix=0, iy=0)
 529   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 530
 531   CubicPolyInterpolator* pCubicInterp = NULL;
 532   if (interpType == Backprojector::INTERP_CUBIC)
 533     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 534
 535   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 536     double curDetPos = detPosColStart;
 537     ImageFileColumn pImCol = v[ix];
 538
 539     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 540       if (interpType == Backprojector::INTERP_NEAREST) {
 541         int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector
 542
 543         if (iDetPos >= 0 && iDetPos < nDet)
 544           *pImCol++ += filteredProj[iDetPos];
 545       } else if (interpType == Backprojector::INTERP_LINEAR) {
 546         double detPosFloor = floor (curDetPos);
 547         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 548         double frac = curDetPos - detPosFloor;  // fraction distance from det
 549         if (iDetPos > 0 && iDetPos < nDet - 1)
 550           *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
 551       } else if (interpType == Backprojector::INTERP_CUBIC) {
 552         double p = iDetCenter + curDetPos;      // position along detector
 553         if (p >= 0 && p < nDet)
 554           *pImCol++  += pCubicInterp->interpolate (p);
 555       }
 556     }   // end for y
 557   }     // end for x
 558
 559   if (interpType == Backprojector::INTERP_CUBIC)
 560     delete pCubicInterp;
 561 }
 562
 563
 564 // CLASS IDENTICATION
 565 //   BackprojectIntDiff
 566 //
 567 // PURPOSE
 568 //   Highly optimized and integer version of BackprojectDiff
 569
 570 void
 571 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 572 {
 573   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
 574   static const int scaleShift = 16;
 575   static const kint32 scale = (1 << scaleShift);
 576   static const kint32 scaleBitmask = scale - 1;
 577   static const kint32 halfScale = scale / 2;
 578   static const double dInvScale = 1. / scale;
 579
 580   const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
 581   const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
 582
 583   // calculate L for first point in image (0, 0)
 584   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 585
 586   double* deltaFilteredProj = NULL;
 587   CubicPolyInterpolator* pCubicInterp = NULL;
 588   if (interpType == Backprojector::INTERP_LINEAR) {
 589     // precalculate scaled difference for linear interpolation
 590     deltaFilteredProj = new double [nDet];
 591     for (int i = 0; i < nDet - 1; i++)
 592       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 593     deltaFilteredProj[nDet - 1] = 0;  // last detector
 594   } else if (interpType == Backprojector::INTERP_CUBIC) {
 595     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 596   }
 597
 598   int iLastDet = nDet - 1;
 599   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 600     kint32 curDetPos = detPosColStart;
 601     ImageFileColumn pImCol = v[ix];
 602
 603     if (interpType == Backprojector::INTERP_NEAREST) {
 604       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 605         const int iDetPos = (curDetPos + halfScale) >> 16;
 606         if (iDetPos >= 0 && iDetPos <= iLastDet)
 607           *pImCol++ += filteredProj[iDetPos];
 608       } // end for iy
 609     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 610       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 611         const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
 612         if (iDetPos >= 0 && iDetPos <= iLastDet)
 613           *pImCol++ += filteredProj[iDetPos];
 614       } // end for iy
 615     } else if (interpType == Backprojector::INTERP_LINEAR) {
 616       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 617         const kint32 iDetPos = curDetPos >> scaleShift;
 618         const kint32 detRemainder = curDetPos & scaleBitmask;
 619         if (iDetPos >= 0 && iDetPos <= iLastDet)
 620           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 621       } // end for iy
 622     } else if (interpType == Backprojector::INTERP_CUBIC) {
 623       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 624         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
 625       }
 626     } // end Cubic
 627   } // end for ix
 628
 629   if (interpType == Backprojector::INTERP_LINEAR)
 630     delete deltaFilteredProj;
 631   else if (interpType == Backprojector::INTERP_CUBIC)
 632     delete pCubicInterp;
 633 }
 634
 635
 636 void
 637 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
 638 {
 639   double beta = view_angle;
 640
 641   CubicPolyInterpolator* pCubicInterp = NULL;
 642   if (interpType == Backprojector::INTERP_CUBIC)
 643     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 644
 645   for (int ix = 0; ix < nx; ix++) {
 646     ImageFileColumn pImCol = v[ix];
 647
 648     for (int iy = 0; iy < ny; iy++) {
 649       double dAngleDiff = beta - phi[ix][iy];
 650       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 651       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 652       double dFLPlusSin = m_dFocalLength + rsin_t;
 653       double gamma =  atan (rcos_t / dFLPlusSin);
 654       double dPos = gamma / detInc;  // position along detector
 655       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
 656
 657       if (interpType == Backprojector::INTERP_NEAREST) {
 658         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 659         if (iDetPos >= 0 && iDetPos < nDet)
 660           pImCol[iy] += filteredProj[iDetPos] / dL2;
 661       } else if (interpType == Backprojector::INTERP_LINEAR) {
 662         double dPosFloor = floor (dPos);
 663         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 664         double frac = dPos - dPosFloor; // fraction distance from det
 665         if (iDetPos >= 0 && iDetPos < nDet - 1)
 666           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
 667       } else if (interpType == Backprojector::INTERP_CUBIC) {
 668         double d = iDetCenter + dPos;           // position along detector
 669         if (d >= 0 && d < nDet)
 670           pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
 671       }
 672     }   // end for y
 673   }     // end for x
 674
 675   if (interpType == Backprojector::INTERP_CUBIC)
 676     delete pCubicInterp;
 677 }
 678
 679 void
 680 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 681 {
 682   double beta = view_angle;
 683
 684   CubicPolyInterpolator* pCubicInterp = NULL;
 685   if (interpType == Backprojector::INTERP_CUBIC)
 686     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 687
 688   for (int ix = 0; ix < nx; ix++) {
 689     ImageFileColumn pImCol = v[ix];
 690
 691     for (int iy = 0; iy < ny; iy++) {
 692       double dAngleDiff = beta - phi[ix][iy];
 693       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 694       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 695
 696       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
 697       double dDetPos =  rcos_t / dU;
 698       // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
 699       dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
 700       double dPos = dDetPos / detInc;  // position along detector array
 701
 702       if (interpType == Backprojector::INTERP_NEAREST) {
 703         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 704         if (iDetPos >= 0 && iDetPos < nDet)
 705           pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
 706       } else if (interpType == Backprojector::INTERP_LINEAR) {
 707         double dPosFloor = floor (dPos);
 708         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 709         double frac = dPos - dPosFloor; // fraction distance from det
 710         if (iDetPos >= 0 && iDetPos < nDet - 1)
 711           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
 712                            / (dU * dU);
 713       } else if (interpType == Backprojector::INTERP_CUBIC) {
 714         double d = iDetCenter + dPos;           // position along detector
 715         if (d >= 0 && d < nDet)
 716           pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
 717       }
 718     }   // end for y
 719   }     // end for x
 720
 721   if (interpType == Backprojector::INTERP_CUBIC)
 722     delete pCubicInterp;
 723 }
 724