Use OpenMP for signal filtering

[ctsim.git] / libctsim / procsignal.cpp

diff --git a/libctsim/procsignal.cpp b/libctsim/procsignal.cpp
index be88ecc31e34cee618746eddd3537a4f6fe85d98..3a28676ff606bbb94c5f4b9753f8b189d2c5268d 100644 (file)
--- a/libctsim/procsignal.cpp
+++ b/libctsim/procsignal.cpp
@@ -7,9 +7,7 @@
 **     Date Started:    Aug 1984
 **
 **  This is part of the CTSim program
-**  Copyright (c) 1983-2001 Kevin Rosenberg
-**
-**  $Id$
+**  Copyright (c) 1983-2009 Kevin Rosenberg
 **
 **  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
@@ -474,13 +472,13 @@ ProcessSignal::convertFilterMethodNameToID (const char* const filterMethodName)
 {
   int fmID = FILTER_METHOD_INVALID;
 
-  for (int i = 0; i < s_iFilterMethodCount; i++)
+  for (int i = 0; i < s_iFilterMethodCount; i++) {
     if (strcasecmp (filterMethodName, s_aszFilterMethodName[i]) == 0) {
       fmID = i;
       break;
     }
-
-    return (fmID);
+  }
+  return (fmID);
 }
 
 const char *
@@ -511,13 +509,13 @@ ProcessSignal::convertFilterGenerationNameToID (const char* const fgName)
 {
   int fgID = FILTER_GENERATION_INVALID;
 
-  for (int i = 0; i < s_iFilterGenerationCount; i++)
+  for (int i = 0; i < s_iFilterGenerationCount; i++) {
     if (strcasecmp (fgName, s_aszFilterGenerationName[i]) == 0) {
       fgID = i;
       break;
     }
-
-    return (fgID);
+  }
+  return (fgID);
 }
 
 const char *
@@ -547,21 +545,34 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
 {
   double* input = new double [m_nSignalPoints];
   int i;
+
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
   for (i = 0; i < m_nSignalPoints; i++)
     input[i] = constInput[i];
 
   if (m_idGeometry == Scanner::GEOMETRY_EQUILINEAR) {
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (int i = 0; i < m_nSignalPoints; i++) {
       int iDetFromCenter = i - (m_nSignalPoints / 2);
       input[i] *= m_dFocalLength / sqrt (m_dFocalLength * m_dFocalLength + iDetFromCenter * iDetFromCenter * m_dSignalInc * m_dSignalInc);
     }
   } else if (m_idGeometry == Scanner::GEOMETRY_EQUIANGULAR) {
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (int i = 0; i < m_nSignalPoints; i++) {
       int iDetFromCenter = i - (m_nSignalPoints / 2);
       input[i] *= m_dFocalLength * cos (iDetFromCenter * m_dSignalInc);
     }
   }
   if (m_idFilterMethod == FILTER_METHOD_CONVOLUTION) {
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (i = 0; i < m_nSignalPoints; i++)
       output[i] = convolve (input, m_dSignalInc, i, m_nSignalPoints);
   } else if (m_idFilterMethod == FILTER_METHOD_FOURIER) {
@@ -573,6 +584,9 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, m_nFilterPoints, FORWARD);
     delete inputSignal;
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
@@ -590,6 +604,9 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
     std::complex<double>* fftSignal = new std::complex<double> [m_nFilterPoints];
     finiteFourierTransform (inputSignal, fftSignal, FORWARD);
     delete inputSignal;
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (i = 0; i < m_nFilterPoints; i++)
       fftSignal[i] *= m_adFilter[i];
     double* inverseFourier = new double [m_nFilterPoints];
@@ -618,6 +635,9 @@ ProcessSignal::filterSignal (const float constInput[], double output[]) const
       m_adComplexFftInput[i][0] = input[i];
 
     fftw_execute (m_complexPlanForward);
+#if HAVE_OPENMP
+    #pragma omp parallel for
+#endif
     for (i = 0; i < m_nFilterPoints; i++) {
       m_adComplexFftSignal[i][0] = m_adFilter[i] * m_adComplexFftOutput[i][0];
       m_adComplexFftSignal[i][1] = m_adFilter[i] * m_adComplexFftOutput[i][1];