projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
r186: *** empty log message ***
[ctsim.git] / libctsim / scanner.cpp
diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index c45d4fd3234c1aa71561728d3bfb952b0840eeba..ef72c5b4cc99a7672a9d3ba2a25eac03a8d26edf 100644 (file)
--- 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.12 2000/08/27 20:32:55 kevin Exp $
+**  $Id: scanner.cpp,v 1.13 2000/08/31 08:38:58 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*);
@@ -101,8 +101,6 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     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 +112,26 @@ 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;
     
     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.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.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.yd2 = m_dYCenter - dHalfDetLen;
     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 +140,26 @@ 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;
     
-    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 +167,20 @@ 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;
+    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 +267,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;
@@ -303,12 +308,12 @@ 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 ();
 
@@ -348,15 +353,15 @@ 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.);
+      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);
+    projectSingleView (phm, detArray, xd1, yd1, xd2, yd2, xs1, ys1, xs2, ys2, viewAngle + 3 * HALFPI);
     detArray.setViewAngle (viewAngle);
       
 #ifdef HAVE_SGP
@@ -366,9 +371,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 +414,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,10 +459,9 @@ 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
 
@@ -470,14 +472,6 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
          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);
-       //      }
 #endif
        if (m_idGeometry == GEOMETRY_EQUIANGULAR)
            dAngle += dAngleSampleInc;
Computed Tomography simulator and reconstructor