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