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.3 2000/06/22 10:17:28 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 Backprojector::Backprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
  30 {
  31   m_fail = false;
  32   m_pBackprojectImplem = NULL;
  33
  34   initBackprojector (proj, im, backprojName, interpName);
  35 }
  36
  37 void
  38 Backprojector::BackprojectView (const double* const viewData, const double viewAngle)
  39 {
  40   if (m_pBackprojectImplem)
  41     m_pBackprojectImplem->BackprojectView (viewData, viewAngle);
  42 }
  43
  44 Backprojector::~Backprojector (void)
  45 {
  46   delete m_pBackprojectImplem;
  47 }
  48
  49 // FUNCTION IDENTIFICATION
  50 //     Backproject* projector = selectBackprojector (...)
  51 //
  52 // PURPOSE
  53 //     Selects a backprojector based on BackprojType
  54 //     and initializes the backprojector
  55
  56 bool
  57 Backprojector::initBackprojector (const Projections& proj, ImageFile& im, const char* const backprojName, const char* const interpName)
  58 {
  59   m_nameBackproject = backprojName;
  60   m_nameInterpolation = interpName;
  61   m_idBackproject = convertBackprojectNameToID (backprojName);
  62   m_idInterpolation = convertInterpolationNameToID (interpName);
  63   m_pBackprojectImplem = NULL;
  64
  65   if (m_idBackproject == BPROJ_INVALID || m_idInterpolation == INTERP_INVALID) {
  66     m_fail = true;
  67     return false;
  68   }
  69
  70   if (m_idBackproject == BPROJ_TRIG)
  71     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTrig (proj, im, m_idInterpolation));
  72   else if (m_idBackproject == BPROJ_TABLE)
  73     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectTable (proj, im, m_idInterpolation));
  74   else if (m_idBackproject == BPROJ_DIFF)
  75     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff (proj, im, m_idInterpolation));
  76   else if (m_idBackproject == BPROJ_DIFF2)
  77     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectDiff2 (proj, im, m_idInterpolation));
  78   else if (m_idBackproject == BPROJ_IDIFF2)
  79     m_pBackprojectImplem = static_cast<Backproject*>(new BackprojectIntDiff2 (proj, im, m_idInterpolation));
  80   else {
  81     m_fail = true;
  82     return false;
  83   }
  84
  85   return true;
  86 }
  87
  88
  89 const Backprojector::BackprojectID
  90 Backprojector::convertBackprojectNameToID (const char* const backprojName)
  91 {
  92   BackprojectID backprojID = BPROJ_INVALID;
  93
  94   if (strcasecmp (backprojName, BPROJ_TRIG_STR) == 0)
  95     backprojID = BPROJ_TRIG;
  96   else if (strcasecmp (backprojName, BPROJ_TABLE_STR) == 0)
  97     backprojID = BPROJ_TABLE;
  98   else if (strcasecmp (backprojName, BPROJ_DIFF_STR) == 0)
  99     backprojID = BPROJ_DIFF;
 100   else if (strcasecmp (backprojName, BPROJ_DIFF2_STR) == 0)
 101     backprojID = BPROJ_DIFF2;
 102   else if (strcasecmp (backprojName, BPROJ_IDIFF2_STR) == 0)
 103     backprojID = BPROJ_IDIFF2;
 104
 105   return (backprojID);
 106 }
 107
 108 const char*
 109 Backprojector::convertBackprojectIDToName (const BackprojectID bprojID)
 110 {
 111   const char *bprojName = "";
 112
 113   if (bprojID == BPROJ_TRIG)
 114     bprojName = BPROJ_TRIG_STR;
 115   else if (bprojID == BPROJ_TABLE)
 116     bprojName = BPROJ_TABLE_STR;
 117   else if (bprojID == BPROJ_DIFF)
 118     bprojName = BPROJ_DIFF_STR;
 119   else if (bprojID == BPROJ_DIFF2)
 120     bprojName = BPROJ_DIFF2_STR;
 121   else if (bprojID == BPROJ_IDIFF2)
 122     bprojName = BPROJ_IDIFF2_STR;
 123
 124   return (bprojName);
 125 }
 126
 127
 128
 129 const Backprojector::InterpolationID
 130 Backprojector::convertInterpolationNameToID (const char* const interpName)
 131 {
 132   InterpolationID interpID = INTERP_INVALID;
 133
 134   if (strcasecmp (interpName, INTERP_NEAREST_STR) == 0)
 135     interpID = INTERP_NEAREST;
 136   else if (strcasecmp (interpName, INTERP_LINEAR_STR) == 0)
 137     interpID = INTERP_LINEAR;
 138 #if HAVE_BSPLINE_INTERP
 139   else if (strcasecmp (interpName, INTERP_BSPLINE_STR) == 0)
 140     interpID = INTERP_BSPLINE;
 141 #endif
 142
 143   return (interpID);
 144 }
 145
 146
 147 /* NAME
 148  *      name_of_interp                  Return name of interpolation method
 149  *
 150  * SYNOPSIS
 151  *      name = name_of_interp (interp_type)
 152  *      char *name                      Name of interpolation method
 153  *      int interp_type                 Method of interpolation
 154  *
 155  * NOTES
 156  *      Returns NULL if interp_type is invalid
 157  */
 158
 159 const char*
 160 Backprojector::convertInterpolationIDToName (const InterpolationID interpID)
 161 {
 162   if (interpID == INTERP_NEAREST)
 163     return (INTERP_NEAREST_STR);
 164   else if (interpID == INTERP_LINEAR)
 165     return (INTERP_LINEAR_STR);
 166 #if HAVE_BSPLINE_INTERP
 167   else if (interpID == INTERP_BSPLINE)
 168     return (INTERP_BSPLINE_STR);
 169 #endif
 170   else
 171     return ("");
 172 }
 173
 174
 175 // CLASS IDENTICATION
 176 //   Backproject
 177 //
 178 // PURPOSE
 179 //   Pure virtual base class for all backprojectors.
 180
 181 Backproject::Backproject (const Projections& proj, ImageFile& im, const Backprojector::InterpolationID interpType)
 182   : proj(proj), im(im), interpType(interpType)
 183 {
 184   detInc = proj.detInc();
 185   nDet = proj.nDet();
 186   iDetCenter = (nDet - 1) / 2;  // index refering to L=0 projection
 187   rotInc = proj.rotInc();
 188
 189   v = im.getArray();
 190   nx = im.nx();
 191   ny = im.ny();
 192   im.arrayDataClear();
 193
 194   xMin = -proj.phmLen() / 2;      // Retangular coords of phantom
 195   xMax = xMin + proj.phmLen();
 196   yMin = -proj.phmLen() / 2;
 197   yMax = yMin + proj.phmLen();
 198
 199   xInc = (xMax - xMin) / nx;    // size of cells
 200   yInc = (yMax - yMin) / ny;
 201
 202   if (interpType != Backprojector::INTERP_NEAREST && interpType != Backprojector::INTERP_LINEAR)
 203     sys_error (ERR_WARNING, "Illegal interpType %d [selectBackprojector]", interpType);
 204 }
 205
 206 Backproject::~Backproject (void)
 207 {}
 208
 209 void
 210 Backproject::ScaleImageByRotIncrement (void)
 211 {
 212   for (int ix = 0; ix < nx; ix++)
 213     for (int iy = 0; iy < ny; iy++)
 214       v[ix][iy] *= rotInc;
 215 }
 216
 217 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double r, double phi, double L, int iDetPos)
 218 {
 219     printf ("r=%f, phi=%f\n", r, phi);
 220     errorIndexOutsideDetector (ix, iy, theta, L, iDetPos);
 221 }
 222
 223 void Backproject::errorIndexOutsideDetector (int ix, int iy, double theta, double L, int iDetPos)
 224 {
 225     printf ("ix=%d, iy=%d\n", ix, iy);
 226     printf ("theta=%f, L=%f, detInc=%f\n", theta, L, detInc);
 227     printf ("proj.ndet=%d, proj.detInc=%.4f, iDetCenter=%d\n", nDet, detInc, iDetCenter);
 228     printf ("xMin=%15.8f, xMax=%15.8f, xInc=%15.8f\n", xMin, xMax, xInc);
 229     printf ("yMin=%15.8f, yMax=%15.8f, yInc=%15.8f\n", yMin, yMax, yInc);
 230     sys_error (ERR_WARNING, "iDetPos index outside bounds: %d [backprojector]", iDetPos);
 231 }
 232
 233
 234 // CLASS IDENTICATION
 235 //   BackprojectTrig
 236 //
 237 // PURPOSE
 238 //   Uses trigometric functions at each point in image for backprojection.
 239
 240 void
 241 BackprojectTrig::BackprojectView (const double* const filteredProj, const double view_angle)
 242 {
 243   double theta = HALFPI + view_angle;   // Add PI/2 to get perpendicular angle to detector
 244   int ix, iy;
 245   double x, y;                  // Rectang coords of center of pixel
 246
 247   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 248     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 249       double r = sqrt (x * x + y * y);   // distance of cell from center
 250       double phi = atan2 (y, x);         // angle of cell from center
 251       double L = r * cos (theta - phi);  // position on detector
 252
 253       if (interpType == Backprojector::INTERP_NEAREST) {
 254         int iDetPos = iDetCenter + nearest<int> (L / detInc); // calc'd index in the filter raysum array
 255
 256         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 257             errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 258         else
 259           v[ix][iy] += rotInc * filteredProj[iDetPos];
 260       } else if (interpType == Backprojector::INTERP_LINEAR) {
 261           double p = L / detInc;        // position along detector
 262           double pFloor = floor (p);
 263           int iDetPos = iDetCenter + static_cast<int>(pFloor);
 264           double frac = p - pFloor;     // fraction distance from det
 265           if (iDetPos < 0 || iDetPos >= nDet - 1)       // check for impossible: index outside of raysum pos
 266             errorIndexOutsideDetector (ix, iy, theta, r, phi, L, iDetPos);
 267           else
 268             v[ix][iy] += rotInc * ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 269       }
 270     }
 271 }
 272
 273
 274 // CLASS IDENTICATION
 275 //   BackprojectTable
 276 //
 277 // PURPOSE
 278 //   Precalculates trigometric function value for each point in image for backprojection.
 279
 280 BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
 281   : Backproject::Backproject (proj, im, interpType)
 282 {
 283   arrayR.initSetSize (nx, ny);
 284   arrayPhi.initSetSize (nx, ny);
 285   r = arrayR.getArray();
 286   phi = arrayPhi.getArray();
 287
 288   double x, y;                  // Rectang coords of center of pixel
 289   int ix, iy;
 290   for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
 291     for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
 292       r[ix][iy] = sqrt (x * x + y * y);
 293       phi[ix][iy] = atan2 (y, x);
 294     }
 295 }
 296
 297 BackprojectTable::~BackprojectTable (void)
 298 {
 299   ScaleImageByRotIncrement();
 300 }
 301
 302 void
 303 BackprojectTable::BackprojectView (const double* const filteredProj, const double view_angle)
 304 {
 305   double theta = HALFPI + view_angle;  // add half PI to view angle to get perpendicular theta angle
 306
 307   for (int ix = 0; ix < nx; ix++) {
 308     ImageFileColumn pImCol = v[ix];
 309
 310     for (int iy = 0; iy < ny; iy++) {
 311       double L = r[ix][iy] * cos (theta - phi[ix][iy]);
 312
 313       if (interpType == Backprojector::INTERP_NEAREST) {
 314         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 315
 316         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 317           errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
 318         else
 319           pImCol[iy] += filteredProj[iDetPos];
 320       } else if (interpType == Backprojector::INTERP_LINEAR) {
 321         double dPos = L / detInc;               // position along detector
 322         double dPosFloor = floor (dPos);
 323         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
 324         double frac = dPos - dPosFloor; // fraction distance from det
 325         if (iDetPos < 0 || iDetPos >= nDet - 1)
 326             errorIndexOutsideDetector (ix, iy, theta, r[ix][iy], phi[ix][iy], L, iDetPos);
 327         else
 328           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 329       }
 330     }   // end for y
 331   }     // end for x
 332 }
 333
 334
 335 // CLASS IDENTICATION
 336 //   BackprojectDiff
 337 //
 338 // PURPOSE
 339 //   Backprojects by precalculating the change in L position for each x & y step in the image.
 340 //   Iterates in x & y direction by adding difference in L position
 341
 342 BackprojectDiff::BackprojectDiff (const Projections& proj, ImageFile& im, Backprojector::InterpolationID interpType)
 343   :  Backproject::Backproject (proj, im, interpType)
 344 {
 345   // calculate center of first pixel v[0][0]
 346   double x = xMin + xInc / 2;
 347   double y = yMin + yInc / 2;
 348   start_r = sqrt (x * x + y * y);
 349   start_phi = atan2 (y, x);
 350
 351   im.arrayDataClear();
 352 }
 353
 354 BackprojectDiff::~BackprojectDiff()
 355 {
 356   ScaleImageByRotIncrement();
 357 }
 358
 359 void
 360 BackprojectDiff::BackprojectView (const double* const filteredProj, const double view_angle)
 361 {
 362   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 363   double det_dx = xInc * sin (theta);
 364   double det_dy = yInc * cos (theta);
 365   double lColStart = start_r * cos (theta - start_phi);  // calculate L for first point in image
 366
 367   for (int ix = 0; ix < nx; ix++, lColStart += det_dx) {
 368     double curDetPos = lColStart;
 369     ImageFileColumn pImCol = v[ix];
 370
 371     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 372 #ifdef DEBUG
 373       printf ("[%2d,%2d]:  %8.5lf  ", ix, iy, curDetPos);
 374 #endif
 375       if (interpType == Backprojector::INTERP_NEAREST) {
 376         int iDetPos = iDetCenter + nearest<int>(curDetPos / detInc);    // calc index in the filtered raysum vector
 377
 378         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 379             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 380         else
 381           pImCol[iy] += filteredProj[iDetPos];
 382       } else if (interpType == Backprojector::INTERP_LINEAR) {
 383         double detPos = curDetPos / detInc;             // position along detector
 384         double detPosFloor = floor (detPos);
 385         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 386         double frac = detPos - detPosFloor;     // fraction distance from det
 387         if (iDetPos < 0 || iDetPos >= nDet - 1)
 388             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 389         else
 390           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 391       }
 392     }   // end for y
 393   }     // end for x
 394 }
 395
 396
 397 // CLASS IDENTICATION
 398 //   BackprojectDiff2
 399 //
 400 // PURPOSE
 401 //   Optimized version of BackprojectDiff
 402
 403 void
 404 BackprojectDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 405 {
 406   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 407
 408   // Distance betw. detectors for an angle given in units of detectors
 409   double det_dx = xInc * sin (theta) / detInc;
 410   double det_dy = yInc * cos (theta) / detInc;
 411
 412   // calculate detPosition for first point in image (ix=0, iy=0)
 413   double detPosColStart = start_r * cos (theta - start_phi) / detInc;
 414
 415 #ifdef DEBUG
 416   printf ("start_r=%8.5f, start_phi=%8.5f, rotInc=%8.5f\n", start_r, start_phi, rotInc);
 417 #endif
 418   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 419     double curDetPos = detPosColStart;
 420     ImageFileColumn pImCol = v[ix];
 421
 422     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 423 #ifdef DEBUG
 424       printf ("[%2d,%2d]: %8.5f %8.5f\n", ix, iy, curDetPos, filteredProj[iDetCenter + nearest<int>(L))]);
 425 #endif
 426       if (interpType == Backprojector::INTERP_NEAREST) {
 427         int iDetPos = iDetCenter + nearest<int> (curDetPos);    // calc index in the filtered raysum vector
 428
 429         if (iDetPos < 0 || iDetPos >= nDet)     // check for impossible: index outside of raysum pos
 430             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 431         else
 432           *pImCol++ += filteredProj[iDetPos];
 433       } else if (interpType == Backprojector::INTERP_LINEAR) {
 434         double detPosFloor = floor (curDetPos);
 435         int iDetPos = iDetCenter + static_cast<int>(detPosFloor);
 436         double frac = curDetPos - detPosFloor;  // fraction distance from det
 437         if (iDetPos < 0 || iDetPos >= nDet - 1)
 438             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 439         else
 440           *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 441       }
 442     }   // end for y
 443   }     // end for x
 444 }
 445
 446 // CLASS IDENTICATION
 447 //   BackprojectIntDiff2
 448 //
 449 // PURPOSE
 450 //   Integer version of BackprojectDiff2
 451
 452 void
 453 BackprojectIntDiff2::BackprojectView (const double* const filteredProj, const double view_angle)
 454 {
 455   double theta = - view_angle;  // add half PI to view angle to get perpendicular theta angle
 456
 457   kint32 scale = 1 << 16;
 458   double dScale = scale;
 459   kint32 halfScale = scale / 2;
 460
 461   kint32 det_dx = nearest<kint32> (xInc * sin (theta) / detInc * scale);
 462   kint32 det_dy = nearest<kint32> (yInc * cos (theta) / detInc * scale);
 463
 464   // calculate L for first point in image (0, 0)
 465   kint32 detPosColStart = nearest<kint32> (start_r * cos (theta - start_phi) / detInc * scale);
 466
 467   for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
 468     kint32 curDetPos = detPosColStart;
 469     ImageFileColumn pImCol = v[ix];
 470
 471     for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
 472       if (interpType == Backprojector::INTERP_NEAREST) {
 473         int detPosNearest = (curDetPos >= 0 ? ((curDetPos + halfScale) / scale) : ((curDetPos - halfScale) / scale));
 474         int iDetPos = iDetCenter + detPosNearest;       // calc index in the filtered raysum vector
 475
 476         if (iDetPos < 0 || iDetPos >= nDet)  // check for impossible: index outside of raysum pos
 477             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 478         else
 479           *pImCol++ += filteredProj[iDetPos];
 480       } else if (interpType == Backprojector::INTERP_LINEAR) {
 481         kint32 detPosFloor = curDetPos / scale;
 482         kint32 detPosRemainder = curDetPos % scale;
 483         if (detPosRemainder < 0) {
 484           detPosFloor--;
 485           detPosRemainder += scale;
 486         }
 487         int iDetPos = iDetCenter + detPosFloor;
 488         double frac = detPosRemainder / dScale;
 489         if (iDetPos < 0 || iDetPos >= nDet - 1)
 490             errorIndexOutsideDetector (ix, iy, theta, curDetPos, iDetPos);
 491         else
 492           *pImCol++ += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
 493       }
 494     }   // end for y
 495   }     // end for x
 496 }