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