X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fscanner.cpp;h=2314595219d76c45801d1a8dbf8a57f91de10f05;hp=7081903f71b7984b17723c905a29c2b810826efc;hb=c953cbb6ffc2fd50e736230f4e6976a025983cff;hpb=c5e7140bd08b8c8f527713e8dc861bcb7ee5f633

diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index 7081903..2314595 100644
--- a/libctsim/scanner.cpp
+++ b/libctsim/scanner.cpp
@@ -9,7 +9,7 @@
 **  This is part of the CTSim program
 **  Copyright (c) 1983-2001 Kevin Rosenberg
 **
-**  $Id: scanner.cpp,v 1.30 2001/02/08 06:25:07 kevin Exp $
+**  $Id: scanner.cpp,v 1.33 2001/03/01 07:30:49 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
@@ -84,8 +84,8 @@ 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 dViewRatio,
-                  double dScanRatio)
+                  const double dFocalLengthRatio, const double dCenterDetectorRatio,
+                  const double dViewRatio, const double dScanRatio)
 {
   m_fail = false;
   m_idGeometry = convertGeometryNameToID (geometryName);
@@ -108,10 +108,14 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
   m_nView    = nView;
   m_nSample  = nSample;
   m_dFocalLengthRatio = dFocalLengthRatio;
+  m_dCenterDetectorRatio = dCenterDetectorRatio;
   m_dViewRatio = dViewRatio;
   m_dScanRatio = dScanRatio;
+
   m_dViewDiameter = phm.getDiameterBoundaryCircle() * m_dViewRatio;
-  m_dFocalLength = (m_dViewDiameter / 2) * dFocalLengthRatio;
+  m_dFocalLength = (m_dViewDiameter / 2) * m_dFocalLengthRatio;
+  m_dCenterDetectorLength = (m_dViewDiameter / 2) * m_dCenterDetectorRatio;
+  m_dSourceDetectorLength = m_dFocalLength + m_dCenterDetectorLength;
   m_dScanDiameter = m_dViewDiameter * m_dScanRatio;
   
   m_dXCenter = phm.xmin() + (phm.xmax() - phm.xmin()) / 2;
@@ -119,21 +123,26 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
   m_rotLen  = rot_anglen;
   m_rotInc  = m_rotLen / m_nView;
   if (m_idGeometry == GEOMETRY_PARALLEL) {
+    m_dFanBeamAngle = 0;
     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 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
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
+      dDetectorArrayEndOffset = m_detInc;
+    }
+
     double dHalfDetLen = m_detLen / 2;
     m_initPos.xs1 = m_dXCenter - dHalfDetLen;
     m_initPos.ys1 = m_dYCenter + m_dFocalLength;
-    m_initPos.xs2 = m_dXCenter + dHalfDetLen;
+    m_initPos.xs2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
     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.yd1 = m_dYCenter - m_dCenterDetectorLength;
+    m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
+    m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUILINEAR) {
   if (m_dScanDiameter / 2 >= m_dFocalLength) {
@@ -141,13 +150,17 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
       m_failMessage = "Invalid geometry: Focal length must be larger than scan length";
       return;
     }
+
     const double dAngle = asin ((m_dScanDiameter / 2) / m_dFocalLength);
-    const double dHalfDetLen = 2 * m_dFocalLength * tan (dAngle);
+    const double dHalfDetLen = m_dSourceDetectorLength * 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
+    double dDetectorArrayEndOffset = 0;
+    if (m_nDet % 2 == 0) { // Adjust for Even number of detectors
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
+      dDetectorArrayEndOffset = m_detInc;
+    }
   
     m_dFanBeamAngle = dAngle * 2;
     m_initPos.angle = 0.0;
@@ -156,9 +169,9 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
     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.yd1 = m_dYCenter - m_dCenterDetectorLength;
+    m_initPos.xd2 = m_dXCenter + dHalfDetLen + dDetectorArrayEndOffset;
+    m_initPos.yd2 = m_dYCenter - m_dCenterDetectorLength;
     m_initPos.angle = 0.0;
   } else if (m_idGeometry == GEOMETRY_EQUIANGULAR) {
     if (m_dScanDiameter / 2 > m_dFocalLength) {
@@ -170,10 +183,15 @@ Scanner::Scanner (const Phantom& phm, const char* const geometryName,
 
     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;
+    double dDetectorArrayEndOffset = 0;
+    if (m_nDet % 2 == 0) { // Adjust for Even number of detectors
+      m_detInc = m_detLen / (m_nDet - 1); // center detector = (nDet/2)
+      dDetectorArrayEndOffset = m_detInc;
+    }
+    double dA1 = acos ((m_dScanDiameter / 2) / m_dCenterDetectorLength);
+    double dAngularScale = 2 * (HALFPI + dAngle - dA1) / m_detLen;
+    m_dAngularDetLen = dAngularScale * (m_detLen + dDetectorArrayEndOffset);
+    m_dAngularDetIncrement = dAngularScale * m_detInc;
     m_initPos.dAngularDet = -m_dAngularDetLen / 2;
     
     m_dFanBeamAngle = dAngle * 2;
@@ -255,7 +273,16 @@ 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)
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
+                             bool bStoreAtViewPosition, const int trace, SGP* pSGP)
+{
+  int iStorageOffset = (bStoreAtViewPosition ? iStartView : 0);
+  collectProjections (proj, phm, iStartView, iNumViews, iStorageOffset, trace, pSGP);
+}
+
+void
+Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView, const int iNumViews, 
+                             int iStorageOffset, const int trace, SGP* pSGP)
 {
   m_trace = trace;
   double start_angle = iStartView * proj.rotInc();
@@ -288,16 +315,14 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   int iView;
   double viewAngle;
   for (iView = 0, viewAngle = start_angle;  iView < iNumViews; iView++, viewAngle += proj.rotInc()) {
-    int iStoragePosition = iView;
-    if (bStoreAtViewPosition)
-      iStoragePosition += iStartView;
-    
+    int iStoragePosition = iView + iStorageOffset;
+
     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
     
 #ifdef HAVE_SGP 
     if (pSGP && m_trace >= Trace::TRACE_PHANTOM) {
       m_pSGP = pSGP;
-      double dWindowSize = dmax (m_detLen, m_dFocalLength * 2) * SQRT2;
+      double dWindowSize = dmax (m_detLen, m_dSourceDetectorLength) * 2;
       double dHalfWindowSize = dWindowSize / 2;
       m_dXMinWin = m_dXCenter - dHalfWindowSize;
       m_dXMaxWin = m_dXCenter + dHalfWindowSize;
@@ -353,7 +378,7 @@ 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, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
+        m_pSGP->drawArc (m_dCenterDetectorLength, viewAngle + 3 * HALFPI - (m_dAngularDetLen/2), viewAngle + 3 * HALFPI + (m_dAngularDetLen/2));
       }
       m_pSGP->setPenWidth (1);
     }
@@ -454,8 +479,8 @@ Scanner::projectSingleView (const Phantom& phm, DetectorArray& detArray, const d
       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);
+          xd = m_dCenterDetectorLength * cos (dAngle);
+          yd = m_dCenterDetectorLength * sin (dAngle);
         }
         
 #ifdef HAVE_SGP