r499: no message
[ctsim.git] / libctsim / scanner.cpp
index 6e2e57ad48bf962af80a9219c9638947f21bc713..7081903f71b7984b17723c905a29c2b810826efc 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.29 2001/02/04 21:28:19 kevin Exp $
+**  $Id: scanner.cpp,v 1.30 2001/02/08 06:25:07 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
@@ -82,10 +82,11 @@ DetectorArray::~DetectorArray (void)
 *   int nSample                Number of rays per detector
 */
 
-Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet, int nView, int nSample, const double rot_anglen, const double dFocalLengthRatio, const double dFieldOfViewRatio)
+Scanner::Scanner (const Phantom& phm, const char* const geometryName, 
+                  int nDet, int nView, int nSample, const double rot_anglen, 
+                  const double dFocalLengthRatio, const double dViewRatio,
+                  double dScanRatio)
 {
-  m_phmLen = phm.maxAxisLength();      // maximal length along an axis
-  
   m_fail = false;
   m_idGeometry = convertGeometryNameToID (geometryName);
   if (m_idGeometry == GEOMETRY_INVALID) {
@@ -107,20 +108,23 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int 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_dViewRatio = dViewRatio;
+  m_dScanRatio = dScanRatio;
+  m_dViewDiameter = phm.getDiameterBoundaryCircle() * m_dViewRatio;
+  m_dFocalLength = (m_dViewDiameter / 2) * dFocalLengthRatio;
+  m_dScanDiameter = m_dViewDiameter * m_dScanRatio;
   
   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_detLen   = m_dScanDiameter;
     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_dFanBeamAngle = 0;
     double dHalfDetLen = m_detLen / 2;
     m_initPos.xs1 = m_dXCenter - dHalfDetLen;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
@@ -132,34 +136,20 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_initPos.yd2 = m_dYCenter - m_dFocalLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
-#if 0
-    double dAngle1 = atan ((m_dFieldOfView / 2) / m_dFocalLength);
-    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 dAngle2 = atan( dHalfSquare / dFocalPastPhm );
-    double dAngle = maxValue<double> (dAngle1, dAngle2);
-    //double dAngle = (m_dFieldOfView / 2) / cos (asin (m_dFieldOfView / 2 / m_dFocalLength));
-#else
-    if (m_dFieldOfView/2 >= m_dFocalLength) {
+  if (m_dScanDiameter / 2 >= m_dFocalLength) {
       m_fail = true;
-      m_failMessage = "Invalid geometry: Focal length must be larger than field of view";
+      m_failMessage = "Invalid geometry: Focal length must be larger than scan length";
       return;
     }
-    double dAngle = asin ((m_dFieldOfView/2) / m_dFocalLength);
-#endif    
-
-    double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
+    const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
+    const double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
     
     m_detLen = dHalfDetLen * 2;
     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_dFanBeamAngle = dAngle * 2;
     m_initPos.angle = 0.0;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
@@ -171,12 +161,12 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_initPos.yd2 = m_dYCenter - m_dFocalLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
-    if (m_dFieldOfView/2 > m_dFocalLength) {
+    if (m_dScanDiameter / 2 > m_dFocalLength) {
       m_fail = true;
-      m_failMessage = "Invalid geometry: Focal length must be larger than field of view";
+      m_failMessage = "Invalid geometry: Focal length must be larger than scan length";
       return;
     }
-    double dAngle = asin ((m_dFieldOfView/2) / m_dFocalLength);
+    const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
 
     m_detLen = 2 * dAngle;
     m_detInc = m_detLen / m_nDet;
@@ -186,6 +176,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName, int nDet,
     m_dAngularDetLen = m_detLen * 2;
     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
     
+    m_dFanBeamAngle = dAngle * 2;
     m_initPos.angle = 0;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;;
@@ -312,14 +303,16 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       m_dXMaxWin = m_dXCenter + dHalfWindowSize;
       m_dYMinWin = m_dYCenter - dHalfWindowSize;
       m_dYMaxWin = m_dYCenter + dHalfWindowSize;
-      double dHalfPhmLen = m_phmLen /  2;
       
       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->drawCircle (m_dViewDiameter / 2);
+
       m_pSGP->moveAbs (0., 0.);
+      m_pSGP->setColor (C_GREEN);
       m_pSGP->drawCircle (m_dFocalLength);
       m_pSGP->setColor (C_BLUE);
       m_pSGP->setTextPointSize (9);
@@ -329,7 +322,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       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 ("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);