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