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-2009 Kevin Rosenberg
  11 **
  12 **  This program is free software; you can redistribute it and/or modify
  13 **  it under the terms of the GNU General Public License (version 2) as
  14 **  published by the Free Software Foundation.
  15 **
  16 **  This program is distributed in the hope that it will be useful,
  17 **  but WITHOUT ANY WARRANTY; without even the implied warranty of
  18 **  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  19 **  GNU General Public License for more details.
  20 **
  21 **  You should have received a copy of the GNU General Public License
  22 **  along with this program; if not, write to the Free Software
  23 **  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  24 ******************************************************************************/
  25
  26 #include "ct.h"
  27
  28
  29 const int Scanner::GEOMETRY_INVALID = -1;
  30 const int Scanner::GEOMETRY_PARALLEL = 0;
  31 const int Scanner::GEOMETRY_EQUIANGULAR = 1;
  32 const int Scanner::GEOMETRY_EQUILINEAR = 2;
  33 const int Scanner::GEOMETRY_LINOGRAM = 3;
  34
  35 const char* const Scanner::s_aszGeometryName[] =
  36 {
  37   "parallel",
  38   "equiangular",
  39   "equilinear",
  40   "linogram",
  41 };
  42
  43 const char* const Scanner::s_aszGeometryTitle[] =
  44 {
  45   "Parallel",
  46   "Equiangular",
  47   "Equilinear",
  48   "Linogram",
  49 };
  50
  51 const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
  52
  53
  54 // NAME
  55 //   DetectorArray       Construct a DetectorArray
  56
  57 DetectorArray::DetectorArray (const int nDet)
  58 {
  59   m_nDet = nDet;
  60   m_detValues = new DetectorValue [m_nDet];
  61 }
  62
  63
  64 // NAME
  65 //   ~DetectorArray             Free memory allocated to a detector array
  66
  67 DetectorArray::~DetectorArray (void)
  68 {
  69   delete [] m_detValues;
  70 }
  71
  72
  73
  74 /* NAME
  75 *   Scanner::Scanner            Construct a user specified detector structure
  76 *
  77 * SYNOPSIS
  78 *   Scanner (phm, nDet, nView, nSample)
  79 *   Phantom& phm                PHANTOM that we are making detector for
  80 *   int geomety                Geometry of detector
  81 *   int nDet                    Number of detector along detector array
  82 *   int nView                   Number of rotated views
  83 *   int nSample         Number of rays per detector
  84 */
  85
  86 Scanner::Scanner (const Phantom& phm, const char* const geometryName,
  87                   int nDet, int nView, int offsetView,
  88                                   int nSample, const double rot_anglen,
  89                   const double dFocalLengthRatio,
  90                                   const double dCenterDetectorRatio,
  91                   const double dViewRatio, const double dScanRatio)
  92 {
  93   m_fail = false;
  94   m_idGeometry = convertGeometryNameToID (geometryName);
  95   if (m_idGeometry == GEOMETRY_INVALID) {
  96     m_fail = true;
  97     m_failMessage = "Invalid geometry name ";
  98     m_failMessage += geometryName;
  99     return;
 100   }
 101
 102   if (nView < 1 || nDet < 1) {
 103     m_fail = true;
 104     m_failMessage = "nView & nDet must be greater than 0";
 105     return;
 106   }
 107   if (nSample < 1)
 108     m_nSample = 1;
 109
 110   m_nDet     = nDet;
 111   m_nView    = nView;
 112   m_iOffsetView = offsetView;
 113   m_nSample  = nSample;
 114   m_dFocalLengthRatio = dFocalLengthRatio;
 115   m_dCenterDetectorRatio = dCenterDetectorRatio;
 116   m_dViewRatio = dViewRatio;
 117   m_dScanRatio = dScanRatio;
 118
 119   m_dViewDiameter = phm.getDiameterBoundaryCircle() * m_dViewRatio;
 120   m_dFocalLength = (m_dViewDiameter / 2) * m_dFocalLengthRatio;
 121   m_dCenterDetectorLength = (m_dViewDiameter / 2) * m_dCenterDetectorRatio;
 122   m_dSourceDetectorLength = m_dFocalLength + m_dCenterDetectorLength;
 123   m_dScanDiameter = m_dViewDiameter * m_dScanRatio;
 124
 125   m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
 126   m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
 127   m_rotLen  = rot_anglen;
 128   m_rotInc  = m_rotLen / m_nView;
 129   if (m_idGeometry == GEOMETRY_PARALLEL) {
 130     m_dFanBeamAngle = 0;
 131     m_detLen   = m_dScanDiameter;
 132     m_detStart = -m_detLen / 2;
 133     m_detInc  = m_detLen / m_nDet;
 134     double dDetectorArrayEndOffset = 0;
 135     // For even number of detectors, make detInc slightly larger so that center lies
 136     // at nDet/2. Also, extend detector array by one detInc so that all of the phantom is scanned
 137     if (isEven (m_nDet)) { // Adjust for Even number of detectors
 138       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
 139       dDetectorArrayEndOffset = m_detInc;
 140     }
 141
 142     double dHalfDetLen = m_detLen / 2;
 143     m_initPos.xs1 = m_dXCenter - dHalfDetLen;
 144     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
 145     m_initPos.xs2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
 146     m_initPos.ys2 = m_dYCenter + m_dFocalLength;
 147     m_initPos.xd1 = m_dXCenter - dHalfDetLen;
 148     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
 149     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
 150     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
 151     m_initPos.angle = m_iOffsetView * m_rotInc;
 152     m_detLen += dDetectorArrayEndOffset;
 153   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
 154   if (m_dScanDiameter / 2 >= m_dFocalLength) {
 155       m_fail = true;
 156       m_failMessage = "Invalid geometry: Focal length must be larger than scan length";
 157       return;
 158     }
 159
 160     const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
 161     const double dHalfDetLen = m_dSourceDetectorLength * tan (dAngle);
 162
 163     m_detLen = dHalfDetLen * 2;
 164     m_detStart = -dHalfDetLen;
 165     m_detInc  = m_detLen / m_nDet;
 166     double dDetectorArrayEndOffset = 0;
 167     if (isEven (m_nDet)) { // Adjust for Even number of detectors
 168       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
 169       dDetectorArrayEndOffset = m_detInc;
 170       m_detLen += dDetectorArrayEndOffset;
 171     }
 172
 173     m_dFanBeamAngle = dAngle * 2;
 174     m_initPos.xs1 = m_dXCenter;
 175     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
 176     m_initPos.xs2 = m_dXCenter;
 177     m_initPos.ys2 = m_dYCenter + m_dFocalLength;
 178     m_initPos.xd1 = m_dXCenter - dHalfDetLen;
 179     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
 180     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
 181     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
 182     m_initPos.angle = m_iOffsetView * m_rotInc;
 183   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 184     if (m_dScanDiameter / 2 > m_dFocalLength) {
 185       m_fail = true;
 186       m_failMessage = "Invalid geometry: Focal length must be larger than scan length";
 187       return;
 188     }
 189     const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
 190
 191     m_detLen = 2 * dAngle;
 192     m_detStart = -dAngle;
 193     m_detInc = m_detLen / m_nDet;
 194     double dDetectorArrayEndOffset = 0;
 195     if (isEven (m_nDet)) { // Adjust for Even number of detectors
 196       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
 197       dDetectorArrayEndOffset = m_detInc;
 198     }
 199     // adjust for center-detector length
 200     double dA1 = acos ((m_dScanDiameter / 2) / m_dCenterDetectorLength);
 201     double dAngularScale = 2 * (HALFPI + dAngle - dA1) / m_detLen;
 202
 203     m_dAngularDetLen = dAngularScale * (m_detLen + dDetectorArrayEndOffset);
 204     m_dAngularDetIncrement = dAngularScale * m_detInc;
 205     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
 206
 207     m_dFanBeamAngle = dAngle * 2;
 208     m_initPos.angle = m_iOffsetView * m_rotInc;
 209     m_initPos.xs1 = m_dXCenter;
 210     m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
 211     m_initPos.xs2 = m_dXCenter;
 212     m_initPos.ys2 = m_dYCenter + m_dFocalLength;
 213     m_detLen += dDetectorArrayEndOffset;
 214   }
 215
 216   // Calculate incrementatal rotation matrix
 217   GRFMTX_2D temp;
 218   xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter);
 219   rot_mtx2 (temp, m_rotInc);
 220   mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
 221   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
 222   mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
 223
 224 }
 225
 226 Scanner::~Scanner (void)
 227 {
 228 }
 229
 230
 231 const char*
 232 Scanner::convertGeometryIDToName (const int geomID)
 233 {
 234   const char *name = "";
 235
 236   if (geomID >= 0 && geomID < s_iGeometryCount)
 237     return (s_aszGeometryName[geomID]);
 238
 239   return (name);
 240 }
 241
 242 const char*
 243 Scanner::convertGeometryIDToTitle (const int geomID)
 244 {
 245   const char *title = "";
 246
 247   if (geomID >= 0 && geomID < s_iGeometryCount)
 248     return (s_aszGeometryName[geomID]);
 249
 250   return (title);
 251 }
 252
 253 int
 254 Scanner::convertGeometryNameToID (const char* const geomName)
 255 {
 256   int id = GEOMETRY_INVALID;
 257
 258   for (int i = 0; i < s_iGeometryCount; i++) {
 259     if (strcasecmp (geomName, s_aszGeometryName[i]) == 0) {
 260       id = i;
 261       break;
 262     }
 263   }
 264   return (id);
 265 }
 266
 267
 268 /* NAME
 269 *   collectProjections          Calculate projections for a Phantom
 270 *
 271 * SYNOPSIS
 272 *   collectProjections (proj, phm, start_view, nView, bStoreViewPos, trace)
 273 *   Projectrions& proj      Projection storage
 274 *   Phantom& phm             Phantom for which we collect projections
 275 *   bool bStoreViewPos      TRUE then storage proj at normal view position
 276 *   int trace                Trace level
 277 */
 278
 279
 280 void
 281 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP)
 282 {
 283   collectProjections (proj, phm, m_startView, proj.nView(), m_iOffsetView, true, trace, pSGP);
 284 }
 285
 286 void
 287 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView,
 288                              const int iNumViews, const int iOffsetView,  bool bStoreAtViewPosition,
 289                              const int trace, SGP* pSGP)
 290 {
 291   int iStorageOffset = (bStoreAtViewPosition ? iStartView : 0);
 292   collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP);
 293 }
 294
 295 static void mtx2_offset_rot (GRFMTX_2D m, double angle, double x, double y) {
 296   GRFMTX_2D temp;
 297   xlat_mtx2 (m, -x, -y);
 298   rot_mtx2 (temp, angle);
 299   mult_mtx2 (m, temp, m);
 300   xlat_mtx2 (temp, x, y);
 301   mult_mtx2 (m, temp, m);
 302 }
 303
 304 void
 305 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView,
 306                              const int iNumViews, const int iOffsetView, int iStorageOffset,
 307                              const int trace, SGP* pSGP)
 308 {
 309   m_trace = trace;
 310   double start_angle = (iStartView + iOffsetView) * proj.rotInc();
 311   int parallel_enabled = 1;
 312   UNUSED(parallel_enabled);
 313
 314 #if HAVE_SGP
 315   if (pSGP && (m_trace >= Trace::TRACE_PHANTOM))
 316     parallel_enabled = 0;
 317 #endif
 318
 319 #if HAVE_OPENMP
 320   #pragma omp parallel for if (parallel_enabled)
 321 #endif
 322   for (int iView = 0;  iView < iNumViews; iView++) {
 323     double viewAngle = start_angle + (iView * proj.rotInc());
 324
 325     // With OpenMP, need to calculate source and detector positions at each view
 326     GRFMTX_2D rotmtx;
 327     mtx2_offset_rot (rotmtx, viewAngle, m_dXCenter, m_dYCenter);
 328     double xd1=0, yd1=0, xd2=0, yd2=0;
 329     if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
 330       xd1 = m_initPos.xd1; yd1 = m_initPos.yd1;
 331       xd2 = m_initPos.xd2; yd2 = m_initPos.yd2;
 332       xform_mtx2 (rotmtx, xd1, yd1);      // rotate detector endpoints
 333       xform_mtx2 (rotmtx, xd2, yd2);      // to initial view_angle
 334     }
 335
 336     double xs1 = m_initPos.xs1, ys1 = m_initPos.ys1;
 337     double xs2 = m_initPos.xs2, ys2 = m_initPos.ys2;
 338     xform_mtx2 (rotmtx, xs1, ys1);      // rotate source endpoints to
 339     xform_mtx2 (rotmtx, xs2, ys2);      // initial view angle
 340
 341     int iStoragePosition = iView + iStorageOffset;
 342     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
 343
 344 #ifdef HAVE_SGP
 345     if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
 346       m_pSGP = pSGP;
 347       double dWindowSize = dmax (m_detLen, m_dSourceDetectorLength) * 2;
 348       double dHalfWindowSize = dWindowSize / 2;
 349       m_dXMinWin = m_dXCenter - dHalfWindowSize;
 350       m_dXMaxWin = m_dXCenter + dHalfWindowSize;
 351       m_dYMinWin = m_dYCenter - dHalfWindowSize;
 352       m_dYMaxWin = m_dYCenter + dHalfWindowSize;
 353
 354       m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
 355       m_pSGP->setRasterOp (RO_COPY);
 356
 357       m_pSGP->setColor (C_RED);
 358       m_pSGP->moveAbs (0., 0.);
 359       m_pSGP->drawCircle (m_dViewDiameter / 2);
 360
 361       m_pSGP->moveAbs (0., 0.);
 362       m_pSGP->setColor (C_GREEN);
 363       m_pSGP->drawCircle (m_dFocalLength);
 364       m_pSGP->setColor (C_BLUE);
 365       m_pSGP->setTextPointSize (9);
 366       phm.draw (*m_pSGP);
 367       m_dTextHeight = m_pSGP->getCharHeight ();
 368
 369       traceShowParam ("Phantom:",       "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str());
 370       traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry));
 371       traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio);
 372       //      traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio);
 373       traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet());
 374       traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView());
 375       traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
 376
 377       m_pSGP->setMarker (SGP::MARKER_BDIAMOND);
 378
 379       m_pSGP->setColor (C_BLACK);
 380       m_pSGP->setPenWidth (2);
 381       if (m_idGeometry == GEOMETRY_PARALLEL) {
 382         m_pSGP->moveAbs (xs1, ys1);
 383         m_pSGP->lineAbs (xs2, ys2);
 384         m_pSGP->moveAbs (xd1, yd1);
 385         m_pSGP->lineAbs (xd2, yd2);
 386       } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
 387         m_pSGP->setPenWidth (4);
 388         m_pSGP->moveAbs (xs1, ys1);
 389         m_pSGP->lineAbs (xs2, ys2);
 390         m_pSGP->setPenWidth (2);
 391         m_pSGP->moveAbs (xd1, yd1);
 392         m_pSGP->lineAbs (xd2, yd2);
 393       } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 394         m_pSGP->setPenWidth (4);
 395         m_pSGP->moveAbs (xs1, ys1);
 396         m_pSGP->lineAbs (xs2, ys2);
 397         m_pSGP->setPenWidth (2);
 398         m_pSGP->moveAbs (0., 0.);
 399         m_pSGP->drawArc (m_dCenterDetectorLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
 400       }
 401       m_pSGP->setPenWidth (1);
 402     }
 403     if (m_trace > Trace::TRACE_CONSOLE)
 404       traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast<double>(m_nView) * 100.);
 405 #endif
 406
 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
 419   } /* for each iView */
 420 }
 421
 422
 423 /* NAME
 424 *    rayview                    Calculate raysums for a view at any angle
 425 *
 426 * SYNOPSIS
 427 *    rayview (phm, detArray, xd1, nSample, yd1, xd2, yd2, xs1, ys1, xs2, ys2)
 428 *    Phantom& phm               Phantom to scan
 429 *    DETARRAY *detArray         Storage of values for detector array
 430 *    Scanner& det               Scanner parameters
 431 *    double xd1, yd1, xd2, yd2  Beginning & ending detector positions
 432 *    double xs1, ys1, xs2, ys2  Beginning & ending source positions
 433 *
 434 * RAY POSITIONING
 435 *         For each detector, have there are a variable number of rays traced.
 436 *     The source of each ray is the center of the source x-ray cell. The
 437 *     detector positions are equally spaced within the cell
 438 *
 439 *         The increments between rays are calculated so that the cells start
 440 *     at the beginning of a detector cell and they end on the endpoint
 441 *     of the cell.  Thus, the last cell starts at (xd2-ddx),(yd2-ddy).
 442 *         The exception to this is if there is only one ray per detector.
 443 *     In that case, the detector position is the center of the detector cell.
 444 */
 445
 446 void
 447 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)
 448 {
 449
 450   double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords
 451   double sdy = (ys2 - ys1) / detArray.nDet();  // between source
 452   double xs_maj = xs1 + (sdx / 2);      // put ray source in center of cell
 453   double ys_maj = ys1 + (sdy / 2);
 454
 455   double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0;
 456   double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0;
 457   if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 458     dAngleInc = m_dAngularDetIncrement;
 459     dAngleSampleInc = dAngleInc / m_nSample;
 460     dAngleSampleOffset = dAngleSampleInc / 2;
 461     dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
 462   } else {
 463     ddx = (xd2 - xd1) / detArray.nDet();  // change in coords
 464     ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
 465     ddx2 = ddx / m_nSample;     // Incr. between rays with detector cell
 466     ddy2 = ddy / m_nSample;  // Doesn't include detector endpoints
 467     ddx2_ofs = ddx2 / 2;    // offset of 1st ray from start of detector cell
 468     ddy2_ofs = ddy2 / 2;
 469
 470     xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
 471     yd_maj = yd1 + ddy2_ofs;
 472   }
 473
 474   DetectorValue* detval = detArray.detValues();
 475
 476   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
 477     for (int d = 0; d < detArray.nDet(); d++)
 478         detval[d] = 0;
 479     detval[ detArray.nDet() / 2 ] = 1;
 480   } else {
 481     for (int d = 0; d < detArray.nDet(); d++) {
 482       double xs = xs_maj;
 483       double ys = ys_maj;
 484       double xd=0, yd=0, dAngle=0;
 485       if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 486         dAngle = dAngleMajor;
 487       } else {
 488         xd = xd_maj;
 489         yd = yd_maj;
 490       }
 491       double sum = 0.0;
 492       for (unsigned int i = 0; i < m_nSample; i++) {
 493         if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
 494           xd = m_dCenterDetectorLength * cos (dAngle);
 495           yd = m_dCenterDetectorLength * sin (dAngle);
 496         }
 497
 498 #ifdef HAVE_SGP
 499         if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
 500           m_pSGP->setColor (C_YELLOW);
 501           m_pSGP->setRasterOp (RO_AND);
 502           m_pSGP->moveAbs (xs, ys);
 503           m_pSGP->lineAbs (xd, yd);
 504         }
 505 #endif
 506
 507         sum += projectSingleLine (phm, xd, yd, xs, ys);
 508
 509 #ifdef HAVE_SGP
 510         //      if (m_trace >= Trace::TRACE_CLIPPING) {
 511         //        traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
 512         //        traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
 513         //      }
 514 #endif
 515         if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 516           dAngle += dAngleSampleInc;
 517         else {
 518           xd += ddx2;
 519           yd += ddy2;
 520         }
 521       } // for each sample in detector
 522
 523       detval[d] = sum / m_nSample;
 524       xs_maj += sdx;
 525       ys_maj += sdy;
 526       if (m_idGeometry == GEOMETRY_EQUIANGULAR)
 527         dAngleMajor += dAngleInc;
 528       else {
 529         xd_maj += ddx;
 530         yd_maj += ddy;
 531       }
 532     } /* for each detector */
 533   } /* if not unit pulse */
 534 }
 535
 536
 537 void
 538 Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int color, ...)
 539 {
 540   va_list arg;
 541   va_start(arg, color);
 542 #ifdef HAVE_SGP
 543   traceShowParamRasterOp (RO_COPY, szLabel, fmt, row, color, arg);
 544 #else
 545   traceShowParamRasterOp (0, szLabel, fmt, row, color, arg);
 546 #endif
 547   va_end(arg);
 548 }
 549
 550 void
 551 Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int color, ...)
 552 {
 553   va_list arg;
 554   va_start(arg, color);
 555 #ifdef HAVE_SGP
 556   traceShowParamRasterOp (RO_XOR, szLabel, fmt, row, color, arg);
 557 #else
 558   traceShowParamRasterOp (0, szLabel, fmt, row, color, arg);
 559 #endif
 560   va_end(arg);
 561 }
 562
 563 void
 564 Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args)
 565 {
 566   char szValue[256];
 567
 568   vsnprintf (szValue, sizeof(szValue), fmt, args);
 569
 570 #ifdef HAVE_SGP
 571   if (m_pSGP) {
 572     m_pSGP->setRasterOp (iRasterOp);
 573     m_pSGP->setTextColor (color, -1);
 574     double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 4;
 575     if (row < 4) {
 576       double dYPos = m_dYMaxWin - (row * m_dTextHeight);
 577       double dXPos = m_dXMinWin;
 578       m_pSGP->moveAbs (dXPos, dYPos);
 579       m_pSGP->drawText (szLabel);
 580       m_pSGP->moveAbs (dXPos + dValueOffset, dYPos);
 581       m_pSGP->drawText (szValue);
 582     } else {
 583       row -= 4;
 584       double dYPos = m_dYMaxWin - (row * m_dTextHeight);
 585       double dXPos = m_dXMinWin + (m_dXMaxWin - m_dXMinWin) * 0.5;
 586       m_pSGP->moveAbs (dXPos, dYPos);
 587       m_pSGP->drawText (szLabel);
 588       m_pSGP->moveAbs (dXPos + dValueOffset, dYPos);
 589       m_pSGP->drawText (szValue);
 590     }
 591   } else
 592 #endif
 593   {
 594     cio_put_str (szLabel);
 595     cio_put_str (szValue);
 596     cio_put_str ("\n");
 597   }
 598 }
 599
 600 void swap_xy_points (double& x1, double& y1, double& x2, double& y2)
 601 {
 602   double temp = x1; x1 = x2; x2 = temp;
 603   temp = y1; y1 = y2; y2 = temp;
 604 }
 605
 606 class WeightedPoint {
 607 public:
 608   int x, y;
 609   double weight;
 610   WeightedPoint (int _x, int _y, double _weight)
 611     : x(_x), y(_y), weight(_weight)
 612   {}
 613 };
 614
 615
 616 /* FUNCTION
 617  *  Name:    projection_pixel_weights
 618  *  Purpose: Returns a vector of WeightedPoint with the length of
 619  *           line that intersects with each pixel
 620  */
 621
 622 void
 623 projection_pixel_weights (std::vector<WeightedPoint>& wp, const int nx, const int ny,
 624                           double x1, double y1, double x2, double y2)
 625 {
 626   double ylen = fabs(y2-y1);
 627   double xlen = fabs(x2-x1);
 628   bool swap_xy = false, invert_slope = false;
 629   double slope;
 630
 631   if (ylen > xlen) {
 632     swap_xy = true;
 633     slope = xlen / ylen;
 634     if (y2 < y1)      // swap start/end so always moving from bottom to top
 635       swap_xy_points (x1, y1, x2, y2);
 636     if (x2 < x1) {
 637       invert_slope = true;
 638     }
 639   } else {
 640 #if DEBUG
 641     if (ylen == xlen)
 642       sys_error(ERR_WARNING, "Slope == 1");
 643 #endif
 644     slope = ylen / xlen;
 645     if (x2 < x1)      // swap start/end so always moving from left to right in image
 646       swap_xy_points (x1, y1, x2, y2);
 647     if (y2 < y1) {
 648       invert_slope = true;
 649     }
 650   }
 651   double angle = atan(fabs(slope));
 652   double minor_dist = sin(angle); // distance along minor axis
 653   double pixel_len = 1 / cos(angle);
 654
 655   int minor_dir = 1;
 656   if (invert_slope) {
 657     minor_dir = -1;
 658     slope = -slope;
 659   }
 660
 661   double x = x1, y = y1;
 662   int ix = floor(x);
 663   int iy = floor(y);
 664   double ydelta = y - iy;
 665   double xdelta = x - ix;
 666
 667   double min_delta;
 668   int *imaj, *imin;
 669   int max_maj, max_min;
 670   if (swap_xy) {
 671     min_delta = xdelta;
 672     imaj = &iy;
 673     imin = &ix;
 674     max_maj = ny;
 675     max_min = nx;
 676   } else {
 677     min_delta = ydelta;
 678     imaj = &ix;
 679     imin = &iy;
 680     max_maj = nx;
 681     max_min = ny;
 682   }
 683
 684 #if DEBUG
 685   sys_error(ERR_TRACE, "m=%6.3f swap_xy=%d invert=%d len=%8.6f min_delta=%.4g minor_dist=%6.3f (%.3f,%.3f)-(%.3f,%.3f)",
 686             slope, swap_xy, invert_slope, pixel_len, min_delta, minor_dist, x1, y1, x2, y2);
 687 #endif
 688
 689   // if position of minor axis is at edge of image, but will be moving into pixel within image
 690   if (*imin == max_min && invert_slope) {
 691     (*imin)--;  // select the pixel within image
 692 #if DEBUG
 693     sys_error(ERR_TRACE, "Moving pixel inside image, adding %f to min_delta", (1+slope));
 694 #endif
 695     min_delta += (1+slope);
 696   }
 697
 698   while (*imaj < max_maj && *imin < max_min && *imin >= 0) {
 699     double next_min_delta = min_delta + slope;
 700
 701     if (((!invert_slope) && (next_min_delta < 1)) ||
 702         (invert_slope && (next_min_delta > 0))) {
 703       // stay within same pixel
 704       double w = pixel_len;
 705       WeightedPoint p (ix, iy, w);
 706       wp.push_back(p);
 707 #if DEBUG
 708       sys_error(ERR_TRACE, "  Full pixel: (%3d,%3d)=%.4g, min_delta=%.4g", ix, iy, w, min_delta);
 709 #endif
 710       min_delta = next_min_delta;
 711     } else {
 712       // Scale partial pixel_len into pixel
 713       double norm_delta = invert_slope ? min_delta : (1 - min_delta);
 714       double p1_line = norm_delta * pixel_len;
 715       WeightedPoint p1 (ix, iy, p1_line);
 716       wp.push_back (p1);
 717 #if DEBUG
 718       sys_error(ERR_TRACE, "  Part pixel: (%3d,%3d)=%.4g, min_delta=%.4g", ix, iy, p1_line, min_delta);
 719 #endif
 720       (*imin) += minor_dir;
 721       min_delta = next_min_delta - minor_dir;
 722     }
 723     (*imaj)++;
 724   }
 725
 726 }
 727
 728
 729 /* NAME
 730 *    projectSingleLine                  INTERNAL: Calculates raysum along a line for a Phantom
 731 *
 732 * SYNOPSIS
 733 *    rsum = phm_ray_attenuation (phm, x1, y1, x2, y2)
 734 *    double rsum                Ray sum of Phantom along given line
 735 *    Phantom& phm;              Phantom from which to calculate raysum
 736 *    double *x1, *y1, *x2, y2   Endpoints of ray path (in Phantom coords)
 737 */
 738
 739 double
 740 Scanner::projectSingleLine (const Phantom& phm, double x1, double y1, double x2, double y2)
 741 {
 742   double rsum = 0.0;
 743
 744   if (phm.isImagefile()) {
 745     // Project through an imagefile
 746
 747     const ImageFile* im = phm.getImagefile();
 748     const ImageFileArray v = im->getArray();
 749
 750     // Get image axis extents
 751     int nx = im->nx(), ny = im->ny();
 752     double xmin=0, xmax=nx, ymin=0, ymax=ny; // default coordinate
 753     if (! im->getAxisExtent (xmin, xmax, ymin, ymax)) {
 754       sys_error(ERR_WARNING, "Axis extent not available [Scanner::projectSingleLine]");
 755     }
 756
 757     // Clip line in image object coordinates
 758     double rect[4];
 759     rect[0] = xmin; rect[1] = ymin;
 760     rect[2] = xmax; rect[3] = ymax;
 761     bool accept = clip_rect (x1, y1, x2, y2, rect);
 762     if (! accept)
 763       return (0.0);
 764
 765     // Convert to pixel coordinates
 766     double xlen = xmax - xmin;
 767     double ylen = ymax - ymin;
 768     double px1 = nx * (x1 - xmin) / xlen;
 769     double px2 = nx * (x2 - xmin) / xlen;
 770     double py1 = ny * (y1 - ymin) / ylen;
 771     double py2 = ny * (y2 - ymin) / ylen;
 772
 773     std::vector<WeightedPoint> wp;
 774     projection_pixel_weights (wp, nx, ny, px1, py1, px2, py2);
 775     for (unsigned int i = 0; i < wp.size(); i++) {
 776       WeightedPoint& p = wp[i];
 777       rsum += v[p.x][p.y] * p.weight;
 778     }
 779   } else {
 780
 781     // Project through each pelem in Phantom
 782     for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
 783       rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
 784   }
 785   return (rsum);
 786 }
 787
 788
 789 /* NAME
 790 *   pelem_ray_attenuation               Calculate raysum of an pelem along one line
 791 *
 792 * SYNOPSIS
 793 *   rsum = pelem_ray_attenuation (pelem, x1, y1, x2, y2)
 794 *   double rsum         Computed raysum
 795 *   PhantomElement& pelem               Pelem to scan
 796 *   double x1, y1, x2, y2       Endpoints of raysum line
 797 */
 798
 799 double
 800 Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
 801 {
 802   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
 803     if (m_trace == Trace::TRACE_CLIPPING)
 804       cio_tone (1000., 0.05);
 805     return (0.0);
 806   }
 807
 808 #ifdef HAVE_SGP
 809   if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) {
 810     m_pSGP->setRasterOp (RO_XOR);
 811     m_pSGP->moveAbs (x1, y1);
 812     m_pSGP->lineAbs (x2, y2);
 813     cio_tone (8000., 0.05);
 814     m_pSGP->moveAbs (x1, y1);
 815     m_pSGP->lineAbs (x2, y2);
 816     m_pSGP->setRasterOp (RO_SET);
 817   }
 818 #endif
 819
 820   double len = lineLength (x1, y1, x2, y2);
 821   return (len * pelem.atten());
 822 }
 823