Add OpenMP for more backprojectors

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 26a5a97c5861f919b666ab03ff1c3563f8366638..2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -371,8 +371,13 @@ BackprojectTrig::BackprojectView (const double* const filteredProj, const double
   if (interpType == Backprojector::INTERP_CUBIC)
     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 
-  double x = xMin + xInc / 2;   // Rectang coords of center of pixel
-  for (int ix = 0; ix < nx; x += xInc, ix++) {
+  double xstart = xMin + xInc / 2;   // Rectang coords of center of pixel
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int ix = 0; ix < nx; ix++) {
+    double x = xstart + (ix * xInc);
+
     double y = yMin + yInc / 2;
     for (int iy = 0; iy < ny; y += yInc, iy++) {
       double r = sqrt (x * x + y * y);   // distance of cell from center
@@ -419,13 +424,19 @@ BackprojectTable::BackprojectTable (const Projections& proj, ImageFile& im, int
   r = arrayR.getArray();
   phi = arrayPhi.getArray();
 
-  double x, y;                  // Rectang coords of center of pixel
-  int ix, iy;
-  for (x = xMin + xInc / 2, ix = 0; ix < nx; x += xInc, ix++)
-    for (y = yMin + yInc / 2, iy = 0; iy < ny; y += yInc, iy++) {
+  double xstart = xMin + xInc / 2;
+
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int ix = 0; ix < nx; ix++) {
+    double x = xstart + (ix * xInc);
+    double y = yMin + yInc / 2;
+    for (int iy = 0; iy < ny; iy++, y += yInc) {
       r[ix][iy] = sqrt (x * x + y * y);
       phi[ix][iy] = atan2 (y, x);
     }
+  }
 }
 
 BackprojectTable::~BackprojectTable ()
@@ -450,6 +461,9 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
   if (interpType == Backprojector::INTERP_CUBIC)
     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
 
@@ -667,6 +681,9 @@ BackprojectEquiangular::BackprojectView (const double* const filteredProj, const
   if (interpType == Backprojector::INTERP_CUBIC)
     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
 
@@ -710,6 +727,9 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
   if (interpType == Backprojector::INTERP_CUBIC)
     pCubicInterp = new CubicPolyInterpolator (filteredProj, nDet);
 
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
   for (int ix = 0; ix < nx; ix++) {
     ImageFileColumn pImCol = v[ix];
 
@@ -719,6 +739,8 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
       double rsin_t = r[ix][iy] * sin (dAngleDiff);
 
       double dU = (m_dFocalLength + rsin_t) / m_dFocalLength;
+      double dU2 = dU * dU;
+
       double dDetPos =  rcos_t / dU;
       // Scale for imaginary detector that passes through origin of phantom, see Kak-Slaney Figure 3.22.
       dDetPos *= m_dSourceDetectorLength / m_dFocalLength;
@@ -727,18 +749,17 @@ BackprojectEquilinear::BackprojectView (const double* const filteredProj, const
       if (interpType == Backprojector::INTERP_NEAREST) {
         int iDetPos = iDetCenter + nearest<int>(dPos);  // calc index in the filtered raysum vector
         if (iDetPos >= 0 && iDetPos < nDet)
-          pImCol[iy] += (filteredProj[iDetPos] / (dU * dU));
+          pImCol[iy] += filteredProj[iDetPos] / dU2;
       } else if (interpType == Backprojector::INTERP_LINEAR) {
         double dPosFloor = floor (dPos);
         int iDetPos = iDetCenter + static_cast<int>(dPosFloor);
         double frac = dPos - dPosFloor; // fraction distance from det
         if (iDetPos >= 0 && iDetPos < nDet - 1)
-          pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos]))
-                           / (dU * dU);
+          pImCol[iy] += (filteredProj[iDetPos] + frac * (filteredProj[iDetPos+1] - filteredProj[iDetPos])) / dU2;
       } else if (interpType == Backprojector::INTERP_CUBIC) {
         double d = iDetCenter + dPos;           // position along detector
         if (d >= 0 && d < nDet)
-          pImCol[iy] += pCubicInterp->interpolate (d) / (dU * dU);
+          pImCol[iy] += pCubicInterp->interpolate (d) / dU2;
       }
     }   // end for y
   }     // end for x