X-Git-Url: http://git.kpe.io/?p=ctsim.git;a=blobdiff_plain;f=libctsim%2Fbackprojectors.cpp;h=2b9eba45cdc4ec7ce4b300030c19a4fafaf5af1e;hp=b492861cbbcad6a5cbcdcb0a22e238a34cbeb2e9;hb=352a0691b9bd67f6c93ea822b353be3c101f4adb;hpb=f13a8c004b8f182b42d9e4df2bcd7c7f030bf1ad diff --git a/libctsim/backprojectors.cpp b/libctsim/backprojectors.cpp index b492861..2b9eba4 100644 --- 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(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(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 (curDetPos); // calc index in the filtered raysum vector - - if (iDetPos >= 0 && iDetPos < nDet) + int iDetPos = nearest (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(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 (yInc * sin (theta) / detInc * scale); // calculate L for first point in image (0, 0) - long detPosColStart = nearest ((start_r * cos (theta - start_phi) / detInc + iDetCenter) * scale); + long detPosColBase = nearest ((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(curDetPos) / scale); - } - } // end Cubic + } else if (interpType == Backprojector::INTERP_CUBIC) { + *pImCol++ += pCubicInterp->interpolate (static_cast(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(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(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