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