Add OpenMP parallelism

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index e4a240824393f9186c5a7c259cf48ef996d503f3..695582aa86d493c6c8ba85615f9c4989bf3acc7f 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -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 (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
     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 (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
     for (int i = 0; i < nDet - 1; i++)
       deltaFilteredProj[i] = (filteredProj[i+1] - filteredProj[i]) * dInvScale;
     deltaFilteredProj[nDet - 1] = 0;  // last detector

diff --git a/libctsim/filter.cpp b/libctsim/filter.cpp
index 397f8a204120926ed04516653b5e1b2ca1d0f654..6382b5609497d5bc5b37e9e08c4ab3a049e54561 100644 (file)
--- a/libctsim/filter.cpp
+++ b/libctsim/filter.cpp
@@ -198,10 +198,14 @@ SignalFilter::~SignalFilter (void)
 void
 SignalFilter::createFrequencyFilter (double* adFilter) const
 {
 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);
     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 {
     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
       if (haveAnalyticSpatial(m_idFilter))
         m_adFilter[i] = spatialResponseAnalytic (x);
       else