diff options
Diffstat (limited to 'include/astra/FanFlatBeamStripKernelProjector2D.inl')
| -rw-r--r-- | include/astra/FanFlatBeamStripKernelProjector2D.inl | 961 |
1 files changed, 961 insertions, 0 deletions
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl new file mode 100644 index 0000000..e3f8b29 --- /dev/null +++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl @@ -0,0 +1,961 @@ +/* +----------------------------------------------------------------------- +Copyright 2012 iMinds-Vision Lab, University of Antwerp + +Contact: astra@ua.ac.be +Website: http://astra.ua.ac.be + + +This file is part of the +All Scale Tomographic Reconstruction Antwerp Toolbox ("ASTRA Toolbox"). + +The ASTRA Toolbox is free software: you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation, either version 3 of the License, or +(at your option) any later version. + +The ASTRA Toolbox is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +$Id$ +*/ + + +using namespace astra; + +//---------------------------------------------------------------------------------------- +// PROJECT ALL +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::project(Policy& p) +{ + ASTRA_ASSERT(m_bIsInitialized); + + // Some variables + float32 theta; + int row, col; + int iAngle, iDetector; + float32 res; + int x1L, x1R; + float32 x2L, x2R; + int iVolumeIndex, iRayIndex; + + CFanFlatProjectionGeometry2D* projgeom = static_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); + + // Other precalculations + float32 PW = m_pVolumeGeometry->getPixelLengthX(); + float32 PH = m_pVolumeGeometry->getPixelLengthY(); + float32 DW = m_pProjectionGeometry->getDetectorWidth(); + float32 inv_PW = 1.0f / PW; + float32 inv_PH = 1.0f / PH; + + // calculate alpha's + float32 alpha; + float32* cos_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + float32* sin_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + for (int i = 0; i < m_pProjectionGeometry->getDetectorCount() + 1; ++i) { + alpha = -atan((i - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW / projgeom->getSourceDetectorDistance()); + cos_alpha[i] = cos(alpha); + sin_alpha[i] = sin(alpha); + } + + // loop angles + for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + + // get values + theta = m_pProjectionGeometry->getProjectionAngle(iAngle); + bool switch_t = false; + if (theta >= 7*PIdiv4) theta -= 2*PI; + if (theta >= 3*PIdiv4) { + theta -= PI; + switch_t = true; + } + + // Precalculate sin, cos, 1/cos + float32 sin_theta = sin(theta); + float32 cos_theta = cos(theta); + + // [-45?,45?] and [135?,225?] + if (theta < PIdiv4) { + + // loop detectors + for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + + iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + + cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + + t_l = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + + t_l = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + } + + float32 inv_cos_theta_left = 1.0f / cos_theta_left; + float32 inv_cos_theta_right = 1.0f / cos_theta_right; + + float32 updateX_left = sin_theta_left * inv_cos_theta_left; + float32 updateX_right = sin_theta_right * inv_cos_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) {S_l = -S_l;} + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) {S_r = -S_r;} + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - sin_theta_left * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_left; + float32 PR = (t_r - sin_theta_right * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_right; + float32 PLimitL = PL + S_l * PH; + float32 PLimitR = PR - S_r * PH; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xR = (PR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + + // for each row + for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected col + for (col = x1L; col <= x1R; ++col) { + + if (col < 0 || col >= m_pVolumeGeometry->getGridColCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end detector loop + + // [45?,135?] and [225?,315?] + // horizontaly + } else { + + // loop detectors + for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + + iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + + t_l = sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + + cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + + t_l = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + } + + float32 inv_sin_theta_left = 1.0f / sin_theta_left; + float32 inv_sin_theta_right = 1.0f / sin_theta_right; + + float32 updateX_left = cos_theta_left * inv_sin_theta_left; + float32 updateX_right = cos_theta_right * inv_sin_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) { S_l = -S_l; } + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) { S_r = -S_r; } + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - cos_theta_left * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_left; + float32 PR = (t_r - cos_theta_right * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_right; + float32 PLimitL = PL - S_l * PW; + float32 PLimitR = PR + S_r * PW; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH; + float32 XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH; + float32 xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH; + float32 xR = (m_pVolumeGeometry->getWindowMaxY() - PR) * inv_PH; + + // for each col + for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected row + for (row = x1L; row <= x1R; ++row) { + + if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end detector loop + + } // end theta switch + + } // end angle loop + + delete[] cos_alpha; + delete[] sin_alpha; +} + + +//---------------------------------------------------------------------------------------- +// PROJECT SINGLE PROJECTION +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + ASTRA_ASSERT(m_bIsInitialized); + + // Some variables + float32 theta; + int row, col; + int iDetector; + float32 res; + int x1L, x1R; + float32 x2L, x2R; + int iVolumeIndex, iRayIndex; + + CFanFlatProjectionGeometry2D* projgeom = static_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); + + // Other precalculations + float32 PW = m_pVolumeGeometry->getPixelLengthX(); + float32 PH = m_pVolumeGeometry->getPixelLengthY(); + float32 DW = m_pProjectionGeometry->getDetectorWidth(); + float32 inv_PW = 1.0f / PW; + float32 inv_PH = 1.0f / PH; + + // calculate alpha's + float32 alpha; + float32* cos_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + float32* sin_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + for (int i = 0; i < m_pProjectionGeometry->getDetectorCount() + 1; ++i) { + alpha = -atan((i - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW / projgeom->getSourceDetectorDistance()); + cos_alpha[i] = cos(alpha); + sin_alpha[i] = sin(alpha); + } + + + // get values + theta = m_pProjectionGeometry->getProjectionAngle(_iProjection); + bool switch_t = false; + if (theta >= 7*PIdiv4) theta -= 2*PI; + if (theta >= 3*PIdiv4) { + theta -= PI; + switch_t = true; + } + + // Precalculate sin, cos, 1/cos + float32 sin_theta = sin(theta); + float32 cos_theta = cos(theta); + + // [-45?,45?] and [135?,225?] + if (theta < PIdiv4) { + + // loop detectors + for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + + iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + + cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + + t_l = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + + t_l = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + } + + float32 inv_cos_theta_left = 1.0f / cos_theta_left; + float32 inv_cos_theta_right = 1.0f / cos_theta_right; + + float32 updateX_left = sin_theta_left * inv_cos_theta_left; + float32 updateX_right = sin_theta_right * inv_cos_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) {S_l = -S_l;} + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) {S_r = -S_r;} + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - sin_theta_left * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_left; + float32 PR = (t_r - sin_theta_right * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_right; + float32 PLimitL = PL + S_l * PH; + float32 PLimitR = PR - S_r * PH; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xR = (PR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + + // for each row + for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected col + for (col = x1L; col <= x1R; ++col) { + + if (col < 0 || col >= m_pVolumeGeometry->getGridColCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end detector loop + + // [45?,135?] and [225?,315?] + // horizontaly + } else { + + // loop detectors + for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + + iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + + t_l = sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector]; + + cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector]; + + t_l = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance(); + } + + float32 inv_sin_theta_left = 1.0f / sin_theta_left; + float32 inv_sin_theta_right = 1.0f / sin_theta_right; + + float32 updateX_left = cos_theta_left * inv_sin_theta_left; + float32 updateX_right = cos_theta_right * inv_sin_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) { S_l = -S_l; } + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) { S_r = -S_r; } + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - cos_theta_left * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_left; + float32 PR = (t_r - cos_theta_right * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_right; + float32 PLimitL = PL - S_l * PW; + float32 PLimitR = PR + S_r * PW; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH; + float32 XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH; + float32 xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH; + float32 xR = (m_pVolumeGeometry->getWindowMaxY() - PR) * inv_PH; + + // for each col + for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected row + for (row = x1L; row <= x1R; ++row) { + + if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end detector loop + + } // end theta switch + + delete[] cos_alpha; + delete[] sin_alpha; +} + + +//---------------------------------------------------------------------------------------- +// PROJECT SINGLE RAY +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + ASTRA_ASSERT(m_bIsInitialized); + + // Some variables + float32 theta; + int row, col; + float32 res; + int x1L, x1R; + float32 x2L, x2R; + int iVolumeIndex, iRayIndex; + + CFanFlatProjectionGeometry2D* projgeom = static_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); + + // Other precalculations + float32 PW = m_pVolumeGeometry->getPixelLengthX(); + float32 PH = m_pVolumeGeometry->getPixelLengthY(); + float32 DW = m_pProjectionGeometry->getDetectorWidth(); + float32 inv_PW = 1.0f / PW; + float32 inv_PH = 1.0f / PH; + + // calculate alpha's + float32 alpha; + float32* cos_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + float32* sin_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1]; + for (int i = 0; i < m_pProjectionGeometry->getDetectorCount() + 1; ++i) { + alpha = -atan((i - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW / projgeom->getSourceDetectorDistance()); + cos_alpha[i] = cos(alpha); + sin_alpha[i] = sin(alpha); + } + + + // get values + theta = m_pProjectionGeometry->getProjectionAngle(_iProjection); + bool switch_t = false; + if (theta >= 7*PIdiv4) theta -= 2*PI; + if (theta >= 3*PIdiv4) { + theta -= PI; + switch_t = true; + } + + // Precalculate sin, cos, 1/cos + float32 sin_theta = sin(theta); + float32 cos_theta = cos(theta); + + // [-45?,45?] and [135?,225?] + if (theta < PIdiv4) { + + iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) { delete[] cos_alpha; delete[] sin_alpha; return; } + + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector]; + + cos_theta_left = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector]; + + t_l = sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[_iDetector] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector]; + sin_theta_right = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector]; + cos_theta_right = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1]; + + t_l = -sin_alpha[_iDetector] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance(); + } + + float32 inv_cos_theta_left = 1.0f / cos_theta_left; + float32 inv_cos_theta_right = 1.0f / cos_theta_right; + + float32 updateX_left = sin_theta_left * inv_cos_theta_left; + float32 updateX_right = sin_theta_right * inv_cos_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) {S_l = -S_l;} + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) {S_r = -S_r;} + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - sin_theta_left * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_left; + float32 PR = (t_r - sin_theta_right * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_right; + float32 PLimitL = PL + S_l * PH; + float32 PLimitR = PR - S_r * PH; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + float32 xR = (PR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; + + // for each row + for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected col + for (col = x1L; col <= x1R; ++col) { + + if (col < 0 || col >= m_pVolumeGeometry->getGridColCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + // [45?,135?] and [225?,315?] + // horizontaly + } else { + + iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) { delete[] cos_alpha; delete[] sin_alpha; return; } + + // get theta_l = alpha_left + theta and theta_r = alpha_right + theta + float32 sin_theta_left, cos_theta_left; + float32 sin_theta_right, cos_theta_right; + float32 t_l, t_r; + if (!switch_t) { + sin_theta_left = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector]; + sin_theta_right = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1]; + + cos_theta_left = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector]; + cos_theta_right = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1]; + + t_l = sin_alpha[_iDetector] * projgeom->getOriginSourceDistance(); + t_r = sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance(); + + } else { + sin_theta_left = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1]; + sin_theta_right = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector]; + + cos_theta_left = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1]; + cos_theta_right = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector]; + + t_l = -sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance(); + t_r = -sin_alpha[_iDetector] * projgeom->getOriginSourceDistance(); + } + + float32 inv_sin_theta_left = 1.0f / sin_theta_left; + float32 inv_sin_theta_right = 1.0f / sin_theta_right; + + float32 updateX_left = cos_theta_left * inv_sin_theta_left; + float32 updateX_right = cos_theta_right * inv_sin_theta_right; + + // Precalculate kernel limits + float32 S_l = -0.5f * updateX_left; + if (S_l > 0) { S_l = -S_l; } + float32 T_l = -S_l; + float32 U_l = 1.0f + S_l; + float32 V_l = 1.0f - S_l; + float32 inv_4T_l = 0.25f / T_l; + + float32 S_r = -0.5f * updateX_right; + if (S_r > 0) { S_r = -S_r; } + float32 T_r = -S_r; + float32 U_r = 1.0f + S_r; + float32 V_r = 1.0f - S_r; + float32 inv_4T_r = 0.25f / T_r; + + // calculate strip extremes (volume coordinates) + float32 PL = (t_l - cos_theta_left * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_left; + float32 PR = (t_r - cos_theta_right * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_right; + float32 PLimitL = PL - S_l * PW; + float32 PLimitR = PR + S_r * PW; + + // calculate strip extremes (pixel coordinates) + float32 XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH; + float32 XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH; + float32 xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH; + float32 xR = (m_pVolumeGeometry->getWindowMaxY() - PR) * inv_PH; + + // for each col + for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { + + // get strip extremes in column indices + x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f); + x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f); + + // get coords w.r.t leftmost column hit by strip + x2L = xL - x1L; + x2R = xR - x1L; + + // update strip extremes for the next row + XLimitL += updateX_left; + XLimitR += updateX_right; + xL += updateX_left; + xR += updateX_right; + + // for each affected row + for (row = x1L; row <= x1R; ++row) { + + if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col); + + // POLICY: PIXEL PRIOR + if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // right + if (x2R >= V_r) res = 1.0f; + else if (x2R > U_r) res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r; + else if (x2R >= T_r) res = x2R; + else if (x2R > S_r) res = (x2R-S_r)*(x2R-S_r) * inv_4T_r; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // left + if (x2L <= S_l) {} + else if (x2L < T_l) res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l; + else if (x2L <= U_l) res -= x2L; + else if (x2L < V_l) res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l; + else { x2L -= 1.0f; x2R -= 1.0f; continue; } + + // POLICY: ADD + p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res); + + // POLICY: PIXEL POSTERIOR + p.pixelPosterior(iVolumeIndex); + + x2L -= 1.0f; + x2R -= 1.0f; + + } // end col loop + + } // end row loop + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } // end theta switch + + delete[] cos_alpha; + delete[] sin_alpha; +} |