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