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.12 2000/08/27 20:32:55 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_dFocalLengthRatio = dFocalLengthRatio;
 111   m_dFieldOfViewRatio = dFieldOfViewRatio;
 112   m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio;
 113   m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio;
 114
 115   m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
 116   m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
 117   if (m_idGeometry == GEOMETRY_PARALLEL) {
 118     m_detLen   = m_dFieldOfView;
 119     m_detInc  = m_detLen / m_nDet;
 120     m_rotLen  = rot_anglen;
 121     m_rotInc  = m_rotLen / m_nView;
 122
 123     double dHalfDetLen = m_detLen / 2;
 124     m_initPos.xs1 = m_dXCenter - m_dFocalLength;
 125     m_initPos.ys1 = m_dYCenter - dHalfDetLen;
 126     m_initPos.xs2 = m_dXCenter - m_dFocalLength;
 127     m_initPos.ys2 = m_dYCenter + dHalfDetLen;
 128     m_initPos.xd1 = m_dXCenter + m_dFocalLength;
 129     m_initPos.yd1 = m_dYCenter - dHalfDetLen;
 130     m_initPos.xd2 = m_dXCenter + m_dFocalLength;
 131     m_initPos.yd2 = m_dYCenter + dHalfDetLen;
 132     m_initPos.angle = 0.0;
 133   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
 134     double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
 135     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
 136     if (dFocalPastPhm <= 0.) {
 137       m_fail = true;
 138       m_failMessage = "Focal Point inside of phantom";
 139       return;
 140     }
 141     double dAngle = atan( dHalfSquare / dFocalPastPhm );
 142     double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
 143
 144     m_detLen = dHalfDetLen * 2;
 145     m_detInc  = m_detLen / m_nDet;
 146     m_rotLen  = rot_anglen;
 147     m_rotInc  = m_rotLen / m_nView;
 148
 149     m_initPos.xs1 = m_dXCenter - m_dFocalLength;
 150     m_initPos.ys1 = m_dYCenter;
 151     m_initPos.xs2 = m_dXCenter - m_dFocalLength;
 152     m_initPos.ys2 = m_dYCenter;
 153     m_initPos.xd1 = m_dXCenter + m_dFocalLength;
 154     m_initPos.yd1 = m_dYCenter - dHalfDetLen;
 155     m_initPos.xd2 = m_dXCenter + m_dFocalLength;
 156     m_initPos.yd2 = m_dYCenter + dHalfDetLen;
 157     m_initPos.angle = 0.0;
 158   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 159     double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
 160     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
 161     if (dFocalPastPhm <= 0.) {
 162       m_fail = true;
 163       m_failMessage = "Focal Point inside of phantom";
 164       return;
 165     }
 166     double dAngle = 2 * atan( dHalfSquare / dFocalPastPhm );
 167     m_detLen = 2 * dAngle;
 168     m_detInc = m_detLen / m_nDet;
 169     m_rotLen = rot_anglen;
 170     m_rotInc = m_rotLen / m_nView;
 171
 172     m_initPos.xs1 = m_dXCenter - m_dFocalLength;
 173     m_initPos.ys1 = m_dYCenter;
 174     m_initPos.xs2 = m_dXCenter - m_dFocalLength;
 175     m_initPos.ys2 = m_dYCenter;
 176     m_initPos.angle = -dAngle;
 177   }
 178
 179   // Calculate incrementatal rotation matrix
 180   GRFMTX_2D temp;
 181   xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter);
 182   rot_mtx2 (temp, m_rotInc);
 183   mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
 184   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
 185   mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
 186
 187 }
 188
 189 Scanner::~Scanner (void)
 190 {
 191 }
 192
 193
 194 const char*
 195 Scanner::convertGeometryIDToName (const int geomID)
 196 {
 197   const char *name = "";
 198
 199   if (geomID >= 0 && geomID < s_iGeometryCount)
 200       return (s_aszGeometryName[geomID]);
 201
 202   return (name);
 203 }
 204
 205 const char*
 206 Scanner::convertGeometryIDToTitle (const int geomID)
 207 {
 208   const char *title = "";
 209
 210   if (geomID >= 0 && geomID < s_iGeometryCount)
 211       return (s_aszGeometryName[geomID]);
 212
 213   return (title);
 214 }
 215
 216 int
 217 Scanner::convertGeometryNameToID (const char* const geomName)
 218 {
 219   int id = GEOMETRY_INVALID;
 220
 221   for (int i = 0; i < s_iGeometryCount; i++)
 222       if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
 223           id = i;
 224           break;
 225       }
 226
 227   return (id);
 228 }
 229
 230
 231 /* NAME
 232  *   collectProjections         Calculate projections for a Phantom
 233  *
 234  * SYNOPSIS
 235  *   collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace)
 236  *   Projectrions& proj      Projection storage
 237  *   Phantom& phm            Phantom for which we collect projections
 238  *   bool bStoreViewPos      TRUE then storage proj at normal view position
 239  *   int trace               Trace level
 240 */
 241
 242
 243 void
 244 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP)
 245 {
 246   collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP);
 247 }
 248
 249 void
 250 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, bool bStoreAtViewPosition, const int trace, SGP* pSGP)
 251 {
 252   m_trace = trace;
 253   double start_angle = iStartView * proj.rotInc();
 254
 255   // Calculate initial rotation matrix
 256   GRFMTX_2D rotmtx_initial, temp;
 257   xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter);
 258   rot_mtx2 (temp, start_angle);
 259   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 260   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
 261   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 262
 263   double xd1=0, yd1=0, xd2=0, yd2=0, dDetAngle=0;
 264   if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 265       dDetAngle = m_initPos.angle + start_angle;
 266   else {
 267       xd1 = m_initPos.xd1;
 268       yd1 = m_initPos.yd1;
 269       xd2 = m_initPos.xd2;
 270       yd2 = m_initPos.yd2;
 271       xform_mtx2 (rotmtx_initial, xd1, yd1);      // rotate detector endpoints
 272       xform_mtx2 (rotmtx_initial, xd2, yd2);      // to initial view_angle
 273   }
 274
 275   double xs1 = m_initPos.xs1;
 276   double ys1 = m_initPos.ys1;
 277   double xs2 = m_initPos.xs2;
 278   double ys2 = m_initPos.ys2;
 279   xform_mtx2 (rotmtx_initial, xs1, ys1);      // rotate source endpoints to
 280   xform_mtx2 (rotmtx_initial, xs2, ys2);      // initial view angle
 281
 282   int iView;
 283   double viewAngle;
 284   for (iView = 0, viewAngle = start_angle;  iView < iNumViews; iView++, viewAngle += proj.rotInc()) {
 285     int iStoragePosition = iView;
 286     if (bStoreAtViewPosition)
 287       iStoragePosition += iStartView;
 288
 289     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
 290
 291 #ifdef HAVE_SGP
 292   if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
 293       m_pSGP = pSGP;
 294       m_pSGP->eraseWindow();
 295       double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2;
 296       double dHalfWindowSize = dWindowSize / 2;
 297       m_dXMinWin = m_dXCenter - dHalfWindowSize;
 298       m_dXMaxWin = m_dXCenter + dHalfWindowSize;
 299       m_dYMinWin = m_dYCenter - dHalfWindowSize;
 300       m_dYMaxWin = m_dYCenter + dHalfWindowSize;
 301       double dHalfPhmLen = m_phmLen /  2;
 302
 303     m_pSGP->eraseWindow ();
 304     m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
 305     m_pSGP->setRasterOp (RO_COPY);
 306     m_pSGP->setColor (C_BLUE);
 307     m_pSGP->moveAbs (0., 0.);
 308     m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen);
 309     m_pSGP->moveAbs (0., 0.);
 310     m_pSGP->drawCircle (m_dFocalLength);
 311     m_pSGP->setColor (C_RED);
 312     phm.draw (*m_pSGP);
 313     m_dTextHeight = m_pSGP->getCharHeight ();
 314
 315     traceShowParam ("Projection Collector", "%s", PROJECTION_TRACE_ROW_TITLE, C_BLACK, " ");
 316     traceShowParam ("________________", "%s", PROJECTION_TRACE_ROW_TITLE2, C_LTGRAY, " ");
 317     traceShowParam ("Phantom:",       "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str());
 318     traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry));
 319     traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio);
 320     traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio);
 321     traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet());
 322     traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView());
 323     traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
 324
 325     m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
 326   }
 327 #endif
 328
 329 #ifdef HAVE_SGP
 330     if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
 331       m_pSGP->setColor (C_BLACK);
 332       m_pSGP->setPenWidth (2);
 333       if (m_idGeometry == GEOMETRY_PARALLEL) {
 334         m_pSGP->moveAbs (xs1, ys1);
 335         m_pSGP->lineAbs (xs2, ys2);
 336         m_pSGP->moveAbs (xd1, yd1);
 337         m_pSGP->lineAbs (xd2, yd2);
 338       } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
 339         m_pSGP->setPenWidth (4);
 340         m_pSGP->moveAbs (xs1, ys1);
 341         m_pSGP->lineAbs (xs2, ys2);
 342         m_pSGP->setPenWidth (2);
 343         m_pSGP->moveAbs (xd1, yd1);
 344         m_pSGP->lineAbs (xd2, yd2);
 345       } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 346         m_pSGP->setPenWidth (4);
 347         m_pSGP->moveAbs (xs1, ys1);
 348         m_pSGP->lineAbs (xs2, ys2);
 349         m_pSGP->setPenWidth (2);
 350         m_pSGP->moveAbs (0., 0.);
 351         m_pSGP->drawArc (m_dFocalLength, start_angle + m_initPos.angle, start_angle - m_initPos.angle);
 352       }
 353       m_pSGP->setPenWidth (1);
 354     }
 355     if (m_trace >= Trace::TRACE_CONSOLE)
 356       traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / m_nView * 100.);
 357 #endif
 358
 359     projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, dDetAngle);
 360     detArray.setViewAngle (viewAngle);
 361
 362 #ifdef HAVE_SGP
 363     if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
 364       //        rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta);
 365     }
 366 #endif
 367     xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
 368     xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
 369     if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 370         dDetAngle += m_detInc;
 371     else {
 372         xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints
 373         xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
 374     }
 375   } /* for each iView */
 376 }
 377
 378
 379 /* NAME
 380  *    rayview                   Calculate raysums for a view at any angle
 381  *
 382  * SYNOPSIS
 383  *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 384  *    Phantom& phm              Phantom to scan
 385  *    DETARRAY *detArray                Storage of values for detector array
 386  *    Scanner& det              Scanner parameters
 387  *    double xd1, yd1, xd2, yd2 Beginning & ending detector positions
 388  *    double xs1, ys1, xs2, ys2 Beginning & ending source positions
 389  *
 390  * RAY POSITIONING
 391  *         For each detector, have there are a variable number of rays traced.
 392  *     The source of each ray is the center of the source x-ray cell. The
 393  *     detector positions are equally spaced within the cell
 394  *
 395  *         The increments between rays are calculated so that the cells start
 396  *     at the beginning of a detector cell and they end on the endpoint
 397  *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 398  *         The exception to this is if there is only one ray per detector.
 399  *     In that case, the detector position is the center of the detector cell.
 400  */
 401
 402 void
 403 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, const double dDetAngle)
 404 {
 405
 406   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords
 407   double sdy = (ys2 - ys1) / detArray.nDet();  // between source
 408   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 409   double ys_maj = ys1 + (sdy / 2);
 410
 411   double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0;
 412   double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0;
 413   if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 414       dAngleInc = m_detInc;
 415       dAngleSampleInc = dAngleInc / m_nSample;
 416       dAngleSampleOffset = dAngleSampleInc / 2;
 417       dAngleMajor = dDetAngle + dAngleSampleOffset;
 418   } else {
 419       ddx = (xd2 - xd1) / detArray.nDet();  // change in coords
 420       ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
 421       ddx2 = ddx / m_nSample;   // Incr. between rays with detector cell
 422       ddy2 = ddy / m_nSample;  // Doesn't include detector endpoints
 423       ddx2_ofs = ddx2 / 2;    // offset of 1st ray from start of detector cell
 424       ddy2_ofs = ddy2 / 2;
 425
 426       xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 427       yd_maj = yd1 + ddy2_ofs;
 428   }
 429
 430   DetectorValue* detval = detArray.detValues();
 431
 432   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 433     for (int d = 0; d < detArray.nDet(); d++)
 434       if (detArray.nDet() / 2 == d && (d % 2) == 1)
 435         detval[d] = 1;
 436       else
 437         detval[d] = 0;
 438   } else {
 439       for (int d = 0; d < detArray.nDet(); d++) {
 440       double xs = xs_maj;
 441       double ys = ys_maj;
 442       double xd=0, yd=0, dAngle=0;
 443       if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 444           dAngle = dAngleMajor;
 445       } else {
 446           xd = xd_maj;
 447           yd = yd_maj;
 448       }
 449       double sum = 0.0;
 450       for (unsigned int i = 0; i < m_nSample; i++) {
 451         if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 452             xd = m_dFocalLength * cos (dAngle);
 453             yd = m_dFocalLength * sin (dAngle);
 454         }
 455
 456 #ifdef HAVE_SGP
 457         if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
 458           m_pSGP->setColor (C_YELLOW);
 459           m_pSGP->setRasterOp (RO_OR_REVERSE);
 460           m_pSGP->moveAbs (xs, ys);
 461           m_pSGP->lineAbs (xd, yd);
 462           m_pSGP->setRasterOp (RO_SET);
 463         }
 464 #endif
 465
 466         sum += projectSingleLine (phm, xd, yd, xs, ys);
 467
 468 #ifdef HAVE_SGP
 469         if (m_trace >= Trace::TRACE_CLIPPING) {
 470           traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
 471           traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
 472         }
 473
 474         //      if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
 475         //        m_pSGP->setColor (C_YELLOW);
 476         //        m_pSGP->setRasterOp (RO_XOR);
 477         //        m_pSGP->moveAbs (xs, ys);
 478         //        m_pSGP->lineAbs (xd, yd);
 479         //        m_pSGP->setRasterOp (RO_SET);
 480         //      }
 481 #endif
 482         if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 483             dAngle += dAngleSampleInc;
 484         else {
 485             xd += ddx2;
 486             yd += ddy2;
 487         }
 488       } // for each sample in detector
 489
 490       detval[d] = sum / m_nSample;
 491       xs_maj += sdx;
 492       ys_maj += sdy;
 493       if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 494           dAngleMajor += dAngleInc;
 495         else {
 496             xd_maj += ddx;
 497             yd_maj += ddy;
 498         }
 499     } /* for each detector */
 500   } /* if not unit pulse */
 501 }
 502
 503
 504 void
 505 Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int color, ...)
 506 {
 507   va_list arg;
 508   va_start(arg, color);
 509   traceShowParamRasterOp (RO_COPY, szLabel, fmt, row, color, arg);
 510   va_end(arg);
 511 }
 512
 513 void
 514 Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int color, ...)
 515 {
 516   va_list arg;
 517   va_start(arg, color);
 518   traceShowParamRasterOp (RO_XOR, szLabel, fmt, row, color, arg);
 519   va_end(arg);
 520 }
 521
 522 void
 523 Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args)
 524 {
 525   char szValue[256];
 526
 527   vsnprintf (szValue, sizeof(szValue), fmt, args);
 528
 529   //  cio_set_cpos (raysum_trace_menu_column, row);
 530   //  cio_set_text_clr (color - 8, 0);
 531   //  cio_set_text_clr (color, 0);
 532
 533   if (m_pSGP) {
 534     m_pSGP->setRasterOp (iRasterOp);
 535     double dYPos = m_dYMaxWin - (row * m_dTextHeight);
 536     m_pSGP->moveAbs (m_dXMinWin, dYPos);
 537     m_pSGP->setTextColor (color, -1);
 538     m_pSGP->drawText (szLabel);
 539     double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 5;
 540     m_pSGP->moveAbs (m_dXMinWin + dValueOffset, dYPos);
 541     m_pSGP->drawText (szValue);
 542   } else {
 543     cio_put_str (szLabel);
 544     cio_put_str (szValue);
 545     cio_put_str ("\n");
 546   }
 547 }
 548
 549
 550
 551 /* NAME
 552  *    projectSingleLine                 INTERNAL: Calculates raysum along a line for a Phantom
 553  *
 554  * SYNOPSIS
 555  *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 556  *    double rsum               Ray sum of Phantom along given line
 557  *    Phantom& phm;             Phantom from which to calculate raysum
 558  *    double *x1, *y1, *x2, y2  Endpoints of ray path (in Phantom coords)
 559  */
 560
 561 double
 562 Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2)
 563 {
 564   // check ray against each pelem in Phantom
 565   double rsum = 0.0;
 566   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 567     rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
 568
 569   return (rsum);
 570 }
 571
 572
 573 /* NAME
 574  *   pelem_ray_attenuation              Calculate raysum of an pelem along one line
 575  *
 576  * SYNOPSIS
 577  *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 578  *   double rsum                Computed raysum
 579  *   PhantomElement& pelem              Pelem to scan
 580  *   double x1, y1, x2, y2      Endpoints of raysum line
 581  */
 582
 583 double
 584 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
 585 {
 586   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 587     if (m_trace == Trace::TRACE_CLIPPING)
 588       cio_tone (1000., 0.05);
 589     return (0.0);
 590   }
 591
 592 #ifdef HAVE_SGP
 593   if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) {
 594     m_pSGP->setRasterOp (RO_XOR);
 595     m_pSGP->moveAbs (x1, y1);
 596     m_pSGP->lineAbs (x2, y2);
 597     cio_tone (8000., 0.05);
 598     m_pSGP->moveAbs (x1, y1);
 599     m_pSGP->lineAbs (x2, y2);
 600     m_pSGP->setRasterOp (RO_SET);
 601   }
 602 #endif
 603
 604   double len = lineLength (x1, y1, x2, y2);
 605   return (len * pelem.atten());
 606 }
 607