r5231: Auto commit for Debian build
[ctsim.git] / libctsim / scanner.cpp
index ff3cd9af366833c74a3e7e3e52641aeeaaa43a2b..aeb51975e82fd7e8436a542aae971d1203617245 100644 (file)
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.34 2001/03/10 23:14:16 kevin Exp $
+**  $Id: scanner.cpp,v 1.44 2003/07/04 21:39:40 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
@@ -32,19 +32,22 @@ const int Scanner::GEOMETRY_INVALID = -1;
 const int Scanner::GEOMETRY_PARALLEL = 0;
 const int Scanner::GEOMETRY_EQUIANGULAR = 1;
 const int Scanner::GEOMETRY_EQUILINEAR = 2;
+const int Scanner::GEOMETRY_LINOGRAM = 3;
 
-const char* Scanner::s_aszGeometryName[] = 
+const char* const Scanner::s_aszGeometryName[] = 
 {
-  {"parallel"},
-  {"equiangular"},
-  {"equilinear"},
+  "parallel",
+  "equiangular",
+  "equilinear",
+  "linogram",
 };
 
-const char* Scanner::s_aszGeometryTitle[] = 
+const char* const Scanner::s_aszGeometryTitle[] = 
 {
-  {"Parallel"},
-  {"Equiangular"},
-  {"Equilinear"},
+  "Parallel",
+  "Equiangular",
+  "Equilinear",
+  "Linogram",
 };
 
 const int Scanner::s_iGeometryCount = sizeof(s_aszGeometryName) / sizeof(const char*);
@@ -83,8 +86,10 @@ DetectorArray::~DetectorArray (void)
 */
 
 Scanner::Scanner (const Phantom& phm, const char* const geometryName, 
-                  int nDet, int nView, int nSample, const double rot_anglen, 
-                  const double dFocalLengthRatio, const double dCenterDetectorRatio,
+                  int nDet, int nView, int offsetView, 
+                                 int nSample, const double rot_anglen, 
+                  const double dFocalLengthRatio, 
+                                 const double dCenterDetectorRatio,
                   const double dViewRatio, const double dScanRatio)
 {
   m_fail = false;
@@ -106,6 +111,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
   
   m_nDet     = nDet;
   m_nView    = nView;
+  m_iOffsetView = offsetView;
   m_nSample  = nSample;
   m_dFocalLengthRatio = dFocalLengthRatio;
   m_dCenterDetectorRatio = dCenterDetectorRatio;
@@ -125,11 +131,12 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
   if (m_idGeometry == GEOMETRY_PARALLEL) {
     m_dFanBeamAngle = 0;
     m_detLen   = m_dScanDiameter;
+    m_detStart = -m_detLen / 2;
     m_detInc  = m_detLen / m_nDet;
     double dDetectorArrayEndOffset = 0;
     // For even number of detectors, make detInc slightly larger so that center lies
     // at nDet/2. Also, extend detector array by one detInc so that all of the phantom is scanned
-    if (m_nDet % 2 == 0) { // Adjust for Even number of detectors
+    if (isEven (m_nDet)) { // Adjust for Even number of detectors
       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
       dDetectorArrayEndOffset = m_detInc;
     }
@@ -143,7 +150,8 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
-    m_initPos.angle = 0.0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
+    m_detLen += dDetectorArrayEndOffset;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
   if (m_dScanDiameter / 2 >= m_dFocalLength) {
       m_fail = true;
@@ -155,15 +163,16 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     const double dHalfDetLen = m_dSourceDetectorLength * tan (dAngle);
     
     m_detLen = dHalfDetLen * 2;
+    m_detStart = -dHalfDetLen;
     m_detInc  = m_detLen / m_nDet;
     double dDetectorArrayEndOffset = 0;
-    if (m_nDet % 2 == 0) { // Adjust for Even number of detectors
+    if (isEven (m_nDet)) { // Adjust for Even number of detectors
       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
       dDetectorArrayEndOffset = m_detInc;
+      m_detLen += dDetectorArrayEndOffset;
     }
   
     m_dFanBeamAngle = dAngle * 2;
-    m_initPos.angle = 0.0;
     m_initPos.xs1 = m_dXCenter;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
     m_initPos.xs2 = m_dXCenter;
@@ -172,7 +181,7 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.yd1 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
     m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
-    m_initPos.angle = 0.0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
     if (m_dScanDiameter / 2 > m_dFocalLength) {
       m_fail = true;
@@ -182,9 +191,10 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
 
     m_detLen = 2 * dAngle;
+    m_detStart = -dAngle;
     m_detInc = m_detLen / m_nDet;
     double dDetectorArrayEndOffset = 0;
-    if (m_nDet % 2 == 0) { // Adjust for Even number of detectors
+    if (isEven (m_nDet)) { // Adjust for Even number of detectors
       m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
       dDetectorArrayEndOffset = m_detInc;
     }
@@ -197,11 +207,12 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
     
     m_dFanBeamAngle = dAngle * 2;
-    m_initPos.angle = 0;
+    m_initPos.angle = m_iOffsetView * m_rotInc;
     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_detLen += dDetectorArrayEndOffset;
   }
   
   // Calculate incrementatal rotation matrix 
@@ -271,23 +282,25 @@ Scanner::convertGeometryNameToID (const char* const geomName)
 void
 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int trace, SGP* pSGP)
 {
-  collectProjections (proj, phm, 0, proj.nView(), true, trace, pSGP);
+  collectProjections (proj, phm, m_startView, proj.nView(), m_iOffsetView, true, trace, pSGP);
 }
 
 void
-Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-                             bool bStoreAtViewPosition, const int trace, SGP* pSGP)
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, 
+                             const int iNumViews, const int iOffsetView,  bool bStoreAtViewPosition, 
+                             const int trace, SGP* pSGP)
 {
   int iStorageOffset = (bStoreAtViewPosition ? iStartView : 0);
-  collectProjections (proj, phm, iStartView, iNumViews, iStorageOffset, trace, pSGP);
+  collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP);
 }
 
 void
-Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
-                             int iStorageOffset, const int trace, SGP* pSGP)
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, 
+                             const int iNumViews, const int iOffsetView, int iStorageOffset, 
+                             const int trace, SGP* pSGP)
 {
   m_trace = trace;
-  double start_angle = iStartView * proj.rotInc();
+  double start_angle = (iStartView + iOffsetView) * proj.rotInc();
   
   // Calculate initial rotation matrix 
   GRFMTX_2D rotmtx_initial, temp;
@@ -342,19 +355,23 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       m_pSGP->setColor (C_GREEN);
       m_pSGP->drawCircle (m_dFocalLength);
       m_pSGP->setColor (C_BLUE);
+#if MSVC
       m_pSGP->setTextPointSize (9);
+#else
+      m_pSGP->setTextPointSize (14);
+#endif
       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 ("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);
       
-      m_pSGP->setMarker (SGP::MARK_BDIAMOND, C_LTGREEN);
+      m_pSGP->setMarker (SGP::MARKER_BDIAMOND);
     }
 #endif
     
@@ -463,10 +480,8 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
   
   if (phm.getComposition() == P_UNIT_PULSE) {  // put unit pulse in center of view
     for (int d = 0; d < detArray.nDet(); d++)
-      if (detArray.nDet() / 2 == d && (d % 2) == 1)
-        detval[d] = 1;
-      else
         detval[d] = 0;
+    detval[ detArray.nDet() / 2 ] = 1;
   } else {
     for (int d = 0; d < detArray.nDet(); d++) {
       double xs = xs_maj;