X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=48f23943a739ec654d361fba8fdcc35f1fad92d5;hp=c4fc4d06a2b8041c246ba1a47aff538412d54c12;hb=9703ace3c451ce079967284bf191783736dbc77f;hpb=1e88cf0f7fa4f690ea9f110e8ed3f2b5338d0a10

diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index c4fc4d0..48f2394 100644
--- a/libctsim/scanner.cpp
+++ b/libctsim/scanner.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.11 2000/08/25 15:59:13 kevin Exp $
+**  $Id: scanner.cpp,v 1.18 2000/12/16 02:31:00 kevin Exp $
 **
 **  This program is free software; you can redistribute it and/or modify
 **  it under the terms of the GNU General Public License (version 2) as
@@ -30,21 +30,21 @@
 
 const int Scanner::GEOMETRY_INVALID = -1;
 const int Scanner::GEOMETRY_PARALLEL = 0;
-const int Scanner::GEOMETRY_EQUILINEAR = 1;
-const int Scanner::GEOMETRY_EQUIANGULAR = 2;
+const int Scanner::GEOMETRY_EQUIANGULAR = 1;
+const int Scanner::GEOMETRY_EQUILINEAR = 2;
 
 const char* Scanner::s_aszGeometryName[] = 
 {
   {"parallel"},
-  {"equilinear"},
   {"equiangular"},
+  {"equilinear"},
 };
 
 const char* Scanner::s_aszGeometryTitle[] = 
 {
   {"Parallel"},
-  {"Equilinear"},
   {"Equiangular"},
+  {"Equilinear"},
 };
 
 const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
@@ -95,39 +95,47 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     return;
   }
 
-  if (nView < 1)
-    nView = 1;
+  if (nView < 1 || nDet < 1) {
+    m_fail = true;
+    m_failMessage = "nView & nDet must be greater than 0";
+    return;
+  }
   if (nSample < 1)
     m_nSample = 1;
-  if (nDet < 1)
-    nDet = 1;
-  //  if ((nDet % 2) == 0)
-  //    ++nDet;		// ensure odd number of detectors
 
   m_nDet     = nDet;
   m_nView    = nView;
   m_nSample  = nSample;
+  m_dFocalLengthRatio = dFocalLengthRatio;
+  m_dFieldOfViewRatio = dFieldOfViewRatio;
   m_dFocalLength = (m_phmLen * SQRT2 / 2) * dFocalLengthRatio;
   m_dFieldOfView = m_phmLen * SQRT2 * dFieldOfViewRatio;
 
+  m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
+  m_dYCenter = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
+  m_rotLen  = rot_anglen;
+  m_rotInc  = m_rotLen / m_nView;
   if (m_idGeometry == GEOMETRY_PARALLEL) {
     m_detLen   = m_dFieldOfView;
     m_detInc  = m_detLen / m_nDet;
-    m_rotLen  = rot_anglen;
-    m_rotInc  = m_rotLen / m_nView;
+    if (m_nDet % 2 == 0) // Adjust for Even number of detectors
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1
     
     double dHalfDetLen = m_detLen / 2;
-    m_initPos.xs1 = -m_dFocalLength;
-    m_initPos.ys1 = -dHalfDetLen;
-    m_initPos.xs2 = -m_dFocalLength;
-    m_initPos.ys2 = dHalfDetLen;
-    m_initPos.xd1 = m_dFocalLength;
-    m_initPos.yd1 = -dHalfDetLen;
-    m_initPos.xd2 = m_dFocalLength;
-    m_initPos.yd2 = dHalfDetLen;
+    m_initPos.xs1 = m_dXCenter - dHalfDetLen;
+    m_initPos.ys1 = m_dYCenter + m_dFocalLength;
+    m_initPos.xs2 = m_dXCenter + dHalfDetLen;
+    m_initPos.ys2 = m_dYCenter + m_dFocalLength;
+    m_initPos.xd1 = m_dXCenter - dHalfDetLen;
+    m_initPos.yd1 = m_dYCenter - m_dFocalLength;
+    m_initPos.xd2 = m_dXCenter + dHalfDetLen;
+    m_initPos.yd2 = m_dYCenter - m_dFocalLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
-    double dHalfSquare = m_phmLen / 2;
+#if 0
+    double dAngle = (m_dFieldOfView / 2) / cos (asin (m_dFieldOfView / 2 / m_dFocalLength));
+#else
+    double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
     if (dFocalPastPhm <= 0.) {
       m_fail = true;
@@ -135,33 +143,60 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
       return;
     }
     double dAngle = atan( dHalfSquare / dFocalPastPhm );
+#endif
     double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
 
     m_detLen = dHalfDetLen * 2;
     m_detInc  = m_detLen / m_nDet;
-    m_rotLen  = rot_anglen;
-    m_rotInc  = m_rotLen / m_nView;
+    if (m_nDet % 2 == 0) // Adjust for Even number of detectors
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1
     
-    m_initPos.xs1 = -m_dFocalLength;
-    m_initPos.ys1 = 0;
-    m_initPos.xs2 = -m_dFocalLength;
-    m_initPos.ys2 = 0;
-    m_initPos.xd1 = m_dFocalLength;
-    m_initPos.yd1 = -dHalfDetLen;
-    m_initPos.xd2 = m_dFocalLength;
-    m_initPos.yd2 = dHalfDetLen;
+    m_initPos.angle = 0.0;
+    m_initPos.xs1 = m_dXCenter;
+    m_initPos.ys1 = m_dYCenter + m_dFocalLength;
+    m_initPos.xs2 = m_dXCenter;
+    m_initPos.ys2 = m_dYCenter + m_dFocalLength;
+    m_initPos.xd1 = m_dXCenter - dHalfDetLen;
+    m_initPos.yd1 = m_dYCenter - m_dFocalLength;
+    m_initPos.xd2 = m_dXCenter + dHalfDetLen;
+    m_initPos.yd2 = m_dYCenter - m_dFocalLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-    double dHalfSquare = m_phmLen / 2;
+#if 0
+    double dAngle = atan ((m_dFieldOfView / 2) / m_dFocalLength);
+#else
+    double dHalfSquare = m_dFieldOfView / SQRT2 / 2;
     double dFocalPastPhm = m_dFocalLength - dHalfSquare;
     if (dFocalPastPhm <= 0.) {
       m_fail = true;
       m_failMessage = "Focal Point inside of phantom";
       return;
     }
-    double dAngle = atan( dHalfSquare / dFocalPastPhm );
-    m_detLen = dAngle * 2;
+    double dAngle =  atan ( dHalfSquare / dFocalPastPhm ); 
+#endif
+    m_detLen = 2 * dAngle;
+    m_detInc = m_detLen / m_nDet;
+    if (m_nDet % 2 == 0) // Adjust for Even number of detectors
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)-1
+    m_dAngularDetIncrement = m_detInc * 2; // Angular Position 2x gamma angle
+    m_dAngularDetLen = m_detLen * 2;
+    m_initPos.dAngularDet = -m_dAngularDetLen / 2;
+
+    m_initPos.angle = 0;
+    m_initPos.xs1 = m_dXCenter;
+    m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
+    m_initPos.xs2 = m_dXCenter;
+    m_initPos.ys2 = m_dYCenter + m_dFocalLength;
   }
+
+  // Calculate incrementatal rotation matrix 
+  GRFMTX_2D temp;
+  xlat_mtx2 (m_rotmtxIncrement, -m_dXCenter, -m_dYCenter);
+  rot_mtx2 (temp, m_rotInc);
+  mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
+  xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
+  mult_mtx2 (m_rotmtxIncrement, temp, m_rotmtxIncrement);
+  
 }
 
 Scanner::~Scanner (void)
@@ -227,119 +262,126 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int tr
 void
 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, bool bStoreAtViewPosition, const int trace, SGP* pSGP)
 {
-  GRFMTX_2D rotmtx_initial, temp;
-  GRFMTX_2D rotmtx_incr;
-
-  double start_angle = iStartView * proj.rotInc();
-  double xcent = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
-  double ycent = phm.ymin() + (phm.ymax() - phm.ymin()) / 2;
-
-  double xd1 = xcent + m_initPos.xd1;
-  double yd1 = ycent + m_initPos.yd1;
-  double xd2 = xcent + m_initPos.xd2;
-  double yd2 = ycent + m_initPos.yd2;
-  double xs1 = xcent + m_initPos.xs1;
-  double ys1 = ycent + m_initPos.ys1;
-  double xs2 = xcent + m_initPos.xs2;
-  double ys2 = ycent + m_initPos.ys2;
- 
   m_trace = trace;
+  double start_angle = iStartView * proj.rotInc();
 
-/* Calculate initial rotation matrix */
-  xlat_mtx2 (rotmtx_initial, -xcent, -ycent);
+  // Calculate initial rotation matrix 
+  GRFMTX_2D rotmtx_initial, temp;
+  xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter);
   rot_mtx2 (temp, start_angle);
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
-  xlat_mtx2 (temp, xcent, ycent);
+  xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 
-  xform_mtx2 (rotmtx_initial, xd1, yd1);	/* rotate detector endpoints */
-  xform_mtx2 (rotmtx_initial, xd2, yd2);      /* to initial view_angle */
-  xform_mtx2 (rotmtx_initial, xs1, ys1);
-  xform_mtx2 (rotmtx_initial, xs2, ys2);
-
-/* Calculate incrementatal rotation matrix */
-  xlat_mtx2 (rotmtx_incr, -xcent, -ycent);
-  rot_mtx2 (temp, proj.rotInc());
-  mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
-  xlat_mtx2 (temp, xcent, ycent);
-  mult_mtx2 (rotmtx_incr, temp, rotmtx_incr);
+  double xd1=0, yd1=0, xd2=0, yd2=0;
+  if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
+      xd1 = m_initPos.xd1;
+      yd1 = m_initPos.yd1;
+      xd2 = m_initPos.xd2;
+      yd2 = m_initPos.yd2;
+      xform_mtx2 (rotmtx_initial, xd1, yd1);      // rotate detector endpoints 
+      xform_mtx2 (rotmtx_initial, xd2, yd2);      // to initial view_angle 
+  }
   
-  int iview;
+  double xs1 = m_initPos.xs1;
+  double ys1 = m_initPos.ys1;
+  double xs2 = m_initPos.xs2;
+  double ys2 = m_initPos.ys2;
+  xform_mtx2 (rotmtx_initial, xs1, ys1);      // rotate source endpoints to
+  xform_mtx2 (rotmtx_initial, xs2, ys2);      // initial view angle
+
+  int iView;
   double viewAngle;
-  for (iview = 0, viewAngle = start_angle;  iview < iNumViews; iview++, viewAngle += proj.rotInc()) {
-    int iStoragePosition = iview;
+  for (iView = 0, viewAngle = start_angle;  iView < iNumViews; iView++, viewAngle += proj.rotInc()) {
+    int iStoragePosition = iView;
     if (bStoreAtViewPosition)
       iStoragePosition += iStartView;
 
     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
 
 #ifdef HAVE_SGP 
-  if (pSGP && m_trace >= TRACE_PHM) {
-    pSGP->eraseWindow();
-    double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2;
-    double dHalfWindowSize = dWindowSize / 2;
-    double dHalfPhmLen = m_phmLen /  2;
-
-    pSGP->setRasterOp (RO_SET);
-    pSGP->eraseWindow ();
-    pSGP->setWindow (xcent - dHalfWindowSize, ycent - dHalfWindowSize, xcent + dHalfWindowSize, ycent + dHalfWindowSize);
-    pSGP->setColor (C_BROWN);
-    pSGP->moveAbs (0., 0.);
-    pSGP->drawRect (xcent - dHalfPhmLen, ycent - dHalfPhmLen, xcent + dHalfPhmLen, ycent + dHalfPhmLen);
-
-#if 0    
-    traceShowParam (pSGP, "X-Ray Simulator", "%s", RAYSUM_TRACE_ROW_TITLE, C_BLACK, " ");
-    traceShowParam (pSGP, "---------------", "%s", RAYSUM_TRACE_ROW_TITLE2, C_BLACK, " ");
-    traceShowParam (pSGP, "Phantom:",       "%s", RAYSUM_TRACE_ROW_PHANT_ID, C_YELLOW, " Herman");
-    traceShowParam (pSGP, "Chomaticity  :", "%s", RAYSUM_TRACE_ROW_CHROMATIC, C_LTGREEN, "Mono");
-    traceShowParam (pSGP, "Scatter      :", "%5.1f", RAYSUM_TRACE_ROW_SCATTER, C_LTGREEN, 0.);
-    traceShowParam (pSGP, "Photon Uncert:", "%5.1f", RAYSUM_TRACE_ROW_PHOT_STAT, C_LTGREEN, 0.);
-    traceShowParam (pSGP, "Num Detectors:", "%5d", RAYSUM_TRACE_ROW_NDET, C_LTRED, proj.nDet());
-    traceShowParam (pSGP, "Num Views    :", "%5d", RAYSUM_TRACE_ROW_NVIEW, C_LTRED, proj.nViews());
-    traceShowParam (pSGP, "Samples / Ray:", "%5d", RAYSUM_TRACE_ROW_SAMPLES, C_LTRED, m_nSample);
-#endif
+  if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
+      m_pSGP = pSGP;
+      m_pSGP->eraseWindow();
+      double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2;
+      double dHalfWindowSize = dWindowSize / 2;
+      m_dXMinWin = m_dXCenter - dHalfWindowSize;
+      m_dXMaxWin = m_dXCenter + dHalfWindowSize;
+      m_dYMinWin = m_dYCenter - dHalfWindowSize;
+      m_dYMaxWin = m_dYCenter + dHalfWindowSize;
+      double dHalfPhmLen = m_phmLen /  2;
+
+    m_pSGP->eraseWindow ();
+    m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
+    m_pSGP->setRasterOp (RO_COPY);
+    m_pSGP->setColor (C_RED);
+    m_pSGP->moveAbs (0., 0.);
+    m_pSGP->drawRect (m_dXCenter - dHalfPhmLen, m_dYCenter - dHalfPhmLen, m_dXCenter + dHalfPhmLen, m_dYCenter + dHalfPhmLen);
+    m_pSGP->moveAbs (0., 0.);
+    m_pSGP->drawCircle (m_dFocalLength);
+    m_pSGP->setColor (C_BLUE);
+    phm.draw (*m_pSGP);
+    m_dTextHeight = m_pSGP->getCharHeight ();
+
+    traceShowParam ("Phantom:",       "%s", PROJECTION_TRACE_ROW_PHANT_ID, C_BLACK, phm.name().c_str());
+    traceShowParam ("Geometry:", "%s", PROJECTION_TRACE_ROW_GEOMETRY, C_BLUE, convertGeometryIDToName(m_idGeometry));
+    traceShowParam ("Focal Length Ratio:", "%.2f", PROJECTION_TRACE_ROW_FOCAL_LENGTH, C_BLUE, m_dFocalLengthRatio);
+    traceShowParam ("Field Of View Ratio:", "%.2f", PROJECTION_TRACE_ROW_FIELD_OF_VIEW, C_BLUE, m_dFieldOfViewRatio);
+    traceShowParam ("Num Detectors:", "%d", PROJECTION_TRACE_ROW_NDET, C_BLUE, proj.nDet());
+    traceShowParam ("Num Views:", "%d", PROJECTION_TRACE_ROW_NVIEW, C_BLUE, proj.nView());
+    traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
     
-    pSGP->setColor (C_RED);
-    phm.draw (*pSGP);
-    
-    pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
+    m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
   }
 #endif
 
 #ifdef HAVE_SGP
-    if (pSGP && m_trace >= TRACE_PHM) {
-      pSGP->setRasterOp (RO_XOR);
-      pSGP->setColor (C_RED);
-      pSGP->moveAbs (xd1, yd1);
-      pSGP->lineAbs (xd2, yd2);
-      pSGP->moveAbs (xs1, ys1);
-      pSGP->lineAbs (xs2, ys2);
-      pSGP->setRasterOp (RO_SET);
+    if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
+      m_pSGP->setColor (C_BLACK);
+      m_pSGP->setPenWidth (2);
+      if (m_idGeometry == GEOMETRY_PARALLEL) {
+	m_pSGP->moveAbs (xs1, ys1);
+	m_pSGP->lineAbs (xs2, ys2);
+	m_pSGP->moveAbs (xd1, yd1);
+	m_pSGP->lineAbs (xd2, yd2);
+      } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {	
+	m_pSGP->setPenWidth (4);
+	m_pSGP->moveAbs (xs1, ys1);
+	m_pSGP->lineAbs (xs2, ys2);
+	m_pSGP->setPenWidth (2);
+	m_pSGP->moveAbs (xd1, yd1);
+	m_pSGP->lineAbs (xd2, yd2);
+      } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+	m_pSGP->setPenWidth (4);
+	m_pSGP->moveAbs (xs1, ys1);
+	m_pSGP->lineAbs (xs2, ys2);
+	m_pSGP->setPenWidth (2);
+	m_pSGP->moveAbs (0., 0.);
+	m_pSGP->drawArc (m_dFocalLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
+      }
+      m_pSGP->setPenWidth (1);
     }
-    if (m_trace >= TRACE_TEXT)
-      traceShowParam (pSGP, "Current View :", "%5d", RAYSUM_TRACE_ROW_CURR_VIEW, C_LTMAGENTA, iview);
+    if (m_trace > Trace::TRACE_CONSOLE)
+      traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / static_cast<double>(m_nView) * 100.);
 #endif
-	    
-    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, pSGP);
+    if (m_trace == Trace::TRACE_CONSOLE)
+	cout << "Current View: " << iView+iStartView << endl;
+
+    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
     detArray.setViewAngle (viewAngle);
       
 #ifdef HAVE_SGP
-    if (pSGP && m_trace >= TRACE_PHM) {
-      //	rs_plot (detArray, xd1, yd1, xcent, ycent, theta);
-      pSGP->setColor (C_RED);
-      pSGP->setRasterOp (RO_XOR);
-      pSGP->moveAbs (xd1, yd1);
-      pSGP->lineAbs (xd2, yd2);
-      pSGP->moveAbs (xs1, ys1);
-      pSGP->lineAbs (xs2, ys2);
-      pSGP->setRasterOp (RO_SET);
+    if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
+      //	rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta);
     }
 #endif
-    xform_mtx2 (rotmtx_incr, xd1, yd1);	 // rotate detector endpoints 
-    xform_mtx2 (rotmtx_incr, xd2, yd2);
-    xform_mtx2 (rotmtx_incr, xs1, ys1);
-    xform_mtx2 (rotmtx_incr, xs2, ys2);
-  } /* for each iview */
+    xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
+    xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
+    if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
+	xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints 
+	xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
+    }
+  } /* for each iView */
 }
 
 
@@ -367,23 +409,33 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
  */
 
 void 
-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, SGP* pSGP)
+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)
 {
-  double ddx = (xd2 - xd1) / detArray.nDet();  // change in coords between detectors
-  double ddy = (yd2 - yd1) / detArray.nDet();
-  double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords between source
-  double sdy = (ys2 - ys1) / detArray.nDet();
-
-  double ddx2 = ddx / m_nSample;	// Incr. between rays with detector cell
-  double ddy2 = ddy / m_nSample;	// Doesn't include detector endpoints 
-  double ddx2_ofs = ddx2 / 2;           // offset of 1st ray from start of detector cell
-  double ddy2_ofs = ddy2 / 2;
-  
-  double xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
-  double yd_maj = yd1 + ddy2_ofs;
+
+  double sdx = (xs2 - xs1) / detArray.nDet();  // change in coords 
+  double sdy = (ys2 - ys1) / detArray.nDet();  // between source
   double xs_maj = xs1 + (sdx / 2);	// put ray source in center of cell 
   double ys_maj = ys1 + (sdy / 2);
 
+  double ddx=0, ddy=0, ddx2=0, ddy2=0, ddx2_ofs=0, ddy2_ofs=0, xd_maj=0, yd_maj=0;
+  double dAngleInc=0, dAngleSampleInc=0, dAngleSampleOffset=0, dAngleMajor=0;
+  if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+      dAngleInc = m_dAngularDetIncrement;
+      dAngleSampleInc = dAngleInc / m_nSample;
+      dAngleSampleOffset = dAngleSampleInc / 2;
+      dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
+  } else {
+      ddx = (xd2 - xd1) / detArray.nDet();  // change in coords 
+      ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
+      ddx2 = ddx / m_nSample;	// Incr. between rays with detector cell
+      ddy2 = ddy / m_nSample;  // Doesn't include detector endpoints 
+      ddx2_ofs = ddx2 / 2;    // offset of 1st ray from start of detector cell
+      ddy2_ofs = ddy2 / 2;
+      
+      xd_maj = xd1 + ddx2_ofs;       // Incr. between detector cells
+      yd_maj = yd1 + ddy2_ofs;
+  }
+
   DetectorValue* detval = detArray.detValues();
 
   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
@@ -393,80 +445,124 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
       else
 	detval[d] = 0;
   } else {
-    for (int d = 0; d < detArray.nDet(); d++) {
-      double xd = xd_maj;
-      double yd = yd_maj;
+      for (int d = 0; d < detArray.nDet(); d++) {
       double xs = xs_maj;
       double ys = ys_maj;
+      double xd=0, yd=0, dAngle=0;
+      if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+	  dAngle = dAngleMajor;
+      } else {
+	  xd = xd_maj;
+	  yd = yd_maj;
+      }
       double sum = 0.0;
       for (unsigned int i = 0; i < m_nSample; i++) {
+	if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
+	    xd = m_dFocalLength * cos (dAngle);
+	    yd = m_dFocalLength * sin (dAngle);
+	}
+
 #ifdef HAVE_SGP
-	if (pSGP && m_trace >= TRACE_RAYS) {
-	  pSGP->setColor (C_LTBLUE);
-	  pSGP->setRasterOp (RO_XOR);
-	  pSGP->moveAbs (xs, ys);
-	  pSGP->lineAbs (xd, yd);
-	  pSGP->setRasterOp (RO_SET);
+	if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
+	  m_pSGP->setColor (C_YELLOW);
+	  m_pSGP->setRasterOp (RO_AND);
+	  m_pSGP->moveAbs (xs, ys);
+	  m_pSGP->lineAbs (xd, yd);
 	}
 #endif
-	sum += projectSingleLine (phm, xd, yd, xs, ys, pSGP);
+
+	sum += projectSingleLine (phm, xd, yd, xs, ys);
 	      
 #ifdef HAVE_SGP
-	if (m_trace >= TRACE_RAYS) {
-	  traceShowParam (pSGP, "Attenuation  :", "%5.2f", RAYSUM_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
-	}
-
-	if (pSGP && m_trace >= TRACE_RAYS) {
-	  pSGP->setColor (C_LTBLUE);
-	  pSGP->setRasterOp (RO_XOR);
-	  pSGP->moveAbs (xs, ys);
-	  pSGP->lineAbs (xd, yd);
-	  pSGP->setRasterOp (RO_SET);
-	}
+	//	if (m_trace >= Trace::TRACE_CLIPPING) {
+	//	  traceShowParam ("Attenuation:", "%s", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, "        ");
+	//	  traceShowParam ("Attenuation:", "%.3f", PROJECTION_TRACE_ROW_ATTEN, C_LTMAGENTA, sum);
+	//	}
 #endif
-	xd += ddx2;
-	yd += ddy2;
-      }
+	if (m_idGeometry == GEOMETRY_EQUIANGULAR)
+	    dAngle += dAngleSampleInc;
+	else {
+	    xd += ddx2;
+	    yd += ddy2;
+	}
+      } // for each sample in detector
 
       detval[d] = sum / m_nSample;
-      xd_maj += ddx;
-      yd_maj += ddy;
       xs_maj += sdx;
       ys_maj += sdy;
+      if (m_idGeometry == GEOMETRY_EQUIANGULAR)
+	  dAngleMajor += dAngleInc;
+	else {
+	    xd_maj += ddx;
+	    yd_maj += ddy;
+	}
     } /* for each detector */
   } /* if not unit pulse */
 }
 
 
 void 
-Scanner::traceShowParam (SGP* pSGP, const char *label, const char *fmt, int row, int color, ...)
+Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int color, ...)
 {  
-  char s[256];
   va_list arg;
+  va_start(arg, color);
+#ifdef HAVE_SGP
+  traceShowParamRasterOp (RO_COPY, szLabel, fmt, row, color, arg);
+#else
+  traceShowParamRasterOp (0, szLabel, fmt, row, color, arg);
+#endif  
+  va_end(arg);
+}
 
+void 
+Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int color, ...)
+{  
+  va_list arg;
   va_start(arg, color);
-  snprintf (s, sizeof(s), label, "%s");
-  string strOut(s);
-  vsnprintf (s, sizeof(s), fmt, arg);
-  strOut += s;
+#ifdef HAVE_SGP
+  traceShowParamRasterOp (RO_XOR, szLabel, fmt, row, color, arg);
+#else
+  traceShowParamRasterOp (0, szLabel, fmt, row, color, arg);
+#endif
+  va_end(arg);
+}
+
+void 
+Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char *fmt, int row, int color, va_list args)
+{  
+  char szValue[256];
+
+#ifdef MSVC
+  vsprintf (szValue, fmt, args);
+#else
+  vsnprintf (szValue, sizeof(szValue), fmt, args);
+#endif
 
   //  cio_set_cpos (raysum_trace_menu_column, row);
   //  cio_set_text_clr (color - 8, 0);
   //  cio_set_text_clr (color, 0);
 
-  if (pSGP) {
-      pSGP->moveAbs (0., row * 0.04);
-      pSGP->setTextColor (color, -1);
-      pSGP->drawText (strOut);
-  } else {
-    cio_put_str (strOut.c_str());
-    cio_put_str ("\n");
-  }
-
-  va_end(arg);
+#ifdef HAVE_SGP
+  if (m_pSGP) {
+    m_pSGP->setRasterOp (iRasterOp);
+    double dYPos = m_dYMaxWin - (row * m_dTextHeight);
+    m_pSGP->moveAbs (m_dXMinWin, dYPos);
+    m_pSGP->setTextColor (color, -1);
+    m_pSGP->drawText (szLabel);
+    double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 4;
+    m_pSGP->moveAbs (m_dXMinWin + dValueOffset, dYPos);
+    m_pSGP->drawText (szValue);
+  } else 
+#endif
+    {
+      cio_put_str (szLabel);
+      cio_put_str (szValue);
+      cio_put_str ("\n");
+    }
 }
 
 
+
 /* NAME
  *    projectSingleLine			INTERNAL: Calculates raysum along a line for a Phantom
  *
@@ -478,12 +574,12 @@ Scanner::traceShowParam (SGP* pSGP, const char *label, const char *fmt, int row,
  */
 
 double 
-Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2, SGP* pSGP)
+Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1, const double x2, const double y2)
 {
   // check ray against each pelem in Phantom 
   double rsum = 0.0;
   for (PElemConstIterator i = phm.listPElem().begin(); i != phm.listPElem().end(); i++)
-    rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2, pSGP);
+    rsum += projectLineAgainstPElem (**i, x1, y1, x2, y2);
 
   return (rsum);
 }
@@ -500,23 +596,23 @@ Scanner::projectSingleLine (const Phantom& phm, const double x1, const double y1
  */
 
 double 
-Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2, SGP* pSGP)
+Scanner::projectLineAgainstPElem (const PhantomElement& pelem, double x1, double y1, double x2, double y2)
 {
   if (! pelem.clipLineWorldCoords (x1, y1, x2, y2)) {
-    if (m_trace == TRACE_CLIPPING)
+    if (m_trace == Trace::TRACE_CLIPPING)
       cio_tone (1000., 0.05);
     return (0.0);
   }
 
 #ifdef HAVE_SGP
-  if (pSGP && m_trace == TRACE_CLIPPING) {
-    pSGP->setRasterOp (RO_XOR);
-    pSGP->moveAbs (x1, y1);
-    pSGP->lineAbs (x2, y2);
+  if (m_pSGP && m_trace == Trace::TRACE_CLIPPING) {
+    m_pSGP->setRasterOp (RO_XOR);
+    m_pSGP->moveAbs (x1, y1);
+    m_pSGP->lineAbs (x2, y2);
     cio_tone (8000., 0.05);
-    pSGP->moveAbs (x1, y1);
-    pSGP->lineAbs (x2, y2);
-    pSGP->setRasterOp (RO_SET);
+    m_pSGP->moveAbs (x1, y1);
+    m_pSGP->lineAbs (x2, y2);
+    m_pSGP->setRasterOp (RO_SET);
   }
 #endif