Merge pull request #183 from wjp/par2d_dd

Add basic implementation of par2d CPU Distance Driven projector

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()