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