Use OpenMP for backprojection

diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index 49e3a217bc8a63cb4fd065fd96d22a224292d114..26a5a97c5861f919b666ab03ff1c3563f8366638 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -459,19 +459,22 @@ BackprojectTable::BackprojectView (const double* const filteredProj, const doubl
       if (interpType == Backprojector::INTERP_NEAREST) {
         int iDetPos = iDetCenter + nearest<int>(L / detInc);    // calc index in the filtered raysum vector
 
-        if (iDetPos >= 0 && iDetPos < nDet)
+        if (iDetPos >= 0 && iDetPos < nDet) {
           pImCol[iy] += filteredProj[iDetPos];
+        }
       } else if (interpType == Backprojector::INTERP_LINEAR) {
         double dPos = L / detInc;               // position along detector
         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)
+        if (iDetPos >= 0 && iDetPos < nDet - 1) {
           pImCol[iy] += ((1-frac) * filteredProj[iDetPos] + frac * filteredProj[iDetPos+1]);
+        }
       } else if (interpType == Backprojector::INTERP_CUBIC) {
         double p = iDetCenter + (L / detInc);   // position along detector
-        if (p >= 0 && p < nDet)
+        if (p >= 0 && p < nDet) {
           pImCol[iy] += pCubicInterp->interpolate (p);
+        }
       }
     }   // end for y
   }     // end for x
@@ -524,7 +527,7 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   double det_dy = yInc * sin (theta) / detInc;
 
   // calculate detPosition for first point in image (ix=0, iy=0)
-  double detPosColStart = iDetCenter + start_r * cos (theta - start_phi) / detInc;
+  double detPosColBase = iDetCenter + start_r * cos (theta - start_phi) / detInc;
 
   CubicPolyInterpolator* pCubicInterp = NULL;
   double* deltaFilteredProj = NULL;
@@ -539,26 +542,32 @@ BackprojectDiff::BackprojectView (const double* const filteredProj, const double
   }
 
   int iLastDet = nDet - 1;
-  for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
-    double curDetPos = detPosColStart;
+
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int ix = 0; ix < nx; ix++) {
+    double detPos = detPosColBase + (ix * det_dx);
     ImageFileColumn pImCol = v[ix];
 
-    for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
+    for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
       if (interpType == Backprojector::INTERP_NEAREST) {
-        int iDetPos = nearest<int> (curDetPos); // calc index in the filtered raysum vector
-
-        if (iDetPos >= 0 && iDetPos < nDet)
+        int iDetPos = nearest<int> (detPos); // calc index in the filtered raysum vector
+        if (iDetPos >= 0 && iDetPos < nDet) {
           *pImCol++ += filteredProj[iDetPos];
+        }
       } else if (interpType == Backprojector::INTERP_LINEAR) {
-        double detPosFloor = floor (curDetPos);
+        double detPosFloor = floor (detPos);
         int iDetPos = static_cast<int>(detPosFloor);
-        double frac = curDetPos - detPosFloor;  // fraction distance from det
-        if (iDetPos >= 0 && iDetPos <= iLastDet)
-            *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
+        double frac = detPos - detPosFloor;  // fraction distance from det
+        if (iDetPos >= 0 && iDetPos <= iLastDet) {
+          *pImCol++ += filteredProj[iDetPos] + (frac * deltaFilteredProj[iDetPos]);
+        }
       } else if (interpType == Backprojector::INTERP_CUBIC) {
-        double p = curDetPos;   // position along detector
-        if (p >= 0 && p < nDet)
+        double p = detPos;   // position along detector
+        if (p >= 0 && p < nDet) {
           *pImCol++  += pCubicInterp->interpolate (p);
+        }
       }
     }   // end for y
   }     // end for x
@@ -594,7 +603,7 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
   const long det_dy = nearest<long> (yInc * sin (theta) / detInc * scale);
 
   // calculate L for first point in image (0, 0)
-  long detPosColStart = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
+  long detPosColBase = nearest<long> ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale);
 
   double* deltaFilteredProj = NULL;
   CubicPolyInterpolator* pCubicInterp = NULL;
@@ -609,41 +618,37 @@ BackprojectIntDiff::BackprojectView (const double* const filteredProj, const dou
   }
 
   int iLastDet = nDet - 1;
-  for (int ix = 0; ix < nx; ix++, detPosColStart += det_dx) {
-    long curDetPos = detPosColStart;
+#if HAVE_OPENMP
+  #pragma omp parallel for
+#endif
+  for (int ix = 0; ix < nx; ix++) {
+    long detPos = detPosColBase + (ix * det_dx);
     ImageFileColumn pImCol = v[ix];
 
-    if (interpType == Backprojector::INTERP_NEAREST) {
-      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        const int iDetPos = (curDetPos + halfScale) >> scaleShift;
-        if (iDetPos >= 0 && iDetPos <= iLastDet)
+    for (int iy = 0; iy < ny; iy++, detPos += det_dy) {
+      if (interpType == Backprojector::INTERP_NEAREST) {
+        const int iDetPos = (detPos + halfScale) >> scaleShift;
+        if (iDetPos >= 0 && iDetPos <= iLastDet) {
           *pImCol++ += filteredProj[iDetPos];
-        else
+        } else
           pImCol++;
-
-      } // end for iy
-    } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
-      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        const int iDetPos = ((curDetPos + halfScale) >> scaleShift) * m_interpFactor;
-        if (iDetPos >= 0 && iDetPos <= iLastDet)
+      } else if (interpType == Backprojector::INTERP_FREQ_PREINTERPOLATION) {
+        const int iDetPos = ((detPos + halfScale) >> scaleShift) * m_interpFactor;
+        if (iDetPos >= 0 && iDetPos <= iLastDet) {
           *pImCol++ += filteredProj[iDetPos];
-        else
+        } else
           pImCol++;
-      } // end for iy
-    } else if (interpType == Backprojector::INTERP_LINEAR) {
-      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        const long iDetPos = curDetPos >> scaleShift;
+      } else if (interpType == Backprojector::INTERP_LINEAR) {
+        const long iDetPos = detPos >> scaleShift;
         if (iDetPos >= 0 && iDetPos <= iLastDet) {
-          const long detRemainder = curDetPos & scaleBitmask;
+          const long detRemainder = detPos & scaleBitmask;
           *pImCol++ += filteredProj[iDetPos] + (detRemainder * deltaFilteredProj[iDetPos]);
         } else
           pImCol++;
-      } // end for iy
-    } else if (interpType == Backprojector::INTERP_CUBIC) {
-      for (int iy = 0; iy < ny; iy++, curDetPos += det_dy) {
-        *pImCol++ += pCubicInterp->interpolate (static_cast<double>(curDetPos) / scale);
-      }
-    } // end Cubic
+      } else if (interpType == Backprojector::INTERP_CUBIC) {
+        *pImCol++ += pCubicInterp->interpolate (static_cast<double>(detPos) / scale);
+      } // end Cubic
+    } // end for iy
   } // end for ix
 
   if (interpType == Backprojector::INTERP_LINEAR)