projects
/
ctsim.git
/ commitdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
| commitdiff |
tree
raw
|
patch
|
inline
| side by side (from parent 1:
747a2ec
)
Add OpenMP for more backprojectors
author
Kevin M. Rosenberg
<kevin@rosenberg.net>
Wed, 21 Mar 2018 23:12:53 +0000
(17:12 -0600)
committer
Kevin M. Rosenberg
<kevin@rosenberg.net>
Wed, 21 Mar 2018 23:12:53 +0000
(17:12 -0600)
libctsim/backprojectors.cpp
patch
|
blob
|
history
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);
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
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();
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);
}
r[ix][iy] = sqrt (x * x + y * y);
phi[ix][iy] = atan2 (y, x);
}
+ }
}
BackprojectTable::~BackprojectTable ()
}
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 (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];
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 (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];
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 (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];
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 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;
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)
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)
} 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)
} 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
}
} // end for y
} // end for x