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.7 2000/07/29 19:50:08 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, 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 wsize = 1.42 * m_phmLen / 2; /* sqrt(2) * radius */
 206
 207     pSGP->setColor (C_LTBLUE);
 208     pSGP->setWindow (xcent - wsize, ycent - wsize, xcent + wsize, ycent + wsize);
 209     pSGP->setColor (C_BROWN);
 210 #if RADIUS
 211     pSGP->drawCircle (m_phmLen / 2);
 212 #else
 213     pSGP->drawRect (xcent - m_phmLen / 2, ycent - m_phmLen / 2,
 214                     xcent + m_phmLen / 2, ycent + m_phmLen / 2);
 215 #endif
 216     pSGP->setColor (C_BROWN);
 217     pSGP->moveAbs (0., 0.);
 218     pSGP->drawCircle (wsize);
 219     //      raysum_trace_menu_column = (crt->xsize * crt->asp) / 8 + 3;
 220     traceShowParam ("X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " ");
 221     traceShowParam ("---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " ");
 222     traceShowParam ("Phantom:",       "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
 223     traceShowParam ("Chomaticity  :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
 224     traceShowParam ("Scatter      :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
 225     traceShowParam ("Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
 226     traceShowParam ("Num Scanners:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
 227     traceShowParam ("Num Views    :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nView());
 228     traceShowParam ("Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
 229
 230     pSGP->setColor (C_LTGREEN);
 231     phm.draw (*pSGP);
 232
 233     pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
 234   }
 235 #endif
 236
 237 /* Calculate initial rotation matrix */
 238   xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
 239   rot_mtx2 (temp, start_angle);
 240   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 241   xlat_mtx2 (temp, xcent, ycent);
 242   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 243
 244   xform_mtx2 (rotmtx_initial, xd1, yd1);        /* rotate detector endpoints */
 245   xform_mtx2 (rotmtx_initial, xd2, yd2);      /* to initial view_angle */
 246   xform_mtx2 (rotmtx_initial, xs1, ys1);
 247   xform_mtx2 (rotmtx_initial, xs2, ys2);
 248
 249 /* Calculate incrementatal rotation matrix */
 250   xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
 251   rot_mtx2 (temp, proj.rotInc());
 252   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 253   xlat_mtx2 (temp, xcent, ycent);
 254   mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
 255
 256   int iview;
 257   double viewAngle;
 258   for (iview = 0, viewAngle = start_angle;  iview < proj.nView(); iview++, viewAngle += proj.rotInc()) {
 259       DetectorArray& detArray = proj.getDetectorArray( iview );
 260
 261 #ifdef HAVE_SGP
 262     if (pSGP && m_trace >= TRACE_PHM) {
 263       pSGP->moveAbs (xd1, yd1);
 264       pSGP->lineAbs (xd2, yd2);
 265       pSGP->moveAbs (xs1, ys1);
 266       pSGP->lineAbs (xs2, ys2);
 267     }
 268 #endif
 269     if (m_trace)
 270       traceShowParam ("Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
 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->moveAbs (xd1, yd1);
 279       pSGP->lineAbs (xd2, yd2);
 280       pSGP->moveAbs (xs1, ys1);
 281       pSGP->lineAbs (xs2, ys2);
 282     }
 283 #endif
 284     xform_mtx2 (rotmtx_incr, xd1, yd1);  // rotate detector endpoints
 285     xform_mtx2 (rotmtx_incr, xd2, yd2);
 286     xform_mtx2 (rotmtx_incr, xs1, ys1);
 287     xform_mtx2 (rotmtx_incr, xs2, ys2);
 288   } /* for each iview */
 289 }
 290
 291
 292 /* NAME
 293  *    rayview                   Calculate raysums for a view at any angle
 294  *
 295  * SYNOPSIS
 296  *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 297  *    Phantom& phm              Phantom to scan
 298  *    DETARRAY *detArray                Storage of values for detector array
 299  *    Scanner& det              Scanner parameters
 300  *    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
 301  *    double xs1, ys1, xs2, ys2 Beginning & ending source positions
 302  *
 303  * RAY POSITIONING
 304  *         For each detector, have there are a variable number of rays traced.
 305  *     The source of each ray is the center of the source x-ray cell. The
 306  *     detector positions are equally spaced within the cell
 307  *
 308  *         The increments between rays are calculated so that the cells start
 309  *     at the beginning of a detector cell and they end on the endpoint
 310  *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 311  *         The exception to this is if there is only one ray per detector.
 312  *     In that case, the detector position is the center of the detector cell.
 313  */
 314
 315 void
 316 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)
 317 {
 318   double ddx = (xd2 - xd1) / detArray.nDet();  // change in coords between detectors
 319   double ddy = (yd2 - yd1) / detArray.nDet();
 320   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords between source
 321   double sdy = (ys2 - ys1) / detArray.nDet();
 322
 323   double ddx2 = ddx / m_nSample;        // Incr. between rays with detector cell
 324   double ddy2 = ddy / m_nSample;        // Doesn't include detector endpoints
 325   double ddx2_ofs = ddx2 / 2;           // offset of 1st ray from start of detector cell
 326   double ddy2_ofs = ddy2 / 2;
 327
 328   double xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 329   double yd_maj = yd1 + ddy2_ofs;
 330   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 331   double ys_maj = ys1 + (sdy / 2);
 332
 333   DetectorValue* detval = detArray.detValues();
 334
 335   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 336     for (int d = 0; d < detArray.nDet(); d++)
 337       if (detArray.nDet() / 2 == d && (d % 2) == 1)
 338         detval[d] = 1;
 339       else
 340         detval[d] = 0;
 341   } else {
 342     for (int d = 0; d < detArray.nDet(); d++) {
 343       double xd = xd_maj;
 344       double yd = yd_maj;
 345       double xs = xs_maj;
 346       double ys = ys_maj;
 347       double sum = 0.0;
 348       for (unsigned int i = 0; i < m_nSample; i++) {
 349 #ifdef HAVE_SGP
 350         if (pSGP && m_trace >= TRACE_RAYS) {
 351           pSGP->moveAbs (xs, ys);
 352           pSGP->lineAbs (xd, yd);
 353         }
 354 #endif
 355         sum += projectSingleLine (phm, xd, yd, xs, ys, pSGP);
 356
 357         if (m_trace >= TRACE_RAYS)
 358           traceShowParam ("Attenuation  :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, "sum");
 359
 360 #ifdef HAVE_SGP
 361         if (pSGP && m_trace >= TRACE_RAYS) {
 362           pSGP->moveAbs (xs, ys);
 363           pSGP->lineAbs (xd, yd);
 364         }
 365 #endif
 366         xd += ddx2;
 367         yd += ddy2;
 368       }
 369
 370       detval[d] = sum / m_nSample;
 371       xd_maj += ddx;
 372       yd_maj += ddy;
 373       xs_maj += sdx;
 374       ys_maj += sdy;
 375     } /* for each detector */
 376   } /* if not unit pulse */
 377 }
 378
 379
 380 void
 381 Scanner::traceShowParam (const char *label, const char *fmt, int row, int color, ...)
 382 {
 383   char s[80];
 384   va_list arg;
 385
 386   va_start(arg, color);
 387   //  cio_set_cpos (raysum_trace_menu_column, row);
 388   snprintf (s, sizeof(s), label, "%s");
 389   //  cio_set_text_clr (color - 8, 0);
 390   cio_put_str (s);
 391   vsnprintf (s, sizeof(s), fmt, arg);
 392   //  cio_set_text_clr (color, 0);
 393   cio_put_str (s);
 394   cout << "\n";
 395   va_end(arg);
 396 }
 397
 398
 399 /* NAME
 400  *    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
 401  *
 402  * SYNOPSIS
 403  *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 404  *    double rsum               Ray sum of Phantom along given line
 405  *    Phantom& phm;             Phantom from which to calculate raysum
 406  *    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
 407  */
 408
 409 double
 410 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2, SGP* pSGP)
 411 {
 412   // check ray against each pelem in Phantom
 413   double rsum = 0.0;
 414   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 415     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2, pSGP);
 416
 417   return (rsum);
 418 }
 419
 420
 421 /* NAME
 422  *   pelem_ray_attenuation              Calculate raysum of an pelem along one line
 423  *
 424  * SYNOPSIS
 425  *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 426  *   double rsum                Computed raysum
 427  *   PhantomElement& pelem              Pelem to scan
 428  *   double x1, y1, x2, y2      Endpoints of raysum line
 429  */
 430
 431 double
 432 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2, SGP* pSGP)
 433 {
 434   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 435     if (m_trace == TRACE_CLIPPING)
 436       cio_tone (1000., 0.05);
 437     return (0.0);
 438   }
 439
 440 #ifdef HAVE_SGP
 441   if (pSGP && m_trace == TRACE_CLIPPING) {
 442     pSGP->moveAbs (x1, y1);
 443     pSGP->lineAbs (x2, y2);
 444     cio_tone (8000., 0.05);
 445     pSGP->moveAbs (x1, y1);
 446     pSGP->lineAbs (x2, y2);
 447   }
 448 #endif
 449
 450   double len = lineLength (x1, y1, x2, y2);
 451   return (len * pelem.atten());
 452 }
 453