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.10 2000/08/03 09:21:12 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)
  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_detLen   = SQRT2 * m_phmLen * ((m_nDet + N_EXTRA_DETECTORS) / static_cast<double>(m_nDet));
 111
 112   m_rotLen   =  rot_anglen;
 113
 114   m_radius  = m_detLen / 2;
 115   m_detInc  = m_detLen / m_nDet;
 116   m_rotInc  = m_rotLen / m_nView;
 117
 118   m_initPos.xd1 = m_detLen / 2;
 119   m_initPos.yd1 = -m_detLen / 2;
 120   m_initPos.xd2 = m_detLen / 2;
 121   m_initPos.yd2 = m_detLen / 2;
 122   m_initPos.xs1 = -m_detLen / 2;
 123   m_initPos.ys1 = -m_detLen / 2;
 124   m_initPos.xs2 = -m_detLen / 2;
 125   m_initPos.ys2 = m_detLen / 2;
 126   m_initPos.angle = 0.0;
 127 }
 128
 129 Scanner::~Scanner (void)
 130 {
 131 }
 132
 133
 134 const char*
 135 Scanner::convertGeometryIDToName (const int geomID)
 136 {
 137   const char *name = "";
 138
 139   if (geomID >= 0 && geomID < s_iGeometryCount)
 140       return (s_aszGeometryName[geomID]);
 141
 142   return (name);
 143 }
 144
 145 const char*
 146 Scanner::convertGeometryIDToTitle (const int geomID)
 147 {
 148   const char *title = "";
 149
 150   if (geomID >= 0 && geomID < s_iGeometryCount)
 151       return (s_aszGeometryName[geomID]);
 152
 153   return (title);
 154 }
 155
 156 int
 157 Scanner::convertGeometryNameToID (const char* const geomName)
 158 {
 159   int id = GEOMETRY_INVALID;
 160
 161   for (int i = 0; i < s_iGeometryCount; i++)
 162       if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
 163           id = i;
 164           break;
 165       }
 166
 167   return (id);
 168 }
 169
 170
 171 /* NAME
 172  *   raysum_collect                             Calculate ray sums for a Phantom
 173  *
 174  * SYNOPSIS
 175  *   rs = raysum_collect (det, phm, start_view, trace, unit_pulse)
 176  *   Scanner& det                       Scanner specifications**
 177  *   RAYSUM *rs                         Calculated ray sum matrix
 178  *   Phantom& phm                       Phantom we are taking ray sums of
 179  *   int trace                          Boolean flag to signal ray sum tracing
 180 */
 181
 182 void
 183 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int start_view = 0, const int trace = TRACE_NONE, SGP* pSGP = NULL)
 184 {
 185   GRFMTX_2D rotmtx_initial, temp;
 186   GRFMTX_2D rotmtx_incr;
 187
 188   double start_angle = start_view * proj.rotInc();
 189   double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
 190   double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
 191
 192   double xd1 = xcent + m_initPos.xd1;
 193   double yd1 = ycent + m_initPos.yd1;
 194   double xd2 = xcent + m_initPos.xd2;
 195   double yd2 = ycent + m_initPos.yd2;
 196   double xs1 = xcent + m_initPos.xs1;
 197   double ys1 = ycent + m_initPos.ys1;
 198   double xs2 = xcent + m_initPos.xs2;
 199   double ys2 = ycent + m_initPos.ys2;
 200
 201   m_trace = trace;
 202
 203 #ifdef HAVE_SGP
 204   if (pSGP && m_trace >= TRACE_PHM) {
 205     double halfPhmLen = m_phmLen / 2;
 206     double wsize = SQRT2 * halfPhmLen;
 207
 208     pSGP->setRasterOp (RO_SET);
 209     pSGP->eraseWindow ();
 210     pSGP->setColor (C_LTBLUE);
 211     pSGP->setWindow (xcent - wsize, ycent - wsize, xcent + wsize, ycent + wsize);
 212     pSGP->setColor (C_BROWN);
 213     pSGP->drawRect (xcent - halfPhmLen, ycent - halfPhmLen, xcent + halfPhmLen, ycent + halfPhmLen);
 214     pSGP->setColor (C_BROWN);
 215     pSGP->moveAbs (0., 0.);
 216     pSGP->drawCircle (wsize);
 217
 218     traceShowParam ("X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " ");
 219     traceShowParam ("---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " ");
 220     traceShowParam ("Phantom:",       "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
 221     traceShowParam ("Chomaticity  :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
 222     traceShowParam ("Scatter      :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
 223     traceShowParam ("Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
 224     traceShowParam ("Num Detectors:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
 225     traceShowParam ("Num Views    :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nView());
 226     traceShowParam ("Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
 227
 228     pSGP->setColor (C_LTGREEN);
 229     phm.draw (*pSGP);
 230
 231     pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
 232   }
 233 #endif
 234
 235 /* Calculate initial rotation matrix */
 236   xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
 237   rot_mtx2 (temp, start_angle);
 238   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 239   xlat_mtx2 (temp, xcent, ycent);
 240   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 241
 242   xform_mtx2 (rotmtx_initial, xd1, yd1);        /* rotate detector endpoints */
 243   xform_mtx2 (rotmtx_initial, xd2, yd2);      /* to initial view_angle */
 244   xform_mtx2 (rotmtx_initial, xs1, ys1);
 245   xform_mtx2 (rotmtx_initial, xs2, ys2);
 246
 247 /* Calculate incrementatal rotation matrix */
 248   xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
 249   rot_mtx2 (temp, proj.rotInc());
 250   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 251   xlat_mtx2 (temp, xcent, ycent);
 252   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 253
 254   int iview;
 255   double viewAngle;
 256   for (iview = 0, viewAngle = start_angle;  iview < proj.nView(); iview++, viewAngle += proj.rotInc()) {
 257       DetectorArray& detArray = proj.getDetectorArray( iview );
 258
 259 #ifdef HAVE_SGP
 260     if (pSGP && m_trace >= TRACE_PHM) {
 261       pSGP->setRasterOp (RO_XOR);
 262       pSGP->setColor (C_RED);
 263       pSGP->moveAbs (xd1, yd1);
 264       pSGP->lineAbs (xd2, yd2);
 265       pSGP->moveAbs (xs1, ys1);
 266       pSGP->lineAbs (xs2, ys2);
 267     }
 268     if (m_trace)
 269       traceShowParam ("Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
 270 #endif
 271
 272     projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, pSGP);
 273     detArray.setViewAngle (viewAngle);
 274
 275 #ifdef HAVE_SGP
 276     if (pSGP && m_trace >= TRACE_PHM) {
 277       //        rs_plot (detArray, xd1, yd1, xcent, ycent, theta);
 278       pSGP->setColor (C_RED);
 279       pSGP->moveAbs (xd1, yd1);
 280       pSGP->lineAbs (xd2, yd2);
 281       pSGP->moveAbs (xs1, ys1);
 282       pSGP->lineAbs (xs2, ys2);
 283     }
 284 #endif
 285     xform_mtx2 (rotmtx_incr, xd1, yd1);  // rotate detector endpoints
 286     xform_mtx2 (rotmtx_incr, xd2, yd2);
 287     xform_mtx2 (rotmtx_incr, xs1, ys1);
 288     xform_mtx2 (rotmtx_incr, xs2, ys2);
 289   } /* for each iview */
 290 }
 291
 292
 293 /* NAME
 294  *    rayview                   Calculate raysums for a view at any angle
 295  *
 296  * SYNOPSIS
 297  *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 298  *    Phantom& phm              Phantom to scan
 299  *    DETARRAY *detArray                Storage of values for detector array
 300  *    Scanner& det              Scanner parameters
 301  *    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
 302  *    double xs1, ys1, xs2, ys2 Beginning & ending source positions
 303  *
 304  * RAY POSITIONING
 305  *         For each detector, have there are a variable number of rays traced.
 306  *     The source of each ray is the center of the source x-ray cell. The
 307  *     detector positions are equally spaced within the cell
 308  *
 309  *         The increments between rays are calculated so that the cells start
 310  *     at the beginning of a detector cell and they end on the endpoint
 311  *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 312  *         The exception to this is if there is only one ray per detector.
 313  *     In that case, the detector position is the center of the detector cell.
 314  */
 315
 316 void
 317 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)
 318 {
 319   double ddx = (xd2 - xd1) / detArray.nDet();  // change in coords between detectors
 320   double ddy = (yd2 - yd1) / detArray.nDet();
 321   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords between source
 322   double sdy = (ys2 - ys1) / detArray.nDet();
 323
 324   double ddx2 = ddx / m_nSample;        // Incr. between rays with detector cell
 325   double ddy2 = ddy / m_nSample;        // Doesn't include detector endpoints
 326   double ddx2_ofs = ddx2 / 2;           // offset of 1st ray from start of detector cell
 327   double ddy2_ofs = ddy2 / 2;
 328
 329   double xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 330   double yd_maj = yd1 + ddy2_ofs;
 331   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 332   double ys_maj = ys1 + (sdy / 2);
 333
 334   DetectorValue* detval = detArray.detValues();
 335
 336   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 337     for (int d = 0; d < detArray.nDet(); d++)
 338       if (detArray.nDet() / 2 == d && (d % 2) == 1)
 339         detval[d] = 1;
 340       else
 341         detval[d] = 0;
 342   } else {
 343     for (int d = 0; d < detArray.nDet(); d++) {
 344       double xd = xd_maj;
 345       double yd = yd_maj;
 346       double xs = xs_maj;
 347       double ys = ys_maj;
 348       double sum = 0.0;
 349       for (unsigned int i = 0; i < m_nSample; i++) {
 350 #ifdef HAVE_SGP
 351         if (pSGP && m_trace >= TRACE_RAYS) {
 352           pSGP->setColor (C_LTBLUE);
 353           pSGP->moveAbs (xs, ys);
 354           pSGP->lineAbs (xd, yd);
 355         }
 356 #endif
 357         sum += projectSingleLine (phm, xd, yd, xs, ys, pSGP);
 358
 359 #ifdef HAVE_SGP
 360         if (m_trace >= TRACE_RAYS)
 361           traceShowParam ("Attenuation  :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, "sum");
 362
 363         if (pSGP && m_trace >= TRACE_RAYS) {
 364           pSGP->setColor (C_LTBLUE);
 365           pSGP->moveAbs (xs, ys);
 366           pSGP->lineAbs (xd, yd);
 367         }
 368 #endif
 369         xd += ddx2;
 370         yd += ddy2;
 371       }
 372
 373       detval[d] = sum / m_nSample;
 374       xd_maj += ddx;
 375       yd_maj += ddy;
 376       xs_maj += sdx;
 377       ys_maj += sdy;
 378     } /* for each detector */
 379   } /* if not unit pulse */
 380 }
 381
 382
 383 void
 384 Scanner::traceShowParam (const char *label, const char *fmt, int row, int color, ...)
 385 {
 386   char s[80];
 387   va_list arg;
 388
 389   va_start(arg, color);
 390   //  cio_set_cpos (raysum_trace_menu_column, row);
 391   snprintf (s, sizeof(s), label, "%s");
 392   //  cio_set_text_clr (color - 8, 0);
 393   cio_put_str (s);
 394   vsnprintf (s, sizeof(s), fmt, arg);
 395   //  cio_set_text_clr (color, 0);
 396   cio_put_str (s);
 397   cout << "\n";
 398   va_end(arg);
 399 }
 400
 401
 402 /* NAME
 403  *    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
 404  *
 405  * SYNOPSIS
 406  *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 407  *    double rsum               Ray sum of Phantom along given line
 408  *    Phantom& phm;             Phantom from which to calculate raysum
 409  *    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
 410  */
 411
 412 double
 413 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2, SGP* pSGP)
 414 {
 415   // check ray against each pelem in Phantom
 416   double rsum = 0.0;
 417   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 418     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2, pSGP);
 419
 420   return (rsum);
 421 }
 422
 423
 424 /* NAME
 425  *   pelem_ray_attenuation              Calculate raysum of an pelem along one line
 426  *
 427  * SYNOPSIS
 428  *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 429  *   double rsum                Computed raysum
 430  *   PhantomElement& pelem              Pelem to scan
 431  *   double x1, y1, x2, y2      Endpoints of raysum line
 432  */
 433
 434 double
 435 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2, SGP* pSGP)
 436 {
 437   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 438     if (m_trace == TRACE_CLIPPING)
 439       cio_tone (1000., 0.05);
 440     return (0.0);
 441   }
 442
 443 #ifdef HAVE_SGP
 444   if (pSGP && m_trace == TRACE_CLIPPING) {
 445     pSGP->moveAbs (x1, y1);
 446     pSGP->lineAbs (x2, y2);
 447     cio_tone (8000., 0.05);
 448     pSGP->moveAbs (x1, y1);
 449     pSGP->lineAbs (x2, y2);
 450   }
 451 #endif
 452
 453   double len = lineLength (x1, y1, x2, y2);
 454   return (len * pelem.atten());
 455 }
 456