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$
  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 #if SIZEOF_LONG == 4
 576   static const int scaleShift = 16;
 577 #elif SIZEOF_LONG == 8
 578   static const int scaleShift = 32;
 579 #endif
 580   static const long scale = (1 << scaleShift);
 581   static const long scaleBitmask = scale - 1;
 582   static const long halfScale = scale / 2;
 583   static const double dInvScale = 1. / scale;
 584
 585   const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
 586   const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
 587
 588   // calculate L for first point in image (0, 0)
 589   long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 590
 591   double* deltaFilteredProj = NULL;
 592   CubicPolyInterpolator* pCubicInterp = NULL;
 593   if (interpType == Backprojector::INTERP_LINEAR) {
 594     // precalculate scaled difference for linear interpolation
 595     deltaFilteredProj = new double [nDet];
 596     for (int i = 0; i < nDet - 1; i++)
 597       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 598     deltaFilteredProj[nDet - 1] = 0;  // last detector
 599   } else if (interpType == Backprojector::INTERP_CUBIC) {
 600     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 601   }
 602
 603   int iLastDet = nDet - 1;
 604   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 605     kint32 curDetPos = detPosColStart;
 606     ImageFileColumn pImCol = v[ix];
 607
 608     if (interpType == Backprojector::INTERP_NEAREST) {
 609       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 610         const int iDetPos = (curDetPos + halfScale) >> scaleShift;
 611         if (iDetPos >= 0 && iDetPos <= iLastDet)
 612           *pImCol++ += filteredProj[iDetPos];
 613       } // end for iy
 614     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 615       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 616         const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
 617         if (iDetPos >= 0 && iDetPos <= iLastDet)
 618           *pImCol++ += filteredProj[iDetPos];
 619       } // end for iy
 620     } else if (interpType == Backprojector::INTERP_LINEAR) {
 621       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 622         const long iDetPos = curDetPos >> scaleShift;
 623         const long detRemainder = curDetPos & scaleBitmask;
 624         if (iDetPos >= 0 && iDetPos <= iLastDet)
 625           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 626       } // end for iy
 627     } else if (interpType == Backprojector::INTERP_CUBIC) {
 628       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 629         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
 630       }
 631     } // end Cubic
 632   } // end for ix
 633
 634   if (interpType == Backprojector::INTERP_LINEAR)
 635     delete deltaFilteredProj;
 636   else if (interpType == Backprojector::INTERP_CUBIC)
 637     delete pCubicInterp;
 638 }
 639
 640
 641 void
 642 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
 643 {
 644   double beta = view_angle;
 645
 646   CubicPolyInterpolator* pCubicInterp = NULL;
 647   if (interpType == Backprojector::INTERP_CUBIC)
 648     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 649
 650   for (int ix = 0; ix < nx; ix++) {
 651     ImageFileColumn pImCol = v[ix];
 652
 653     for (int iy = 0; iy < ny; iy++) {
 654       double dAngleDiff = beta - phi[ix][iy];
 655       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 656       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 657       double dFLPlusSin = m_dFocalLength + rsin_t;
 658       double gamma =  atan (rcos_t / dFLPlusSin);
 659       double dPos = gamma / detInc;  // position along detector
 660       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
 661
 662       if (interpType == Backprojector::INTERP_NEAREST) {
 663         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 664         if (iDetPos >= 0 && iDetPos < nDet)
 665           pImCol[iy] += filteredProj[iDetPos] / dL2;
 666       } else if (interpType == Backprojector::INTERP_LINEAR) {
 667         double dPosFloor = floor (dPos);
 668         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 669         double frac = dPos - dPosFloor; // fraction distance from det
 670         if (iDetPos >= 0 && iDetPos < nDet - 1)
 671           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
 672       } else if (interpType == Backprojector::INTERP_CUBIC) {
 673         double d = iDetCenter + dPos;           // position along detector
 674         if (d >= 0 && d < nDet)
 675           pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
 676       }
 677     }   // end for y
 678   }     // end for x
 679
 680   if (interpType == Backprojector::INTERP_CUBIC)
 681     delete pCubicInterp;
 682 }
 683
 684 void
 685 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 686 {
 687   double beta = view_angle;
 688
 689   CubicPolyInterpolator* pCubicInterp = NULL;
 690   if (interpType == Backprojector::INTERP_CUBIC)
 691     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 692
 693   for (int ix = 0; ix < nx; ix++) {
 694     ImageFileColumn pImCol = v[ix];
 695
 696     for (int iy = 0; iy < ny; iy++) {
 697       double dAngleDiff = beta - phi[ix][iy];
 698       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 699       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 700
 701       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
 702       double dDetPos =  rcos_t / dU;
 703       // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
 704       dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
 705       double dPos = dDetPos / detInc;  // position along detector array
 706
 707       if (interpType == Backprojector::INTERP_NEAREST) {
 708         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 709         if (iDetPos >= 0 && iDetPos < nDet)
 710           pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
 711       } else if (interpType == Backprojector::INTERP_LINEAR) {
 712         double dPosFloor = floor (dPos);
 713         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 714         double frac = dPos - dPosFloor; // fraction distance from det
 715         if (iDetPos >= 0 && iDetPos < nDet - 1)
 716           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
 717                            / (dU * dU);
 718       } else if (interpType == Backprojector::INTERP_CUBIC) {
 719         double d = iDetCenter + dPos;           // position along detector
 720         if (d >= 0 && d < nDet)
 721           pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
 722       }
 723     }   // end for y
 724   }     // end for x
 725
 726   if (interpType == Backprojector::INTERP_CUBIC)
 727     delete pCubicInterp;
 728 }
 729