diff options
26 files changed, 764 insertions, 716 deletions
diff --git a/build/linux/Makefile.in b/build/linux/Makefile.in index 7c4089a..078a1a2 100644 --- a/build/linux/Makefile.in +++ b/build/linux/Makefile.in @@ -168,6 +168,7 @@ BASE_OBJECTS=\ src/Globals.lo \ src/Logging.lo \ src/ParallelBeamBlobKernelProjector2D.lo \ + src/ParallelBeamDistanceDrivenProjector2D.lo \ src/ParallelBeamLinearKernelProjector2D.lo \ src/ParallelBeamLineKernelProjector2D.lo \ src/ParallelBeamStripKernelProjector2D.lo \ diff --git a/build/msvc/gen.py b/build/msvc/gen.py index 42d0e0e..2200ecc 100644 --- a/build/msvc/gen.py +++ b/build/msvc/gen.py @@ -249,6 +249,7 @@ P_astra["filters"]["Projectors\\source"] = [ "src\\FanFlatBeamLineKernelProjector2D.cpp", "src\\FanFlatBeamStripKernelProjector2D.cpp", "src\\ParallelBeamBlobKernelProjector2D.cpp", +"src\\ParallelBeamDistanceDrivenProjector2D.cpp", "src\\ParallelBeamLinearKernelProjector2D.cpp", "src\\ParallelBeamLineKernelProjector2D.cpp", "src\\ParallelBeamStripKernelProjector2D.cpp", @@ -393,6 +394,7 @@ P_astra["filters"]["Projectors\\headers"] = [ "include\\astra\\FanFlatBeamLineKernelProjector2D.h", "include\\astra\\FanFlatBeamStripKernelProjector2D.h", "include\\astra\\ParallelBeamBlobKernelProjector2D.h", +"include\\astra\\ParallelBeamDistanceDrivenProjector2D.h", "include\\astra\\ParallelBeamLinearKernelProjector2D.h", "include\\astra\\ParallelBeamLineKernelProjector2D.h", "include\\astra\\ParallelBeamStripKernelProjector2D.h", diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.h b/include/astra/FanFlatBeamLineKernelProjector2D.h index 1bfaf07..4492bdf 100644 --- a/include/astra/FanFlatBeamLineKernelProjector2D.h +++ b/include/astra/FanFlatBeamLineKernelProjector2D.h @@ -134,14 +134,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.h b/include/astra/FanFlatBeamStripKernelProjector2D.h index a6a303c..4942f6c 100644 --- a/include/astra/FanFlatBeamStripKernelProjector2D.h +++ b/include/astra/FanFlatBeamStripKernelProjector2D.h @@ -132,14 +132,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.h b/include/astra/ParallelBeamBlobKernelProjector2D.h index 46b0b52..ffae707 100644 --- a/include/astra/ParallelBeamBlobKernelProjector2D.h +++ b/include/astra/ParallelBeamBlobKernelProjector2D.h @@ -160,14 +160,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.h b/include/astra/ParallelBeamDistanceDrivenProjector2D.h new file mode 100644 index 0000000..6d1d633 --- /dev/null +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.h @@ -0,0 +1,197 @@ +/* +----------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#ifndef _INC_ASTRA_PARALLELDISTANCEDRIVENPROJECTOR +#define _INC_ASTRA_PARALLELDISTANCEDRIVENPROJECTOR + +#include "ParallelProjectionGeometry2D.h" +#include "ParallelVecProjectionGeometry2D.h" +#include "Float32Data2D.h" +#include "Projector2D.h" + +#include "Float32ProjectionData2D.h" +#include "Float32VolumeData2D.h" + +namespace astra +{ + +/** This class implements a "distance driven" two-dimensional projector. + * + * Reference: + * De Man, Bruno, and Samit Basu. “Distance-Driven Projection and Backprojection in Three Dimensions.” Physics in Medicine and Biology 49, no. 11 (2004): 2463. + * + * \par XML Configuration + * \astra_xml_item{ProjectionGeometry, xml node, The geometry of the projection.} + * \astra_xml_item{VolumeGeometry, xml node, The geometry of the volume.} + * + * \par MATLAB example + * \astra_code{ + * cfg = astra_struct('distance_driven');\n + * cfg.ProjectionGeometry = proj_geom;\n + * cfg.VolumeGeometry = vol_geom;\n + * proj_id = astra_mex_projector('create'\, cfg);\n + * } + */ +class _AstraExport CParallelBeamDistanceDrivenProjector2D : public CProjector2D { + +protected: + + /** Initial clearing. Only to be used by constructors. + */ + virtual void _clear(); + + /** Check the values of this object. If everything is ok, the object can be set to the initialized state. + * The following statements are then guaranteed to hold: + * - no NULL pointers + * - all sub-objects are initialized properly + * - blobvalues are ok + */ + virtual bool _check(); + +public: + + // type of the projector, needed to register with CProjectorFactory + static std::string type; + + /** Default constructor. + */ + CParallelBeamDistanceDrivenProjector2D(); + + /** Constructor. + * + * @param _pProjectionGeometry Information class about the geometry of the projection. Will be HARDCOPIED. + * @param _pReconstructionGeometry Information class about the geometry of the reconstruction volume. Will be HARDCOPIED. + */ + CParallelBeamDistanceDrivenProjector2D(CParallelProjectionGeometry2D* _pProjectionGeometry, + CVolumeGeometry2D* _pReconstructionGeometry); + + /** Destructor, is virtual to show that we are aware subclass destructor are called. + */ + ~CParallelBeamDistanceDrivenProjector2D(); + + /** Initialize the projector with a config object. + * + * @param _cfg Configuration Object + * @return initialization successful? + */ + virtual bool initialize(const Config& _cfg); + + /** Initialize the projector. + * + * @param _pProjectionGeometry Information class about the geometry of the projection. Will be HARDCOPIED. + * @param _pReconstructionGeometry Information class about the geometry of the reconstruction volume. Will be HARDCOPIED. + * @return initialization successful? + */ + virtual bool initialize(CParallelProjectionGeometry2D* _pProjectionGeometry, + CVolumeGeometry2D* _pVolumeGeometry); + + /** Clear this class. + */ + virtual void clear(); + + /** Returns the number of weights required for storage of all weights of one projection. + * + * @param _iProjectionIndex Index of the projection (zero-based). + * @return Size of buffer (given in SPixelWeight elements) needed to store weighted pixels. + */ + virtual int getProjectionWeightsCount(int _iProjectionIndex); + + /** Compute the pixel weights for a single ray, from the source to a detector pixel. + * + * @param _iProjectionIndex Index of the projection + * @param _iDetectorIndex Index of the detector pixel + * @param _pWeightedPixels Pointer to a pre-allocated array, consisting of _iMaxPixelCount elements + * of type SPixelWeight. On return, this array contains a list of the index + * and weight for all pixels on the ray. + * @param _iMaxPixelCount Maximum number of pixels (and corresponding weights) that can be stored in _pWeightedPixels. + * This number MUST be greater than the total number of pixels on the ray. + * @param _iStoredPixelCount On return, this variable contains the total number of pixels on the + * ray (that have been stored in the list _pWeightedPixels). + */ + virtual void computeSingleRayWeights(int _iProjectionIndex, + int _iDetectorIndex, + SPixelWeight* _pWeightedPixels, + int _iMaxPixelCount, + int& _iStoredPixelCount); + + + /** Policy-based projection of all rays. This function will calculate each non-zero projection + * weight and use this value for a task provided by the policy object. + * + * @param _policy Policy object. Should contain prior, addWeight and posterior function. + */ + template <typename Policy> + void project(Policy& _policy); + + /** Policy-based projection of all rays of a single projection. This function will calculate + * each non-zero projection weight and use this value for a task provided by the policy object. + * + * @param _iProjection Which projection should be projected? + * @param _policy Policy object. Should contain prior, addWeight and posterior function. + */ + template <typename Policy> + void projectSingleProjection(int _iProjection, Policy& _policy); + + /** Policy-based projection of a single ray. This function will calculate each non-zero + * projection weight and use this value for a task provided by the policy object. + * + * @param _iProjection Which projection should be projected? + * @param _iDetector Which detector should be projected? + * @param _policy Policy object. Should contain prior, addWeight and posterior function. + */ + template <typename Policy> + void projectSingleRay(int _iProjection, int _iDetector, Policy& _policy); + + /** Return the type of this projector. + * + * @return identification type of this projector + */ + 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 CParallelBeamDistanceDrivenProjector2D::getType() +{ + return type; +} + + +} // namespace astra + + +#endif + diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl new file mode 100644 index 0000000..aedcee9 --- /dev/null +++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl @@ -0,0 +1,220 @@ +/* +----------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#define policy_weight(p,rayindex,volindex,weight) do { if (p.pixelPrior(volindex)) { p.addWeight(rayindex, volindex, weight); p.pixelPosterior(volindex); } } while (false) + +template <typename Policy> +void CParallelBeamDistanceDrivenProjector2D::project(Policy& p) +{ + projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(), + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamDistanceDrivenProjector2D::projectSingleProjection(int _iProjection, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + 0, m_pProjectionGeometry->getDetectorCount(), p); +} + +template <typename Policy> +void CParallelBeamDistanceDrivenProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p) +{ + projectBlock_internal(_iProjection, _iProjection + 1, + _iDetector, _iDetector + 1, p); +} + + + + + +template <typename Policy> +void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p) +{ + // get vector geometry + const CParallelVecProjectionGeometry2D* pVecProjectionGeometry; + if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) { + pVecProjectionGeometry = dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry(); + } else { + pVecProjectionGeometry = dynamic_cast<CParallelVecProjectionGeometry2D*>(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]; + + 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) { + + const float32 RxOverRy = proj->fRayX/proj->fRayY; + const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 deltac = -pixelLengthY * 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; + + // loop rows + for (int row = 0; row < rowCount; ++row, c+= deltac) { + + // horizontal extent of ray in center of this row: + // [ c - deltad , c + deltad ] + + // |-gapBegin-*---|------|----*-gapEnd-| + // * = ray extent intercepts; c - deltad and c + deltad + // | = pixel column edges + + const int colBegin = (int)floor(c - deltad); + const int colEnd = (int)ceil(c + deltad); + + 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, + 2.0f * deltad * lengthPerRow); + } else { + + if (colBegin >= 0) { + const float gapBegin = (c - deltad) - (float32)colBegin; + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapBegin) * 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); + } + + iVolumeIndex = row * colCount + colEnd - 1; + if (colEnd <= colCount) { + const float gapEnd = (float32)colEnd - (c + deltad); + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapEnd) * lengthPerRow); + } + } + + } + + } else { + + const float32 RyOverRx = proj->fRayY/proj->fRayX; + const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY(); + const float32 deltar = -pixelLengthX * 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; + + // loop columns + for (int col = 0; col < colCount; ++col, r+= deltar) { + + // vertical extent of ray in center of this column: + // [ r - deltad , r + deltad ] + + const int rowBegin = (int)floor(r - deltad); + const int rowEnd = (int)ceil(r + deltad); + + 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, + 2.0f * deltad * lengthPerCol); + } else { + if (rowBegin >= 0) { + const float gapBegin = (r - deltad) - (float32)rowBegin; + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapBegin) * 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); + } + + iVolumeIndex = (rowEnd - 1) * colCount + col; + if (rowEnd <= rowCount) { + const float gapEnd = (float32)rowEnd - (r + deltad); + policy_weight(p, iRayIndex, iVolumeIndex, + (1.0f - gapEnd) * lengthPerCol); + } + } + + } + + } + + // POLICY: RAY POSTERIOR + p.rayPosterior(iRayIndex); + + } + } + + if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry)) + delete pVecProjectionGeometry; +} diff --git a/include/astra/ParallelBeamLineKernelProjector2D.h b/include/astra/ParallelBeamLineKernelProjector2D.h index f709e1d..4238919 100644 --- a/include/astra/ParallelBeamLineKernelProjector2D.h +++ b/include/astra/ParallelBeamLineKernelProjector2D.h @@ -132,14 +132,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h index a09dcd1..67d940e 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.h +++ b/include/astra/ParallelBeamLinearKernelProjector2D.h @@ -135,14 +135,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl index 485eef6..10d4892 100644 --- a/include/astra/ParallelBeamLinearKernelProjector2D.inl +++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl @@ -156,7 +156,7 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY); - bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); + const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY); if (vertical) { RxOverRy = proj->fRayX/proj->fRayY; lengthPerRow = detSize * m_pVolumeGeometry->getPixelLengthX() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY); diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h index 153a454..597f8dc 100644 --- a/include/astra/ParallelBeamStripKernelProjector2D.h +++ b/include/astra/ParallelBeamStripKernelProjector2D.h @@ -133,14 +133,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/include/astra/Projector2D.h b/include/astra/Projector2D.h index 8e0b7a8..9d9130f 100644 --- a/include/astra/Projector2D.h +++ b/include/astra/Projector2D.h @@ -149,14 +149,6 @@ public: SPixelWeight* _pfWeightedPixels, int* _piRayStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol) = 0; - /** Returns the number of weights required for storage of all weights of one projection ray. * * @param _iProjectionIndex Index of the projection (zero-based). diff --git a/include/astra/Projector2DImpl.inl b/include/astra/Projector2DImpl.inl index b7edb1f..1308141 100644 --- a/include/astra/Projector2DImpl.inl +++ b/include/astra/Projector2DImpl.inl @@ -26,6 +26,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. */ +#include "ParallelBeamDistanceDrivenProjector2D.inl" #include "ParallelBeamLinearKernelProjector2D.inl" #include "ParallelBeamLineKernelProjector2D.inl" #include "ParallelBeamStripKernelProjector2D.inl" diff --git a/include/astra/ProjectorTypelist.h b/include/astra/ProjectorTypelist.h index 1a65aad..063693d 100644 --- a/include/astra/ProjectorTypelist.h +++ b/include/astra/ProjectorTypelist.h @@ -35,6 +35,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>. #include "Projector2D.h" #include "ParallelBeamLineKernelProjector2D.h" #include "ParallelBeamLinearKernelProjector2D.h" +#include "ParallelBeamDistanceDrivenProjector2D.h" #include "ParallelBeamBlobKernelProjector2D.h" #include "ParallelBeamStripKernelProjector2D.h" #include "SparseMatrixProjector2D.h" @@ -46,9 +47,10 @@ namespace astra{ #ifdef ASTRA_CUDA - typedef TYPELIST_8( + typedef TYPELIST_9( CFanFlatBeamLineKernelProjector2D, CFanFlatBeamStripKernelProjector2D, + CParallelBeamDistanceDrivenProjector2D, CParallelBeamLinearKernelProjector2D, CParallelBeamLineKernelProjector2D, CParallelBeamBlobKernelProjector2D, @@ -59,9 +61,10 @@ namespace astra{ #else - typedef TYPELIST_7( + typedef TYPELIST_8( CFanFlatBeamLineKernelProjector2D, CFanFlatBeamStripKernelProjector2D, + CParallelBeamDistanceDrivenProjector2D, CParallelBeamLinearKernelProjector2D, CParallelBeamLineKernelProjector2D, CParallelBeamBlobKernelProjector2D, diff --git a/include/astra/SparseMatrixProjector2D.h b/include/astra/SparseMatrixProjector2D.h index db8c4e2..1fc44c5 100644 --- a/include/astra/SparseMatrixProjector2D.h +++ b/include/astra/SparseMatrixProjector2D.h @@ -133,14 +133,6 @@ public: int _iMaxPixelCount, int& _iStoredPixelCount); - /** Create a list of detectors that are influenced by point [_iRow, _iCol]. - * - * @param _iRow row of the point - * @param _iCol column of the point - * @return list of SDetector2D structs - */ - virtual std::vector<SDetector2D> projectPoint(int _iRow, int _iCol); - /** Policy-based projection of all rays. This function will calculate each non-zero projection * weight and use this value for a task provided by the policy object. * diff --git a/matlab/mex/astra_mex_projector3d_c.cpp b/matlab/mex/astra_mex_projector3d_c.cpp index 5581d6c..e760e3d 100644 --- a/matlab/mex/astra_mex_projector3d_c.cpp +++ b/matlab/mex/astra_mex_projector3d_c.cpp @@ -179,175 +179,6 @@ void astra_mex_projector3d_get_volume_geometry(int nlhs, mxArray* plhs[], int nr } //----------------------------------------------------------------------------------------- -/** -* [weights] = astra_mex_projector3d('weights_single_ray', pid, projection_index, detector_index); -*/ -void astra_mex_projector_weights_single_ray(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) -{ - //// step1: get input - //if (nrhs < 4) { - // mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); - // return; - //} - //int iPid = (int)(mxGetScalar(prhs[1])); - //int iProjectionIndex = (int)(mxGetScalar(prhs[2])); - //int iDetectorIndex = (int)(mxGetScalar(prhs[3])); - - //// step2: get projector - //CProjector3D* pProjector; - //if (!(pProjector = CProjector3DManager::getSingleton().get(iPid))) { - // mexErrMsgTxt("Projector not found.\n"); - // return; - //} - // - //// step3: create output vars - //int iStoredPixelCount; - //int iMaxPixelCount = pProjector->getProjectionWeightsCount(iProjectionIndex); - //SWeightedPixel* pPixelsWeights = new SWeightedPixel3D[iMaxPixelCount]; - // - //// step4: perform operation - //pProjector->computeSingleRayWeights(iProjectionIndex, - // iDetectorIndex, - // pPixelsWeights, - // iMaxPixelCount, - // iStoredPixelCount); - - //// step5: return output - //if (1 <= nlhs) { - // mwSize dims[2]; - // dims[0] = iStoredPixelCount; - // dims[1] = 2; - - // plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL); - // double* out = mxGetPr(plhs[0]); - - // for (int col = 0; col < iStoredPixelCount; col++) { - // out[col] = pPixelsWeights[col].m_iIndex; - // out[iStoredPixelCount+col] = pPixelsWeights[col].m_fWeight; - // //cout << pPixelsWeights[col].m_iIndex << " " << pPixelsWeights[col].m_fWeight <<endl; - // } - //} - // - //// garbage collection - //delete[] pPixelsWeights; -} - -////----------------------------------------------------------------------------------------- -///** -//* [weights] = astra_mex_projector('weights_projection', pid, projection_index); -//*/ -//void astra_mex_projector_weights_projection(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) -//{ -// // step1: get input -// if (nrhs < 3) { -// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); -// return; -// } -// int iPid = (int)(mxGetScalar(prhs[1])); -// int iProjectionIndex = (int)(mxGetScalar(prhs[2])); -// -// // step2: get projector -// CProjector2D* pProjector; -// if (!(pProjector = CProjectorManager::getSingleton().get(iPid))) { -// mexErrMsgTxt("Projector not found.\n"); -// return; -// } -// -// // step3: create output vars -// SWeightedPixel2D* pPixelsWheights = new SWeightedPixel2D[pProjector->getProjectionWeightsCount(iProjectionIndex)]; -// int* piRayStoredPixelCount = new int[pProjector->getProjectionGeometry()->getDetectorCount()]; -// -// // step4: perform operation -// pProjector->computeProjectionRayWeights(iProjectionIndex, pPixelsWheights, piRayStoredPixelCount); -// -// // step5: return output -// if (1 <= nlhs) { -// // get basic values -// int iMatrixSize = pProjector->getVolumeGeometry()->getWindowLengthX() * -// pProjector->getVolumeGeometry()->getWindowLengthY(); -// int iDetectorCount = pProjector->getProjectionGeometry()->getDetectorCount(); -// int iTotalStoredPixelCount = 0; -// for (int i = 0; i < iDetectorCount; i++) { -// iTotalStoredPixelCount += piRayStoredPixelCount[i]; -// } -// -// // create matlab sparse matrix -// plhs[0] = mxCreateSparse(iMatrixSize, // number of rows (#pixels) -// iDetectorCount, // number of columns (#detectors) -// iTotalStoredPixelCount, // number of non-zero elements -// mxREAL); // element type -// double* values = mxGetPr(plhs[0]); -// mwIndex* rows = mxGetIr(plhs[0]); -// mwIndex* cols = mxGetJc(plhs[0]); -// -// int currentBase = 0; -// int currentIndex = 0; -// for (int i = 0; i < iDetectorCount; i++) { -// for (int j = 0; j < piRayStoredPixelCount[i]; j++) { -// values[currentIndex + j] = pPixelsWheights[currentBase + j].m_fWeight; -// rows[currentIndex + j] = pPixelsWheights[currentBase + j].m_iIndex; -// } -// -// currentBase += pProjector->getProjectionWeightsCount(iProjectionIndex) / pProjector->getProjectionGeometry()->getDetectorCount(); -// currentIndex += piRayStoredPixelCount[i]; -// } -// cols[0] = piRayStoredPixelCount[0]; -// for (int j = 1; j < iDetectorCount; j++) { -// cols[j] = cols[j-1] + piRayStoredPixelCount[j]; -// } -// cols[iDetectorCount] = iTotalStoredPixelCount; -// } -// -//} -// -////----------------------------------------------------------------------------------------- -///** -//* output = astra_mex_projector('splat', pid, x, y); -//*/ -//void astra_mex_projector_splat(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) -//{ -// // step1: get input -// if (nrhs < 4) { -// mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); -// return; -// } -// int iPid = (int)(mxGetScalar(prhs[1])); -// int iX = (int)(mxGetScalar(prhs[2])); -// int iY = (int)(mxGetScalar(prhs[3])); -// -// // step2: get projector -// CProjector2D* pProjector; -// if (!(pProjector = CProjectorManager::getSingleton().get(iPid))) { -// mexErrMsgTxt("Projector not found.\n"); -// return; -// } -// -// // step3: perform action -// vector<SDetector2D> detinfo = pProjector->projectPoint(iX, iY); -// -// // step4: output -// if (nlhs <= 1) { -// plhs[0] = mxCreateDoubleMatrix(detinfo.size(), // # rows -// 2, // # cols -// mxREAL); // datatype 32-bits -// double* out = mxGetPr(plhs[0]); -// -// // fill up output -// int i = 0; -// for (int x = 0; x < detinfo.size() ; x++) { -// out[i] = detinfo[x].m_iAngleIndex; -// i++; -// } -// for (int x = 0; x < detinfo.size() ; x++) { -// out[i] = detinfo[x].m_iDetectorIndex; -// i++; -// } -// } -// -// -//} - -//----------------------------------------------------------------------------------------- /** result = astra_mex_projector3d('is_cuda', id); * * Return is the specified projector is a cuda projector. @@ -395,8 +226,7 @@ static void printHelp() { mexPrintf("Please specify a mode of operation.\n"); mexPrintf("Valid modes: create, delete, clear, get_projection_geometry,\n"); - mexPrintf(" get_volume_geometry, weights_single_ray, weights_projection\n"); - mexPrintf(" splat, is_cuda\n"); + mexPrintf(" get_volume_geometry, is_cuda\n"); } //----------------------------------------------------------------------------------------- @@ -428,12 +258,6 @@ void mexFunction(int nlhs, mxArray* plhs[], astra_mex_projector3d_get_projection_geometry(nlhs, plhs, nrhs, prhs); } else if (sMode == "get_volume_geometry") { astra_mex_projector3d_get_volume_geometry(nlhs, plhs, nrhs, prhs); - } else if (sMode == "weights_single_ray") { - astra_mex_projector_weights_single_ray(nlhs, plhs, nrhs, prhs); - //} else if (sMode == "weights_projection") { - // astra_mex_projector_weights_projection(nlhs, plhs, nrhs, prhs); - //} else if (sMode == "splat") { - // astra_mex_projector_splat(nlhs, plhs, nrhs, prhs); } else if (sMode == "is_cuda") { astra_mex_projector3d_is_cuda(nlhs, plhs, nrhs, prhs); } else if (sMode == "help") { diff --git a/matlab/mex/astra_mex_projector_c.cpp b/matlab/mex/astra_mex_projector_c.cpp index 27c4b20..16da126 100644 --- a/matlab/mex/astra_mex_projector_c.cpp +++ b/matlab/mex/astra_mex_projector_c.cpp @@ -337,59 +337,6 @@ void astra_mex_projector_weights_projection(int nlhs, mxArray* plhs[], int nrhs, } //----------------------------------------------------------------------------------------- -/** hit_detectors = astra_mex_projector('splat', id, col, row); - * - * Create a list of detector indices which have a nonzero contribution to the projection matrix for a pixel [row,col]. - * id: identifier of the projector object as stored in the astra-library. - * col: column of the pixel - * row: row of the pixel - * hit_detectors: list of detector indices [angle_index, det_index] that are hit - */ -void astra_mex_projector_splat(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) -{ - // step1: get input - if (nrhs < 4) { - mexErrMsgTxt("Not enough arguments. See the help document for a detailed argument list. \n"); - return; - } - int iPid = (int)(mxGetScalar(prhs[1])); - int iX = (int)(mxGetScalar(prhs[2])); - int iY = (int)(mxGetScalar(prhs[3])); - - // step2: get projector - CProjector2D* pProjector = CProjector2DManager::getSingleton().get(iPid); - if (!pProjector || !pProjector->isInitialized()) { - mexErrMsgTxt("Projector not initialized.\n"); - return; - } - - // step3: perform action - vector<SDetector2D> detinfo = pProjector->projectPoint(iX, iY); - - // step4: output - if (nlhs <= 1) { - plhs[0] = mxCreateDoubleMatrix(detinfo.size(), // # rows - 2, // # cols - mxREAL); // datatype 32-bits - double* out = mxGetPr(plhs[0]); - - // fill up output - int i = 0; - int x; - for (x = 0; x < detinfo.size(); ++x) { - out[i] = detinfo[x].m_iAngleIndex; - ++i; - } - for (x = 0; x < detinfo.size(); ++x) { - out[i] = detinfo[x].m_iDetectorIndex; - ++i; - } - } - - -} - -//----------------------------------------------------------------------------------------- /** matrix_id = astra_mex_projector('matrix', id); * * Create an explicit projection matrix for this projector. @@ -467,7 +414,7 @@ static void printHelp() mexPrintf("Please specify a mode of operation.\n"); mexPrintf("Valid modes: create, delete, clear, info, projection_geometry,\n"); mexPrintf(" volume_geometry, weights_single_ray, weights_projection\n"); - mexPrintf(" splat, matrix, is_cuda\n"); + mexPrintf(" matrix, is_cuda\n"); } @@ -506,8 +453,6 @@ void mexFunction(int nlhs, mxArray* plhs[], astra_mex_projector_weights_single_ray(nlhs, plhs, nrhs, prhs); } else if (sMode == "weights_projection") { astra_mex_projector_weights_projection(nlhs, plhs, nrhs, prhs); - } else if (sMode == "splat") { - astra_mex_projector_splat(nlhs, plhs, nrhs, prhs); } else if (sMode == "matrix") { astra_mex_projector_matrix(nlhs, plhs, nrhs, prhs); } else if (sMode == "is_cuda") { diff --git a/src/FanFlatBeamLineKernelProjector2D.cpp b/src/FanFlatBeamLineKernelProjector2D.cpp index 6475ab2..a092504 100644 --- a/src/FanFlatBeamLineKernelProjector2D.cpp +++ b/src/FanFlatBeamLineKernelProjector2D.cpp @@ -160,14 +160,6 @@ void CFanFlatBeamLineKernelProjector2D::computeSingleRayWeights(int _iProjection } //---------------------------------------------------------------------------------------- -// Splat a single point -std::vector<SDetector2D> CFanFlatBeamLineKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - std::vector<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> CParallelBeamBlobKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - // float32 x = m_pVolumeGeometry->pixelColToCenterX(_iCol); - // float32 y = m_pVolumeGeometry->pixelRowToCenterY(_iRow); - - // std::vector<SDetector2D> 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<SDetector2D>(); -} diff --git a/src/ParallelBeamDistanceDrivenProjector2D.cpp b/src/ParallelBeamDistanceDrivenProjector2D.cpp new file mode 100644 index 0000000..531bef9 --- /dev/null +++ b/src/ParallelBeamDistanceDrivenProjector2D.cpp @@ -0,0 +1,161 @@ +/* +----------------------------------------------------------------------- +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 <http://www.gnu.org/licenses/>. + +----------------------------------------------------------------------- +*/ + +#include "astra/ParallelBeamDistanceDrivenProjector2D.h" + +#include <cmath> +#include <algorithm> + +#include "astra/DataProjectorPolicies.h" + +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<CParallelProjectionGeometry2D*>(m_pProjectionGeometry) || dynamic_cast<CParallelVecProjectionGeometry2D*>(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 = 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; +} + +//---------------------------------------------------------------------------------------- +// 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(); +} + +} 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<SDetector2D> CParallelBeamLineKernelProjector2D::projectPoint(int _iRow, int _iCol) -{ - std::vector<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> 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<SDetector2D> ret; - - const CSparseMatrix* pMatrix = dynamic_cast<CSparseMatrixProjectionGeometry2D*>(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; -} - -//---------------------------------------------------------------------------------------- diff --git a/tests/python/test_line2d.py b/tests/python/test_line2d.py index de68033..7b46458 100644 --- a/tests/python/test_line2d.py +++ b/tests/python/test_line2d.py @@ -58,6 +58,53 @@ def intersect_line_rectangle_interval(src, det, xmin, xmax, ymin, ymax, f): c = intersect_line_rectangle_feather(src, det, xmin, xmax, ymin, ymax, f) return (a,b,c) + +# x-coord of intersection of the line (src, det) with the horizontal line at y +def intersect_line_horizontal(src, det, y): + EPS = 1e-5 + + if np.abs(src[1] - det[1]) < EPS: + return np.nan + + t = (y - src[1]) / (det[1] - src[1]) + + return src[0] + t * (det[0] - src[0]) + + +# length of the intersection of the strip with boundaries edge1, edge2 with the horizontal +# segment at y, with horizontal extent x_seg +def intersect_ray_horizontal_segment(edge1, edge2, y, x_seg): + e1 = intersect_line_horizontal(edge1[0], edge1[1], y) + e2 = intersect_line_horizontal(edge2[0], edge2[1], y) + + if not (np.isfinite(e1) and np.isfinite(e2)): + return np.nan + + (e1, e2) = np.sort([e1, e2]) + (x1, x2) = np.sort(x_seg) + l = np.max([e1, x1]) + r = np.min([e2, x2]) + #print(edge1, edge2, y, x_seg, r-l) + return np.max([r-l, 0.0]) + +def intersect_ray_vertical_segment(edge1, edge2, x, y_seg): + # mirror edge1 and edge2 + edge1 = [ (a[1], a[0]) for a in edge1 ] + edge2 = [ (a[1], a[0]) for a in edge2 ] + return intersect_ray_horizontal_segment(edge1, edge2, x, y_seg) + + + + + + +# LINE GENERATORS +# --------------- +# +# Per ray these yield three lines, at respectively the center and two edges of the detector pixel. +# Each line is given by two points on the line. +# ( ( (p0x, p0y), (q0x, q0y) ), ( (p1x, p1y), (q1x, q1y) ), ( (p2x, p2y), (q2x, q2y) ) ) + def gen_lines_fanflat(proj_geom): angles = proj_geom['ProjectionAngles'] for theta in angles: @@ -76,7 +123,9 @@ def gen_lines_fanflat(proj_geom): detb= detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (src, detb + i * detu) + yield ((src, detb + i * detu), + (src, detb + (i - 0.5) * detu), + (src, detb + (i + 0.5) * detu)) def gen_lines_fanflat_vec(proj_geom): v = proj_geom['Vectors'] @@ -87,7 +136,9 @@ def gen_lines_fanflat_vec(proj_geom): detb = detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (src, detb + i * detu) + yield ((src, detb + i * detu), + (src, detb + (i - 0.5) * detu), + (src, detb + (i + 0.5) * detu)) def gen_lines_parallel(proj_geom): angles = proj_geom['ProjectionAngles'] @@ -106,7 +157,10 @@ def gen_lines_parallel(proj_geom): detb= detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (detb + i * detu - ray, detb + i * detu) + yield ((detb + i * detu - ray, detb + i * detu), + (detb + (i - 0.5) * detu - ray, detb + (i - 0.5) * detu), + (detb + (i + 0.5) * detu - ray, detb + (i + 0.5) * detu)) + def gen_lines_parallel_vec(proj_geom): v = proj_geom['Vectors'] @@ -118,7 +172,9 @@ def gen_lines_parallel_vec(proj_geom): detb = detc + (0.5 - 0.5*proj_geom['DetectorCount']) * detu for i in range(proj_geom['DetectorCount']): - yield (detb + i * detu - ray, detb + i * detu) + yield ((detb + i * detu - ray, detb + i * detu), + (detb + (i - 0.5) * detu - ray, detb + (i - 0.5) * detu), + (detb + (i + 0.5) * detu - ray, detb + (i + 0.5) * detu)) def gen_lines(proj_geom): @@ -179,8 +235,8 @@ def gen_random_geometry_parallel_vec(): nloops = 50 seed = 123 -class TestLineKernel(unittest.TestCase): - def single_test(self, type): +class Test2DKernel(unittest.TestCase): + def single_test(self, type, proj_type): shape = np.random.randint(*range2d, size=2) # these rectangles are biased, but that shouldn't matter rect_min = [ np.random.randint(0, a) for a in shape ] @@ -201,16 +257,16 @@ class TestLineKernel(unittest.TestCase): if type == 'parallel': pg = gen_random_geometry_parallel() - projector_id = astra.create_projector('line', pg, vg) + projector_id = astra.create_projector(proj_type, pg, vg) elif type == 'parallel_vec': pg = gen_random_geometry_parallel_vec() - projector_id = astra.create_projector('line', pg, vg) + projector_id = astra.create_projector(proj_type, pg, vg) elif type == 'fanflat': pg = gen_random_geometry_fanflat() - projector_id = astra.create_projector('line_fanflat', pg, vg) + projector_id = astra.create_projector(proj_type + '_fanflat', pg, vg) elif type == 'fanflat_vec': pg = gen_random_geometry_fanflat_vec() - projector_id = astra.create_projector('line_fanflat', pg, vg) + projector_id = astra.create_projector(proj_type + '_fanflat', pg, vg) data = np.zeros((shape[1], shape[0]), dtype=np.float32) @@ -225,84 +281,129 @@ class TestLineKernel(unittest.TestCase): astra.projector.delete(projector_id) - a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) - b = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) - c = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) - - i = 0 - #print( origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], origin[1] + (-0.5 * shape[1] + rect_min[1]) * pixsize[1], origin[1] + (-0.5 * shape[1] + rect_max[1]) * pixsize[1]) - for src, det in gen_lines(pg): - #print(src,det) - - # NB: Flipped y-axis here, since that is how astra interprets 2D volumes - # We compute line intersections with slightly bigger (cw) and - # smaller (aw) rectangles, and see if the kernel falls - # between these two values. - (aw,bw,cw) = intersect_line_rectangle_interval(src, det, - origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], - origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], - origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1], - origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1], - 1e-3) - a[i] = aw - b[i] = bw - c[i] = cw - i += 1 - # Add weight for pixel / voxel size - try: - detweight = pg['DetectorWidth'] - except KeyError: - detweight = np.sqrt(pg['Vectors'][0,4]*pg['Vectors'][0,4] + pg['Vectors'][0,5]*pg['Vectors'][0,5] ) - a *= detweight - b *= detweight - c *= detweight - a = a.reshape(astra.functions.geom_size(pg)) - b = b.reshape(astra.functions.geom_size(pg)) - c = c.reshape(astra.functions.geom_size(pg)) - - # Check if sinogram lies between a and c - y = np.min(sinogram-a) - z = np.min(c-sinogram) - x = np.max(np.abs(sinogram-b)) # ideally this is small, but can be large - # due to discontinuities in line kernel - self.assertFalse(z < 0 or y < 0) - if z < 0 or y < 0: - print(y,z,x) - pylab.gray() - pylab.imshow(data) - pylab.figure() - pylab.imshow(sinogram) - pylab.figure() - pylab.imshow(b) - pylab.figure() - pylab.imshow(a) - pylab.figure() - pylab.imshow(c) - pylab.figure() - pylab.imshow(sinogram-a) - pylab.figure() - pylab.imshow(c-sinogram) - pylab.show() + if proj_type == 'line': + + a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) + b = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) + c = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) + + i = 0 + #print( origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], origin[1] + (-0.5 * shape[1] + rect_min[1]) * pixsize[1], origin[1] + (-0.5 * shape[1] + rect_max[1]) * pixsize[1]) + for center, edge1, edge2 in gen_lines(pg): + (src, det) = center + #print(src,det) + + # NB: Flipped y-axis here, since that is how astra interprets 2D volumes + # We compute line intersections with slightly bigger (cw) and + # smaller (aw) rectangles, and see if the kernel falls + # between these two values. + (aw,bw,cw) = intersect_line_rectangle_interval(src, det, + origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0], + origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], + origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1], + origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1], + 1e-3) + a[i] = aw + b[i] = bw + c[i] = cw + i += 1 + # Add weight for pixel / voxel size + try: + detweight = pg['DetectorWidth'] + except KeyError: + detweight = np.sqrt(pg['Vectors'][0,4]*pg['Vectors'][0,4] + pg['Vectors'][0,5]*pg['Vectors'][0,5] ) + a *= detweight + b *= detweight + c *= detweight + a = a.reshape(astra.functions.geom_size(pg)) + b = b.reshape(astra.functions.geom_size(pg)) + c = c.reshape(astra.functions.geom_size(pg)) + + # Check if sinogram lies between a and c + y = np.min(sinogram-a) + z = np.min(c-sinogram) + x = np.max(np.abs(sinogram-b)) # ideally this is small, but can be large + # due to discontinuities in line kernel + self.assertFalse(z < 0 or y < 0) + if z < 0 or y < 0: + print(y,z,x) + pylab.gray() + pylab.imshow(data) + pylab.figure() + pylab.imshow(sinogram) + pylab.figure() + pylab.imshow(b) + pylab.figure() + pylab.imshow(a) + pylab.figure() + pylab.imshow(c) + pylab.figure() + pylab.imshow(sinogram-a) + pylab.figure() + pylab.imshow(c-sinogram) + pylab.show() + elif proj_type == 'distance_driven': + a = np.zeros(np.prod(astra.functions.geom_size(pg)), dtype=np.float32) + i = 0 + for center, edge1, edge2 in gen_lines(pg): + (xd, yd) = center[1] - center[0] + l = 0.0 + if np.abs(xd) > np.abs(yd): # horizontal ray + y_seg = (origin[1] + (+0.5 * shape[1] - rect_max[1]) * pixsize[1], + origin[1] + (+0.5 * shape[1] - rect_min[1]) * pixsize[1]) + for j in range(rect_min[0], rect_max[0]): + x = origin[0] + (-0.5 * shape[0] + j + 0.5) * pixsize[0] + l += intersect_ray_vertical_segment(edge1, edge2, x, y_seg) * pixsize[0] + else: + x_seg = (origin[0] + (-0.5 * shape[0] + rect_max[0]) * pixsize[0], + origin[0] + (-0.5 * shape[0] + rect_min[0]) * pixsize[0]) + for j in range(rect_min[1], rect_max[1]): + y = origin[1] + (+0.5 * shape[1] - j - 0.5) * pixsize[1] + l += intersect_ray_horizontal_segment(edge1, edge2, y, x_seg) * pixsize[1] + a[i] = l + i += 1 + a = a.reshape(astra.functions.geom_size(pg)) + x = np.max(np.abs(sinogram-a)) + if x > 2e-3: + pylab.gray() + pylab.imshow(data) + pylab.figure() + pylab.imshow(sinogram) + pylab.figure() + pylab.imshow(a) + pylab.figure() + pylab.imshow(sinogram-a) + pylab.show() + self.assertFalse(x > 2e-3) + def test_par(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('parallel') + self.single_test('parallel', 'line') + def test_par_dd(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel', 'distance_driven') def test_fan(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('fanflat') + self.single_test('fanflat', 'line') def test_parvec(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('parallel_vec') + self.single_test('parallel_vec', 'line') + def test_parvec_dd(self): + np.random.seed(seed) + for _ in range(nloops): + self.single_test('parallel_vec', 'distance_driven') def test_fanvec(self): np.random.seed(seed) for _ in range(nloops): - self.single_test('fanflat_vec') + self.single_test('fanflat_vec', 'line') + - if __name__ == '__main__': unittest.main() |