Add OpenMP parallelism
authorKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 22 Mar 2018 17:53:59 +0000 (11:53 -0600)
committerKevin M. Rosenberg <kevin@rosenberg.net>
Thu, 22 Mar 2018 17:53:59 +0000 (11:53 -0600)
libctsim/backprojectors.cpp
libctsim/filter.cpp

index e4a2408..695582a 100644 (file)
@@ -551,6 +551,9 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   if (interpType == Backprojector::INTERP_LINEAR) {
     // precalculate scaled difference for linear interpolation
     deltaFilteredProj = new double [nDet];
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
     for (int i = 0; i < nDet - 1; i++)
       deltaFilteredProj[i] = filteredProj[i+1] - filteredProj[i];
     deltaFilteredProj[nDet - 1] = 0;  // last detector
@@ -627,6 +630,9 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
   if (interpType == Backprojector::INTERP_LINEAR) {
     // precalculate scaled difference for linear interpolation
     deltaFilteredProj = new double [nDet];
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (int i = 0; i < nDet - 1; i++)
       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
     deltaFilteredProj[nDet - 1] = 0;  // last detector
index 397f8a2..6382b56 100644 (file)
@@ -198,10 +198,14 @@ SignalFilter::~SignalFilter (void)
 void
 SignalFilter::createFrequencyFilter (double* adFilter) const
 {
-  double x;
-  int i;
-  for (x = m_dFilterMin, i = 0; i < m_nFilterPoints; x += m_dFilterInc, i++)
+  double xstart = m_dFilterMin;
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int i = 0; i < m_nFilterPoints; i++) {
+    double x = xstart + (i * m_dFilterInc);
     adFilter[i] = frequencyResponse (x);
+  }
 }
 
 
@@ -218,8 +222,12 @@ SignalFilter::createSpatialFilter (double* adFilter) const
     for (int i = 1; i <= sidelen; i++ )
       m_adFilter [center + i] = m_adFilter [center - i] = c / (4 * (i * i) - 1);
   } else {
-    double x = m_dFilterMin;
-    for (int i = 0; i < m_nFilterPoints; i++, x += m_dFilterInc) {
+    double xstart = m_dFilterMin;
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
+    for (int i = 0; i < m_nFilterPoints; i++) {
+      double x = xstart + (i * m_dFilterInc);
       if (haveAnalyticSpatial(m_idFilter))
         m_adFilter[i] = spatialResponseAnalytic (x);
       else