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