summaryrefslogtreecommitdiffstats
path: root/include
diff options
context:
space:
mode:
Diffstat (limited to 'include')
-rw-r--r--include/astra/ParallelBeamDistanceDrivenProjector2D.h205
-rw-r--r--include/astra/ParallelBeamDistanceDrivenProjector2D.inl227
-rw-r--r--include/astra/Projector2DImpl.inl1
-rw-r--r--include/astra/ProjectorTypelist.h7
4 files changed, 438 insertions, 2 deletions
diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.h b/include/astra/ParallelBeamDistanceDrivenProjector2D.h
new file mode 100644
index 0000000..67dd21b
--- /dev/null
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.h
@@ -0,0 +1,205 @@
+/*
+-----------------------------------------------------------------------
+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);
+
+ /** 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.
+ *
+ * @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..6bf3b56
--- /dev/null
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
@@ -0,0 +1,227 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2018, imec Vision Lab, University of Antwerp
+ 2014-2018, CWI, Amsterdam
+
+Contact: astra@astra-toolbox.com
+Website: http://www.astra-toolbox.com/
+
+This file is part of the ASTRA Toolbox.
+
+
+The ASTRA Toolbox is free software: you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation, either version 3 of the License, or
+(at your option) any later version.
+
+The ASTRA Toolbox is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with the ASTRA Toolbox. If not, see <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];
+
+ float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
+
+ const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
+
+ const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+ const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
+
+ // loop detectors
+ for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
+
+ const int iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
+
+ // POLICY: RAY PRIOR
+ if (!p.rayPrior(iRayIndex)) continue;
+
+ const float32 Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+ const float32 Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+
+ if (vertical && true) {
+
+ const float32 RxOverRy = proj->fRayX/proj->fRayY;
+ // TODO: Determine det/pixel scaling factors
+ const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
+ const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+ const float32 deltad = fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX);
+
+ // calculate c for row 0
+ float32 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX + 0.5f;
+
+ // loop rows
+ for (int row = 0; row < rowCount; ++row, c+= deltac) {
+
+ // horizontal extent of ray in center of this row:
+ // [ c - deltad/2 , c + deltad/2 ]
+
+ // |-gapBegin-*---|------|----*-gapEnd-|
+ // * = ray extent intercepts; c - deltad/2 and c + deltad/2
+ // | = pixel column edges
+
+ const int colBegin = (int)floor(c - deltad/2.0f);
+ const int colEnd = (int)ceil(c + deltad/2.0f);
+
+ // TODO: Optimize volume edge checks
+
+ int iVolumeIndex = row * colCount + colBegin;
+ if (colBegin + 1 == colEnd) {
+
+ if (colBegin >= 0 && colBegin < colCount)
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ deltad * lengthPerRow);
+ } else {
+ const float gapBegin = (c - deltad/2.0f) - (float32)colBegin;
+ const float gapEnd = (float32)colEnd - (c + deltad/2.0f);
+ float tot = 1.0f - gapBegin;
+ if (colBegin >= 0 && colBegin < colCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ (1.0f - gapBegin) * lengthPerRow);
+ }
+ iVolumeIndex++;
+
+ for (int col = colBegin + 1; col + 1 < colEnd; ++col, ++iVolumeIndex) {
+ tot += 1.0f;
+ if (col >= 0 && col < colCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow);
+ }
+ }
+ assert(iVolumeIndex == row * colCount + colEnd - 1);
+ tot += 1.0f - gapEnd;
+ if (colEnd > 0 && colEnd <= colCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ (1.0f - gapEnd) * lengthPerRow);
+ }
+ assert(fabs(tot - deltad) < 0.0001);
+ }
+
+ }
+
+ } else if (!vertical && true) {
+
+ const float32 RyOverRx = proj->fRayY/proj->fRayX;
+ // TODO: Determine det/pixel scaling factors
+ const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
+ const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
+ const float32 deltad = fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY);
+
+ // calculate r for col 0
+ float32 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY + 0.5f;
+
+ // loop columns
+ for (int col = 0; col < colCount; ++col, r+= deltar) {
+
+ // vertical extent of ray in center of this column:
+ // [ r - deltad/2 , r + deltad/2 ]
+
+ const int rowBegin = (int)floor(r - deltad/2.0f);
+ const int rowEnd = (int)ceil(r + deltad/2.0f);
+
+ // TODO: Optimize volume edge checks
+
+ int iVolumeIndex = rowBegin * colCount + col;
+ if (rowBegin + 1 == rowEnd) {
+
+ if (rowBegin >= 0 && rowBegin < rowCount)
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ deltad * lengthPerCol);
+ } else {
+ const float gapBegin = (r - deltad/2.0f) - (float32)rowBegin;
+ const float gapEnd = (float32)rowEnd - (r + deltad/2.0f);
+ float tot = 1.0f - gapBegin;
+
+ if (rowBegin >= 0 && rowBegin < rowCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ (1.0f - gapBegin) * lengthPerCol);
+ }
+ iVolumeIndex += colCount;
+
+ for (int row = rowBegin + 1; row + 1 < rowEnd; ++row, iVolumeIndex += colCount) {
+ tot += 1.0f;
+ if (row >= 0 && row < rowCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol);
+ }
+ }
+ assert(iVolumeIndex == (rowEnd - 1) * colCount + col);
+ tot += 1.0f - gapEnd;
+ if (rowEnd > 0 && rowEnd <= rowCount) {
+ policy_weight(p, iRayIndex, iVolumeIndex,
+ (1.0f - gapEnd) * lengthPerCol);
+ }
+ assert(fabs(tot - deltad) < 0.0001);
+ }
+
+ }
+
+ }
+
+ // POLICY: RAY POSTERIOR
+ p.rayPosterior(iRayIndex);
+
+ }
+ }
+
+ if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+ delete pVecProjectionGeometry;
+}
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,
generated by cgit v1.2.3 (git 2.25.1) at 2024-12-27 10:03:45 +0000