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 --- src/ParallelBeamDistanceDrivenProjector2D.cpp | 174 ++++++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 src/ParallelBeamDistanceDrivenProjector2D.cpp (limited to 'src') diff --git a/src/ParallelBeamDistanceDrivenProjector2D.cpp b/src/ParallelBeamDistanceDrivenProjector2D.cpp new file mode 100644 index 0000000..a4ac450 --- /dev/null +++ b/src/ParallelBeamDistanceDrivenProjector2D.cpp @@ -0,0 +1,174 @@ +/* +----------------------------------------------------------------------- +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 . + +----------------------------------------------------------------------- +*/ + +#include "astra/ParallelBeamDistanceDrivenProjector2D.h" + +#include +#include + +#include "astra/DataProjectorPolicies.h" + +using namespace std; +using namespace astra; + +#include "astra/ParallelBeamDistanceDrivenProjector2D.inl" + +// type of the projector, needed to register with CProjectorFactory +std::string CParallelBeamDistanceDrivenProjector2D::type = "distance_driven"; + +//---------------------------------------------------------------------------------------- +// default constructor +CParallelBeamDistanceDrivenProjector2D::CParallelBeamDistanceDrivenProjector2D() +{ + _clear(); +} + + +//---------------------------------------------------------------------------------------- +// constructor +CParallelBeamDistanceDrivenProjector2D::CParallelBeamDistanceDrivenProjector2D(CParallelProjectionGeometry2D* _pProjectionGeometry, + CVolumeGeometry2D* _pReconstructionGeometry) + +{ + _clear(); + initialize(_pProjectionGeometry, _pReconstructionGeometry); +} + +//---------------------------------------------------------------------------------------- +// destructor +CParallelBeamDistanceDrivenProjector2D::~CParallelBeamDistanceDrivenProjector2D() +{ + clear(); +} + +//--------------------------------------------------------------------------------------- +// Clear - CParallelBeamDistanceDrivenProjector2D +void CParallelBeamDistanceDrivenProjector2D::_clear() +{ + CProjector2D::_clear(); + m_bIsInitialized = false; +} + +//--------------------------------------------------------------------------------------- +// Clear - Public +void CParallelBeamDistanceDrivenProjector2D::clear() +{ + CProjector2D::clear(); + m_bIsInitialized = false; +} + +//--------------------------------------------------------------------------------------- +// Check +bool CParallelBeamDistanceDrivenProjector2D::_check() +{ + // check base class + ASTRA_CONFIG_CHECK(CProjector2D::_check(), "ParallelBeamDistanceDrivenProjector2D", "Error in Projector2D initialization"); + + ASTRA_CONFIG_CHECK(dynamic_cast(m_pProjectionGeometry) || dynamic_cast(m_pProjectionGeometry), "ParallelBeamDistanceDrivenProjector2D", "Unsupported projection geometry"); + + ASTRA_CONFIG_CHECK(abs(m_pVolumeGeometry->getPixelLengthX() / m_pVolumeGeometry->getPixelLengthY()) - 1 < eps, "ParallelBeamDistanceDrivenProjector2D", "Pixel height must equal pixel width."); + + // success + return true; +} + +//--------------------------------------------------------------------------------------- +// Initialize, use a Config object +bool CParallelBeamDistanceDrivenProjector2D::initialize(const Config& _cfg) +{ + ASTRA_ASSERT(_cfg.self); + + // if already initialized, clear first + if (m_bIsInitialized) { + clear(); + } + + // initialization of parent class + if (!CProjector2D::initialize(_cfg)) { + return false; + } + + // success + m_bIsInitialized = _check(); + return m_bIsInitialized; +} + +//--------------------------------------------------------------------------------------- +// Initialize +bool CParallelBeamDistanceDrivenProjector2D::initialize(CParallelProjectionGeometry2D* _pProjectionGeometry, + CVolumeGeometry2D* _pVolumeGeometry) +{ + // if already initialized, clear first + if (m_bIsInitialized) { + clear(); + } + + // hardcopy geometries + m_pProjectionGeometry = _pProjectionGeometry->clone(); + m_pVolumeGeometry = _pVolumeGeometry->clone(); + + // success + m_bIsInitialized = _check(); + return m_bIsInitialized; +} + +//---------------------------------------------------------------------------------------- +// Get maximum amount of weights on a single ray +int CParallelBeamDistanceDrivenProjector2D::getProjectionWeightsCount(int _iProjectionIndex) +{ + int maxDim = max(m_pVolumeGeometry->getGridRowCount(), m_pVolumeGeometry->getGridColCount()); + int scale = m_pProjectionGeometry->getDetectorWidth() / min(m_pVolumeGeometry->getPixelLengthX(), m_pVolumeGeometry->getPixelLengthY()); + return maxDim * scale * 10 + 1; +} + +//---------------------------------------------------------------------------------------- +// Single Ray Weights +void CParallelBeamDistanceDrivenProjector2D::computeSingleRayWeights(int _iProjectionIndex, + int _iDetectorIndex, + SPixelWeight* _pWeightedPixels, + int _iMaxPixelCount, + int& _iStoredPixelCount) +{ + ASTRA_ASSERT(m_bIsInitialized); + StorePixelWeightsPolicy p(_pWeightedPixels, _iMaxPixelCount); + projectSingleRay(_iProjectionIndex, _iDetectorIndex, p); + _iStoredPixelCount = p.getStoredPixelCount(); +} + +//---------------------------------------------------------------------------------------- +// Splat a single point +std::vector CParallelBeamDistanceDrivenProjector2D::projectPoint(int _iRow, int _iCol) +{ + // TODO + assert(false); + + std::vector res; + + // return result vector + return res; + +} -- 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. --- src/FanFlatBeamLineKernelProjector2D.cpp | 8 ---- src/FanFlatBeamStripKernelProjector2D.cpp | 63 -------------------------- src/ParallelBeamBlobKernelProjector2D.cpp | 45 ------------------- src/ParallelBeamDistanceDrivenProjector2D.cpp | 19 ++------ src/ParallelBeamLineKernelProjector2D.cpp | 64 --------------------------- src/ParallelBeamLinearKernelProjector2D.cpp | 61 ------------------------- src/ParallelBeamStripKernelProjector2D.cpp | 61 ------------------------- src/SparseMatrixProjector2D.cpp | 41 ----------------- 8 files changed, 3 insertions(+), 359 deletions(-) (limited to 'src') diff --git a/src/FanFlatBeamLineKernelProjector2D.cpp b/src/FanFlatBeamLineKernelProjector2D.cpp index 6475ab2..a092504 100644 --- a/src/FanFlatBeamLineKernelProjector2D.cpp +++ b/src/FanFlatBeamLineKernelProjector2D.cpp @@ -159,14 +159,6 @@ void CFanFlatBeamLineKernelProjector2D::computeSingleRayWeights(int _iProjection _iStoredPixelCount = p.getStoredPixelCount(); } -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CFanFlatBeamLineKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - std::vector res; - return res; -} - //---------------------------------------------------------------------------------------- //Result is always in [-PI/2; PI/2] float32 CFanFlatBeamLineKernelProjector2D::angleBetweenVectors(float32 _fAX, float32 _fAY, float32 _fBX, float32 _fBY) diff --git a/src/FanFlatBeamStripKernelProjector2D.cpp b/src/FanFlatBeamStripKernelProjector2D.cpp index 7d25a55..c974b82 100644 --- a/src/FanFlatBeamStripKernelProjector2D.cpp +++ b/src/FanFlatBeamStripKernelProjector2D.cpp @@ -157,66 +157,3 @@ void CFanFlatBeamStripKernelProjector2D::computeSingleRayWeights(int _iProjectio projectSingleRay(_iProjectionIndex, _iDetectorIndex, p); _iStoredPixelCount = p.getStoredPixelCount(); } - -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CFanFlatBeamStripKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - //float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - //float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - //float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - //float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - //float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - //float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - //float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - //float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - - std::vector res; - //// loop projectors and detectors - //for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // // get projection angle - // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - // if (theta >= 7*PIdiv4) theta -= 2*PI; - // bool inverse = false; - // if (theta >= 3*PIdiv4) { - // theta -= PI; - // inverse = true; - // } - - // // calculate distance from the center of the voxel to the ray though the origin - // float32 tUL = xUL * cos(theta) + yUL * sin(theta); - // float32 tUR = xUR * cos(theta) + yUR * sin(theta); - // float32 tLL = xLL * cos(theta) + yLL * sin(theta); - // float32 tLR = xLR * cos(theta) + yLR * sin(theta); - // if (inverse) { - // tUL *= -1.0f; - // tUR *= -1.0f; - // tLL *= -1.0f; - // tLR *= -1.0f; - // } - // float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - // float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // // calculate the offset on the detectorarray (in indices) - // int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - // int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // // add affected detectors to the list - // for (int i = dmin; i <= dmax; ++i) { - // if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - // SDetector2D det; - // det.m_iAngleIndex = iProjection; - // det.m_iDetectorIndex = i; - // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - // res.push_back(det); - // } - // } - //} - - //// return result vector - return res; - -} - -//---------------------------------------------------------------------------------------- diff --git a/src/ParallelBeamBlobKernelProjector2D.cpp b/src/ParallelBeamBlobKernelProjector2D.cpp index da21024..0fb658f 100644 --- a/src/ParallelBeamBlobKernelProjector2D.cpp +++ b/src/ParallelBeamBlobKernelProjector2D.cpp @@ -212,48 +212,3 @@ void CParallelBeamBlobKernelProjector2D::computeSingleRayWeights(int _iProjectio projectSingleRay(_iProjectionIndex, _iDetectorIndex, p); _iStoredPixelCount = p.getStoredPixelCount(); } -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CParallelBeamBlobKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - // float32 x = m_pVolumeGeometry->pixelColToCenterX(_iCol); - // float32 y = m_pVolumeGeometry->pixelRowToCenterY(_iRow); - - // std::vector res; - // // loop projectors and detectors - // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // // get projection angle - // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - // if (theta >= 7*PIdiv4) theta -= 2*PI; - // bool inverse = false; - // if (theta >= 3*PIdiv4) { - // theta -= PI; - // inverse = true; - // } - - // // calculate distance from the center of the voxel to the ray though the origin - // float32 t = x * cos(theta) + y * sin(theta); - // if (inverse) t *= -1.0f; - - // // calculate the offset on the detectorarray (in indices) - // float32 d = m_pProjectionGeometry->detectorOffsetToIndexFloat(t); - // int dmin = (int)ceil(d - m_fBlobSize); - // int dmax = (int)floor(d + m_fBlobSize); - - // // add affected detectors to the list - // for (int i = dmin; i <= dmax; ++i) { - // if (d >= 0 && d < m_pProjectionGeometry->getDetectorCount()) { - // SDetector2D det; - // det.m_iAngleIndex = iProjection; - // det.m_iDetectorIndex = i; - // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - // res.push_back(det); - // } - // } - // } - - // // return result vector - // return res; - return std::vector(); -} diff --git a/src/ParallelBeamDistanceDrivenProjector2D.cpp b/src/ParallelBeamDistanceDrivenProjector2D.cpp index a4ac450..531bef9 100644 --- a/src/ParallelBeamDistanceDrivenProjector2D.cpp +++ b/src/ParallelBeamDistanceDrivenProjector2D.cpp @@ -32,8 +32,7 @@ along with the ASTRA Toolbox. If not, see . #include "astra/DataProjectorPolicies.h" -using namespace std; -using namespace astra; +namespace astra { #include "astra/ParallelBeamDistanceDrivenProjector2D.inl" @@ -140,8 +139,8 @@ bool CParallelBeamDistanceDrivenProjector2D::initialize(CParallelProjectionGeome // Get maximum amount of weights on a single ray int CParallelBeamDistanceDrivenProjector2D::getProjectionWeightsCount(int _iProjectionIndex) { - int maxDim = max(m_pVolumeGeometry->getGridRowCount(), m_pVolumeGeometry->getGridColCount()); - int scale = m_pProjectionGeometry->getDetectorWidth() / min(m_pVolumeGeometry->getPixelLengthX(), m_pVolumeGeometry->getPixelLengthY()); + int maxDim = std::max(m_pVolumeGeometry->getGridRowCount(), m_pVolumeGeometry->getGridColCount()); + int scale = m_pProjectionGeometry->getDetectorWidth() / std::min(m_pVolumeGeometry->getPixelLengthX(), m_pVolumeGeometry->getPixelLengthY()); return maxDim * scale * 10 + 1; } @@ -159,16 +158,4 @@ void CParallelBeamDistanceDrivenProjector2D::computeSingleRayWeights(int _iProje _iStoredPixelCount = p.getStoredPixelCount(); } -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CParallelBeamDistanceDrivenProjector2D::projectPoint(int _iRow, int _iCol) -{ - // TODO - assert(false); - - std::vector res; - - // return result vector - return res; - } diff --git a/src/ParallelBeamLineKernelProjector2D.cpp b/src/ParallelBeamLineKernelProjector2D.cpp index 3654fb0..1765df3 100644 --- a/src/ParallelBeamLineKernelProjector2D.cpp +++ b/src/ParallelBeamLineKernelProjector2D.cpp @@ -157,67 +157,3 @@ void CParallelBeamLineKernelProjector2D::computeSingleRayWeights(int _iProjectio _iStoredPixelCount = p.getStoredPixelCount(); } -//---------------------------------------------------------------------------------------- -// Project Point -std::vector CParallelBeamLineKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - std::vector res; - return res; - - // float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - - // std::vector res; - // // loop projectors and detectors - // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // // get projection angle - // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - // if (theta >= 7*PIdiv4) theta -= 2*PI; - // bool inverse = false; - // if (theta >= 3*PIdiv4) { - // theta -= PI; - // inverse = true; - // } - - // // calculate distance from the center of the voxel to the ray though the origin - // float32 tUL = xUL * cos(theta) + yUL * sin(theta); - // float32 tUR = xUR * cos(theta) + yUR * sin(theta); - // float32 tLL = xLL * cos(theta) + yLL * sin(theta); - // float32 tLR = xLR * cos(theta) + yLR * sin(theta); - // if (inverse) { - // tUL *= -1.0f; - // tUR *= -1.0f; - // tLL *= -1.0f; - // tLR *= -1.0f; - // } - // float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - // float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // // calculate the offset on the detectorarray (in indices) - // int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - // int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // // add affected detectors to the list - // for (int i = dmin; i <= dmax; ++i) { - // if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - // SDetector2D det; - // det.m_iAngleIndex = iProjection; - // det.m_iDetectorIndex = i; - // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - // res.push_back(det); - // } - // } - // } - - // // return result vector - // return res; -} - -//---------------------------------------------------------------------------------------- diff --git a/src/ParallelBeamLinearKernelProjector2D.cpp b/src/ParallelBeamLinearKernelProjector2D.cpp index 5ca1082..120b1a2 100644 --- a/src/ParallelBeamLinearKernelProjector2D.cpp +++ b/src/ParallelBeamLinearKernelProjector2D.cpp @@ -158,64 +158,3 @@ void CParallelBeamLinearKernelProjector2D::computeSingleRayWeights(int _iProject projectSingleRay(_iProjectionIndex, _iDetectorIndex, p); _iStoredPixelCount = p.getStoredPixelCount(); } - -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CParallelBeamLinearKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 1.5f; - float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 1.5f; - float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 1.5f; - float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 1.5f; - float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 1.5f; - float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 1.5f; - float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 1.5f; - float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 1.5f; - - std::vector res; - // loop projectors and detectors - for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // get projection angle - float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - if (theta >= 7*PIdiv4) theta -= 2*PI; - bool inverse = false; - if (theta >= 3*PIdiv4) { - theta -= PI; - inverse = true; - } - - // calculate distance from the center of the voxel to the ray though the origin - float32 tUL = xUL * cos(theta) + yUL * sin(theta); - float32 tUR = xUR * cos(theta) + yUR * sin(theta); - float32 tLL = xLL * cos(theta) + yLL * sin(theta); - float32 tLR = xLR * cos(theta) + yLR * sin(theta); - if (inverse) { - tUL *= -1.0f; - tUR *= -1.0f; - tLL *= -1.0f; - tLR *= -1.0f; - } - float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // calculate the offset on the detectorarray (in indices) - int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // add affected detectors to the list - for (int i = dmin; i <= dmax; ++i) { - if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - SDetector2D det; - det.m_iAngleIndex = iProjection; - det.m_iDetectorIndex = i; - det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - res.push_back(det); - } - } - } - - // return result vector - return res; - -} diff --git a/src/ParallelBeamStripKernelProjector2D.cpp b/src/ParallelBeamStripKernelProjector2D.cpp index 4b10b68..849168d 100644 --- a/src/ParallelBeamStripKernelProjector2D.cpp +++ b/src/ParallelBeamStripKernelProjector2D.cpp @@ -160,64 +160,3 @@ void CParallelBeamStripKernelProjector2D::computeSingleRayWeights(int _iProjecti } //---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CParallelBeamStripKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - // float32 xUL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yUL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xUR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yUR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) - m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xLL = m_pVolumeGeometry->pixelColToCenterX(_iCol) - m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yLL = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - // float32 xLR = m_pVolumeGeometry->pixelColToCenterX(_iCol) + m_pVolumeGeometry->getPixelLengthX() * 0.5f; - // float32 yLR = m_pVolumeGeometry->pixelRowToCenterY(_iRow) + m_pVolumeGeometry->getPixelLengthY() * 0.5f; - - std::vector res; - // // loop projectors and detectors - // for (int iProjection = 0; iProjection < m_pProjectionGeometry->getProjectionAngleCount(); ++iProjection) { - - // // get projection angle - // float32 theta = m_pProjectionGeometry->getProjectionAngle(iProjection); - // if (theta >= 7*PIdiv4) theta -= 2*PI; - // bool inverse = false; - // if (theta >= 3*PIdiv4) { - // theta -= PI; - // inverse = true; - // } - - // // calculate distance from the center of the voxel to the ray though the origin - // float32 tUL = xUL * cos(theta) + yUL * sin(theta); - // float32 tUR = xUR * cos(theta) + yUR * sin(theta); - // float32 tLL = xLL * cos(theta) + yLL * sin(theta); - // float32 tLR = xLR * cos(theta) + yLR * sin(theta); - // if (inverse) { - // tUL *= -1.0f; - // tUR *= -1.0f; - // tLL *= -1.0f; - // tLR *= -1.0f; - // } - // float32 tMin = min(tUL, min(tUR, min(tLL,tLR))); - // float32 tMax = max(tUL, max(tUR, max(tLL,tLR))); - - // // calculate the offset on the detectorarray (in indices) - // int dmin = (int)floor(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMin)); - // int dmax = (int)ceil(m_pProjectionGeometry->detectorOffsetToIndexFloat(tMax)); - - // // add affected detectors to the list - // for (int i = dmin; i <= dmax; ++i) { - // if (i >= 0 && i < m_pProjectionGeometry->getDetectorCount()) { - // SDetector2D det; - // det.m_iAngleIndex = iProjection; - // det.m_iDetectorIndex = i; - // det.m_iIndex = iProjection * getProjectionGeometry()->getDetectorCount() + i; - // res.push_back(det); - // } - // } - // } - - // return result vector - return res; - -} - -//---------------------------------------------------------------------------------------- diff --git a/src/SparseMatrixProjector2D.cpp b/src/SparseMatrixProjector2D.cpp index 6fded6a..46302f3 100644 --- a/src/SparseMatrixProjector2D.cpp +++ b/src/SparseMatrixProjector2D.cpp @@ -174,44 +174,3 @@ void CSparseMatrixProjector2D::computeSingleRayWeights(int _iProjectionIndex, _iStoredPixelCount = p.getStoredPixelCount(); } -//---------------------------------------------------------------------------------------- -// Splat a single point -std::vector CSparseMatrixProjector2D::projectPoint(int _iRow, int _iCol) -{ - unsigned int iVolumeIndex = _iCol * m_pVolumeGeometry->getGridRowCount() + _iRow; - - // NOTE: This is very slow currently because we don't have the - // sparse matrix stored in an appropriate form for this function. - std::vector ret; - - const CSparseMatrix* pMatrix = dynamic_cast(m_pProjectionGeometry)->getMatrix(); - - for (int iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) - { - for (int iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) - { - int iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector; - const unsigned int* piColIndices; - const float32* pfValues; - unsigned int iSize; - - pMatrix->getRowData(iRayIndex, iSize, pfValues, piColIndices); - - for (unsigned int i = 0; i < iSize; ++i) { - if (piColIndices[i] == iVolumeIndex) { - SDetector2D s; - s.m_iIndex = iRayIndex; - s.m_iAngleIndex = iAngle; - s.m_iDetectorIndex = iDetector; - ret.push_back(s); - break; - } else if (piColIndices[i] > iVolumeIndex) { - break; - } - } - } - } - return ret; -} - -//---------------------------------------------------------------------------------------- -- cgit v1.2.3