r261: Use explicit std:: namespace
[ctsim.git] / libctsim / scanner.cpp
index c45d4fd3234c1aa71561728d3bfb952b0840eeba..59bcf943fa1eefd17fdbd1a74b8dae5c5bdf4330 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (C) 1983-2000 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.12 2000/08/27 20:32:55 kevin Exp $
+**  $Id: scanner.cpp,v 1.20 2000/12/16 06:12:47 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
 
 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,14 +95,13 @@ 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;
@@ -114,23 +113,28 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
 
   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_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter - dHalfDetLen;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys2 = m_dYCenter + dHalfDetLen;
-    m_initPos.xd1 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd1 = m_dYCenter - dHalfDetLen;
-    m_initPos.xd2 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd2 = m_dYCenter + 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) {
+#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.) {
@@ -139,23 +143,28 @@ 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_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys2 = m_dYCenter;
-    m_initPos.xd1 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd1 = m_dYCenter - dHalfDetLen;
-    m_initPos.xd2 = m_dXCenter + m_dFocalLength;
-    m_initPos.yd2 = m_dYCenter + 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) {
+#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.) {
@@ -163,17 +172,21 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
       m_failMessage = "Focal Point inside of phantom";
       return;
     }
-    double dAngle = 2 * atan( dHalfSquare / dFocalPastPhm );
+    double dAngle =  atan ( dHalfSquare / dFocalPastPhm ); 
+#endif
     m_detLen = 2 * dAngle;
     m_detInc = m_detLen / m_nDet;
-    m_rotLen = rot_anglen;
-    m_rotInc = m_rotLen / m_nView;
-
-    m_initPos.xs1 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys1 = m_dYCenter;
-    m_initPos.xs2 = m_dXCenter - m_dFocalLength;
-    m_initPos.ys2 = m_dYCenter;
-    m_initPos.angle = -dAngle;
+    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 
@@ -260,10 +273,8 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
   mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
 
-  double xd1=0, yd1=0, xd2=0, yd2=0, dDetAngle=0;
-  if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-      dDetAngle = m_initPos.angle + start_angle;
-  else {
+  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;
@@ -292,7 +303,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
       m_pSGP = pSGP;
       m_pSGP->eraseWindow();
-      double dWindowSize = max(m_detLen, m_dFocalLength * 2) * SQRT2;
+      double dWindowSize = dmax (m_detLen, m_dFocalLength * 2) * SQRT2;
       double dHalfWindowSize = dWindowSize / 2;
       m_dXMinWin = m_dXCenter - dHalfWindowSize;
       m_dXMaxWin = m_dXCenter + dHalfWindowSize;
@@ -303,17 +314,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
     m_pSGP->eraseWindow ();
     m_pSGP->setWindow (m_dXMinWin, m_dYMinWin, m_dXMaxWin, m_dYMaxWin);
     m_pSGP->setRasterOp (RO_COPY);
-    m_pSGP->setColor (C_BLUE);
+    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_RED);
+    m_pSGP->setColor (C_BLUE);
     phm.draw (*m_pSGP);
     m_dTextHeight = m_pSGP->getCharHeight ();
 
-    traceShowParam ("Projection Collector", "%s", PROJECTION_TRACE_ROW_TITLE, C_BLACK, " ");
-    traceShowParam ("________________", "%s", PROJECTION_TRACE_ROW_TITLE2, C_LTGRAY, " ");
     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);
@@ -348,15 +357,17 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
        m_pSGP->lineAbs (xs2, ys2);
        m_pSGP->setPenWidth (2);
        m_pSGP->moveAbs (0., 0.);
-       m_pSGP->drawArc (m_dFocalLength, start_angle + m_initPos.angle, start_angle - m_initPos.angle);
+       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::TRACE_CONSOLE)
-      traceShowParam ("Current View:", "%d (%.0f%%)", PROJECTION_TRACE_ROW_CURR_VIEW, C_RED, iView + iStartView, (iView + iStartView) / m_nView * 100.);
+    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, dDetAngle);
+    if (m_trace == Trace::TRACE_CONSOLE)
+               std::cout << "Current View: " << iView+iStartView << std::endl;
+
+    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
     detArray.setViewAngle (viewAngle);
       
 #ifdef HAVE_SGP
@@ -366,9 +377,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
 #endif
     xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
     xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
-    if (m_idGeometry == GEOMETRY_EQUIANGULAR)
-       dDetAngle += m_detInc;
-    else {
+    if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
        xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints 
        xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
     }
@@ -411,10 +420,10 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
   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_detInc;
+      dAngleInc = m_dAngularDetIncrement;
       dAngleSampleInc = dAngleInc / m_nSample;
       dAngleSampleOffset = dAngleSampleInc / 2;
-      dAngleMajor = dDetAngle + dAngleSampleOffset;
+      dAngleMajor = dDetAngle - (m_dAngularDetLen/2) + dAngleSampleOffset;
   } else {
       ddx = (xd2 - xd1) / detArray.nDet();  // change in coords 
       ddy = (yd2 - yd1) / detArray.nDet();  // between detectors
@@ -456,27 +465,18 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
 #ifdef HAVE_SGP
        if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
          m_pSGP->setColor (C_YELLOW);
-         m_pSGP->setRasterOp (RO_OR_REVERSE);
+         m_pSGP->setRasterOp (RO_AND);
          m_pSGP->moveAbs (xs, ys);
          m_pSGP->lineAbs (xd, yd);
-         m_pSGP->setRasterOp (RO_SET);
        }
 #endif
 
        sum += projectSingleLine (phm, xd, yd, xs, ys);
              
 #ifdef HAVE_SGP
-       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);
-       }
-
-       //      if (m_pSGP && m_trace >= Trace::TRACE_PROJECTIONS) {
-       //        m_pSGP->setColor (C_YELLOW);
-       //        m_pSGP->setRasterOp (RO_XOR);
-       //        m_pSGP->moveAbs (xs, ys);
-       //        m_pSGP->lineAbs (xd, yd);
-       //        m_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
        if (m_idGeometry == GEOMETRY_EQUIANGULAR)
@@ -506,7 +506,11 @@ Scanner::traceShowParam (const char *szLabel, const char *fmt, int row, int colo
 {  
   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);
 }
 
@@ -515,7 +519,11 @@ Scanner::traceShowParamXOR (const char *szLabel, const char *fmt, int row, int c
 {  
   va_list arg;
   va_start(arg, color);
+#ifdef HAVE_SGP
   traceShowParamRasterOp (RO_XOR, szLabel, fmt, row, color, arg);
+#else
+  traceShowParamRasterOp (0, szLabel, fmt, row, color, arg);
+#endif
   va_end(arg);
 }
 
@@ -524,26 +532,33 @@ Scanner::traceShowParamRasterOp (int iRasterOp, const char *szLabel, const char
 {  
   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);
 
+#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) / 5;
+    double dValueOffset = (m_dXMaxWin - m_dXMinWin) / 4;
     m_pSGP->moveAbs (m_dXMinWin + dValueOffset, dYPos);
     m_pSGP->drawText (szValue);
-  } else {
-    cio_put_str (szLabel);
-    cio_put_str (szValue);
-    cio_put_str ("\n");
-  }
+  } else 
+#endif
+    {
+      cio_put_str (szLabel);
+      cio_put_str (szValue);
+      cio_put_str ("\n");
+    }
 }