From: Kevin M. Rosenberg <kevin@rosenberg.net>
Date: Wed, 21 Mar 2018 04:31:14 +0000 (-0600)
Subject: Use OpenMP for scanner
X-Git-Tag: v6.0.2~14
X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=commitdiff_plain;h=9ac3b88884957e2c07bf365c2503c6c1fbeaa60e

Use OpenMP for scanner
---

diff --git a/include/ctsupport.h b/include/ctsupport.h
index b509eee..6d1c1dc 100644
--- a/include/ctsupport.h
+++ b/include/ctsupport.h
@@ -257,6 +257,7 @@ void scale_mtx2 (GRFMTX_2D m, const double sx, const double sy);
 void rot_mtx2 (GRFMTX_2D m, const double theta);
 void mult_mtx2 (const GRFMTX_2D m1, const GRFMTX_2D m2, GRFMTX_2D result);
 void xform_mtx2 (const GRFMTX_2D m, double& x, double& y);
+void copy_mtx2 (GRFMTX_2D to, const GRFMTX_2D from);
 void rotate2d (double x[], double y[], int pts, double angle);
 void xlat2d (double x[], double y[], int pts, double xoffset, double yoffset);
 void scale2d (double x[], double y[], int pts, double xfact, double yfact);
diff --git a/libctsim/scanner.cpp b/libctsim/scanner.cpp
index 8974506..1567bd8 100644
--- 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 */
 }
 
diff --git a/libctsupport/xform.cpp b/libctsupport/xform.cpp
index cf35c4c..a8f0956 100644
--- a/libctsupport/xform.cpp
+++ b/libctsupport/xform.cpp
@@ -119,6 +119,13 @@ rot_mtx2 (GRFMTX_2D m, const double theta)
   m[1][0] = -s;  m[1][1] = c;
 }
 
+void
+copy_mtx2 (GRFMTX_2D to, const GRFMTX_2D from) {
+  for (int r = 0; r < 3; r++)
+    for (int c = 0; c < 3; c++)
+      to[r][c] = from[r][c];
+}
+
 void
 mult_mtx2 (const GRFMTX_2D m1, const GRFMTX_2D m2, GRFMTX_2D result)
 {
@@ -131,10 +138,7 @@ mult_mtx2 (const GRFMTX_2D m1, const GRFMTX_2D m2, GRFMTX_2D result)
             temp[row][col] += m1[row][calc] * m2[calc][col];
     }
   }
-
-  for (int r = 0; r < 3; r++)
-    for (int col = 0; col < 3; col++)
-      result[r][col] = temp[r][col];
+  copy_mtx2 (result, temp);
 }
 
 void