libctsim/scanner.cpp

   1 /*****************************************************************************
   2 ** FILE IDENTIFICATION
   3 **
   4 **   Name:          scanner.cpp
   5 **   Purpose:       Classes for CT scanner
   6 **   Programmer:    Kevin Rosenberg
   7 **   Date Started:  1984
   8 **
   9 **  This is part of the CTSim program
  10 **  Copyright (C) 1983-2000 Kevin Rosenberg
  11 **
  12 **  $Id: scanner.cpp,v 1.11 2000/08/25 15:59:13 kevin Exp $
  13 **
  14 **  This program is free software; you can redistribute it and/or modify
  15 **  it under the terms of the GNU General Public License (version 2) as
  16 **  published by the Free Software Foundation.
  17 **
  18 **  This program is distributed in the hope that it will be useful,
  19 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 **  GNU General Public License for more details.
  22 **
  23 **  You should have received a copy of the GNU General Public License
  24 **  along with this program; if not, write to the Free Software
  25 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  26 ******************************************************************************/
  27
  28 #include "ct.h"
  29
  30
  31 const int Scanner::GEOMETRY_INVALID = -1;
  32 const int Scanner::GEOMETRY_PARALLEL = 0;
  33 const int Scanner::GEOMETRY_EQUILINEAR = 1;
  34 const int Scanner::GEOMETRY_EQUIANGULAR = 2;
  35
  36 const char* Scanner::s_aszGeometryName[] =
  37 {
  38   {"parallel"},
  39   {"equilinear"},
  40   {"equiangular"},
  41 };
  42
  43 const char* Scanner::s_aszGeometryTitle[] =
  44 {
  45   {"Parallel"},
  46   {"Equilinear"},
  47   {"Equiangular"},
  48 };
  49
  50 const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
  51
  52
  53 // NAME
  54 //   DetectorArray       Construct a DetectorArray
  55
  56 DetectorArray::DetectorArray (const int nDet)
  57 {
  58   m_nDet = nDet;
  59   m_detValues = new DetectorValue [m_nDet];
  60 }
  61
  62
  63 // NAME
  64 //   ~DetectorArray             Free memory allocated to a detector array
  65
  66 DetectorArray::~DetectorArray (void)
  67 {
  68   delete [] m_detValues;
  69 }
  70
  71
  72
  73 /* NAME
  74  *   Scanner::Scanner           Construct a user specified detector structure
  75  *
  76  * SYNOPSIS
  77  *   Scanner (phm, nDet, nView, nSample)
  78  *   Phantom& phm               PHANTOM that we are making detector for
  79  *   int geomety                Geometry of detector
  80  *   int nDet                   Number of detector along detector array
  81  *   int nView                  Number of rotated views
  82  *   int nSample                Number of rays per detector
  83  */
  84
  85 Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen, const double dFocalLengthRatio, const double dFieldOfViewRatio)
  86 {
  87   m_phmLen = phm.maxAxisLength();      // maximal length along an axis
  88
  89   m_fail = false;
  90   m_idGeometry = convertGeometryNameToID (geometryName);
  91   if (m_idGeometry == GEOMETRY_INVALID) {
  92     m_fail = true;
  93     m_failMessage = "Invalid geometry name ";
  94     m_failMessage += geometryName;
  95     return;
  96   }
  97
  98   if (nView < 1)
  99     nView = 1;
 100   if (nSample < 1)
 101     m_nSample = 1;
 102   if (nDet < 1)
 103     nDet = 1;
 104   //  if ((nDet % 2) == 0)
 105   //    ++nDet;         // ensure odd number of detectors
 106
 107   m_nDet     = nDet;
 108   m_nView    = nView;
 109   m_nSample  = nSample;
 110   m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio;
 111   m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio;
 112
 113   if (m_idGeometry == GEOMETRY_PARALLEL) {
 114     m_detLen   = m_dFieldOfView;
 115     m_detInc  = m_detLen / m_nDet;
 116     m_rotLen  = rot_anglen;
 117     m_rotInc  = m_rotLen / m_nView;
 118
 119     double dHalfDetLen = m_detLen / 2;
 120     m_initPos.xs1 = -m_dFocalLength;
 121     m_initPos.ys1 = -dHalfDetLen;
 122     m_initPos.xs2 = -m_dFocalLength;
 123     m_initPos.ys2 = dHalfDetLen;
 124     m_initPos.xd1 = m_dFocalLength;
 125     m_initPos.yd1 = -dHalfDetLen;
 126     m_initPos.xd2 = m_dFocalLength;
 127     m_initPos.yd2 = dHalfDetLen;
 128     m_initPos.angle = 0.0;
 129   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
 130     double dHalfSquare = m_phmLen / 2;
 131     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
 132     if (dFocalPastPhm <= 0.) {
 133       m_fail = true;
 134       m_failMessage = "Focal Point inside of phantom";
 135       return;
 136     }
 137     double dAngle = atan( dHalfSquare / dFocalPastPhm );
 138     double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
 139
 140     m_detLen = dHalfDetLen * 2;
 141     m_detInc  = m_detLen / m_nDet;
 142     m_rotLen  = rot_anglen;
 143     m_rotInc  = m_rotLen / m_nView;
 144
 145     m_initPos.xs1 = -m_dFocalLength;
 146     m_initPos.ys1 = 0;
 147     m_initPos.xs2 = -m_dFocalLength;
 148     m_initPos.ys2 = 0;
 149     m_initPos.xd1 = m_dFocalLength;
 150     m_initPos.yd1 = -dHalfDetLen;
 151     m_initPos.xd2 = m_dFocalLength;
 152     m_initPos.yd2 = dHalfDetLen;
 153     m_initPos.angle = 0.0;
 154   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 155     double dHalfSquare = m_phmLen / 2;
 156     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
 157     if (dFocalPastPhm <= 0.) {
 158       m_fail = true;
 159       m_failMessage = "Focal Point inside of phantom";
 160       return;
 161     }
 162     double dAngle = atan( dHalfSquare / dFocalPastPhm );
 163     m_detLen = dAngle * 2;
 164   }
 165 }
 166
 167 Scanner::~Scanner (void)
 168 {
 169 }
 170
 171
 172 const char*
 173 Scanner::convertGeometryIDToName (const int geomID)
 174 {
 175   const char *name = "";
 176
 177   if (geomID >= 0 && geomID < s_iGeometryCount)
 178       return (s_aszGeometryName[geomID]);
 179
 180   return (name);
 181 }
 182
 183 const char*
 184 Scanner::convertGeometryIDToTitle (const int geomID)
 185 {
 186   const char *title = "";
 187
 188   if (geomID >= 0 && geomID < s_iGeometryCount)
 189       return (s_aszGeometryName[geomID]);
 190
 191   return (title);
 192 }
 193
 194 int
 195 Scanner::convertGeometryNameToID (const char* const geomName)
 196 {
 197   int id = GEOMETRY_INVALID;
 198
 199   for (int i = 0; i < s_iGeometryCount; i++)
 200       if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
 201           id = i;
 202           break;
 203       }
 204
 205   return (id);
 206 }
 207
 208
 209 /* NAME
 210  *   collectProjections         Calculate projections for a Phantom
 211  *
 212  * SYNOPSIS
 213  *   collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace)
 214  *   Projectrions& proj      Projection storage
 215  *   Phantom& phm            Phantom for which we collect projections
 216  *   bool bStoreViewPos      TRUE then storage proj at normal view position
 217  *   int trace               Trace level
 218 */
 219
 220
 221 void
 222 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP)
 223 {
 224   collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP);
 225 }
 226
 227 void
 228 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, bool bStoreAtViewPosition, const int trace, SGP* pSGP)
 229 {
 230   GRFMTX_2D rotmtx_initial, temp;
 231   GRFMTX_2D rotmtx_incr;
 232
 233   double start_angle = iStartView * proj.rotInc();
 234   double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
 235   double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
 236
 237   double xd1 = xcent + m_initPos.xd1;
 238   double yd1 = ycent + m_initPos.yd1;
 239   double xd2 = xcent + m_initPos.xd2;
 240   double yd2 = ycent + m_initPos.yd2;
 241   double xs1 = xcent + m_initPos.xs1;
 242   double ys1 = ycent + m_initPos.ys1;
 243   double xs2 = xcent + m_initPos.xs2;
 244   double ys2 = ycent + m_initPos.ys2;
 245
 246   m_trace = trace;
 247
 248 /* Calculate initial rotation matrix */
 249   xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
 250   rot_mtx2 (temp, start_angle);
 251   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 252   xlat_mtx2 (temp, xcent, ycent);
 253   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 254
 255   xform_mtx2 (rotmtx_initial, xd1, yd1);        /* rotate detector endpoints */
 256   xform_mtx2 (rotmtx_initial, xd2, yd2);      /* to initial view_angle */
 257   xform_mtx2 (rotmtx_initial, xs1, ys1);
 258   xform_mtx2 (rotmtx_initial, xs2, ys2);
 259
 260 /* Calculate incrementatal rotation matrix */
 261   xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
 262   rot_mtx2 (temp, proj.rotInc());
 263   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 264   xlat_mtx2 (temp, xcent, ycent);
 265   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 266
 267   int iview;
 268   double viewAngle;
 269   for (iview = 0, viewAngle = start_angle;  iview < iNumViews; iview++, viewAngle += proj.rotInc()) {
 270     int iStoragePosition = iview;
 271     if (bStoreAtViewPosition)
 272       iStoragePosition += iStartView;
 273
 274     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
 275
 276 #ifdef HAVE_SGP
 277   if (pSGP && m_trace >= TRACE_PHM) {
 278     pSGP->eraseWindow();
 279     double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2;
 280     double dHalfWindowSize = dWindowSize / 2;
 281     double dHalfPhmLen = m_phmLen /  2;
 282
 283     pSGP->setRasterOp (RO_SET);
 284     pSGP->eraseWindow ();
 285     pSGP->setWindow (xcent - dHalfWindowSize, ycent - dHalfWindowSize, xcent + dHalfWindowSize, ycent + dHalfWindowSize);
 286     pSGP->setColor (C_BROWN);
 287     pSGP->moveAbs (0., 0.);
 288     pSGP->drawRect (xcent - dHalfPhmLen, ycent - dHalfPhmLen, xcent + dHalfPhmLen, ycent + dHalfPhmLen);
 289
 290 #if 0
 291     traceShowParam (pSGP, "X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " ");
 292     traceShowParam (pSGP, "---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " ");
 293     traceShowParam (pSGP, "Phantom:",       "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
 294     traceShowParam (pSGP, "Chomaticity  :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
 295     traceShowParam (pSGP, "Scatter      :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
 296     traceShowParam (pSGP, "Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
 297     traceShowParam (pSGP, "Num Detectors:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
 298     traceShowParam (pSGP, "Num Views    :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nViews());
 299     traceShowParam (pSGP, "Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
 300 #endif
 301
 302     pSGP->setColor (C_RED);
 303     phm.draw (*pSGP);
 304
 305     pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
 306   }
 307 #endif
 308
 309 #ifdef HAVE_SGP
 310     if (pSGP && m_trace >= TRACE_PHM) {
 311       pSGP->setRasterOp (RO_XOR);
 312       pSGP->setColor (C_RED);
 313       pSGP->moveAbs (xd1, yd1);
 314       pSGP->lineAbs (xd2, yd2);
 315       pSGP->moveAbs (xs1, ys1);
 316       pSGP->lineAbs (xs2, ys2);
 317       pSGP->setRasterOp (RO_SET);
 318     }
 319     if (m_trace >= TRACE_TEXT)
 320       traceShowParam (pSGP, "Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
 321 #endif
 322
 323     projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, pSGP);
 324     detArray.setViewAngle (viewAngle);
 325
 326 #ifdef HAVE_SGP
 327     if (pSGP && m_trace >= TRACE_PHM) {
 328       //        rs_plot (detArray, xd1, yd1, xcent, ycent, theta);
 329       pSGP->setColor (C_RED);
 330       pSGP->setRasterOp (RO_XOR);
 331       pSGP->moveAbs (xd1, yd1);
 332       pSGP->lineAbs (xd2, yd2);
 333       pSGP->moveAbs (xs1, ys1);
 334       pSGP->lineAbs (xs2, ys2);
 335       pSGP->setRasterOp (RO_SET);
 336     }
 337 #endif
 338     xform_mtx2 (rotmtx_incr, xd1, yd1);  // rotate detector endpoints
 339     xform_mtx2 (rotmtx_incr, xd2, yd2);
 340     xform_mtx2 (rotmtx_incr, xs1, ys1);
 341     xform_mtx2 (rotmtx_incr, xs2, ys2);
 342   } /* for each iview */
 343 }
 344
 345
 346 /* NAME
 347  *    rayview                   Calculate raysums for a view at any angle
 348  *
 349  * SYNOPSIS
 350  *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 351  *    Phantom& phm              Phantom to scan
 352  *    DETARRAY *detArray                Storage of values for detector array
 353  *    Scanner& det              Scanner parameters
 354  *    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
 355  *    double xs1, ys1, xs2, ys2 Beginning & ending source positions
 356  *
 357  * RAY POSITIONING
 358  *         For each detector, have there are a variable number of rays traced.
 359  *     The source of each ray is the center of the source x-ray cell. The
 360  *     detector positions are equally spaced within the cell
 361  *
 362  *         The increments between rays are calculated so that the cells start
 363  *     at the beginning of a detector cell and they end on the endpoint
 364  *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 365  *         The exception to this is if there is only one ray per detector.
 366  *     In that case, the detector position is the center of the detector cell.
 367  */
 368
 369 void
 370 Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const double xd1, const double yd1, const double xd2, const double yd2, const double xs1, const double ys1, const double xs2, const double ys2, SGP* pSGP)
 371 {
 372   double ddx = (xd2 - xd1) / detArray.nDet();  // change in coords between detectors
 373   double ddy = (yd2 - yd1) / detArray.nDet();
 374   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords between source
 375   double sdy = (ys2 - ys1) / detArray.nDet();
 376
 377   double ddx2 = ddx / m_nSample;        // Incr. between rays with detector cell
 378   double ddy2 = ddy / m_nSample;        // Doesn't include detector endpoints
 379   double ddx2_ofs = ddx2 / 2;           // offset of 1st ray from start of detector cell
 380   double ddy2_ofs = ddy2 / 2;
 381
 382   double xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 383   double yd_maj = yd1 + ddy2_ofs;
 384   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 385   double ys_maj = ys1 + (sdy / 2);
 386
 387   DetectorValue* detval = detArray.detValues();
 388
 389   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 390     for (int d = 0; d < detArray.nDet(); d++)
 391       if (detArray.nDet() / 2 == d && (d % 2) == 1)
 392         detval[d] = 1;
 393       else
 394         detval[d] = 0;
 395   } else {
 396     for (int d = 0; d < detArray.nDet(); d++) {
 397       double xd = xd_maj;
 398       double yd = yd_maj;
 399       double xs = xs_maj;
 400       double ys = ys_maj;
 401       double sum = 0.0;
 402       for (unsigned int i = 0; i < m_nSample; i++) {
 403 #ifdef HAVE_SGP
 404         if (pSGP && m_trace >= TRACE_RAYS) {
 405           pSGP->setColor (C_LTBLUE);
 406           pSGP->setRasterOp (RO_XOR);
 407           pSGP->moveAbs (xs, ys);
 408           pSGP->lineAbs (xd, yd);
 409           pSGP->setRasterOp (RO_SET);
 410         }
 411 #endif
 412         sum += projectSingleLine (phm, xd, yd, xs, ys, pSGP);
 413
 414 #ifdef HAVE_SGP
 415         if (m_trace >= TRACE_RAYS) {
 416           traceShowParam (pSGP, "Attenuation  :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
 417         }
 418
 419         if (pSGP && m_trace >= TRACE_RAYS) {
 420           pSGP->setColor (C_LTBLUE);
 421           pSGP->setRasterOp (RO_XOR);
 422           pSGP->moveAbs (xs, ys);
 423           pSGP->lineAbs (xd, yd);
 424           pSGP->setRasterOp (RO_SET);
 425         }
 426 #endif
 427         xd += ddx2;
 428         yd += ddy2;
 429       }
 430
 431       detval[d] = sum / m_nSample;
 432       xd_maj += ddx;
 433       yd_maj += ddy;
 434       xs_maj += sdx;
 435       ys_maj += sdy;
 436     } /* for each detector */
 437   } /* if not unit pulse */
 438 }
 439
 440
 441 void
 442 Scanner::traceShowParam (SGP* pSGP, const char *label, const char *fmt, int row, int color, ...)
 443 {
 444   char s[256];
 445   va_list arg;
 446
 447   va_start(arg, color);
 448   snprintf (s, sizeof(s), label, "%s");
 449   string strOut(s);
 450   vsnprintf (s, sizeof(s), fmt, arg);
 451   strOut += s;
 452
 453   //  cio_set_cpos (raysum_trace_menu_column, row);
 454   //  cio_set_text_clr (color - 8, 0);
 455   //  cio_set_text_clr (color, 0);
 456
 457   if (pSGP) {
 458       pSGP->moveAbs (0., row * 0.04);
 459       pSGP->setTextColor (color, -1);
 460       pSGP->drawText (strOut);
 461   } else {
 462     cio_put_str (strOut.c_str());
 463     cio_put_str ("\n");
 464   }
 465
 466   va_end(arg);
 467 }
 468
 469
 470 /* NAME
 471  *    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
 472  *
 473  * SYNOPSIS
 474  *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 475  *    double rsum               Ray sum of Phantom along given line
 476  *    Phantom& phm;             Phantom from which to calculate raysum
 477  *    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
 478  */
 479
 480 double
 481 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2, SGP* pSGP)
 482 {
 483   // check ray against each pelem in Phantom
 484   double rsum = 0.0;
 485   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 486     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2, pSGP);
 487
 488   return (rsum);
 489 }
 490
 491
 492 /* NAME
 493  *   pelem_ray_attenuation              Calculate raysum of an pelem along one line
 494  *
 495  * SYNOPSIS
 496  *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 497  *   double rsum                Computed raysum
 498  *   PhantomElement& pelem              Pelem to scan
 499  *   double x1, y1, x2, y2      Endpoints of raysum line
 500  */
 501
 502 double
 503 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2, SGP* pSGP)
 504 {
 505   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 506     if (m_trace == TRACE_CLIPPING)
 507       cio_tone (1000., 0.05);
 508     return (0.0);
 509   }
 510
 511 #ifdef HAVE_SGP
 512   if (pSGP && m_trace == TRACE_CLIPPING) {
 513     pSGP->setRasterOp (RO_XOR);
 514     pSGP->moveAbs (x1, y1);
 515     pSGP->lineAbs (x2, y2);
 516     cio_tone (8000., 0.05);
 517     pSGP->moveAbs (x1, y1);
 518     pSGP->lineAbs (x2, y2);
 519     pSGP->setRasterOp (RO_SET);
 520   }
 521 #endif
 522
 523   double len = lineLength (x1, y1, x2, y2);
 524   return (len * pelem.atten());
 525 }
 526