From b9563df34df98480d35a53caa2b81d998440f9c5 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 24 Jan 2019 13:26:30 +0100
Subject: Add basic implementation of par2d CPU Distance Driven projector

---
 .../astra/ParallelBeamDistanceDrivenProjector2D.h  | 205 +++++++++++++++++++
 .../ParallelBeamDistanceDrivenProjector2D.inl      | 227 +++++++++++++++++++++
 include/astra/Projector2DImpl.inl                  |   1 +
 include/astra/ProjectorTypelist.h                  |   7 +-
 4 files changed, 438 insertions(+), 2 deletions(-)
 create mode 100644 include/astra/ParallelBeamDistanceDrivenProjector2D.h
 create mode 100644 include/astra/ParallelBeamDistanceDrivenProjector2D.inl

(limited to 'include')

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,
-- 
cgit v1.2.3


From e2d5375ecc5c8fc45b796e0bd9d948cd729abcc9 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 24 Jan 2019 13:48:49 +0100
Subject: Remove largely unimplemented CProjector2D::projectPoint method

This includes the astra_mex_projector('splat') matlab function.
---
 include/astra/FanFlatBeamLineKernelProjector2D.h        | 8 --------
 include/astra/FanFlatBeamStripKernelProjector2D.h       | 8 --------
 include/astra/ParallelBeamBlobKernelProjector2D.h       | 8 --------
 include/astra/ParallelBeamDistanceDrivenProjector2D.h   | 8 --------
 include/astra/ParallelBeamDistanceDrivenProjector2D.inl | 2 --
 include/astra/ParallelBeamLineKernelProjector2D.h       | 8 --------
 include/astra/ParallelBeamLinearKernelProjector2D.h     | 8 --------
 include/astra/ParallelBeamLinearKernelProjector2D.inl   | 2 +-
 include/astra/ParallelBeamStripKernelProjector2D.h      | 8 --------
 include/astra/Projector2D.h                             | 8 --------
 include/astra/SparseMatrixProjector2D.h                 | 8 --------
 11 files changed, 1 insertion(+), 75 deletions(-)

(limited to 'include')

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
index 67dd21b..6d1d633 100644
--- a/include/astra/ParallelBeamDistanceDrivenProjector2D.h
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.h
@@ -138,14 +138,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.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
index 6bf3b56..7b45ed1 100644
--- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
@@ -81,8 +81,6 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 
 		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
 
-		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
-
 		const bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
 
 		const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
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/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.
 	 *
-- 
cgit v1.2.3


From 37643b309520bcd12afcee430210632db305d21f Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Thu, 24 Jan 2019 14:18:11 +0100
Subject: Some basic optimizations

---
 .../ParallelBeamDistanceDrivenProjector2D.inl      | 87 ++++++++++------------
 1 file changed, 41 insertions(+), 46 deletions(-)

(limited to 'include')

diff --git a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
index 7b45ed1..aedcee9 100644
--- a/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
+++ b/include/astra/ParallelBeamDistanceDrivenProjector2D.inl
@@ -97,13 +97,12 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 			const float32 Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
 			const float32 Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
 
-			if (vertical && true) {
+			if (vertical) {
 
 				const float32 RxOverRy = proj->fRayX/proj->fRayY;
-				// TODO: Determine det/pixel scaling factors
 				const float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
 				const float32 deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
-				const float32 deltad = fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX);
+				const float32 deltad = 0.5f * fabs((proj->fDetUX - proj->fDetUY * RxOverRy) * inv_pixelLengthX);
 
 				// calculate c for row 0
 				float32 c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX + 0.5f;
@@ -112,57 +111,55 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 				for (int row = 0; row < rowCount; ++row, c+= deltac) {
 
 					// horizontal extent of ray in center of this row:
-					// [ c - deltad/2 , c + deltad/2 ]
+					// [ c - deltad , c + deltad ]
 
 					// |-gapBegin-*---|------|----*-gapEnd-|
-					// * = ray extent intercepts; c - deltad/2 and c + deltad/2
+					// * = ray extent intercepts; c - deltad and c + deltad
 					// | = pixel column edges
 
-					const int colBegin = (int)floor(c - deltad/2.0f);
-					const int colEnd = (int)ceil(c + deltad/2.0f);
+					const int colBegin = (int)floor(c - deltad);
+					const int colEnd = (int)ceil(c + deltad);
 
-					// TODO: Optimize volume edge checks
+					if (colBegin >= colCount || colEnd <= 0)
+						continue;
 
 					int iVolumeIndex = row * colCount + colBegin;
 					if (colBegin + 1 == colEnd) {
 
 						if (colBegin >= 0 && colBegin < colCount)
 							policy_weight(p, iRayIndex, iVolumeIndex,
-							              deltad * lengthPerRow);
+							              2.0f * deltad * lengthPerRow);
 					} else {
-						const float gapBegin = (c - deltad/2.0f) - (float32)colBegin;
-						const float gapEnd = (float32)colEnd - (c + deltad/2.0f);
-						float tot = 1.0f - gapBegin;
-						if (colBegin >= 0 && colBegin < colCount) {
+
+						if (colBegin >= 0) {
+							const float gapBegin = (c - deltad) - (float32)colBegin;
 							policy_weight(p, iRayIndex, iVolumeIndex,
 							              (1.0f - gapBegin) * lengthPerRow);
 						}
-						iVolumeIndex++;
 
-						for (int col = colBegin + 1; col + 1 < colEnd; ++col, ++iVolumeIndex) {
-							tot += 1.0f;
-							if (col >= 0 && col < colCount) {
-								policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow);
-							}
+						const int clippedMColBegin = std::max(colBegin + 1, 0);
+						const int clippedMColEnd = std::min(colEnd - 1, colCount);
+						iVolumeIndex = row * colCount + clippedMColBegin;
+						for (int col = clippedMColBegin; col < clippedMColEnd; ++col, ++iVolumeIndex) {
+							policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow);
 						}
-						assert(iVolumeIndex == row * colCount + colEnd - 1);
-						tot += 1.0f - gapEnd;
-						if (colEnd > 0 && colEnd <= colCount) {
+
+						iVolumeIndex = row * colCount + colEnd - 1;
+						if (colEnd <= colCount) {
+							const float gapEnd = (float32)colEnd - (c + deltad);
 							policy_weight(p, iRayIndex, iVolumeIndex,
 						 	             (1.0f - gapEnd) * lengthPerRow);
 						}
-						assert(fabs(tot - deltad) < 0.0001);
 					}
 
 				}
 				
-			} else if (!vertical && true) {
+			} else {
 
 				const float32 RyOverRx = proj->fRayY/proj->fRayX;
-				// TODO: Determine det/pixel scaling factors
 				const float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * m_pVolumeGeometry->getPixelLengthY();
 				const float32 deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
-				const float32 deltad = fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY);
+				const float32 deltad = 0.5f * fabs((proj->fDetUY - proj->fDetUX * RyOverRx) * inv_pixelLengthY);
 
 				// calculate r for col 0
 				float32 r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY + 0.5f;
@@ -171,43 +168,41 @@ void CParallelBeamDistanceDrivenProjector2D::projectBlock_internal(int _iProjFro
 				for (int col = 0; col < colCount; ++col, r+= deltar) {
 
 					// vertical extent of ray in center of this column:
-					// [ r - deltad/2 , r + deltad/2 ]
+					// [ r - deltad , r + deltad ]
 
-					const int rowBegin = (int)floor(r - deltad/2.0f);
-					const int rowEnd = (int)ceil(r + deltad/2.0f);
+					const int rowBegin = (int)floor(r - deltad);
+					const int rowEnd = (int)ceil(r + deltad);
 
-					// TODO: Optimize volume edge checks
+					if (rowBegin >= rowCount || rowEnd <= 0)
+						continue;
 
 					int iVolumeIndex = rowBegin * colCount + col;
 					if (rowBegin + 1 == rowEnd) {
 
 						if (rowBegin >= 0 && rowBegin < rowCount)
 							policy_weight(p, iRayIndex, iVolumeIndex,
-							              deltad * lengthPerCol);
+							              2.0f * deltad * lengthPerCol);
 					} else {
-						const float gapBegin = (r - deltad/2.0f) - (float32)rowBegin;
-						const float gapEnd = (float32)rowEnd - (r + deltad/2.0f);
-						float tot = 1.0f - gapBegin;
-
-						if (rowBegin >= 0 && rowBegin < rowCount) {
+						if (rowBegin >= 0) {
+							const float gapBegin = (r - deltad) - (float32)rowBegin;
 							policy_weight(p, iRayIndex, iVolumeIndex,
 							              (1.0f - gapBegin) * lengthPerCol);
 						}
-						iVolumeIndex += colCount;
 
-						for (int row = rowBegin + 1; row + 1 < rowEnd; ++row, iVolumeIndex += colCount) {
-							tot += 1.0f;
-							if (row >= 0 && row < rowCount) {
-								policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol);
-							}
+						const int clippedMRowBegin = std::max(rowBegin + 1, 0);
+						const int clippedMRowEnd = std::min(rowEnd - 1, rowCount);
+						iVolumeIndex = clippedMRowBegin * colCount + col;
+
+						for (int row = clippedMRowBegin; row < clippedMRowEnd; ++row, iVolumeIndex += colCount) {
+							policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol);
 						}
-						assert(iVolumeIndex == (rowEnd - 1) * colCount + col);
-						tot += 1.0f - gapEnd;
-						if (rowEnd > 0 && rowEnd <= rowCount) {
+
+						iVolumeIndex = (rowEnd - 1) * colCount + col;
+						if (rowEnd <= rowCount) {
+							const float gapEnd = (float32)rowEnd - (r + deltad);
 							policy_weight(p, iRayIndex, iVolumeIndex,
 						 	             (1.0f - gapEnd) * lengthPerCol);
 						}
-						assert(fabs(tot - deltad) < 0.0001);
 					}
 
 				}
-- 
cgit v1.2.3