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
 274   if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
 275     rotScale = PI / proj.nView(); // scale by number of PI rotations
 276   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
 277     rotScale =  (2 * PI) / proj.nView(); // 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 = iDetCenter + start_r * cos (theta - start_phi) / detInc;
 530
 531   CubicPolyInterpolator* pCubicInterp = NULL;
 532   double* deltaFilteredProj = NULL;
 533   if (interpType == Backprojector::INTERP_LINEAR) {
 534     // precalculate scaled difference for linear interpolation
 535     deltaFilteredProj = new double [nDet];
 536     for (int i = 0; i < nDet - 1; i++)
 537       deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
 538     deltaFilteredProj[nDet - 1] = 0;  // last detector
 539   } else if (interpType == Backprojector::INTERP_CUBIC) {
 540     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 541   }
 542
 543   int iLastDet = nDet - 1;
 544   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 545     double curDetPos = detPosColStart;
 546     ImageFileColumn pImCol = v[ix];
 547
 548     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 549       if (interpType == Backprojector::INTERP_NEAREST) {
 550         int iDetPos = nearest<int> (curDetPos); // calc index in the filtered raysum vector
 551
 552         if (iDetPos >= 0 && iDetPos < nDet)
 553           *pImCol++ += filteredProj[iDetPos];
 554       } else if (interpType == Backprojector::INTERP_LINEAR) {
 555         double detPosFloor = floor (curDetPos);
 556         int iDetPos = static_cast<int>(detPosFloor);
 557         double frac = curDetPos - detPosFloor;  // fraction distance from det
 558         if (iDetPos >= 0 && iDetPos <= iLastDet)
 559             *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
 560       } else if (interpType == Backprojector::INTERP_CUBIC) {
 561         double p = curDetPos;   // position along detector
 562         if (p >= 0 && p < nDet)
 563           *pImCol++  += pCubicInterp->interpolate (p);
 564       }
 565     }   // end for y
 566   }     // end for x
 567
 568   if (interpType == Backprojector::INTERP_LINEAR)
 569     delete deltaFilteredProj;
 570   else if (interpType == Backprojector::INTERP_CUBIC)
 571     delete pCubicInterp;
 572 }
 573
 574
 575 // CLASS IDENTICATION
 576 //   BackprojectIntDiff
 577 //
 578 // PURPOSE
 579 //   Highly optimized and integer version of BackprojectDiff
 580
 581 void
 582 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 583 {
 584   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
 585 #if SIZEOF_LONG == 4
 586   static const int scaleShift = 16;
 587 #elif SIZEOF_LONG == 8
 588   static const int scaleShift = 32;
 589 #endif
 590   static const long scale = (1L << scaleShift);
 591   static const long scaleBitmask = scale - 1;
 592   static const long halfScale = scale / 2;
 593   static const double dInvScale = 1. / scale;
 594
 595   const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
 596   const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
 597
 598   // calculate L for first point in image (0, 0)
 599   long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 600
 601   double* deltaFilteredProj = NULL;
 602   CubicPolyInterpolator* pCubicInterp = NULL;
 603   if (interpType == Backprojector::INTERP_LINEAR) {
 604     // precalculate scaled difference for linear interpolation
 605     deltaFilteredProj = new double [nDet];
 606     for (int i = 0; i < nDet - 1; i++)
 607       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 608     deltaFilteredProj[nDet - 1] = 0;  // last detector
 609   } else if (interpType == Backprojector::INTERP_CUBIC) {
 610     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 611   }
 612
 613   int iLastDet = nDet - 1;
 614   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 615     long curDetPos = detPosColStart;
 616     ImageFileColumn pImCol = v[ix];
 617
 618     if (interpType == Backprojector::INTERP_NEAREST) {
 619       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 620         const int iDetPos = (curDetPos + halfScale) >> scaleShift;
 621         if (iDetPos >= 0 && iDetPos <= iLastDet)
 622           *pImCol++ += filteredProj[iDetPos];
 623         else
 624           pImCol++;
 625
 626       } // end for iy
 627     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 628       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 629         const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
 630         if (iDetPos >= 0 && iDetPos <= iLastDet)
 631           *pImCol++ += filteredProj[iDetPos];
 632         else
 633           pImCol++;
 634       } // end for iy
 635     } else if (interpType == Backprojector::INTERP_LINEAR) {
 636       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 637         const long iDetPos = curDetPos >> scaleShift;
 638         if (iDetPos >= 0 && iDetPos <= iLastDet) {
 639           const long detRemainder = curDetPos & scaleBitmask;
 640           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 641         } else
 642           pImCol++;
 643       } // end for iy
 644     } else if (interpType == Backprojector::INTERP_CUBIC) {
 645       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 646         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
 647       }
 648     } // end Cubic
 649   } // end for ix
 650
 651   if (interpType == Backprojector::INTERP_LINEAR)
 652     delete deltaFilteredProj;
 653   else if (interpType == Backprojector::INTERP_CUBIC)
 654     delete pCubicInterp;
 655 }
 656
 657
 658 void
 659 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
 660 {
 661   double beta = view_angle;
 662
 663   CubicPolyInterpolator* pCubicInterp = NULL;
 664   if (interpType == Backprojector::INTERP_CUBIC)
 665     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 666
 667   for (int ix = 0; ix < nx; ix++) {
 668     ImageFileColumn pImCol = v[ix];
 669
 670     for (int iy = 0; iy < ny; iy++) {
 671       double dAngleDiff = beta - phi[ix][iy];
 672       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 673       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 674       double dFLPlusSin = m_dFocalLength + rsin_t;
 675       double gamma =  atan (rcos_t / dFLPlusSin);
 676       double dPos = gamma / detInc;  // position along detector
 677       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
 678
 679       if (interpType == Backprojector::INTERP_NEAREST) {
 680         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 681         if (iDetPos >= 0 && iDetPos < nDet)
 682           pImCol[iy] += filteredProj[iDetPos] / dL2;
 683       } else if (interpType == Backprojector::INTERP_LINEAR) {
 684         double dPosFloor = floor (dPos);
 685         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 686         double frac = dPos - dPosFloor; // fraction distance from det
 687         if (iDetPos >= 0 && iDetPos < nDet - 1)
 688           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
 689       } else if (interpType == Backprojector::INTERP_CUBIC) {
 690         double d = iDetCenter + dPos;           // position along detector
 691         if (d >= 0 && d < nDet)
 692           pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
 693       }
 694     }   // end for y
 695   }     // end for x
 696
 697   if (interpType == Backprojector::INTERP_CUBIC)
 698     delete pCubicInterp;
 699 }
 700
 701 void
 702 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 703 {
 704   double beta = view_angle;
 705
 706   CubicPolyInterpolator* pCubicInterp = NULL;
 707   if (interpType == Backprojector::INTERP_CUBIC)
 708     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 709
 710   for (int ix = 0; ix < nx; ix++) {
 711     ImageFileColumn pImCol = v[ix];
 712
 713     for (int iy = 0; iy < ny; iy++) {
 714       double dAngleDiff = beta - phi[ix][iy];
 715       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 716       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 717
 718       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
 719       double dDetPos =  rcos_t / dU;
 720       // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
 721       dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
 722       double dPos = dDetPos / detInc;  // position along detector array
 723
 724       if (interpType == Backprojector::INTERP_NEAREST) {
 725         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 726         if (iDetPos >= 0 && iDetPos < nDet)
 727           pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
 728       } else if (interpType == Backprojector::INTERP_LINEAR) {
 729         double dPosFloor = floor (dPos);
 730         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 731         double frac = dPos - dPosFloor; // fraction distance from det
 732         if (iDetPos >= 0 && iDetPos < nDet - 1)
 733           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
 734                            / (dU * dU);
 735       } else if (interpType == Backprojector::INTERP_CUBIC) {
 736         double d = iDetCenter + dPos;           // position along detector
 737         if (d >= 0 && d < nDet)
 738           pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
 739       }
 740     }   // end for y
 741   }     // end for x
 742
 743   if (interpType == Backprojector::INTERP_CUBIC)
 744     delete pCubicInterp;
 745 }
 746