projects / ctsim.git / blobdiff
? search:
summary | shortlog | log | commit | commitdiff | tree
raw | inline | side by side
Add OpenMP for more backprojectors
[ctsim.git] / libctsim / backprojectors.cpp
diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp
index b492861cbbcad6a5cbcdcb0a22e238a34cbeb2e9..2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e 100644 (file)
--- a/libctsim/backprojectors.cpp
+++ b/libctsim/backprojectors.cpp
@@ -185,13 +185,13 @@ Backprojector::convertBackprojectNameToID (const char* const backprojName)
 {
   int backprojID = BPROJ_INVALID;
 
-  for (int i = 0; i < s_iBackprojectCount; i++)
+  for (int i = 0; i < s_iBackprojectCount; i++) {
     if (strcasecmp (backprojName, s_aszBackprojectName[i]) == 0) {
       backprojID = i;
       break;
     }
-
-    return (backprojID);
+  }
+  return (backprojID);
 }
 
 const char*
@@ -222,13 +222,13 @@ Backprojector::convertInterpNameToID (const char* const interpName)
 {
   int interpID = INTERP_INVALID;
 
-  for (int i = 0; i < s_iInterpCount; i++)
+  for (int i = 0; i < s_iInterpCount; i++) {
     if (strcasecmp (interpName, s_aszInterpName[i]) == 0) {
       interpID = i;
       break;
     }
-
-    return (interpID);
+  }
+  return (interpID);
 }
 
 const char*
@@ -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];
 
@@ -459,19 +473,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 +541,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 +556,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 +617,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 +632,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)
@@ -662,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];
 
@@ -705,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];
 
@@ -714,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;
@@ -722,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
Computed Tomography simulator and reconstructor