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.25 2001/02/09 01:54:20 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_DIFF2 = 3;
  34 const int Backprojector::BPROJ_IDIFF2 = 4;
  35 const int Backprojector::BPROJ_IDIFF3 = 5;
  36
  37 const char* Backprojector::s_aszBackprojectName[] =
  38 {
  39   {"trig"},
  40   {"table"},
  41   {"diff"},
  42   {"diff2"},
  43   {"idiff2"},
  44   {"idiff3"},
  45 };
  46
  47 const char* Backprojector::s_aszBackprojectTitle[] =
  48 {
  49   {"Direct Trigometric"},
  50   {"Trigometric Table"},
  51   {"Difference Iteration"},
  52   {"Difference Iteration Optimized"},
  53   {"Integer Difference Iteration Optimized"},
  54   {"Integer Difference Iteration Highly-Optimized"},
  55 };
  56
  57 const int Backprojector::s_iBackprojectCount = sizeof(s_aszBackprojectName) / sizeof(const char*);
  58
  59 const int Backprojector::INTERP_INVALID = -1;
  60 const int Backprojector::INTERP_NEAREST = 0;
  61 const int Backprojector::INTERP_LINEAR = 1;
  62 const int Backprojector::INTERP_CUBIC = 2;
  63 const int Backprojector::INTERP_FREQ_PREINTERPOLATION = 3;
  64 #if HAVE_BSPLINE_INTERP
  65 const int Backprojector::INTERP_BSPLINE = 4;
  66 const int Backprojector::INTERP_1BSPLINE = 5;
  67 const int Backprojector::INTERP_2BSPLINE = 6;
  68 const int Backprojector::INTERP_3BSPLINE = 7;
  69 #endif
  70
  71 const char* Backprojector::s_aszInterpName[] =
  72 {
  73   {"nearest"},
  74   {"linear"},
  75   {"cubic"},
  76   {"freq_preinterpolationj"},
  77 #if HAVE_BSPLINE_INTERP
  78   {"bspline"},
  79   {"1bspline"},
  80   {"2bspline"},
  81   {"3bspline"},
  82 #endif
  83 };
  84
  85 const char* Backprojector::s_aszInterpTitle[] =
  86 {
  87   {"Nearest"},
  88   {"Linear"},
  89   {"Cubic"},
  90   {"Frequency Preinterpolation"},
  91 #if HAVE_BSPLINE_INTERP
  92   {"B-Spline"},
  93   {"B-Spline 1st Order"},
  94   {"B-Spline 2nd Order"},
  95   {"B-Spline 3rd Order"},
  96 #endif
  97 };
  98
  99 const int Backprojector::s_iInterpCount = sizeof(s_aszInterpName) / sizeof(const char*);
 100
 101
 102
 103 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 104 {
 105   m_fail = false;
 106   m_pBackprojectImplem = NULL;
 107
 108   initBackprojector (proj, im, backprojName, interpName, interpFactor);
 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 Backprojector::~Backprojector ()
 119 {
 120   delete m_pBackprojectImplem;
 121 }
 122
 123 // FUNCTION IDENTIFICATION
 124 //     Backproject* projector = selectBackprojector (...)
 125 //
 126 // PURPOSE
 127 //     Selects a backprojector based on BackprojType
 128 //     and initializes the backprojector
 129
 130 bool
 131 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName, const int interpFactor)
 132 {
 133   m_nameBackproject = backprojName;
 134   m_nameInterpolation = interpName;
 135   m_pBackprojectImplem = NULL;
 136   m_idBackproject = convertBackprojectNameToID (backprojName);
 137   if (m_idBackproject == BPROJ_INVALID) {
 138     m_fail = true;
 139     m_failMessage = "Invalid backprojection name ";
 140     m_failMessage += backprojName;
 141   }
 142   m_idInterpolation = convertInterpNameToID (interpName);
 143   if (m_idInterpolation == INTERP_INVALID) {
 144     m_fail = true;
 145     m_failMessage = "Invalid interpolation name ";
 146     m_failMessage += interpName;
 147   }
 148
 149   if (m_fail || m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
 150     m_fail = true;
 151     return false;
 152   }
 153
 154   if (proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
 155     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquilinear(proj, im, m_idInterpolation, interpFactor));
 156   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR)
 157     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectEquiangular(proj, im, m_idInterpolation, interpFactor));
 158   else if (proj.geometry() == Scanner::GEOMETRY_PARALLEL) {
 159     if (m_idBackproject == BPROJ_TRIG)
 160       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation, interpFactor));
 161     else if (m_idBackproject == BPROJ_TABLE)
 162       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation, interpFactor));
 163     else if (m_idBackproject == BPROJ_DIFF)
 164       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation, interpFactor));
 165     else if (m_idBackproject == BPROJ_DIFF2)
 166       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation, interpFactor));
 167     else if (m_idBackproject == BPROJ_IDIFF2)
 168       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation, interpFactor));
 169     else if (m_idBackproject == BPROJ_IDIFF3)
 170       m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff3 (proj, im, m_idInterpolation, interpFactor));
 171   } else {
 172     m_fail = true;
 173     m_failMessage = "Unable to select a backprojection method [Backprojector::initBackprojector]";
 174     return false;
 175   }
 176
 177   return true;
 178 }
 179
 180
 181 int
 182 Backprojector::convertBackprojectNameToID (const char* const backprojName)
 183 {
 184   int backprojID = BPROJ_INVALID;
 185
 186   for (int i = 0; i < s_iBackprojectCount; i++)
 187     if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
 188       backprojID = i;
 189       break;
 190     }
 191
 192     return (backprojID);
 193 }
 194
 195 const char*
 196 Backprojector::convertBackprojectIDToName (int bprojID)
 197 {
 198   static const char *bprojName = "";
 199
 200   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 201     return (s_aszBackprojectName[bprojID]);
 202
 203   return (bprojName);
 204 }
 205
 206 const char*
 207 Backprojector::convertBackprojectIDToTitle (const int bprojID)
 208 {
 209   static const char *bprojTitle = "";
 210
 211   if (bprojID >= 0 && bprojID < s_iBackprojectCount)
 212     return (s_aszBackprojectTitle[bprojID]);
 213
 214   return (bprojTitle);
 215 }
 216
 217
 218 int
 219 Backprojector::convertInterpNameToID (const char* const interpName)
 220 {
 221   int interpID = INTERP_INVALID;
 222
 223   for (int i = 0; i < s_iInterpCount; i++)
 224     if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
 225       interpID = i;
 226       break;
 227     }
 228
 229     return (interpID);
 230 }
 231
 232 const char*
 233 Backprojector::convertInterpIDToName (const int interpID)
 234 {
 235   static const char *interpName = "";
 236
 237   if (interpID >= 0 && interpID < s_iInterpCount)
 238     return (s_aszInterpName[interpID]);
 239
 240   return (interpName);
 241 }
 242
 243 const char*
 244 Backprojector::convertInterpIDToTitle (const int interpID)
 245 {
 246   static const char *interpTitle = "";
 247
 248   if (interpID >= 0 && interpID < s_iInterpCount)
 249     return (s_aszInterpTitle[interpID]);
 250
 251   return (interpTitle);
 252 }
 253
 254
 255
 256 // CLASS IDENTICATION
 257 //   Backproject
 258 //
 259 // PURPOSE
 260 //   Pure virtual base class for all backprojectors.
 261
 262 Backproject::Backproject (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
 263 : proj(proj), im(im), interpType(interpType), m_interpFactor(interpFactor)
 264 {
 265   detInc = proj.detInc();
 266   nDet = proj.nDet();
 267   iDetCenter = (nDet - 1) / 2;  // index refering to L=0 projection
 268   rotScale = proj.rotInc();
 269
 270   if (proj.geometry() == Scanner::GEOMETRY_PARALLEL)
 271     rotScale /= (proj.nView() * proj.rotInc() / PI); // scale by number of PI rotations
 272   else if (proj.geometry() == Scanner::GEOMETRY_EQUIANGULAR || proj.geometry() == Scanner::GEOMETRY_EQUILINEAR)
 273     rotScale /= (proj.nView() * proj.rotInc() / (2 * PI)); // scale by number of 2PI rotations
 274   else
 275     sys_error (ERR_SEVERE, "Invalid geometry type %d [Backproject::Backproject]", proj.geometry());
 276
 277   v = im.getArray();
 278   nx = im.nx();
 279   ny = im.ny();
 280   im.arrayDataClear();
 281
 282   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
 283   xMax = xMin + proj.phmLen();
 284   yMin = -proj.phmLen() / 2;
 285   yMax = yMin + proj.phmLen();
 286
 287   xInc = (xMax - xMin) / nx;    // size of cells
 288   yInc = (yMax - yMin) / ny;
 289
 290   m_dFocalLength = proj.focalLength();
 291 }
 292
 293 Backproject::~Backproject ()
 294 {}
 295
 296 void
 297 Backproject::ScaleImageByRotIncrement ()
 298 {
 299   for (int ix = 0; ix < nx; ix++)
 300     for (int iy = 0; iy < ny; iy++)
 301       v[ix][iy] *= rotScale;
 302 }
 303
 304 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 305 {
 306   sys_error (ERR_WARNING, "r=%f, phi=%f", r, phi);
 307   errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 308 }
 309
 310 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
 311 {
 312 #if 1
 313   std::ostringstream os;
 314   os << "ix=" << ix << ", iy=" << iy << ", theta=" << theta << ", L=" << L << ", detinc=" << detInc << "\n";
 315   os << "ndet=" << nDet << ", detInc=" << detInc << ", iDetCenter=" << iDetCenter << "\n";
 316   os << "xMin=" << xMin << ", xMax=" << xMax << ", xInc=" << xInc << "\n";
 317   os << "yMin=" << yMin << ", yMax=" << yMax << ", yInc=" << yInc << "\n";
 318   os << "iDetPos index outside bounds: " << iDetPos << " [backprojector]";;
 319
 320   sys_error (ERR_WARNING, os.str().c_str());
 321 #endif
 322 }
 323
 324
 325 // CLASS IDENTICATION
 326 //   BackprojectTrig
 327 //
 328 // PURPOSE
 329 //   Uses trigometric functions at each point in image for backprojection.
 330
 331 void
 332 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 333 {
 334   double theta = view_angle;
 335
 336   CubicInterpolator* pCubicInterp = NULL;
 337   if (interpType == Backprojector::INTERP_CUBIC)
 338     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 339
 340   double x = xMin + xInc / 2;   // Rectang coords of center of pixel
 341   for (int ix = 0; ix < nx; x += xInc, ix++) {
 342     double y = yMin + yInc / 2;
 343     for (int iy = 0; iy < ny; y += yInc, iy++) {
 344       double r = sqrt (x * x + y * y);   // distance of cell from center
 345       double phi = atan2 (y, x);         // angle of cell from center
 346       double L = r * cos (theta - phi);  // position on detector
 347
 348       if (interpType == Backprojector::INTERP_NEAREST) {
 349         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 350
 351         if (iDetPos >= 0 && iDetPos < nDet)
 352           v[ix][iy] += rotScale * filteredProj[iDetPos];
 353       } else if (interpType == Backprojector::INTERP_LINEAR) {
 354         double p = L / detInc;  // position along detector
 355         double pFloor = floor (p);
 356         int iDetPos = iDetCenter + static_cast<int>(pFloor);
 357         double frac = p - pFloor;       // fraction distance from det
 358         if (iDetPos >= 0 && iDetPos < nDet - 1)
 359           v[ix][iy] += rotScale * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 360       } else if (interpType = Backprojector::INTERP_CUBIC) {
 361         double p = iDetCenter + (L / detInc);   // position along detector
 362         if (p >= 0 && p < nDet)
 363           v[ix][iy] += rotScale * pCubicInterp->interpolate (p);
 364       }
 365     }
 366   }
 367
 368   if (interpType == Backprojector::INTERP_CUBIC)
 369     delete pCubicInterp;
 370 }
 371
 372
 373 // CLASS IDENTICATION
 374 //   BackprojectTable
 375 //
 376 // PURPOSE
 377 //   Precalculates trigometric function value for each point in image for backprojection.
 378
 379 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
 380 : Backproject (proj, im, interpType, interpFactor)
 381 {
 382   arrayR.initSetSize (im.nx(), im.ny());
 383   arrayPhi.initSetSize (im.nx(), im.ny());
 384   r = arrayR.getArray();
 385   phi = arrayPhi.getArray();
 386
 387   double x, y;                  // Rectang coords of center of pixel
 388   int ix, iy;
 389   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 390     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 391       r[ix][iy] = sqrt (x * x + y * y);
 392       phi[ix][iy] = atan2 (y, x);
 393     }
 394 }
 395
 396 BackprojectTable::~BackprojectTable ()
 397 {
 398   ScaleImageByRotIncrement();
 399 }
 400
 401 void
 402 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 403 {
 404   double theta = view_angle;
 405
 406   CubicInterpolator* pCubicInterp = NULL;
 407   if (interpType == Backprojector::INTERP_CUBIC)
 408     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 409
 410   for (int ix = 0; ix < nx; ix++) {
 411     ImageFileColumn pImCol = v[ix];
 412
 413     for (int iy = 0; iy < ny; iy++) {
 414       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 415
 416       if (interpType == Backprojector::INTERP_NEAREST) {
 417         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 418
 419         if (iDetPos >= 0 && iDetPos < nDet)
 420           pImCol[iy] += filteredProj[iDetPos];
 421       } else if (interpType == Backprojector::INTERP_LINEAR) {
 422         double dPos = L / detInc;               // position along detector
 423         double dPosFloor = floor (dPos);
 424         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 425         double frac = dPos - dPosFloor; // fraction distance from det
 426         if (iDetPos >= 0 && iDetPos < nDet - 1)
 427           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 428       } else if (interpType = Backprojector::INTERP_CUBIC) {
 429         double p = iDetCenter + (L / detInc);   // position along detector
 430         if (p >= 0 && p < nDet)
 431           pImCol[iy] += pCubicInterp->interpolate (p);
 432       }
 433     }   // end for y
 434   }     // end for x
 435
 436   if (interpType == Backprojector::INTERP_CUBIC)
 437     delete pCubicInterp;
 438 }
 439
 440
 441 // CLASS IDENTICATION
 442 //   BackprojectDiff
 443 //
 444 // PURPOSE
 445 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 446 //   Iterates in x & y direction by adding difference in L position
 447
 448 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, int interpType, const int interpFactor)
 449 :  Backproject (proj, im, interpType, interpFactor)
 450 {
 451   // calculate center of first pixel v[0][0]
 452   double x = xMin + xInc / 2;
 453   double y = yMin + yInc / 2;
 454   start_r = sqrt (x * x + y * y);
 455   start_phi = atan2 (y, x);
 456
 457   im.arrayDataClear();
 458 }
 459
 460 BackprojectDiff::~BackprojectDiff()
 461 {
 462   ScaleImageByRotIncrement();
 463 }
 464
 465 void
 466 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 467 {
 468   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
 469   double det_dx = xInc * cos (theta);
 470   double det_dy = yInc * sin (theta);
 471   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
 472
 473   CubicInterpolator* pCubicInterp = NULL;
 474   if (interpType == Backprojector::INTERP_CUBIC)
 475     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 476
 477   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
 478     double curDetPos = lColStart;
 479     ImageFileColumn pImCol = v[ix];
 480
 481     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 482 #ifdef DEBUG
 483       printf ("[%2d,%2d]:  %8.5f  ", ix, iy, curDetPos);
 484 #endif
 485       if (interpType == Backprojector::INTERP_NEAREST) {
 486         int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);    // calc index in the filtered raysum vector
 487
 488         if (iDetPos >= 0 && iDetPos < nDet)
 489           pImCol[iy] += filteredProj[iDetPos];
 490       } else if (interpType == Backprojector::INTERP_LINEAR) {
 491         double detPos = curDetPos / detInc;             // position along detector
 492         double detPosFloor = floor (detPos);
 493         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 494         double frac = detPos - detPosFloor;     // fraction distance from det
 495         if (iDetPos >= 0 && iDetPos < nDet - 1)
 496           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 497       } else if (interpType = Backprojector::INTERP_CUBIC) {
 498         double p = iDetCenter + (curDetPos / detInc);   // position along detector
 499         if (p >= 0 && p < nDet)
 500           pImCol[iy] += pCubicInterp->interpolate (p);
 501       }
 502     }   // end for y
 503   }     // end for x
 504
 505   if (interpType == Backprojector::INTERP_CUBIC)
 506     delete pCubicInterp;
 507 }
 508
 509
 510 // CLASS IDENTICATION
 511 //   BackprojectDiff2
 512 //
 513 // PURPOSE
 514 //   Optimized version of BackprojectDiff
 515
 516 void
 517 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 518 {
 519   double theta = view_angle;
 520
 521   CubicInterpolator* pCubicInterp = NULL;
 522   if (interpType == Backprojector::INTERP_CUBIC)
 523     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 524
 525   // Distance betw. detectors for an angle given in units of detectors
 526   double det_dx = xInc * cos (theta) / detInc;
 527   double det_dy = yInc * sin (theta) / detInc;
 528
 529   // calculate detPosition for first point in image (ix=0, iy=0)
 530   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 531
 532 #ifdef DEBUG
 533   printf ("start_r=%8.5f, start_phi=%8.5f, rotScale=%8.5f\n", start_r, start_phi, rotScale);
 534 #endif
 535   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 536     double curDetPos = detPosColStart;
 537     ImageFileColumn pImCol = v[ix];
 538
 539     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 540 #ifdef DEBUG
 541       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(curDetPos)]);
 542 #endif
 543       if (interpType == Backprojector::INTERP_NEAREST) {
 544         int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector
 545
 546         if (iDetPos >= 0 && iDetPos < nDet)
 547           *pImCol++ += filteredProj[iDetPos];
 548       } else if (interpType == Backprojector::INTERP_LINEAR) {
 549         double detPosFloor = floor (curDetPos);
 550         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 551         double frac = curDetPos - detPosFloor;  // fraction distance from det
 552         if (iDetPos > 0 && iDetPos < nDet - 1)
 553           *pImCol++ += filteredProj[iDetPos] + (frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]));
 554       } else if (interpType = Backprojector::INTERP_CUBIC) {
 555         double p = iDetCenter + curDetPos;      // position along detector
 556         if (p >= 0 && p < nDet)
 557           *pImCol++  += pCubicInterp->interpolate (p);
 558       }
 559     }   // end for y
 560   }     // end for x
 561
 562   if (interpType == Backprojector::INTERP_CUBIC)
 563     delete pCubicInterp;
 564 }
 565
 566 // CLASS IDENTICATION
 567 //   BackprojectIntDiff2
 568 //
 569 // PURPOSE
 570 //   Integer version of BackprojectDiff2
 571
 572 void
 573 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 574 {
 575   double theta = view_angle;
 576
 577   CubicInterpolator* pCubicInterp = NULL;
 578   if (interpType == Backprojector::INTERP_CUBIC)
 579     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 580
 581   static const kint32 scale = 1 << 16;
 582   static const double dScale = scale;
 583   static const kint32 halfScale = scale / 2;
 584
 585   const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
 586   const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
 587
 588   // calculate L for first point in image (0, 0)
 589   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
 590
 591   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 592     kint32 curDetPos = detPosColStart;
 593     ImageFileColumn pImCol = v[ix];
 594
 595     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 596       if (interpType == Backprojector::INTERP_NEAREST) {
 597         int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
 598         int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector
 599
 600         if (iDetPos < 0 || iDetPos >= nDet)  // check for index outside of raysum pos
 601           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 602         else
 603           *pImCol++ += filteredProj[iDetPos];
 604       } else if (interpType == Backprojector::INTERP_LINEAR) {
 605         kint32 detPosFloor = curDetPos / scale;
 606         kint32 detPosRemainder = curDetPos % scale;
 607         if (detPosRemainder < 0) {
 608           detPosFloor--;
 609           detPosRemainder += scale;
 610         }
 611         int iDetPos = iDetCenter + detPosFloor;
 612         double frac = detPosRemainder / dScale;
 613         if (iDetPos < 0 || iDetPos >= nDet - 1)
 614           errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 615         else
 616           *pImCol++ += ((1.-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 617       } else if (interpType = Backprojector::INTERP_CUBIC) {
 618         double p = iDetCenter + (static_cast<double>(curDetPos) / scale);       // position along detector
 619         if (p >= 0 && p < nDet)
 620           *pImCol++  += pCubicInterp->interpolate (p);
 621       }
 622     }   // end for y
 623   }     // end for x
 624
 625   if (interpType == Backprojector::INTERP_CUBIC)
 626     delete pCubicInterp;
 627 }
 628
 629 // CLASS IDENTICATION
 630 //   BackprojectIntDiff3
 631 //
 632 // PURPOSE
 633 //   Highly optimized version of BackprojectIntDiff2
 634
 635 void
 636 BackprojectIntDiff3::BackprojectView (const double* const filteredProj, const double view_angle)
 637 {
 638   double theta = view_angle;  // add half PI to view angle to get perpendicular theta angle
 639   static const int scaleShift = 16;
 640   static const kint32 scale = (1 << scaleShift);
 641   static const kint32 scaleBitmask = scale - 1;
 642   static const kint32 halfScale = scale / 2;
 643   static const double dInvScale = 1. / scale;
 644
 645   const kint32 det_dx = nearest<kint32> (xInc * cos (theta) / detInc * scale);
 646   const kint32 det_dy = nearest<kint32> (yInc * sin (theta) / detInc * scale);
 647
 648   // calculate L for first point in image (0, 0)
 649   kint32 detPosColStart = nearest<kint32> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 650
 651   double* deltaFilteredProj = NULL;
 652   CubicInterpolator* pCubicInterp = NULL;
 653   if (interpType == Backprojector::INTERP_LINEAR) {
 654     // precalculate scaled difference for linear interpolation
 655     deltaFilteredProj = new double [nDet];
 656     for (int i = 0; i < nDet - 1; i++)
 657       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
 658     deltaFilteredProj[nDet - 1] = 0;  // last detector
 659   } else if (interpType == Backprojector::INTERP_CUBIC) {
 660     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 661   }
 662
 663   int iLastDet = nDet - 1;
 664   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 665     kint32 curDetPos = detPosColStart;
 666     ImageFileColumn pImCol = v[ix];
 667
 668     if (interpType == Backprojector::INTERP_NEAREST) {
 669       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 670         const int iDetPos = (curDetPos + halfScale) >> 16;
 671         if (iDetPos >= 0 && iDetPos <= iLastDet)
 672           *pImCol++ += filteredProj[iDetPos];
 673       } // end for iy
 674     } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
 675       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 676         const int iDetPos = ((curDetPos + halfScale) >> 16) * m_interpFactor;
 677         if (iDetPos >= 0 && iDetPos <= iLastDet)
 678           *pImCol++ += filteredProj[iDetPos];
 679       } // end for iy
 680     } else if (interpType == Backprojector::INTERP_LINEAR) {
 681       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 682         const kint32 iDetPos = curDetPos >> scaleShift;
 683         const kint32 detRemainder = curDetPos & scaleBitmask;
 684         if (iDetPos >= 0 && iDetPos <= iLastDet)
 685           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
 686       } // end for iy
 687     } else if (interpType = Backprojector::INTERP_CUBIC) {
 688       for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 689         *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / 65536);
 690       }
 691     } // end cubic
 692   } // end for ix
 693
 694   if (interpType == Backprojector::INTERP_LINEAR)
 695     delete deltaFilteredProj;
 696   else if (interpType == Backprojector::INTERP_CUBIC)
 697     delete pCubicInterp;
 698 }
 699
 700
 701 void
 702 BackprojectEquiangular::BackprojectView (const double* const filteredProj, const double view_angle)
 703 {
 704   double beta = view_angle;
 705
 706   CubicInterpolator* pCubicInterp = NULL;
 707   if (interpType == Backprojector::INTERP_CUBIC)
 708     pCubicInterp = new CubicInterpolator (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       double dFLPlusSin = m_dFocalLength + rsin_t;
 718       double gamma =  atan (rcos_t / dFLPlusSin);
 719       double dPos = gamma / detInc;  // position along detector
 720       double dL2 = dFLPlusSin * dFLPlusSin + (rcos_t * rcos_t);
 721
 722       if (interpType == Backprojector::INTERP_NEAREST) {
 723         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 724         if (iDetPos >= 0 && iDetPos < nDet)
 725           pImCol[iy] += filteredProj[iDetPos] / dL2;
 726       } else if (interpType == Backprojector::INTERP_LINEAR) {
 727         double dPosFloor = floor (dPos);
 728         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 729         double frac = dPos - dPosFloor; // fraction distance from det
 730         if (iDetPos >= 0 && iDetPos < nDet - 1)
 731           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dL2;
 732       } else if (interpType == Backprojector::INTERP_CUBIC) {
 733         double d = iDetCenter + dPos;           // position along detector
 734         if (d >= 0 && d < nDet)
 735           pImCol[iy] += pCubicInterp->interpolate (d) / dL2;
 736       }
 737     }   // end for y
 738   }     // end for x
 739
 740   if (interpType == Backprojector::INTERP_CUBIC)
 741     delete pCubicInterp;
 742 }
 743
 744 void
 745 BackprojectEquilinear::BackprojectView (const double* const filteredProj, const double view_angle)
 746 {
 747   double beta = view_angle;
 748
 749   CubicInterpolator* pCubicInterp = NULL;
 750   if (interpType == Backprojector::INTERP_CUBIC)
 751     pCubicInterp = new CubicInterpolator (filteredProj, nDet);
 752
 753   for (int ix = 0; ix < nx; ix++) {
 754     ImageFileColumn pImCol = v[ix];
 755
 756     for (int iy = 0; iy < ny; iy++) {
 757       double dAngleDiff = beta - phi[ix][iy];
 758       double rcos_t = r[ix][iy] * cos (dAngleDiff);
 759       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 760
 761       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
 762       double dDetPos =  rcos_t / dU;
 763       // double to scale for imaginary detector that passes through origin
 764       // of phantom, see Kak-Slaney Figure 3.22. This assumes that the detector is also
 765       // located focal-length away from the origin.
 766       dDetPos *= 2;
 767       double dPos = dDetPos / detInc;  // position along detector array
 768
 769       if (interpType == Backprojector::INTERP_NEAREST) {
 770         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
 771         if (iDetPos >= 0 && iDetPos < nDet)
 772           pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
 773       } else if (interpType == Backprojector::INTERP_LINEAR) {
 774         double dPosFloor = floor (dPos);
 775         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 776         double frac = dPos - dPosFloor; // fraction distance from det
 777         if (iDetPos >= 0 && iDetPos < nDet - 1)
 778           pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
 779                            / (dU * dU);
 780       } else if (interpType == Backprojector::INTERP_CUBIC) {
 781         double d = iDetCenter + dPos;           // position along detector
 782         if (d >= 0 && d < nDet)
 783           pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
 784       }
 785     }   // end for y
 786   }     // end for x
 787
 788   if (interpType == Backprojector::INTERP_CUBIC)
 789     delete pCubicInterp;
 790 }
 791