Use OpenMP for scanner

[ctsim.git] / libctsim / scanner.cpp

diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index 89745061240059b7aa3aaeece7dad440afcad907..1567bd8b0f2b32d20577c0b270b78961111ab1e3 100644 (file)
--- a/libctsim/scanner.cpp
+++ b/libctsim/scanner.cpp
@@ -292,6 +292,15 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   collectProjections (proj, phm, iStartView, iNumViews, iOffsetView, iStorageOffset, trace, pSGP);
 }
 
+static void mtx2_offset_rot (GRFMTX_2D m, double angle, double x, double y) {
+  GRFMTX_2D temp;
+  xlat_mtx2 (m, -x, -y);
+  rot_mtx2 (temp, angle);
+  mult_mtx2 (m, temp, m);
+  xlat_mtx2 (temp, x, y);
+  mult_mtx2 (m, temp, m);
+}
+
 void
 Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iStartView,
                              const int iNumViews, const int iOffsetView, int iStorageOffset,
@@ -300,36 +309,48 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
   m_trace = trace;
   double start_angle = (iStartView + iOffsetView) * proj.rotInc();
 
-  // Calculate initial rotation matrix
-  GRFMTX_2D rotmtx_initial, temp;
-  xlat_mtx2 (rotmtx_initial, -m_dXCenter, -m_dYCenter);
-  rot_mtx2 (temp, start_angle);
-  mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
-  xlat_mtx2 (temp, m_dXCenter, m_dYCenter);
-  mult_mtx2 (rotmtx_initial, temp, rotmtx_initial);
-
+#ifndef HAVE_OPENMP
+  // Without OpenMP, precalculate source and detector at view 0, then rotate incrementally each time
+  GRFMTX_2D rotmtx_initial;
+  mtx2_offset_rot (rotmtx_initial, start_angle, m_dXCenter, m_dYCenter);
   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;
-    yd2 = m_initPos.yd2;
+    xd1 = m_initPos.xd1; yd1 = m_initPos.yd1;
+    xd2 = m_initPos.xd2; yd2 = m_initPos.yd2;
     xform_mtx2 (rotmtx_initial, xd1, yd1);      // rotate detector endpoints
     xform_mtx2 (rotmtx_initial, xd2, yd2);      // to initial view_angle
   }
 
-  double xs1 = m_initPos.xs1;
-  double ys1 = m_initPos.ys1;
-  double xs2 = m_initPos.xs2;
-  double ys2 = m_initPos.ys2;
+  double xs1 = m_initPos.xs1, ys1 = m_initPos.ys1;
+  double xs2 = m_initPos.xs2, ys2 = m_initPos.ys2;
   xform_mtx2 (rotmtx_initial, xs1, ys1);      // rotate source endpoints to
   xform_mtx2 (rotmtx_initial, xs2, ys2);      // initial view angle
-
-  int iView;
-  double viewAngle;
-  for (iView = 0, viewAngle = start_angle;  iView < iNumViews; iView++, viewAngle += proj.rotInc()) {
+#else
+  #pragma omp parallel for
+#endif
+  
+  for (int iView = 0;  iView < iNumViews; iView++) {
+    double viewAngle = start_angle + (iView * proj.rotInc());
+    
+#ifdef HAVE_OPENMP
+    // With OpenMP, need to calculate source and detector positions at each view
+    GRFMTX_2D rotmtx;
+    mtx2_offset_rot (rotmtx, viewAngle, m_dXCenter, m_dYCenter);
+    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; yd2 = m_initPos.yd2;
+      xform_mtx2 (rotmtx, xd1, yd1);      // rotate detector endpoints
+      xform_mtx2 (rotmtx, xd2, yd2);      // to initial view_angle
+    }
+    
+    double xs1 = m_initPos.xs1, ys1 = m_initPos.ys1;
+    double xs2 = m_initPos.xs2, ys2 = m_initPos.ys2;
+    xform_mtx2 (rotmtx, xs1, ys1);      // rotate source endpoints to
+    xform_mtx2 (rotmtx, xs2, ys2);      // initial view angle
+#endif
+    
     int iStoragePosition = iView + iStorageOffset;
-
     DetectorArray& detArray = proj.getDetectorArray( iStoragePosition );
 
 #ifdef HAVE_SGP
@@ -366,11 +387,8 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       traceShowParam ("Samples / Ray:", "%d", PROJECTION_TRACE_ROW_SAMPLES, C_BLUE, m_nSample);
 
       m_pSGP->setMarker (SGP::MARKER_BDIAMOND);
-    }
-#endif
 
-#ifdef HAVE_SGP
-    if (m_pSGP && m_trace >= Trace::TRACE_PHANTOM) {
+      
       m_pSGP->setColor (C_BLACK);
       m_pSGP->setPenWidth (2);
       if (m_idGeometry == GEOMETRY_PARALLEL) {
@@ -398,6 +416,7 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
     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
+
     if (m_trace == Trace::TRACE_CONSOLE)
       std::cout << "Current View: " << iView+iStartView << std::endl;
 
@@ -409,12 +428,16 @@ Scanner::collectProjections (Projections& proj, const Phantom& phm, const int iS
       //        rs_plot (detArray, xd1, yd1, dXCenter, dYCenter, theta);
     }
 #endif
+
+#ifndef HAVE_OPENMP
+    // Without OpenMP, incrementally rotate source and detectors 
     xform_mtx2 (m_rotmtxIncrement, xs1, ys1);
     xform_mtx2 (m_rotmtxIncrement, xs2, ys2);
     if (m_idGeometry != GEOMETRY_EQUIANGULAR) {
       xform_mtx2 (m_rotmtxIncrement, xd1, yd1);  // rotate detector endpoints
       xform_mtx2 (m_rotmtxIncrement, xd2, yd2);
     }
+#endif    
   } /* for each iView */
 }