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-2009 Kevin Rosenberg
  10 **
  11 **  This program is free software; you can redistribute it and/or modify
  12 **  it under the terms of the GNU General Public License (version 2) as
  13 **  published by the Free Software Foundation.
  14 **
  15 **  This program is distributed in the hope that it will be useful,
  16 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  17 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18 **  GNU General Public License for more details.
  19 **
  20 **  You should have received a copy of the GNU General Public License
  21 **  along with this program; if not, write to the Free Software
  22 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  23 ******************************************************************************/
  24
  25 #include "ct.h"
  26 #include "interpolator.h"
  27
  28 const int Backprojector::BPROJ_INVALID = -1;
  29 const int Backprojector::BPROJ_TRIG = 0;
  30 const int Backprojector::BPROJ_TABLE = 1;
  31 const int Backprojector::BPROJ_DIFF = 2;
  32 const int Backprojector::BPROJ_IDIFF = 3;
  33
  34 const char* const Backprojector::s_aszBackprojectName[] =
  35 {
  36   "trig",
  37   "table",
  38   "diff",
  39   "idiff",
  40 };
  41
  42 const char* const Backprojector::s_aszBackprojectTitle[] =
  43 {
  44   "Direct Trigometric",
  45   "Trigometric Table",
  46   "Difference Iteration",
  47   "Integer Difference Iteration",
  48 };
  49
  50 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
  51
  52 const int Backprojector::INTERP_INVALID = -1;
  53 const int Backprojector::INTERP_NEAREST = 0;
  54 const int Backprojector::INTERP_LINEAR = 1;
  55 const int Backprojector::INTERP_CUBIC = 2;
  56 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
  57 #if HAVE_BSPLINE_INTERP
  58 const int Backprojector::INTERP_BSPLINE = 4;
  59 const int Backprojector::INTERP_1BSPLINE = 5;
  60 const int Backprojector::INTERP_2BSPLINE = 6;
  61 const int Backprojector::INTERP_3BSPLINE = 7;
  62 #endif
  63
  64 const char* const Backprojector::s_aszInterpName[] =
  65 {
  66   "nearest",
  67   "linear",
  68   "cubic",
  69 #if HAVE_FREQ_PREINTERP
  70   "freq_preinterpolationj",
  71 #endif
  72 #if HAVE_BSPLINE_INTERP
  73   "bspline",
  74   "1bspline",
  75   "2bspline",
  76   "3bspline",
  77 #endif
  78 };
  79
  80 const char* const Backprojector::s_aszInterpTitle[] =
  81 {
  82   "Nearest",
  83   "Linear",
  84   "Cubic",
  85 #if HAVE_FREQ_PREINTERP
  86   "Frequency Preinterpolation",
  87 #endif
  88 #if HAVE_BSPLINE_INTERP
  89   "B-Spline",
  90   "B-Spline 1st Order",
  91   "B-Spline 2nd Order",
  92   "B-Spline 3rd Order",
  93 #endif
  94 };
  95
  96 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
  97
  98
  99
 100 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName,
 101                               const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
 102 {
 103   m_fail = false;
 104   m_pBackprojectImplem = NULL;
 105
 106   initBackprojector (proj, im, backprojName, interpName, interpFactor, pROI);
 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,
 137                                   const char* const interpName, const int interpFactor, const ReconstructionROI* pROI)
 138 {
 139   m_nameBackproject = backprojName;
 140   m_nameInterpolation = interpName;
 141   m_pBackprojectImplem = NULL;
 142   m_idBackproject = convertBackprojectNameToID (backprojName);
 143   if (m_idBackproject == BPROJ_INVALID) {
 144     m_fail = true;
 145     m_failMessage = "Invalid backprojection name ";
 146     m_failMessage += backprojName;
 147   }
 148   m_idInterpolation = convertInterpNameToID (interpName);
 149   if (m_idInterpolation == INTERP_INVALID) {
 150     m_fail = true;
 151     m_failMessage = "Invalid interpolation name ";
 152     m_failMessage += interpName;
 153   }
 154
 155   if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
 156     m_fail = true;
 157     return false;
 158   }
 159
 160   if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
 161     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor, pROI));
 162   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
 163     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor, pROI));
 164   else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
 165     if (m_idBackproject == BPROJ_TRIG)
 166       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor, pROI));
 167     else if (m_idBackproject == BPROJ_TABLE)
 168       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor, pROI));
 169     else if (m_idBackproject == BPROJ_DIFF)
 170       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor, pROI));
 171     else if (m_idBackproject == BPROJ_IDIFF)
 172       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff (proj, im, m_idInterpolation, interpFactor, pROI));
 173   } else {
 174     m_fail = true;
 175     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
 176     return false;
 177   }
 178
 179   return true;
 180 }
 181
 182
 183 int
 184 Backprojector::convertBackprojectNameToID (const char* const backprojName)
 185 {
 186   int backprojID = BPROJ_INVALID;
 187
 188   for (int i = 0; i < s_iBackprojectCount; i++) {
 189     if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
 190       backprojID = i;
 191       break;
 192     }
 193   }
 194   return (backprojID);
 195 }
 196
 197 const char*
 198 Backprojector::convertBackprojectIDToName (int bprojID)
 199 {
 200   static const char *bprojName = "";
 201
 202   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 203     return (s_aszBackprojectName[bprojID]);
 204
 205   return (bprojName);
 206 }
 207
 208 const char*
 209 Backprojector::convertBackprojectIDToTitle (const int bprojID)
 210 {
 211   static const char *bprojTitle = "";
 212
 213   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 214     return (s_aszBackprojectTitle[bprojID]);
 215
 216   return (bprojTitle);
 217 }
 218
 219
 220 int
 221 Backprojector::convertInterpNameToID (const char* const interpName)
 222 {
 223   int interpID = INTERP_INVALID;
 224
 225   for (int i = 0; i < s_iInterpCount; i++) {
 226     if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
 227       interpID = i;
 228       break;
 229     }
 230   }
 231   return (interpID);
 232 }
 233
 234 const char*
 235 Backprojector::convertInterpIDToName (const int interpID)
 236 {
 237   static const char *interpName = "";
 238
 239   if (interpID >= 0 && interpID < s_iInterpCount)
 240     return (s_aszInterpName[interpID]);
 241
 242   return (interpName);
 243 }
 244
 245 const char*
 246 Backprojector::convertInterpIDToTitle (const int interpID)
 247 {
 248   static const char *interpTitle = "";
 249
 250   if (interpID >= 0 && interpID < s_iInterpCount)
 251     return (s_aszInterpTitle[interpID]);
 252
 253   return (interpTitle);
 254 }
 255
 256
 257
 258 // CLASS IDENTICATION
 259 //   Backproject
 260 //
 261 // PURPOSE
 262 //   Pure virtual base class for all backprojectors.
 263
 264 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor,
 265                           const ReconstructionROI* pROI)
 266 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor), m_bPostProcessingDone(false)
 267 {
 268   detInc = proj.detInc();
 269   nDet = proj.nDet();
 270   iDetCenter = (nDet - 1) / 2;  // index refering to L=0 projection
 271
 272   if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
 273     rotScale = PI / proj.nView(); // scale by number of PI rotations
 274   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
 275     rotScale =  (2 * PI) / proj.nView(); // scale by number of 2PI rotations
 276   else
 277     sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
 278
 279   v = im.getArray();
 280   nx = im.nx();
 281   ny = im.ny();
 282   im.arrayDataClear();
 283
 284   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
 285   xMax = xMin + proj.phmLen();
 286   yMin = -proj.phmLen() / 2;
 287   yMax = yMin + proj.phmLen();
 288
 289   if (pROI) {
 290     if (pROI->m_dXMin > xMin)
 291       xMin = pROI->m_dXMin;
 292     if (pROI->m_dXMax < xMax)
 293       xMax = pROI->m_dXMax;
 294     if (pROI->m_dYMin > yMin)
 295       yMin = pROI->m_dYMin;
 296     if (pROI->m_dYMax < yMax)
 297       yMax = pROI->m_dYMax;
 298
 299     if (xMin > xMax) {
 300       double temp = xMin;
 301       xMin = xMax;
 302       xMax = temp;
 303     }
 304     if (yMin > yMax) {
 305       double temp = yMin;
 306       yMin = yMax;
 307       yMax = temp;
 308     }
 309   }
 310
 311   xInc = (xMax - xMin) / nx;    // size of cells
 312   yInc = (yMax - yMin) / ny;
 313
 314   im.setAxisIncrement (xInc, yInc);
 315   im.setAxisExtent (xMin, xMax, yMin, yMax);
 316
 317   m_dFocalLength = proj.focalLength();
 318   m_dSourceDetectorLength = proj.sourceDetectorLength();
 319 }
 320
 321 Backproject::~Backproject ()
 322 {}
 323
 324 void
 325 Backproject::PostProcessing()
 326 {
 327   m_bPostProcessingDone = true;
 328 }
 329
 330 void
 331 Backproject::ScaleImageByRotIncrement ()
 332 {
 333 #if HAVE_OPENMP
 334   #pragma omp parallel for
 335 #endif
 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 xstart = xMin + xInc / 2;   // Rectang coords of center of pixel
 378 #if HAVE_OPENMP
 379   #pragma omp parallel for
 380 #endif
 381   for (int ix = 0; ix < nx; ix++) {
 382     double x = xstart + (ix * xInc);
 383
 384     double y = yMin + yInc / 2;
 385     for (int iy = 0; iy < ny; y += yInc, iy++) {
 386       double r = sqrt (x * x + y * y);   // distance of cell from center
 387       double phi = atan2 (y, x);         // angle of cell from center
 388       double L = r * cos (theta - phi);  // position on detector
 389
 390       if (interpType == Backprojector::INTERP_NEAREST) {
 391         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 392
 393         if (iDetPos >= 0 && iDetPos < nDet)
 394           v[ix][iy] += rotScale * filteredProj[iDetPos];
 395       } else if (interpType == Backprojector::INTERP_LINEAR) {
 396         double p = L / detInc;  // position along detector
 397         double pFloor = floor (p);
 398         int iDetPos = iDetCenter + static_cast<int>(pFloor);
 399         double frac = p - pFloor;       // fraction distance from det
 400         if (iDetPos >= 0 && iDetPos < nDet - 1)
 401           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 402       } else if (interpType == Backprojector::INTERP_CUBIC) {
 403         double p = iDetCenter + (L / detInc);   // position along detector
 404         if (p >= 0 && p < nDet)
 405           v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
 406       }
 407     }
 408   }
 409
 410   if (interpType == Backprojector::INTERP_CUBIC)
 411     delete pCubicInterp;
 412 }
 413
 414
 415 // CLASS IDENTICATION
 416 //   BackprojectTable
 417 //
 418 // PURPOSE
 419 //   Precalculates trigometric function value for each point in image for backprojection.
 420
 421 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType,
 422                                     const int interpFactor, const ReconstructionROI* pROI)
 423 : Backproject (proj, im, interpType, interpFactor, pROI)
 424 {
 425   arrayR.initSetSize (im.nx(), im.ny());
 426   arrayPhi.initSetSize (im.nx(), im.ny());
 427   r = arrayR.getArray();
 428   phi = arrayPhi.getArray();
 429
 430   double xstart = xMin + xInc / 2;
 431
 432 #if HAVE_OPENMP
 433   #pragma omp parallel for
 434 #endif
 435   for (int ix = 0; ix < nx; ix++) {
 436     double x = xstart + (ix * xInc);
 437     double y = yMin + yInc / 2;
 438     for (int iy = 0; iy < ny; iy++, y += yInc) {
 439       r[ix][iy] = sqrt (x * x + y * y);
 440       phi[ix][iy] = atan2 (y, x);
 441     }
 442   }
 443 }
 444
 445 BackprojectTable::~BackprojectTable ()
 446 {
 447 }
 448
 449 void
 450 BackprojectTable::PostProcessing()
 451 {
 452   if (! m_bPostProcessingDone) {
 453     ScaleImageByRotIncrement();
 454     m_bPostProcessingDone = true;
 455   }
 456 }
 457
 458 void
 459 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 460 {
 461   double theta = view_angle;
 462
 463   CubicPolyInterpolator* pCubicInterp = NULL;
 464   if (interpType == Backprojector::INTERP_CUBIC)
 465     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 466
 467 #if HAVE_OPENMP
 468   #pragma omp parallel for
 469 #endif
 470   for (int ix = 0; ix < nx; ix++) {
 471     ImageFileColumn pImCol = v[ix];
 472
 473     for (int iy = 0; iy < ny; iy++) {
 474       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 475
 476       if (interpType == Backprojector::INTERP_NEAREST) {
 477         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 478
 479         if (iDetPos >= 0 && iDetPos < nDet) {
 480           pImCol[iy] += filteredProj[iDetPos];
 481         }
 482       } else if (interpType == Backprojector::INTERP_LINEAR) {
 483         double dPos = L / detInc;               // position along detector
 484         double dPosFloor = floor (dPos);
 485         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 486         double frac = dPos - dPosFloor; // fraction distance from det
 487         if (iDetPos >= 0 && iDetPos < nDet - 1) {
 488           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 489         }
 490       } else if (interpType == Backprojector::INTERP_CUBIC) {
 491         double p = iDetCenter + (L / detInc);   // position along detector
 492         if (p >= 0 && p < nDet) {
 493           pImCol[iy] += pCubicInterp->interpolate (p);
 494         }
 495       }
 496     }   // end for y
 497   }     // end for x
 498
 499   if (interpType == Backprojector::INTERP_CUBIC)
 500     delete pCubicInterp;
 501 }
 502
 503
 504 // CLASS IDENTICATION
 505 //   BackprojectDiff
 506 //
 507 // PURPOSE
 508 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 509 //   Iterates in x & y direction by adding difference in L position
 510
 511 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType,
 512                                   const int interpFactor, const ReconstructionROI* pROI)
 513 :  Backproject (proj, im, interpType, interpFactor, pROI)
 514 {
 515   // calculate center of first pixel v[0][0]
 516   double x = xMin + xInc / 2;
 517   double y = yMin + yInc / 2;
 518   start_r = sqrt (x * x + y * y);
 519   start_phi = atan2 (y, x);
 520
 521   im.arrayDataClear();
 522 }
 523
 524 BackprojectDiff::~BackprojectDiff ()
 525 {
 526 }
 527
 528 void
 529 BackprojectDiff::PostProcessing()
 530 {
 531   if (! m_bPostProcessingDone) {
 532     ScaleImageByRotIncrement();
 533     m_bPostProcessingDone = true;
 534   }
 535 }
 536
 537 void
 538 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 539 {
 540   double theta = view_angle;
 541
 542   // Distance between detectors for an angle given in units of detectors
 543   double det_dx = xInc * cos (theta) / detInc;
 544   double det_dy = yInc * sin (theta) / detInc;
 545
 546   // calculate detPosition for first point in image (ix=0, iy=0)
 547   double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc;
 548
 549   CubicPolyInterpolator* pCubicInterp = NULL;
 550   double* deltaFilteredProj = NULL;
 551   if (interpType == Backprojector::INTERP_LINEAR) {
 552     // precalculate scaled difference for linear interpolation
 553     deltaFilteredProj = new double [nDet];
 554 #if HAVE_OPENMP
 555   #pragma omp parallel for
 556 #endif
 557     for (int i = 0; i < nDet - 1; i++)
 558       deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
 559     deltaFilteredProj[nDet - 1] = 0;  // last detector
 560   } else if (interpType == Backprojector::INTERP_CUBIC) {
 561     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 562   }
 563
 564   int iLastDet = nDet - 1;
 565
 566 #if HAVE_OPENMP
 567   #pragma omp parallel for
 568 #endif
 569   for (int ix = 0; ix < nx; ix++) {
 570     double detPos = detPosColBase + (ix * det_dx);
 571     ImageFileColumn pImCol = v[ix];
 572
 573     for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
 574       if (interpType == Backprojector::INTERP_NEAREST) {
 575         int iDetPos = nearest<int> (detPos); // calc index in the filtered raysum vector
 576         if (iDetPos >= 0 && iDetPos < nDet) {
 577           *pImCol++ += filteredProj[iDetPos];
 578         }
 579       } else if (interpType == Backprojector::INTERP_LINEAR) {
 580         double detPosFloor = floor (detPos);
 581         int iDetPos = static_cast<int>(detPosFloor);
 582         double frac = detPos - detPosFloor;  // fraction distance from det
 583         if (iDetPos >= 0 && iDetPos <= iLastDet) {
 584           *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
 585         }
 586       } else if (interpType == Backprojector::INTERP_CUBIC) {
 587         double p = detPos;   // position along detector
 588         if (p >= 0 && p < nDet) {
 589           *pImCol++  += pCubicInterp->interpolate (p);
 590         }
 591       }
 592     }   // end for y
 593   }     // end for x
 594
 595   if (interpType == Backprojector::INTERP_LINEAR)
 596     delete deltaFilteredProj;
 597   else if (interpType == Backprojector::INTERP_CUBIC)
 598     delete pCubicInterp;
 599 }
 600
 601
 602 // CLASS IDENTICATION
 603 //   BackprojectIntDiff
 604 //
 605 // PURPOSE
 606 //   Highly optimized and integer version of BackprojectDiff
 607
 608 void
 609 BackprojectIntDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 610 {
 611   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
 612 #if SIZEOF_LONG == 4
 613   static const int scaleShift = 16;
 614 #elif SIZEOF_LONG == 8
 615   static const int scaleShift = 32;
 616 #endif
 617   static const long scale = (1L << scaleShift);
 618   static const long scaleBitmask = scale - 1;
 619   static const long halfScale = scale / 2;
 620   static const double dInvScale = 1. / scale;
 621
 622   const long det_dx = nearest<long> (xInc * cos (theta) / detInc * scale);
 623   const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
 624
 625   // calculate L for first point in image (0, 0)
 626   long detPosColBase = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 627
 628   double* deltaFilteredProj = NULL;
 629   CubicPolyInterpolator* pCubicInterp = NULL;
 630   if (interpType == Backprojector::INTERP_LINEAR) {
 631     // precalculate scaled difference for linear interpolation
 632     deltaFilteredProj = new double [nDet];
 633 #if HAVE_OPENMP
 634     #pragma omp parallel for
 635 #endif
 636     for (int i = 0; i < nDet - 1; i++)
 637       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 638     deltaFilteredProj[nDet - 1] = 0;  // last detector
 639   } else if (interpType == Backprojector::INTERP_CUBIC) {
 640     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 641   }
 642
 643   int iLastDet = nDet - 1;
 644 #if HAVE_OPENMP
 645   #pragma omp parallel for
 646 #endif
 647   for (int ix = 0; ix < nx; ix++) {
 648     long detPos = detPosColBase + (ix * det_dx);
 649     ImageFileColumn pImCol = v[ix];
 650
 651     for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
 652       if (interpType == Backprojector::INTERP_NEAREST) {
 653         const int iDetPos = (detPos + halfScale) >> scaleShift;
 654         if (iDetPos >= 0 && iDetPos <= iLastDet) {
 655           *pImCol++ += filteredProj[iDetPos];
 656         } else
 657           pImCol++;
 658       } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 659         const int iDetPos = ((detPos + halfScale) >> scaleShift) * m_interpFactor;
 660         if (iDetPos >= 0 && iDetPos <= iLastDet) {
 661           *pImCol++ += filteredProj[iDetPos];
 662         } else
 663           pImCol++;
 664       } else if (interpType == Backprojector::INTERP_LINEAR) {
 665         const long iDetPos = detPos >> scaleShift;
 666         if (iDetPos >= 0 && iDetPos <= iLastDet) {
 667           const long detRemainder = detPos & scaleBitmask;
 668           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 669         } else
 670           pImCol++;
 671       } else if (interpType == Backprojector::INTERP_CUBIC) {
 672         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(detPos) / scale);
 673       } // end Cubic
 674     } // end for iy
 675   } // end for ix
 676
 677   if (interpType == Backprojector::INTERP_LINEAR)
 678     delete deltaFilteredProj;
 679   else if (interpType == Backprojector::INTERP_CUBIC)
 680     delete pCubicInterp;
 681 }
 682
 683
 684 void
 685 BackprojectEquiangular::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 #if HAVE_OPENMP
 694   #pragma omp parallel for
 695 #endif
 696   for (int ix = 0; ix < nx; ix++) {
 697     ImageFileColumn pImCol = v[ix];
 698
 699     for (int iy = 0; iy < ny; iy++) {
 700       double dAngleDiff = beta - phi[ix][iy];
 701       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 702       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 703       double dFLPlusSin = m_dFocalLength + rsin_t;
 704       double gamma =  atan (rcos_t / dFLPlusSin);
 705       double dPos = gamma / detInc;  // position along detector
 706       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
 707
 708       if (interpType == Backprojector::INTERP_NEAREST) {
 709         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 710         if (iDetPos >= 0 && iDetPos < nDet)
 711           pImCol[iy] += filteredProj[iDetPos] / dL2;
 712       } else if (interpType == Backprojector::INTERP_LINEAR) {
 713         double dPosFloor = floor (dPos);
 714         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 715         double frac = dPos - dPosFloor; // fraction distance from det
 716         if (iDetPos >= 0 && iDetPos < nDet - 1)
 717           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
 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) / dL2;
 722       }
 723     }   // end for y
 724   }     // end for x
 725
 726   if (interpType == Backprojector::INTERP_CUBIC)
 727     delete pCubicInterp;
 728 }
 729
 730 void
 731 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 732 {
 733   double beta = view_angle;
 734
 735   CubicPolyInterpolator* pCubicInterp = NULL;
 736   if (interpType == Backprojector::INTERP_CUBIC)
 737     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 738
 739 #if HAVE_OPENMP
 740   #pragma omp parallel for
 741 #endif
 742   for (int ix = 0; ix < nx; ix++) {
 743     ImageFileColumn pImCol = v[ix];
 744
 745     for (int iy = 0; iy < ny; iy++) {
 746       double dAngleDiff = beta - phi[ix][iy];
 747       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 748       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 749
 750       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
 751       double dU2 = dU * dU;
 752
 753       double dDetPos =  rcos_t / dU;
 754       // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
 755       dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
 756       double dPos = dDetPos / detInc;  // position along detector array
 757
 758       if (interpType == Backprojector::INTERP_NEAREST) {
 759         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 760         if (iDetPos >= 0 && iDetPos < nDet)
 761           pImCol[iy] += filteredProj[iDetPos] / dU2;
 762       } else if (interpType == Backprojector::INTERP_LINEAR) {
 763         double dPosFloor = floor (dPos);
 764         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 765         double frac = dPos - dPosFloor; // fraction distance from det
 766         if (iDetPos >= 0 && iDetPos < nDet - 1)
 767           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2;
 768       } else if (interpType == Backprojector::INTERP_CUBIC) {
 769         double d = iDetCenter + dPos;           // position along detector
 770         if (d >= 0 && d < nDet)
 771           pImCol[iy] += pCubicInterp->interpolate (d) / dU2;
 772       }
 773     }   // end for y
 774   }     // end for x
 775
 776   if (interpType == Backprojector::INTERP_CUBIC)
 777     delete pCubicInterp;
 778 }
 779