diff options
| -rw-r--r-- | include/astra/FanFlatBeamLineKernelProjector2D.h | 7 | ||||
| -rw-r--r-- | include/astra/FanFlatBeamLineKernelProjector2D.inl | 494 | ||||
| -rw-r--r-- | include/astra/FanFlatBeamStripKernelProjector2D.h | 6 | ||||
| -rw-r--r-- | include/astra/FanFlatBeamStripKernelProjector2D.inl | 643 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLineKernelProjector2D.h | 7 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLineKernelProjector2D.inl | 490 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLinearKernelProjector2D.h | 9 | ||||
| -rw-r--r-- | include/astra/ParallelBeamLinearKernelProjector2D.inl | 280 | ||||
| -rw-r--r-- | include/astra/ParallelBeamStripKernelProjector2D.h | 7 | ||||
| -rw-r--r-- | include/astra/ParallelBeamStripKernelProjector2D.inl | 489 |
10 files changed, 158 insertions, 2274 deletions
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.h b/include/astra/FanFlatBeamLineKernelProjector2D.h index 154504a..8ead007 100644 --- a/include/astra/FanFlatBeamLineKernelProjector2D.h +++ b/include/astra/FanFlatBeamLineKernelProjector2D.h @@ -178,6 +178,13 @@ public: float32 angleBetweenVectors(float32 _fAX, float32 _fAY, float32 _fBX, float32 _fBY); +protected: + /** Internal policy-based projection of a range of angles and range. + * (_i*From is inclusive, _i*To exclusive) */ + template <typename Policy> + void projectBlock_internal(int _iProjFrom, int _iProjTo, + int _iDetFrom, int _iDetTo, Policy& _policy); + }; //---------------------------------------------------------------------------------------- diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl index e686532..7ea2749 100644 --- a/include/astra/FanFlatBeamLineKernelProjector2D.inl +++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl @@ -29,11 +29,32 @@ $Id$ using namespace astra; -//---------------------------------------------------------------------------------------- -// PROJECT ALL template <typename Policy> void CFanFlatBeamLineKernelProjector2D::project(Policy& p) { + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CFanFlatBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} + +//---------------------------------------------------------------------------------------- +// PROJECT BLOCK +template <typename Policy> +void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) +{ // variables float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2; float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY; @@ -47,7 +68,7 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p) const SFanProjection * proj = 0; // loop angles - for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { // get theta if (pProjectionGeometry) { @@ -69,7 +90,7 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p) } // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; @@ -273,468 +294,3 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p) } // end loop detector } // end loop angles } - - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE PROJECTION -template <typename Policy> -void CFanFlatBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) -{ - // variables - float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2; - float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY; - int iVolumeIndex, iRayIndex, row, col, iDetector, x1; - bool switch_t; - - const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); - const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry); - - float32 old_theta, theta, alpha; - const SFanProjection * proj = 0; - - //get theta - if (pProjectionGeometry) { - old_theta = pProjectionGeometry->getProjectionAngle(_iProjection); - } - else if (pVecProjectionGeometry) { - proj = &pVecProjectionGeometry->getProjectionVectors()[_iProjection]; - old_theta = atan2(-proj->fSrcX, proj->fSrcY); - if (old_theta < 0) old_theta += 2*PI; - } else { - assert(false); - } - - switch_t = false; - if (old_theta >= 7*PIdiv4) old_theta -= 2*PI; - if (old_theta >= 3*PIdiv4) { - old_theta -= PI; - switch_t = true; - } - - // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { - - iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) continue; - - if (pProjectionGeometry) { - t = -pProjectionGeometry->indexToDetectorOffset(iDetector); - alpha = atan(t / pProjectionGeometry->getSourceDetectorDistance()); - t = sin(alpha) * pProjectionGeometry->getOriginSourceDistance(); - } - else if (pVecProjectionGeometry) { - float32 detX = proj->fDetSX + proj->fDetUX*(0.5f + iDetector); - float32 detY = proj->fDetSY + proj->fDetUY*(0.5f + iDetector); - alpha = angleBetweenVectors(-proj->fSrcX, -proj->fSrcY, detX - proj->fSrcX, detY - proj->fSrcY); - t = sin(alpha) * sqrt(proj->fSrcX*proj->fSrcX + proj->fSrcY*proj->fSrcY); - } else { - assert(false); - } - - if (switch_t) t = -t; - theta = old_theta + alpha; - - // precalculate sin, cos, 1/cos - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_sin_theta = 1.0f / sin_theta; - inv_cos_theta = 1.0f / cos_theta; - - // precalculate kernel limits - lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - updatePerRow = sin_theta * inv_cos_theta; - inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); - - // precalculate kernel limits - lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - updatePerCol = cos_theta * inv_sin_theta; - inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); - - // precalculate S and T - S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow); - T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol); - - // vertically - if (old_theta <= PIdiv4) { - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - // get coords - x1 = int((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerRow; - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; - - // left - if (x2 < 0.5f-S) { - I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridColCount()*/) {//x1 is always less than or equal to gridColCount because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+S) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // right - else if (x2 <= 1.0f) { - I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridColCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } - - // horizontally - else { - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY; - - // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - // get coords - x1 = int((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerCol; - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; - - // up - if (x2 < 0.5f-T) { - I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridRowCount()*/) {//x1 is always less than or equal to gridRowCount because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+T) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // down - else if (x2 <= 1.0f) { - I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridRowCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } // end loop col - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - - } // end loop detector - -} - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE RAY -template <typename Policy> -void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) -{ - // variables - float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2; - float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY; - int iVolumeIndex, iRayIndex, row, col, x1; - bool switch_t; - - const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry); - const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry); - - iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) return; - - float32 old_theta, theta, alpha; - if (pProjectionGeometry) { - old_theta = pProjectionGeometry->getProjectionAngle(_iProjection); - t = -pProjectionGeometry->indexToDetectorOffset(_iDetector); - alpha = atan(t / pProjectionGeometry->getSourceDetectorDistance()); - t = sin(alpha) * pProjectionGeometry->getOriginSourceDistance(); - } - else if (pVecProjectionGeometry) { - const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[_iProjection]; - old_theta = atan2(-proj->fSrcX, proj->fSrcY); - if (old_theta < 0) old_theta += 2*PI; - - float32 detX = proj->fDetSX + proj->fDetUX*(0.5f + _iDetector); - float32 detY = proj->fDetSY + proj->fDetUY*(0.5f + _iDetector); - alpha = angleBetweenVectors(-proj->fSrcX, -proj->fSrcY, detX - proj->fSrcX, detY - proj->fSrcY); - t = sin(alpha) * sqrt(proj->fSrcX*proj->fSrcX + proj->fSrcY*proj->fSrcY); - } else { - assert(false); - } - - switch_t = false; - if (old_theta >= 7*PIdiv4) old_theta -= 2*PI; - if (old_theta >= 3*PIdiv4) { - old_theta -= PI; - switch_t = true; - } - if (switch_t) t = -t; - theta = old_theta + alpha; - - // precalculate sin, cos, 1/cos - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_sin_theta = 1.0f / sin_theta; - inv_cos_theta = 1.0f / cos_theta; - - // precalculate kernel limits - lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - updatePerRow = sin_theta * inv_cos_theta; - inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); - - // precalculate kernel limits - lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - updatePerCol = cos_theta * inv_sin_theta; - inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); - - // precalculate S and T - S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow); - T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol); - - // vertically - if (old_theta <= PIdiv4) { - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - // get coords - x1 = int((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerRow; - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; - - // left - if (x2 < 0.5f-S) { - I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridColCount()*/) {//x1 is always less than or equal to gridColCount because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+S) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // right - else if (x2 <= 1.0f) { - I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridColCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } - - // horizontally - else { - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY; - - // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - // get coords - x1 = int((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerCol; - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; - - // up - if (x2 < 0.5f-T) { - I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridRowCount()*/) {//x1 is always less than or equal to gridRowCount because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+T) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // down - else if (x2 <= 1.0f) { - I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridRowCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } // end loop col - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); -} diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.h b/include/astra/FanFlatBeamStripKernelProjector2D.h index 7fed4c8..6b95427 100644 --- a/include/astra/FanFlatBeamStripKernelProjector2D.h +++ b/include/astra/FanFlatBeamStripKernelProjector2D.h @@ -174,6 +174,12 @@ public: */ virtual std::string getType(); +protected: + /** Internal policy-based projection of a range of angles and range. + * (_i*From is inclusive, _i*To exclusive) */ + template <typename Policy> + void projectBlock_internal(int _iProjFrom, int _iProjTo, + int _iDetFrom, int _iDetTo, Policy& _policy); }; //---------------------------------------------------------------------------------------- diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl index e3f8b29..f95a065 100644 --- a/include/astra/FanFlatBeamStripKernelProjector2D.inl +++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl @@ -29,11 +29,32 @@ $Id$ using namespace astra; -//---------------------------------------------------------------------------------------- -// PROJECT ALL template <typename Policy> void CFanFlatBeamStripKernelProjector2D::project(Policy& p) { + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} + +//---------------------------------------------------------------------------------------- +// PROJECT BLOCK +template <typename Policy> +void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) +{ ASTRA_ASSERT(m_bIsInitialized); // Some variables @@ -65,7 +86,7 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p) } // loop angles - for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { // get values theta = m_pProjectionGeometry->getProjectionAngle(iAngle); @@ -84,7 +105,7 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p) if (theta < PIdiv4) { // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; @@ -345,617 +366,3 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p) 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; -} diff --git a/include/astra/ParallelBeamLineKernelProjector2D.h b/include/astra/ParallelBeamLineKernelProjector2D.h index 9e2cc5a..64b5ef3 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.h +++ b/include/astra/ParallelBeamLineKernelProjector2D.h @@ -173,6 +173,13 @@ public: */ virtual std::string getType(); + +protected: + /** Internal policy-based projection of a range of angles and range. + * (_i*From is inclusive, _i*To exclusive) */ + template <typename Policy> + void projectBlock_internal(int _iProjFrom, int _iProjTo, + int _iDetFrom, int _iDetTo, Policy& _policy); }; inline std::string CParallelBeamLineKernelProjector2D::getType() diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl index 08bbe5f..149ac46 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.inl +++ b/include/astra/ParallelBeamLineKernelProjector2D.inl @@ -26,12 +26,32 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. $Id$ */ +template <typename Policy> +void CParallelBeamLineKernelProjector2D::project(Policy& p) +{ + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} //---------------------------------------------------------------------------------------- -// PROJECT ALL +// PROJECT BLOCK template <typename Policy> -void CParallelBeamLineKernelProjector2D::project(Policy& p) +void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) { // variables float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2; @@ -40,7 +60,7 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p) bool switch_t; // loop angles - for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { // get theta theta = m_pProjectionGeometry->getProjectionAngle(iAngle); @@ -72,7 +92,7 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p) T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol); // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; @@ -267,465 +287,3 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p) } - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE PROJECTION -template <typename Policy> -void CParallelBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) -{ - // variables - float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2; - float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY; - int iVolumeIndex, iRayIndex, row, col, iDetector, x1; - bool switch_t; - - // get theta - theta = m_pProjectionGeometry->getProjectionAngle(_iProjection); - switch_t = false; - if (theta >= 7*PIdiv4) theta -= 2*PI; - if (theta >= 3*PIdiv4) { - theta -= PI; - switch_t = true; - } - - // precalculate sin, cos, 1/cos - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_sin_theta = 1.0f / sin_theta; - inv_cos_theta = 1.0f / cos_theta; - - // precalculate kernel limits - lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - updatePerRow = sin_theta * inv_cos_theta; - inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); - - // precalculate kernel limits - lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - updatePerCol = cos_theta * inv_sin_theta; - inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); - - // precalculate S and T - S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow); - T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol); - - // 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 t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) t = -t; - - // vertically - if (theta <= PIdiv4) { - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX; - - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - x1 = nextx1; - x2 = nextx2; - - nextx2 += updatePerRow; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; - - // left - if (x2 < 0.5f-S) { - I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+S) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // right - else { - I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } - - // horizontally - else if (PIdiv4 <= theta && theta <= 3*PIdiv4) { - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY; - - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; - - // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - x1 = nextx1; - x2 = nextx2; - - nextx2 += updatePerCol; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; - - // up - if (x2 < 0.5f-T) { - I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+T) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // down - else { - I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } // end loop col - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - - } // end loop detector -} - - - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE RAY -template <typename Policy> -void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) -{ - int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) return; - - // variables - float32 t, I, P, x, x2; - int iVolumeIndex, row, col, x1; - - // get theta - float32 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; - } - - // get t - t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector); - if (switch_t) t = -t; - - // vertically - if (theta <= PIdiv4) { - - float32 sin_theta = sin(theta); - float32 inv_cos_theta = 1.0f / cos(theta); - - // precalculate kernel limits - float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - float32 updatePerRow = sin_theta * inv_cos_theta; - float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); - float32 S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow); - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX; - - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - x1 = nextx1; - x2 = nextx2; - - nextx2 += updatePerRow; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue; - - // left - if (x2 < 0.5f-S) { - I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+S) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // right - else { - I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } - - // horizontally - else if (PIdiv4 <= theta && theta <= 3*PIdiv4) { - - float32 cos_theta = cos(theta); - float32 inv_sin_theta = 1.0f / sin(theta); - - // precalculate kernel limits - float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - float32 updatePerCol = cos_theta * inv_sin_theta; - float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); - float32 T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol); - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY; - - // get coords - int nextx1 = int((x > 0.0f) ? x : x-1.0f); - float nextx2 = x - nextx1; - - // for each col - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - x1 = nextx1; - x2 = nextx2; - - nextx2 += updatePerCol; - while (nextx2 >= 1.0f) { - nextx2 -= 1.0f; - nextx1++; - } - while (nextx2 < 0.0f) { - nextx2 += 1.0f; - nextx1--; - } - - if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue; - - // up - if (x2 < 0.5f-T) { - I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // center - else if (x2 <= 0.5f+T) { - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - - // down - else { - I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol; - - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, I); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - } // end loop col - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - -} diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h index ac5899e..0a011f1 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.h +++ b/include/astra/ParallelBeamLinearKernelProjector2D.h @@ -177,6 +177,15 @@ public: */ virtual std::string getType(); + +protected: + /** Internal policy-based projection of a range of angles and range. + * (_i*From is inclusive, _i*To exclusive) */ + template <typename Policy> + void projectBlock_internal(int _iProjFrom, int _iProjTo, + int _iDetFrom, int _iDetTo, Policy& _policy); + + }; //---------------------------------------------------------------------------------------- diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl index 66b1cb2..d545f30 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.inl +++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl @@ -27,11 +27,31 @@ $Id$ */ +template <typename Policy> +void CParallelBeamLinearKernelProjector2D::project(Policy& p) +{ + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamLinearKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} //---------------------------------------------------------------------------------------- -// PROJECT ALL +// PROJECT BLOCK template <typename Policy> -void CParallelBeamLinearKernelProjector2D::project(Policy& p) +void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) { // variables float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, t; @@ -43,7 +63,7 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p) float32 P,x,x2; // loop angles - for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { // get theta theta = m_pProjectionGeometry->getProjectionAngle(iAngle); @@ -71,7 +91,7 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p) inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; @@ -162,255 +182,3 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p) } - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE PROJECTION -template <typename Policy> -void CParallelBeamLinearKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) -{ - // variables - float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, t; - float32 lengthPerRow, updatePerRow, inv_pixelLengthX; - float32 lengthPerCol, updatePerCol, inv_pixelLengthY; - bool switch_t; - int iDetector, iVolumeIndex, iRayIndex; - int row, col, x1; - float32 P,x,x2; - - // get theta - theta = m_pProjectionGeometry->getProjectionAngle(_iProjection); - switch_t = false; - if (theta >= 7*PIdiv4) theta -= 2*PI; - if (theta >= 3*PIdiv4) { - theta -= PI; - switch_t = true; - } - - // precalculate sin, cos, 1/cos - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_cos_theta = 1.0f / cos_theta; - inv_sin_theta = 1.0f / sin_theta; - - // precalculate kernel limits - lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - updatePerRow = sin_theta * inv_cos_theta; - inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX(); - - // precalculate kernel limits - lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - updatePerCol = cos_theta * inv_sin_theta; - inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY(); - - // 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 t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) { - t = -t; - } - - // vertically - if (theta <= PIdiv4) { - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = m_pVolumeGeometry->coordXToColF(P) - 0.5f; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - // get coords - x1 = (int)((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerRow; - - // add weights - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (x2) * lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - - // horizontally - else if (PIdiv4 <= theta && theta <= 3*PIdiv4) { - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = m_pVolumeGeometry->coordYToRowF(P) - 0.5f; - - // for each row - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - // get coords - x1 = (int)((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerCol; - - // add weights - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, x2 * lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - - } // end loop detector -} - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE RAY -template <typename Policy> -void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) -{ - int iVolumeIndex, iRayIndex; - - iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) return; - - // get theta - float32 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; - } - - // get t - float32 t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector); - if (switch_t) { - t = -t; - } - - // vertically - if (theta <= PIdiv4) { - - // precalculate sin, 1/cos - float32 sin_theta = sin(theta); - float32 inv_cos_theta = 1.0f / cos(theta); - - // precalculate kernel limits - float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta; - float32 updatePerRow = sin_theta * inv_cos_theta; - - int row, x1; - float32 P,x,x2; - - // calculate x for row 0 - P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta; - x = m_pVolumeGeometry->coordXToColF(P) - 0.5f; - - // for each row - for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) { - - // get coords - x1 = (int)((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerRow; - - // add weights - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (x2) * lengthPerRow); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - - // horizontally - else if (PIdiv4 <= theta && theta <= 3*PIdiv4) { - - // precalculate cos 1/sin - float32 cos_theta = cos(theta); - float32 inv_sin_theta = 1.0f / sin(theta); - - // precalculate kernel limits - float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta; - float32 updatePerCol = cos_theta * inv_sin_theta; - - int col, x1; - float32 P,x,x2; - - // calculate point P - P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta; - x = m_pVolumeGeometry->coordYToRowF(P) - 0.5f; - - // for each row - for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) { - - // get coords - x1 = (int)((x > 0.0f) ? x : x-1.0f); - x2 = x - x1; - x += updatePerCol; - - // add weights - if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) { - iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col); - // POLICY: PIXEL PRIOR + ADD + POSTERIOR - if (p.pixelPrior(iVolumeIndex)) { - p.addWeight(iRayIndex, iVolumeIndex, x2 * lengthPerCol); - p.pixelPosterior(iVolumeIndex); - } - } - } - } - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); -} diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h index de056de..e3808ca 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.h +++ b/include/astra/ParallelBeamStripKernelProjector2D.h @@ -176,6 +176,13 @@ protected: */ virtual std::string getType(); + + /** Internal policy-based projection of a range of angles and range. + * (_i*From is inclusive, _i*To exclusive) */ + template <typename Policy> + void projectBlock_internal(int _iProjFrom, int _iProjTo, + int _iDetFrom, int _iDetTo, Policy& _policy); + }; //---------------------------------------------------------------------------------------- diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl index c55fa8e..a9168ee 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.inl +++ b/include/astra/ParallelBeamStripKernelProjector2D.inl @@ -26,11 +26,31 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. $Id$ */ +template <typename Policy> +void CParallelBeamStripKernelProjector2D::project(Policy& p) +{ + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} //---------------------------------------------------------------------------------------- -// PROJECT ALL +// PROJECT BLOCK template <typename Policy> -void CParallelBeamStripKernelProjector2D::project(Policy& p) +void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) { ASTRA_ASSERT(m_bIsInitialized); @@ -52,7 +72,7 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p) float32 S, T, U, V, inv_4T; // loop angles - for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) { + for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { // get values theta = m_pProjectionGeometry->getProjectionAngle(iAngle); @@ -94,7 +114,7 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p) updateX = sin_theta * inv_cos_theta; // loop detectors - for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) { + for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; @@ -276,464 +296,3 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p) } -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE PROJECTION -template <typename Policy> -void CParallelBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p) -{ - // Some variables - float32 theta, t; - int row, col, iDetector; - float32 res; - float32 PL, PLimitL, PLimitR; - float32 xL, xR, XLimitL, XLimitR; - int x1L,x1R; - float32 x2L, x2R, updateX; - int iVolumeIndex, iRayIndex; - - float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta; - float32 fabs_sin_theta, fabs_cos_theta, fabs_inv_sin_theta, fabs_inv_cos_theta; - float32 PW, PH, DW, inv_PW, inv_PH; - float32 S, T, U, V, inv_4T; - - // 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 - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_cos_theta = 1.0f / cos_theta; - inv_sin_theta = 1.0f / sin_theta; - - fabs_sin_theta = (sin_theta < 0.0f) ? -sin_theta : sin_theta; - fabs_cos_theta = (cos_theta < 0.0f) ? -cos_theta : cos_theta; - fabs_inv_cos_theta = (inv_cos_theta < 0.0f) ? -inv_cos_theta : inv_cos_theta; - fabs_inv_sin_theta = (inv_sin_theta < 0.0f) ? -inv_sin_theta : inv_sin_theta; - - // Other precalculations - PW = m_pVolumeGeometry->getPixelLengthX(); - PH = m_pVolumeGeometry->getPixelLengthY(); - DW = m_pProjectionGeometry->getDetectorWidth(); - inv_PW = 1.0f / PW; - inv_PH = 1.0f / PH; - - // [-45?,45?] and [135?,225?] - if (theta < PIdiv4) { - - // Precalculate kernel limits - S = -0.5f * fabs_sin_theta * fabs_inv_cos_theta; - T = -S; - U = 1.0f + S; - V = 1.0f - S; - inv_4T = 0.25f / T; - - updateX = sin_theta * inv_cos_theta; - - // 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 t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) t = -t; - - // calculate left strip extremes (volume coordinates) - PL = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0) - DW*0.5f) * inv_cos_theta; - PLimitL = PL - 0.5f * fabs_sin_theta * fabs_inv_cos_theta * PH; - PLimitR = PLimitL + DW * inv_cos_theta + PH * fabs_sin_theta * fabs_inv_cos_theta; - - // calculate strip extremes (pixel coordinates) - XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - xR = xL + (DW * inv_cos_theta) * 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; - XLimitR += updateX; - xL += updateX; - xR += updateX; - - // 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) res = 1.0f; - else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T; - else if (x2R >= T) res = x2R; - else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T; - else { x2L -= 1.0f; x2R -= 1.0f; continue; } - - // left - if (x2L <= S) {} // - 0.0f - else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T; - else if (x2L <= U) res -= x2L; - else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T; - 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 { - - // Precalculate kernel limits - S = -0.5f * fabs_cos_theta * fabs_inv_sin_theta; - T = -S; - U = 1.0f + S; - V = 1.0f - S; - inv_4T = 0.25f / T; - - updateX = cos_theta * inv_sin_theta; - - // 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 t - t = m_pProjectionGeometry->indexToDetectorOffset(iDetector); - if (switch_t) t = -t; - - // calculate left strip extremes (volume coordinates) - PL = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0) + DW*0.5f) * inv_sin_theta; - PLimitL = PL + 0.5f * fabs_cos_theta * fabs_inv_sin_theta * PW; - PLimitR = PLimitL - DW * inv_sin_theta - PH * fabs_cos_theta * fabs_inv_sin_theta; - - // calculate strip extremes (pixel coordinates) - XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH; - XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH; - xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH; - xR = xL + (DW * fabs_inv_sin_theta) * 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; - XLimitR += updateX; - xL += updateX; - xR += updateX; - - // for each affected col - 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) res = 1.0f; - else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T; - else if (x2R >= T) res = x2R; - else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T; - else { x2L -= 1.0f; x2R -= 1.0f; continue; } - - // left - if (x2L <= S) {} // - 0.0f - else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T; - else if (x2L <= U) res -= x2L; - else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T; - 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 row loop - - } // end col loop - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - - } // end detector loop - - } // end theta switch -} - - -//---------------------------------------------------------------------------------------- -// PROJECT SINGLE RAY -template <typename Policy> -void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) -{ - int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector; - - // POLICY: RAY PRIOR - if (!p.rayPrior(iRayIndex)) return; - - // Some variables - float32 theta, t; - int row, col; - float32 res; - float32 PL, PLimitL, PLimitR; - float32 xL, xR, XLimitL, XLimitR; - int x1L,x1R; - float32 x2L, x2R, updateX; - int iVolumeIndex; - - float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta; - float32 fabs_sin_theta, fabs_cos_theta, fabs_inv_sin_theta, fabs_inv_cos_theta; - float32 PW, PH, DW, inv_PW, inv_PH; - float32 S, T, U, V, inv_4T; - - // 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 - sin_theta = sin(theta); - cos_theta = cos(theta); - inv_cos_theta = 1.0f / cos_theta; - inv_sin_theta = 1.0f / sin_theta; - - fabs_sin_theta = (sin_theta < 0.0f) ? -sin_theta : sin_theta; - fabs_cos_theta = (cos_theta < 0.0f) ? -cos_theta : cos_theta; - fabs_inv_cos_theta = (inv_cos_theta < 0.0f) ? -inv_cos_theta : inv_cos_theta; - fabs_inv_sin_theta = (inv_sin_theta < 0.0f) ? -inv_sin_theta : inv_sin_theta; - - // Other precalculations - PW = m_pVolumeGeometry->getPixelLengthX(); - PH = m_pVolumeGeometry->getPixelLengthY(); - DW = m_pProjectionGeometry->getDetectorWidth(); - inv_PW = 1.0f / PW; - inv_PH = 1.0f / PH; - - // [-45?,45?] and [135?,225?] - if (theta < PIdiv4) { - - // Precalculate kernel limits - S = -0.5f * fabs_sin_theta * fabs_inv_cos_theta; - T = -S; - U = 1.0f + S; - V = 1.0f - S; - inv_4T = 0.25f / T; - - updateX = sin_theta * inv_cos_theta; - - // get t - t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector); - if (switch_t) t = -t; - - // calculate left strip extremes (volume coordinates) - PL = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0) - DW*0.5f) * inv_cos_theta; - PLimitL = PL - 0.5f * fabs_sin_theta * fabs_inv_cos_theta * PH; - PLimitR = PLimitL + DW * inv_cos_theta + PH * fabs_sin_theta * fabs_inv_cos_theta; - - // calculate strip extremes (pixel coordinates) - XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW; - xR = xL + (DW * inv_cos_theta) * 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; - XLimitR += updateX; - xL += updateX; - xR += updateX; - - // 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) res = 1.0f; - else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T; - else if (x2R >= T) res = x2R; - else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T; - else { x2L -= 1.0f; x2R -= 1.0f; continue; } - - // left - if (x2L <= S) {} // - 0.0f - else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T; - else if (x2L <= U) res -= x2L; - else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T; - 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 - - // [45?,135?] and [225?,315?] - // horizontaly - } else { - - // Precalculate kernel limits - S = -0.5f * fabs_cos_theta * fabs_inv_sin_theta; - T = -S; - U = 1.0f + S; - V = 1.0f - S; - inv_4T = 0.25f / T; - - updateX = cos_theta * inv_sin_theta; - - // get t - t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector); - if (switch_t) t = -t; - - // calculate left strip extremes (volume coordinates) - PL = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0) + DW*0.5f) * inv_sin_theta; - PLimitL = PL + 0.5f * fabs_cos_theta * fabs_inv_sin_theta * PW; - PLimitR = PLimitL - DW * inv_sin_theta - PH * fabs_cos_theta * fabs_inv_sin_theta; - - // calculate strip extremes (pixel coordinates) - XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH; - XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH; - xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH; - xR = xL + (DW * fabs_inv_sin_theta) * 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; - XLimitR += updateX; - xL += updateX; - xR += updateX; - - // for each affected col - 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) res = 1.0f; - else if (x2R > U) res = x2R - (x2R-U)*(x2R-U)*inv_4T; - else if (x2R >= T) res = x2R; - else if (x2R > S) res = (x2R-S)*(x2R-S) * inv_4T; - else { x2L -= 1.0f; x2R -= 1.0f; continue; } - - // left - if (x2L <= S) {} // - 0.0f - else if (x2L < T) res -= (x2L-S)*(x2L-S) * inv_4T; - else if (x2L <= U) res -= x2L; - else if (x2L < V) res -= x2L - (x2L-U)*(x2L-U)*inv_4T; - 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 row loop - - } // end col loop - - } // end theta switch - - // POLICY: RAY POSTERIOR - p.rayPosterior(iRayIndex); - -} |