From b9563df34df98480d35a53caa2b81d998440f9c5 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 24 Jan 2019 13:26:30 +0100 Subject: Add basic implementation of par2d CPU Distance Driven projector --- .../ParallelBeamDistanceDrivenProjector2D.inl | 227 +++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 include/astra/ParallelBeamDistanceDrivenProjector2D.inl (limited to 'include/astra/ParallelBeamDistanceDrivenProjector2D.inl') diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl new file mode 100644 index 0000000..6bf3b56 --- /dev/null +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl @@ -0,0 +1,227 @@ +/* +----------------------------------------------------------------------- +Copyright: 2010-2018, imec Vision Lab, University of Antwerp + 2014-2018, CWI, Amsterdam + +Contact: astra@astra-toolbox.com +Website: http://www.astra-toolbox.com/ + +This file is part of the 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 . + +----------------------------------------------------------------------- +*/ + +#define policy_weight(p,rayindex,volindex,weight) do { if (p.pixelPrior(volindex)) { p.addWeight(rayindex, volindex, weight); p.pixelPosterior(volindex); } } while (false) + +template +void CParallelBeamDistanceDrivenProjector2D::project(Policy& p) +{ + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template +void CParallelBeamDistanceDrivenProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template +void CParallelBeamDistanceDrivenProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} + + + + + +template +void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) +{ + // get vector geometry + const CParallelVecProjectionGeometry2D* pVecProjectionGeometry; + if (dynamic_cast(m_pProjectionGeometry)) { + pVecProjectionGeometry = dynamic_cast(m_pProjectionGeometry)->toVectorGeometry(); + } else { + pVecProjectionGeometry = dynamic_cast(m_pProjectionGeometry); + } + + // precomputations + const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX(); + const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY(); + const float32 inv_pixelLengthX = 1.0f / pixelLengthX; + const float32 inv_pixelLengthY = 1.0f / pixelLengthY; + const int colCount = m_pVolumeGeometry->getGridColCount(); + const int rowCount = m_pVolumeGeometry->getGridRowCount(); + + // Performance note: + // This is not a very well optimizated version of the distance driven + // projector. The CPU projector model in ASTRA requires ray-driven iteration, + // which limits re-use of intermediate computations. + + // loop angles + for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) { + + const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; + + float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); + + const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); + + const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; + const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f; + + // loop detectors + for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) { + + const int iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; + + // POLICY: RAY PRIOR + if (!p.rayPrior(iRayIndex)) continue; + + const float32 Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; + const float32 Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; + + if (vertical && true) { + + const float32 RxOverRy = proj->fRayX/proj->fRayY; + // TODO: Determine det/pixel scaling factors + const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; + const float32 deltad = fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX); + + // calculate c for row 0 + float32 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX + 0.5f; + + // loop rows + for (int row = 0; row < rowCount; ++row, c+= deltac) { + + // horizontal extent of ray in center of this row: + // [ c - deltad/2 , c + deltad/2 ] + + // |-gapBegin-*---|------|----*-gapEnd-| + // * = ray extent intercepts; c - deltad/2 and c + deltad/2 + // | = pixel column edges + + const int colBegin = (int)floor(c - deltad/2.0f); + const int colEnd = (int)ceil(c + deltad/2.0f); + + // TODO: Optimize volume edge checks + + int iVolumeIndex = row * colCount + colBegin; + if (colBegin + 1 == colEnd) { + + if (colBegin >= 0 && colBegin < colCount) + policy_weight(p, iRayIndex, iVolumeIndex, + deltad * lengthPerRow); + } else { + const float gapBegin = (c - deltad/2.0f) - (float32)colBegin; + const float gapEnd = (float32)colEnd - (c + deltad/2.0f); + float tot = 1.0f - gapBegin; + if (colBegin >= 0 && colBegin < colCount) { + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapBegin) * lengthPerRow); + } + iVolumeIndex++; + + for (int col = colBegin + 1; col + 1 < colEnd; ++col, ++iVolumeIndex) { + tot += 1.0f; + if (col >= 0 && col < colCount) { + policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow); + } + } + assert(iVolumeIndex == row * colCount + colEnd - 1); + tot += 1.0f - gapEnd; + if (colEnd > 0 && colEnd <= colCount) { + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapEnd) * lengthPerRow); + } + assert(fabs(tot - deltad) < 0.0001); + } + + } + + } else if (!vertical && true) { + + const float32 RyOverRx = proj->fRayY/proj->fRayX; + // TODO: Determine det/pixel scaling factors + const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; + const float32 deltad = fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY); + + // calculate r for col 0 + float32 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY + 0.5f; + + // loop columns + for (int col = 0; col < colCount; ++col, r+= deltar) { + + // vertical extent of ray in center of this column: + // [ r - deltad/2 , r + deltad/2 ] + + const int rowBegin = (int)floor(r - deltad/2.0f); + const int rowEnd = (int)ceil(r + deltad/2.0f); + + // TODO: Optimize volume edge checks + + int iVolumeIndex = rowBegin * colCount + col; + if (rowBegin + 1 == rowEnd) { + + if (rowBegin >= 0 && rowBegin < rowCount) + policy_weight(p, iRayIndex, iVolumeIndex, + deltad * lengthPerCol); + } else { + const float gapBegin = (r - deltad/2.0f) - (float32)rowBegin; + const float gapEnd = (float32)rowEnd - (r + deltad/2.0f); + float tot = 1.0f - gapBegin; + + if (rowBegin >= 0 && rowBegin < rowCount) { + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapBegin) * lengthPerCol); + } + iVolumeIndex += colCount; + + for (int row = rowBegin + 1; row + 1 < rowEnd; ++row, iVolumeIndex += colCount) { + tot += 1.0f; + if (row >= 0 && row < rowCount) { + policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol); + } + } + assert(iVolumeIndex == (rowEnd - 1) * colCount + col); + tot += 1.0f - gapEnd; + if (rowEnd > 0 && rowEnd <= rowCount) { + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapEnd) * lengthPerCol); + } + assert(fabs(tot - deltad) < 0.0001); + } + + } + + } + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } + } + + if (dynamic_cast(m_pProjectionGeometry)) + delete pVecProjectionGeometry; +} -- cgit v1.2.3 From e2d5375ecc5c8fc45b796e0bd9d948cd729abcc9 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 24 Jan 2019 13:48:49 +0100 Subject: Remove largely unimplemented CProjector2D::projectPoint method This includes the astra_mex_projector('splat') matlab function. --- include/astra/ParallelBeamDistanceDrivenProjector2D.inl | 2 -- 1 file changed, 2 deletions(-) (limited to 'include/astra/ParallelBeamDistanceDrivenProjector2D.inl') diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl index 6bf3b56..7b45ed1 100644 --- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl @@ -81,8 +81,6 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle]; - float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); - const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f; -- cgit v1.2.3 From 37643b309520bcd12afcee430210632db305d21f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Thu, 24 Jan 2019 14:18:11 +0100 Subject: Some basic optimizations --- .../ParallelBeamDistanceDrivenProjector2D.inl | 87 ++++++++++------------ 1 file changed, 41 insertions(+), 46 deletions(-) (limited to 'include/astra/ParallelBeamDistanceDrivenProjector2D.inl') diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl index 7b45ed1..aedcee9 100644 --- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl @@ -97,13 +97,12 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro const float32 Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX; const float32 Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY; - if (vertical && true) { + if (vertical) { const float32 RxOverRy = proj->fRayX/proj->fRayY; - // TODO: Determine det/pixel scaling factors const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX; - const float32 deltad = fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX); + const float32 deltad = 0.5f * fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX); // calculate c for row 0 float32 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX + 0.5f; @@ -112,57 +111,55 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro for (int row = 0; row < rowCount; ++row, c+= deltac) { // horizontal extent of ray in center of this row: - // [ c - deltad/2 , c + deltad/2 ] + // [ c - deltad , c + deltad ] // |-gapBegin-*---|------|----*-gapEnd-| - // * = ray extent intercepts; c - deltad/2 and c + deltad/2 + // * = ray extent intercepts; c - deltad and c + deltad // | = pixel column edges - const int colBegin = (int)floor(c - deltad/2.0f); - const int colEnd = (int)ceil(c + deltad/2.0f); + const int colBegin = (int)floor(c - deltad); + const int colEnd = (int)ceil(c + deltad); - // TODO: Optimize volume edge checks + if (colBegin >= colCount || colEnd <= 0) + continue; int iVolumeIndex = row * colCount + colBegin; if (colBegin + 1 == colEnd) { if (colBegin >= 0 && colBegin < colCount) policy_weight(p, iRayIndex, iVolumeIndex, - deltad * lengthPerRow); + 2.0f * deltad * lengthPerRow); } else { - const float gapBegin = (c - deltad/2.0f) - (float32)colBegin; - const float gapEnd = (float32)colEnd - (c + deltad/2.0f); - float tot = 1.0f - gapBegin; - if (colBegin >= 0 && colBegin < colCount) { + + if (colBegin >= 0) { + const float gapBegin = (c - deltad) - (float32)colBegin; policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - gapBegin) * lengthPerRow); } - iVolumeIndex++; - for (int col = colBegin + 1; col + 1 < colEnd; ++col, ++iVolumeIndex) { - tot += 1.0f; - if (col >= 0 && col < colCount) { - policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow); - } + const int clippedMColBegin = std::max(colBegin + 1, 0); + const int clippedMColEnd = std::min(colEnd - 1, colCount); + iVolumeIndex = row * colCount + clippedMColBegin; + for (int col = clippedMColBegin; col < clippedMColEnd; ++col, ++iVolumeIndex) { + policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow); } - assert(iVolumeIndex == row * colCount + colEnd - 1); - tot += 1.0f - gapEnd; - if (colEnd > 0 && colEnd <= colCount) { + + iVolumeIndex = row * colCount + colEnd - 1; + if (colEnd <= colCount) { + const float gapEnd = (float32)colEnd - (c + deltad); policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - gapEnd) * lengthPerRow); } - assert(fabs(tot - deltad) < 0.0001); } } - } else if (!vertical && true) { + } else { const float32 RyOverRx = proj->fRayY/proj->fRayX; - // TODO: Determine det/pixel scaling factors const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY; - const float32 deltad = fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY); + const float32 deltad = 0.5f * fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY); // calculate r for col 0 float32 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY + 0.5f; @@ -171,43 +168,41 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro for (int col = 0; col < colCount; ++col, r+= deltar) { // vertical extent of ray in center of this column: - // [ r - deltad/2 , r + deltad/2 ] + // [ r - deltad , r + deltad ] - const int rowBegin = (int)floor(r - deltad/2.0f); - const int rowEnd = (int)ceil(r + deltad/2.0f); + const int rowBegin = (int)floor(r - deltad); + const int rowEnd = (int)ceil(r + deltad); - // TODO: Optimize volume edge checks + if (rowBegin >= rowCount || rowEnd <= 0) + continue; int iVolumeIndex = rowBegin * colCount + col; if (rowBegin + 1 == rowEnd) { if (rowBegin >= 0 && rowBegin < rowCount) policy_weight(p, iRayIndex, iVolumeIndex, - deltad * lengthPerCol); + 2.0f * deltad * lengthPerCol); } else { - const float gapBegin = (r - deltad/2.0f) - (float32)rowBegin; - const float gapEnd = (float32)rowEnd - (r + deltad/2.0f); - float tot = 1.0f - gapBegin; - - if (rowBegin >= 0 && rowBegin < rowCount) { + if (rowBegin >= 0) { + const float gapBegin = (r - deltad) - (float32)rowBegin; policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - gapBegin) * lengthPerCol); } - iVolumeIndex += colCount; - for (int row = rowBegin + 1; row + 1 < rowEnd; ++row, iVolumeIndex += colCount) { - tot += 1.0f; - if (row >= 0 && row < rowCount) { - policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol); - } + const int clippedMRowBegin = std::max(rowBegin + 1, 0); + const int clippedMRowEnd = std::min(rowEnd - 1, rowCount); + iVolumeIndex = clippedMRowBegin * colCount + col; + + for (int row = clippedMRowBegin; row < clippedMRowEnd; ++row, iVolumeIndex += colCount) { + policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol); } - assert(iVolumeIndex == (rowEnd - 1) * colCount + col); - tot += 1.0f - gapEnd; - if (rowEnd > 0 && rowEnd <= rowCount) { + + iVolumeIndex = (rowEnd - 1) * colCount + col; + if (rowEnd <= rowCount) { + const float gapEnd = (float32)rowEnd - (r + deltad); policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - gapEnd) * lengthPerCol); } - assert(fabs(tot - deltad) < 0.0001); } } -- cgit v1.2.3