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