Merge branch 'parallel_vec'

16 files changed, 1128 insertions, 912 deletions

diff --git a/include/astra/CudaFilteredBackProjectionAlgorithm.h b/include/astra/CudaFilteredBackProjectionAlgorithm.h
index 164b51b..55191ef 100644
--- a/include/astra/CudaFilteredBackProjectionAlgorithm.h
+++ b/include/astra/CudaFilteredBackProjectionAlgorithm.h
@@ -25,30 +25,26 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 -----------------------------------------------------------------------
 */
 
-#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM2_H
-#define CUDAFILTEREDBACKPROJECTIONALGORITHM2_H
+#ifndef CUDAFILTEREDBACKPROJECTIONALGORITHM_H
+#define CUDAFILTEREDBACKPROJECTIONALGORITHM_H
 
 #ifdef ASTRA_CUDA
 
 #include <astra/Float32ProjectionData2D.h>
 #include <astra/Float32VolumeData2D.h>
-#include <astra/ReconstructionAlgorithm2D.h>
+#include <astra/CudaReconstructionAlgorithm2D.h>
 
 #include "../../cuda/2d/astra.h"
 
 namespace astra
 {
 
-class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CReconstructionAlgorithm2D
+class _AstraExport CCudaFilteredBackProjectionAlgorithm : public CCudaReconstructionAlgorithm2D
 {
 public:
 	static std::string type;
 
 private:
-	CFloat32ProjectionData2D * m_pSinogram;
-	CFloat32VolumeData2D * m_pReconstruction;
-	int m_iGPUIndex;
-	int m_iPixelSuperSampling;
 	E_FBPFILTER m_eFilter;
 	float * m_pfFilter;
 	int m_iFilterWidth;	// number of elements per projection direction in filter
@@ -65,11 +61,6 @@ public:
 	virtual bool initialize(const Config& _cfg);
 	bool initialize(CFloat32ProjectionData2D * _pSinogram, CFloat32VolumeData2D * _pReconstruction, E_FBPFILTER _eFilter, const float * _pfFilter = NULL, int _iFilterWidth = 0, int _iGPUIndex = -1, float _fFilterParameter = -1.0f);
 
-	virtual void run(int _iNrIterations = 0);
-
-	static int calcIdealRealFilterWidth(int _iDetectorCount);
-	static int calcIdealFourierFilterWidth(int _iDetectorCount);
-	
 	/** Get a description of the class.
 	 *
 	 * @return description string
@@ -79,12 +70,7 @@ public:
 protected:
 	bool check();
 
-	AstraFBP* m_pFBP;
-
-	bool m_bAstraFBPInit;
-
-	void initializeFromProjector();
-	virtual bool requiresProjector() const { return false; }
+	virtual void initCUDAAlgorithm();
 };
 
 // inline functions
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl
index d967844..eb73de8 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl
@@ -25,6 +25,7 @@ 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 CFanFlatBeamLineKernelProjector2D::project(Policy& p)
@@ -48,246 +49,166 @@ void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _
 }
 
 //----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+// PROJECT BLOCK - vector projection geometry
 template <typename Policy>
 void CFanFlatBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
 {
-	// variables
-	float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2;
-	float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY;
-	int iVolumeIndex, iRayIndex, row, col, iAngle, iDetector, x1;
-	bool switch_t;
+	// get vector geometry
+	const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry;
+	if (dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry)) {
+		pVecProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry)->toVectorGeometry();
+	} else {
+		pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(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();
+	const int detCount = pVecProjectionGeometry->getDetectorCount();
+	const float32 Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+	const float32 Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
-	const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
-	const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry);
+	// loop angles
+	for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 
-	float32 old_theta, theta, alpha;
-	const SFanProjection * proj = 0;
+		// variables
+		float32 Dx, Dy, Rx, Ry, S, T, weight, c, r, deltac, deltar, offset, RxOverRy, RyOverRx;
+		float32 lengthPerRow, lengthPerCol, invTminSTimesLengthPerRow, invTminSTimesLengthPerCol;
+		int iVolumeIndex, iRayIndex, row, col, iDetector;
 
-	// loop angles
-	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
-		// get theta
-		if (pProjectionGeometry) {
-			old_theta = pProjectionGeometry->getProjectionAngle(iAngle);
-		}
-		else if (pVecProjectionGeometry) {
-			proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
-			old_theta = atan2(-proj->fSrcX, proj->fSrcY);
-			if (old_theta < 0) old_theta += 2*PI;
-		} else {
-			assert(false);
-		}
-
-		switch_t = false;
-		if (old_theta >= 7*PIdiv4) old_theta -= 2*PI;
-		if (old_theta >= 3*PIdiv4) {
-			old_theta -= PI;
-			switch_t = true;
-		}
+		const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
 
 		// loop detectors
 		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
-			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
+			iRayIndex = iAngle * detCount + iDetector;
 
 			// POLICY: RAY PRIOR
 			if (!p.rayPrior(iRayIndex)) continue;
+	
+			Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+			Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
 
-			// get values
-			if (pProjectionGeometry) {
-				t = -pProjectionGeometry->indexToDetectorOffset(iDetector);
-				alpha = atan(t / pProjectionGeometry->getSourceDetectorDistance());
-				t = sin(alpha) * pProjectionGeometry->getOriginSourceDistance();
-			}
-			else if (pVecProjectionGeometry) {
-				float32 detX = proj->fDetSX + proj->fDetUX*(0.5f + iDetector);
-				float32 detY = proj->fDetSY + proj->fDetUY*(0.5f + iDetector);
-				alpha = angleBetweenVectors(-proj->fSrcX, -proj->fSrcY, detX - proj->fSrcX, detY - proj->fSrcY);
-				t = sin(alpha) * sqrt(proj->fSrcX*proj->fSrcX + proj->fSrcY*proj->fSrcY);
-			} else {
-				assert(false);
-			}
-
-			if (switch_t) t = -t;
-			theta = old_theta + alpha;
-
-			// precalculate sin, cos, 1/cos
-			sin_theta = sin(theta);
-			cos_theta = cos(theta);
-			inv_sin_theta = 1.0f / sin_theta; 
-			inv_cos_theta = 1.0f / cos_theta; 
-
-			// precalculate kernel limits
-			lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-			updatePerRow = sin_theta * inv_cos_theta;
-			inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-
-			// precalculate kernel limits
-			lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
-			updatePerCol = cos_theta * inv_sin_theta;
-			inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
+			Rx = proj->fSrcX - Dx;
+			Ry = proj->fSrcY - Dy;
 
-			// precalculate S and T
-			S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
-			T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
+			bool vertical = fabs(Rx) < fabs(Ry);
+			bool isin = false;
 
 			// vertically
-			if (old_theta <= PIdiv4) {
-			
-				// calculate x for row 0
-				P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-				x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
+			if (vertical) {
+				RxOverRy = Rx/Ry;
+				lengthPerRow = detSize * pixelLengthX * sqrt(Rx*Rx + Ry*Ry) / abs(Ry);
+				deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+				S = 0.5f - 0.5f*fabs(RxOverRy);
+				T = 0.5f + 0.5f*fabs(RxOverRy);
+				invTminSTimesLengthPerRow = lengthPerRow / (T - S);
+
+				// calculate c for row 0
+				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
 				// for each row
-				for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-					
-					// get coords
-					x1 = int((x > 0.0f) ? x : x-1.0f);
-					x2 = x - x1; 
-					x += updatePerRow;
+				for (row = 0; row < rowCount; ++row, c += deltac) {
 
-					if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
+					col = int(floor(c+0.5f));
+					if (col < -1 || col > colCount) { if (!isin) continue; else break; }
+					offset = c - float32(col);
 
 					// left
-					if (x2 < 0.5f-S) {
-						I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
-						if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridColCount()*/) {//x1 is always less than or equal to gridColCount because of the "continue" in the beginning of the for-loop
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-					}
+					if (offset < -S) {
+						weight = (offset + T) * invTminSTimesLengthPerRow;
+
+						iVolumeIndex = row * colCount + col - 1;
+						if (col > 0) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight); }
 
-					// center
-					else if (x2 <= 0.5f+S) {
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}					
+						iVolumeIndex++;
+						if (col >= 0 && col < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
 					// right
-					else  if (x2 <= 1.0f) {
-						I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-						if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridColCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
+					else if (S < offset) {
+						weight = (offset - S) * invTminSTimesLengthPerRow;
+
+						iVolumeIndex = row * colCount + col;
+						if (col >= 0 && col < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight); }
+
+						iVolumeIndex++;
+						if (col + 1 < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); } 
+					}
+
+					// centre
+					else if (col >= 0 && col < colCount) {
+						iVolumeIndex = row * colCount + col;
+						policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow);
 					}
+					isin = true;
 				}
 			}
 
 			// horizontally
-			//else if (PIdiv4 <= old_theta && old_theta <= 3*PIdiv4) {
 			else {
+				RyOverRx = Ry/Rx;
+				lengthPerCol = detSize * pixelLengthY * sqrt(Rx*Rx + Ry*Ry) / abs(Rx);
+				deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
+				S = 0.5f - 0.5f*fabs(RyOverRx);
+				T = 0.5f + 0.5f*fabs(RyOverRx);
+				invTminSTimesLengthPerCol = lengthPerCol / (T - S);
 
-				// calculate point P
-				P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-				x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY;
+				// calculate r for col 0
+				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
 
 				// for each col
-				for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-				
-					// get coords
-					x1 = int((x > 0.0f) ? x : x-1.0f);
-					x2 = x - x1; 
-					x += updatePerCol;
+				for (col = 0; col < colCount; ++col, r += deltar) {
 
-					if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
+					row = int(floor(r+0.5f));
+					if (row < -1 || row > rowCount) { if (!isin) continue; else break; }
+					offset = r - float32(row);
 
 					// up
-					if (x2 < 0.5f-T) {
-						I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-						if (x1-1 >= 0 /*&& x1-1 < m_pVolumeGeometry->getGridRowCount()*/) {//x1 is always less than or equal to gridRowCount because of the "continue" in the beginning of the for-loop
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-					}
+					if (offset < -S) {
+						weight = (offset + T) * invTminSTimesLengthPerCol;
+
+						iVolumeIndex = (row-1) * colCount + col;
+						if (row > 0) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight); }
 
-					// center
-					else if (x2 <= 0.5f+T) {
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}					
+						iVolumeIndex += colCount;
+						if (row >= 0 && row < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
 					// down
-					else  if (x2 <= 1.0f) {
-						I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-						if (/*x1+1 >= 0 &&*/ x1+1 < m_pVolumeGeometry->getGridRowCount()) {//x1 is always greater than or equal to -1 because of the "continue" in the beginning of the for-loop
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
+					else if (S < offset) {
+						weight = (offset - S) * invTminSTimesLengthPerCol;
+
+						iVolumeIndex = row * colCount + col;
+						if (row >= 0 && row < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight); }
+
+						iVolumeIndex += colCount;
+						if (row + 1 < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
+					}
+
+					// centre
+					else if (row >= 0 && row < rowCount) {
+						iVolumeIndex = row * colCount + col;
+						policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol);
 					}
+					isin = true;
 				}
-			} // end loop col
+			}
 	
 			// POLICY: RAY POSTERIOR
 			p.rayPosterior(iRayIndex);
 	
 		} // end loop detector
 	} // end loop angles
+
+	// Delete created vec geometry if required
+	if (dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry))
+		delete pVecProjectionGeometry;
+
 }
diff --git a/include/astra/FanFlatProjectionGeometry2D.h b/include/astra/FanFlatProjectionGeometry2D.h
index e9a0535..c45eb78 100644
--- a/include/astra/FanFlatProjectionGeometry2D.h
+++ b/include/astra/FanFlatProjectionGeometry2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #define _INC_ASTRA_FANFLATPROJECTIONGEOMETRY2D
 
 #include "ProjectionGeometry2D.h"
+#include "FanFlatVecProjectionGeometry2D.h"
 
 #include <cmath>
 
@@ -189,6 +190,10 @@ public:
 	 * @return a unit vector describing the direction
 	 */
 	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex);
+
+	/** Create a vector geom
+	*/
+	CFanFlatVecProjectionGeometry2D* toVectorGeometry();	
 };
 
 
diff --git a/include/astra/GeometryUtil2D.h b/include/astra/GeometryUtil2D.h
index 6434d3c..914e40d 100644
--- a/include/astra/GeometryUtil2D.h
+++ b/include/astra/GeometryUtil2D.h
@@ -30,17 +30,77 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 
 namespace astra {
 
+struct SParProjection {
+	// the ray direction
+	float fRayX, fRayY;
+
+	// the start of the (linear) detector
+	float fDetSX, fDetSY;
+
+	// the length of a single detector pixel
+	float fDetUX, fDetUY;
+
+
+	void translate(double dx, double dy) {
+		fDetSX += dx;
+		fDetSY += dy;
+	}
+	void scale(double factor) {
+		fRayX *= factor;
+		fRayY *= factor;
+		fDetSX *= factor;
+		fDetSY *= factor;
+		fDetUX *= factor;
+		fDetUY *= factor;
+	}
+};
+
+
 struct SFanProjection {
-        // the source
-        float fSrcX, fSrcY;
+	// the source
+	float fSrcX, fSrcY;
+
+	// the start of the (linear) detector
+	float fDetSX, fDetSY;
 
-        // the start of the (linear) detector
-        float fDetSX, fDetSY;
+	// the length of a single detector pixel
+	float fDetUX, fDetUY;
 
-        // the length of a single detector pixel
-        float fDetUX, fDetUY;
+	void translate(double dx, double dy) {
+		fSrcX += dx;
+		fSrcY += dy;
+		fDetSX += dx;
+		fDetSY += dy;
+	}
+	void scale(double factor) {
+		fSrcX *= factor;
+		fSrcY *= factor;
+		fDetSX *= factor;
+		fDetSY *= factor;
+		fDetUX *= factor;
+		fDetUY *= factor;
+	}
 };
 
+
+
+SParProjection* genParProjections(unsigned int iProjAngles,
+                                  unsigned int iProjDets,
+                                  double fDetSize,
+                                  const float *pfAngles,
+                                  const float *pfExtraOffsets);
+
+SFanProjection* genFanProjections(unsigned int iProjAngles,
+                                  unsigned int iProjDets,
+                                  double fOriginSource, double fOriginDetector,
+                                  double fDetSize,
+                                  const float *pfAngles);
+
+bool getParParameters(const SParProjection &proj, unsigned int iProjDets, float &fAngle, float &fDetSize, float &fOffset);
+
+bool getFanParameters(const SFanProjection &proj, unsigned int iProjDets, float &fAngle, float &fOriginSource, float &fOriginDetector, float &fDetSize, float &fOffset);
+
+
 }
 
 #endif
diff --git a/include/astra/Globals.h b/include/astra/Globals.h
index 128422c..8d0d619 100644
--- a/include/astra/Globals.h
+++ b/include/astra/Globals.h
@@ -127,7 +127,7 @@ namespace astra {
 	const float32 PI32 = 3.14159265358979323846264338328f;
 	const float32 PIdiv2 = PI / 2;
 	const float32 PIdiv4 = PI / 4;
-	const float32 eps = 1e-7f;
+	const float32 eps = 1e-6f;
 	
 	extern _AstraExport bool running_in_matlab;
 }
diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.h b/include/astra/ParallelBeamBlobKernelProjector2D.h
index 529cc10..12bee5f 100644
--- a/include/astra/ParallelBeamBlobKernelProjector2D.h
+++ b/include/astra/ParallelBeamBlobKernelProjector2D.h
@@ -214,7 +214,12 @@ protected:
 	float32 m_fBlobSampleRate; //< At which interval are the inserted blob values evaluated?
 	int m_iBlobSampleCount; //< Number of evaluated blob samples
 	float32* m_pfBlobValues; //< Evaluated blob values
-	float32* m_pfBlobValuesNeg; //< Evaluated blob values
+
+	/** 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);
 
 };
 
diff --git a/include/astra/ParallelBeamBlobKernelProjector2D.inl b/include/astra/ParallelBeamBlobKernelProjector2D.inl
index c5040b7..ccd2166 100644
--- a/include/astra/ParallelBeamBlobKernelProjector2D.inl
+++ b/include/astra/ParallelBeamBlobKernelProjector2D.inl
@@ -26,186 +26,204 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 */
 
 
-
-//----------------------------------------------------------------------------------------
-// PROJECT ALL 
 template <typename Policy>
 void CParallelBeamBlobKernelProjector2D::project(Policy& p)
 {
-	for (int iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
-		for (int iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-			projectSingleRay(iAngle, iDetector, p);
-		}
-	}
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+						  0, m_pProjectionGeometry->getDetectorCount(), p);
 }
 
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
 template <typename Policy>
 void CParallelBeamBlobKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
 {
-	for (int iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		projectSingleRay(_iProjection, iDetector, p);
-	}
+	projectBlock_internal(_iProjection, _iProjection + 1,
+						  0, m_pProjectionGeometry->getDetectorCount(), p);
 }
 
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
 template <typename Policy>
 void CParallelBeamBlobKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
 {
-	ASTRA_ASSERT(m_bIsInitialized);
+	projectBlock_internal(_iProjection, _iProjection + 1,
+						  _iDetector, _iDetector + 1, p);
+}
 
-	int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
+//----------------------------------------------------------------------------------------
+// PROJECT BLOCK - vector projection geometry
+// 
+// Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+//
+// For each angle/detector pair:
+// 
+// Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+// let R=(Rx,Ry) denote the direction of the ray (vector).
+// 
+// For mainly vertical rays (|Rx|<=|Ry|), 
+// let E=(Ex,Ey) denote the centre of the most upper left pixel:
+//    E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+// and let F=(Fx,Fy) denote a vector to the next pixel 
+//    F = (PixelLengthX, 0)
+// 
+// The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+//    { Dx + a*Rx = Ex + b*Fx
+//    { Dy + a*Ry = Ey + b*Fy
+// Solving for (a,b) results in:
+//    a = (Ey + b*Fy - Dy)/Ry
+//      = (Ey - Dy)/Ry
+//    b = (Dx + a*Rx - Ex)/Fx
+//      = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+//
+// Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates. 
+//    c = b 
+//
+// The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+//    E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+// expressed in x-value pixel coordinates is
+//    c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+// And thus:
+//    deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+//                    = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+//                    = [Ey' - Ey]*(Rx/Ry)/Fx
+//                    = [Ey' - Ey]*(Rx/Ry)/Fx
+//                    = -PixelLengthY*(Rx/Ry)/Fx.
+// 
+// Given c on a certain row, its pixel directly on its left (col), and the distance (offset) to it, can be found: 
+//    col = floor(c)
+//    offset = c - col
+//
+// The index of this pixel is
+//    volumeIndex = row * colCount + col
+//
+//
+// Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+//
+//    E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+//    F = (0, -PixelLengthX)
+//
+//    a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+//    b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+//    r = b
+//    deltar = PixelLengthX*(Ry/Rx)/Fy.
+//    row = floor(r+1/2)
+//    offset = r - row
+//
+template <typename Policy>
+void CParallelBeamBlobKernelProjector2D::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);
+	}
 
-	// POLICY: RAY PRIOR
-	if (!p.rayPrior(iRayIndex)) return;
+	// precomputations
+	const float32 pixelLengthX = m_pVolumeGeometry->getPixelLengthX();
+	const float32 pixelLengthY = m_pVolumeGeometry->getPixelLengthY();	
+	const float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
+	const float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
+	const int colCount = m_pVolumeGeometry->getGridColCount();
+	const int rowCount = m_pVolumeGeometry->getGridRowCount();
+
+	// loop angles
+	for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+		// variables
+		float32 Dx, Dy, Ex, Ey, c, r, deltac, deltar, offset, invBlobExtent, RxOverRy, RyOverRx;
+		int iVolumeIndex, iRayIndex, row, col, iDetector;
+		int col_left, col_right, row_top, row_bottom, index;
+
+		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+		bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
+		if (vertical) {
+			RxOverRy = proj->fRayX/proj->fRayY;
+			deltac = -m_pVolumeGeometry->getPixelLengthY() * (proj->fRayX/proj->fRayY) * inv_pixelLengthX;
+			invBlobExtent = m_pVolumeGeometry->getPixelLengthY() / abs(m_fBlobSize * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / proj->fRayY);
+		} else {
+			RyOverRx = proj->fRayY/proj->fRayX;
+			deltar = -m_pVolumeGeometry->getPixelLengthX() * (proj->fRayY/proj->fRayX) * inv_pixelLengthY;
+			invBlobExtent = m_pVolumeGeometry->getPixelLengthX() / abs(m_fBlobSize * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / proj->fRayX);
+		}
 
-	// get values
-	float32 t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector);
-	float32 theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	if (theta >= 7*PIdiv4) theta -= 2*PI;
+		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
-	bool flip = false;
+		// loop detectors
+		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
+			
+			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
-	if (theta >= 3*PIdiv4) {
-		theta -= PI;
-		t = -t;
-		flip = true;
-	}
+			// POLICY: RAY PRIOR
+			if (!p.rayPrior(iRayIndex)) continue;
+	
+			Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+			Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
 
+			// vertically
+			if (vertical) {
 
-	if (theta <= PIdiv4) { // -pi/4 <= theta <= pi/4
+				// calculate c for row 0
+				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
-		// precalculate sin, cos, 1/cos
-		float32 sin_theta = sin(theta);
-		float32 cos_theta = cos(theta);
-		float32 inv_cos_theta = 1.0f / cos_theta; 
+				// loop rows
+				for (row = 0; row < rowCount; ++row, c += deltac) {
 
-		// precalculate other stuff
-		float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-		float32 updatePerRow = sin_theta * lengthPerRow;
-		float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-		float32 pixelLengthX_over_blobSize = m_pVolumeGeometry->getPixelLengthX() / m_fBlobSize;
-		
-		// some variables
-		int row, col, xmin, xmax;
-		float32 P, x, d;
+					col_left = int(c - 0.5f - m_fBlobSize);
+					col_right = int(c + 0.5f + m_fBlobSize);
 
-		// calculate P and x for row 0
-		P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-		x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f;
+					if (col_left < 0) col_left = 0; 
+					if (col_right > colCount-1) col_right = colCount-1; 
 
-		// for each row
-		for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-			
-			// calculate extent
-			xmin = (int)ceil((P - m_fBlobSize - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f);
-			xmax = (int)floor((P + m_fBlobSize - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX - 0.5f);
-	
-			// add pixels
-			for (col = xmin; col <= xmax; col++) {
-				if (col >= 0 && col < m_pVolumeGeometry->getGridColCount()) {
-					//d = abs(x - col) * pixelLengthX_over_blobSize;
-					//index = (int)(d*m_iBlobSampleCount+0.5f);
-					//float32 fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)] * lengthPerRow;
-
-					float32 fWeight;
-					int index;
-					if ((x >= col) ^ flip) {
-						d = abs(x - col) * pixelLengthX_over_blobSize * cos_theta;
-						index = (int)(d*m_iBlobSampleCount+0.5f);
-						fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)];
-					} else {
-						d = abs(x - col) * pixelLengthX_over_blobSize * cos_theta;
-						index = (int)(d*m_iBlobSampleCount+0.5f);
-						fWeight = m_pfBlobValuesNeg[min(index,m_iBlobSampleCount-1)];
-					}
+					// loop columns
+					for (col = col_left; col <= col_right; ++col) {
 
-					int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
-					// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-					if (p.pixelPrior(iVolumeIndex)) {
-						p.addWeight(iRayIndex, iVolumeIndex, fWeight);
-						p.pixelPosterior(iVolumeIndex);
+						iVolumeIndex = row * colCount + col;
+						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
+						if (p.pixelPrior(iVolumeIndex)) {
+							offset = abs(c - float32(col)) * invBlobExtent;
+							index = (int)(offset*m_iBlobSampleCount+0.5f);
+							p.addWeight(iRayIndex, iVolumeIndex, m_pfBlobValues[min(index,m_iBlobSampleCount-1)]);
+							p.pixelPosterior(iVolumeIndex);
+						}
 					}
 				}
 			}
 
-			// update P and x
-			P += updatePerRow;
-			x += updatePerRow * inv_pixelLengthX;
-		}
+			// horizontally
+			else {
 
-	} else { // pi/4 < theta < 3pi/4
-
-		// precalculate sin cos
-		float32 sin_90_theta = sin(PIdiv2-theta);
-		float32 cos_90_theta = cos(PIdiv2-theta);
-		float32 inv_cos_90_theta = 1.0f / cos_90_theta; 
-
-		// precalculate other stuff
-		float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_cos_90_theta;
-		float32 updatePerCol = sin_90_theta * lengthPerCol;
-		float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
-		float32 pixelLengthY_over_blobSize = m_pVolumeGeometry->getPixelLengthY() / m_fBlobSize;
-
-		// some variables
-		int row, col, xmin, xmax;
-		float32 P,x, d;
-
-		// calculate P and x for col 0
-		P = (sin_90_theta * m_pVolumeGeometry->pixelColToCenterX(0) - t) * inv_cos_90_theta;
-		x = (P - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f;
-
-		// for each col
-		for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
-			// calculate extent
-			xmin = (int)ceil((P - m_fBlobSize - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f);
-			xmax = (int)floor((P + m_fBlobSize - m_pVolumeGeometry->getWindowMinY()) * inv_pixelLengthY - 0.5f);
-
-			// add pixels
-			for (row = xmin; row <= xmax; row++) {
-				if (row >= 0 && row < m_pVolumeGeometry->getGridRowCount()) {
-					//d = abs(x - row) * pixelLengthY_over_blobSize;
-					//int index = (int)(d*m_iBlobSampleCount+0.5f);
-					//float32 fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)] * lengthPerCol;
-
-					float32 fWeight;
-					int index;
-					if ((x <= row) ^ flip) {
-						d = abs(x - row) * pixelLengthY_over_blobSize * cos_90_theta;
-						index = (int)(d*m_iBlobSampleCount+0.5f);
-						fWeight = m_pfBlobValues[min(index,m_iBlobSampleCount-1)];
-					} else {
-						d = abs(x - row) * pixelLengthY_over_blobSize * cos_90_theta;
-						index = (int)(d*m_iBlobSampleCount+0.5f);
-						fWeight = m_pfBlobValuesNeg[min(index,m_iBlobSampleCount-1)];
-					}
-
-
-					int iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
-					// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-					if (p.pixelPrior(iVolumeIndex)) {
-						p.addWeight(iRayIndex, iVolumeIndex, fWeight);
-						p.pixelPosterior(iVolumeIndex);
-					}
-				}
-			}
+				// calculate r for col 0
+				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
 
-			// update P and x
-			P += updatePerCol;
-			x += updatePerCol * inv_pixelLengthY;
-		}
+				// loop columns
+				for (col = 0; col < colCount; ++col, r += deltar) {
 
-	}
+					row_top = int(r - 0.5f - m_fBlobSize);
+					row_bottom = int(r + 0.5f + m_fBlobSize);
 
-	// POLICY: RAY POSTERIOR
-	p.rayPosterior(iRayIndex);
+					if (row_top < 0) row_top = 0; 
+					if (row_bottom > rowCount-1) row_bottom = rowCount-1; 
 
+					// loop rows
+					for (row = row_top; row <= row_bottom; ++row) {
 
+						iVolumeIndex = row * colCount + col;
+						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
+						if (p.pixelPrior(iVolumeIndex)) {
+							offset = abs(r - float32(row)) * invBlobExtent;
+							index = (int)(offset*m_iBlobSampleCount+0.5f);
+							p.addWeight(iRayIndex, iVolumeIndex, m_pfBlobValues[min(index,m_iBlobSampleCount-1)]);
+							p.pixelPosterior(iVolumeIndex);
+						}
+					}
+				}
+			}
+	
+			// POLICY: RAY POSTERIOR
+			p.rayPosterior(iRayIndex);
+	
+		} // end loop detector
+	} // end loop angles
 
 }
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.h b/include/astra/ParallelBeamLineKernelProjector2D.h
index 0f25d83..e0b7b46 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.h
+++ b/include/astra/ParallelBeamLineKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #define _INC_ASTRA_PARALLELBEAMLINEKERNELPROJECTOR
 
 #include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
 #include "Float32Data2D.h"
 #include "Projector2D.h"
 
@@ -179,6 +180,7 @@ protected:
 	template <typename Policy>
 	void projectBlock_internal(int _iProjFrom, int _iProjTo,
 	                           int _iDetFrom, int _iDetTo, Policy& _policy);
+
 };
 
 inline std::string CParallelBeamLineKernelProjector2D::getType() 
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index dcd11bd..d07f989 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -24,12 +24,13 @@ 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 CParallelBeamLineKernelProjector2D::project(Policy& p)
 {
 	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
-	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+		                  0, m_pProjectionGeometry->getDetectorCount(), p);
 }
 
 template <typename Policy>
@@ -48,235 +49,245 @@ void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int
 
 
 //----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+/* PROJECT BLOCK - vector projection geometry
+   
+   Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+  
+   For each angle/detector pair:
+   
+   Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+   let R=(Rx,Ry) denote the direction of the ray (vector).
+   
+   For mainly vertical rays (|Rx|<=|Ry|), 
+   let E=(Ex,Ey) denote the centre of the most upper left pixel:
+      E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+   and let F=(Fx,Fy) denote a vector to the next pixel 
+      F = (PixelLengthX, 0)
+   
+   The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+      { Dx + a*Rx = Ex + b*Fx
+      { Dy + a*Ry = Ey + b*Fy
+   Solving for (a,b) results in:
+      a = (Ey + b*Fy - Dy)/Ry
+        = (Ey - Dy)/Ry
+      b = (Dx + a*Rx - Ex)/Fx
+        = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+  
+   Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates. 
+      c = b 
+  
+   The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+      E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+   expressed in x-value pixel coordinates is
+      c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+   And thus:
+      deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+                      = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+                      = [Ey' - Ey]*(Rx/Ry)/Fx
+                      = [Ey' - Ey]*(Rx/Ry)/Fx
+                      = -PixelLengthY*(Rx/Ry)/Fx.
+   
+   Given c on a certain row, its closest pixel (col), and the distance (offset) to it, can be found: 
+      col = floor(c+1/2)
+      offset = c - col
+  
+   The index of this pixel is
+      volumeIndex = row * colCount + col
+  
+   The projection kernel is defined by
+  
+           _____       LengthPerRow
+         /|  |  |\
+        / |  |  | \
+     __/  |  |  |  \__ 0
+      -T -S  0  S  T
+  
+   with S = 1/2 - 1/2*|Rx/Ry|, T = 1/2 + 1/2*|Rx/Ry|, and LengthPerRow = pixelLengthX * sqrt(Rx^2+Ry^2) / |Ry|
+  
+   And thus
+                              { (offset+T)/(T-S) * LengthPerRow    if  -T <= offset < S
+      W_(rayIndex,volIndex) = { LengthPerRow                       if  -S <= offset <= S
+                              { (offset-S)/(T-S) * LengthPerRow    if   S < offset <= T
+  
+   If -T <= offset < S, the weight for the pixel directly to the left is
+      W_(rayIndex,volIndex-1) = LengthPerRow - (offset+T)/(T-S) * LengthPerRow,
+   and if S < offset <= T, the weight for the pixel directly to the right is
+      W_(rayIndex,volIndex+1) = LengthPerRow - (offset-S)/(T-S) * LengthPerRow.
+  
+  
+   Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+  
+      E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+      F = (0, -PixelLengthX)
+  
+      a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+      b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+      r = b
+      deltar = PixelLengthX*(Ry/Rx)/Fy.
+      row = floor(r+1/2)
+      offset = r - row
+      S = 1/2 - 1/2*|Ry/Rx|
+      T = 1/2 + 1/2*|Ry/Rx|
+      LengthPerCol = pixelLengthY * sqrt(Rx^2+Ry^2) / |Rx|
+  
+                              { (offset+T)/(T-S) * LengthPerCol    if  -T <= offset < S
+      W_(rayIndex,volIndex) = { LengthPerCol                       if  -S <= offset <= S
+                              { (offset-S)/(T-S) * LengthPerCol    if   S < offset <= T
+  
+      W_(rayIndex,volIndex-colcount) = LengthPerCol - (offset+T)/(T-S) * LengthPerCol
+      W_(rayIndex,volIndex+colcount) = LengthPerCol - (offset-S)/(T-S) * LengthPerCol
+*/
 template <typename Policy>
 void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
 {
-	// variables
-	float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, S, T, t, I, P, x, x2;
-	float32 lengthPerRow, updatePerRow, inv_pixelLengthX, lengthPerCol, updatePerCol, inv_pixelLengthY;
-	int iVolumeIndex, iRayIndex, row, col, iAngle, iDetector, x1;
-	bool switch_t;
+	// 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();
 
 	// loop angles
-	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
-		// get theta
-		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
-		switch_t = false;
-		if (theta >= 7*PIdiv4) theta -= 2*PI;
-		if (theta >= 3*PIdiv4) {
-			theta -= PI;
-			switch_t = true;
+	for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+		// variables
+		float32 Dx, Dy, Ex, Ey, S, T, weight, c, r, deltac, deltar, offset;
+		float32 RxOverRy, RyOverRx, lengthPerRow, lengthPerCol, invTminSTimesLengthPerRow, invTminSTimesLengthPerCol;
+		int iVolumeIndex, iRayIndex, row, col;
+
+		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
+
+		bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
+		if (vertical) {
+			RxOverRy = proj->fRayX/proj->fRayY;
+			lengthPerRow = detSize * pixelLengthX * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayY);
+			deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+			S = 0.5f - 0.5f*fabs(RxOverRy);
+			T = 0.5f + 0.5f*fabs(RxOverRy);
+			invTminSTimesLengthPerRow = lengthPerRow / (T - S);
+		} else {
+			RyOverRx = proj->fRayY/proj->fRayX;
+			lengthPerCol = detSize * pixelLengthY * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
+			deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
+			S = 0.5f - 0.5f*fabs(RyOverRx);
+			T = 0.5f + 0.5f*fabs(RyOverRx);
+			invTminSTimesLengthPerCol = lengthPerCol / (T - S);
 		}
 
-		// precalculate sin, cos, 1/cos
-		sin_theta = sin(theta);
-		cos_theta = cos(theta);
-		inv_sin_theta = 1.0f / sin_theta; 
-		inv_cos_theta = 1.0f / cos_theta; 
-
-		// precalculate kernel limits
-		lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-		updatePerRow = sin_theta * inv_cos_theta;
-		inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-
-		// precalculate kernel limits
-		lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
-		updatePerCol = cos_theta * inv_sin_theta;
-		inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
-
-		// precalculate S and T
-		S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
-		T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
+		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
 		// loop detectors
-		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-			
+		for (int iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
+
 			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
 			// POLICY: RAY PRIOR
 			if (!p.rayPrior(iRayIndex)) continue;
-	
-			// get t
-			t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
-			if (switch_t) t = -t;
 
-			// vertically
-			if (theta <= PIdiv4) {
+			Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+			Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+
+			bool isin = false;
 			
-				// calculate x for row 0
-				P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-				x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
+			// vertically
+			if (vertical) {
 
-				// get coords
-				int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-				float nextx2 = x - nextx1;
+				// calculate c for row 0
+				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
-				// for each row
-				for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
+				// loop rows
+				for (row = 0; row < rowCount; ++row, c += deltac) {
 
-					x1 = nextx1;
-					x2 = nextx2;
+					col = int(floor(c+0.5f));
+					if (col < -1 || col > colCount) { if (!isin) continue; else break; }
+					offset = c - float32(col);
 
-					nextx2 += updatePerRow;
-					while (nextx2 >= 1.0f) {
-						nextx2 -= 1.0f;
-						nextx1++;
-					}
-					while (nextx2 < 0.0f) {
-						nextx2 += 1.0f;
-						nextx1--;
-					}
+					// left
+					if (offset < -S) {
+						weight = (offset + T) * invTminSTimesLengthPerRow;
 
-					if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
+						iVolumeIndex = row * colCount + col - 1;
+						if (col > 0) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight); }
 
-					// left
-					if (x2 < 0.5f-S) {
-						I = (0.5f - S + x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
-						if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1-1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
+						iVolumeIndex++;
+						if (col >= 0 && col < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
-					// center
-					else if (x2 <= 0.5f+S) {
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}					
+					// right
+					else if (S < offset) {
+						weight = (offset - S) * invTminSTimesLengthPerRow;
+
+						iVolumeIndex = row * colCount + col;
+						if (col >= 0 && col < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow-weight); }
+
+						iVolumeIndex++;
+						if (col + 1 < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
-					// right
-					else  {
-						I = (1.5f - S - x2) / (1.0f - 2.0f*S) * lengthPerRow;
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-						if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerRow-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
+					// centre
+					else if (col >= 0 && col < colCount) {
+						iVolumeIndex = row * colCount + col;
+						policy_weight(p, iRayIndex, iVolumeIndex, lengthPerRow);
 					}
+					isin = true;
 				}
 			}
 
 			// horizontally
-			else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
+			else {
 
-				// calculate point P
-				P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-				x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY;
+				// calculate r for col 0
+				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
 
-				// get coords
-				int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-				float nextx2 = x - nextx1;
+				// loop columns
+				for (col = 0; col < colCount; ++col, r += deltar) {
 
-				// for each col
-				for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
+					row = int(floor(r+0.5f));
+					if (row < -1 || row > rowCount) { if (!isin) continue; else break; }
+					offset = r - float32(row);
 
-					x1 = nextx1;
-					x2 = nextx2;
+					// up
+					if (offset < -S) {
+						weight = (offset + T) * invTminSTimesLengthPerCol;
 
-					nextx2 += updatePerCol;
-					while (nextx2 >= 1.0f) {
-						nextx2 -= 1.0f;
-						nextx1++;
-					}
-					while (nextx2 < 0.0f) {
-						nextx2 += 1.0f;
-						nextx1--;
+						iVolumeIndex = (row-1) * colCount + col;
+						if (row > 0) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight); }
+
+						iVolumeIndex += colCount;
+						if (row >= 0 && row < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
-					if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
+					// down
+					else if (S < offset) {
+						weight = (offset - S) * invTminSTimesLengthPerCol;
 
-					// up
-					if (x2 < 0.5f-T) {
-						I = (0.5f - T + x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-						if (x1-1 >= 0 && x1-1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1-1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-					}
+						iVolumeIndex = row * colCount + col;
+						if (row >= 0 && row < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol-weight); }
 
-					// center
-					else if (x2 <= 0.5f+T) {
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}					
+						iVolumeIndex += colCount;
+						if (row + 1 < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, weight); }
 					}
 
-					// down
-					else  {
-						I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-						if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
-						if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) {
-							iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
-							// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-							if (p.pixelPrior(iVolumeIndex)) {
-								p.addWeight(iRayIndex, iVolumeIndex, lengthPerCol-I);
-								p.pixelPosterior(iVolumeIndex);
-							}
-						}
+					// centre
+					else if (row >= 0 && row < rowCount) {
+						iVolumeIndex = row * colCount + col;
+						policy_weight(p, iRayIndex, iVolumeIndex, lengthPerCol);
 					}
+					isin = true;
 				}
-			} // end loop col
+			}
 	
 			// POLICY: RAY POSTERIOR
 			p.rayPosterior(iRayIndex);
@@ -284,5 +295,8 @@ void CParallelBeamLineKernelProjector2D::projectBlock_internal(int _iProjFrom, i
 		} // end loop detector
 	} // end loop angles
 
-}
+	// Delete created vec geometry if required
+	if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+		delete pVecProjectionGeometry;
 
+}
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h
index 7b40628..3e81fa3 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.h
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #define _INC_ASTRA_PARALLELLINEARKERNELPROJECTOR
 
 #include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
 #include "Float32Data2D.h"
 #include "Projector2D.h"
 
@@ -184,7 +185,6 @@ protected:
 	void projectBlock_internal(int _iProjFrom, int _iProjTo,
 	                           int _iDetFrom, int _iDetTo, Policy& _policy);
 
-
 };
 
 //----------------------------------------------------------------------------------------
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index 5dd4781..61e4973 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -25,12 +25,13 @@ 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 CParallelBeamLinearKernelProjector2D::project(Policy& p)
 {
 	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
-	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+		                  0, m_pProjectionGeometry->getDetectorCount(), p);
 }
 
 template <typename Policy>
@@ -47,45 +48,127 @@ void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, in
 	                      _iDetector, _iDetector + 1, p);
 }
 
+
+
 //----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+/* PROJECT BLOCK - vector projection geometry
+   
+   Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+  
+   For each angle/detector pair:
+   
+   Let D=(Dx,Dy) denote the centre of the detector (point) in volume coordinates, and
+   let R=(Rx,Ry) denote the direction of the ray (vector).
+   
+   For mainly vertical rays (|Rx|<=|Ry|), 
+   let E=(Ex,Ey) denote the centre of the most upper left pixel:
+      E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+   and let F=(Fx,Fy) denote a vector to the next pixel 
+      F = (PixelLengthX, 0)
+   
+   The intersection of the ray (D+aR) with the centre line of the upper row of pixels (E+bF) is
+      { Dx + a*Rx = Ex + b*Fx
+      { Dy + a*Ry = Ey + b*Fy
+   Solving for (a,b) results in:
+      a = (Ey + b*Fy - Dy)/Ry
+        = (Ey - Dy)/Ry
+      b = (Dx + a*Rx - Ex)/Fx
+        = (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+  
+   Define c as the x-value of the intersection of the ray with the upper row in pixel coordinates. 
+      c = b 
+  
+   The intersection of the ray (D+aR) with the centre line of the second row of pixels (E'+bF) with
+      E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+   expressed in x-value pixel coordinates is
+      c' = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx.
+   And thus:
+      deltac = c' - c = (Dx + (Ey' - Dy)*Rx/Ry - Ex)/Fx - (Dx + (Ey - Dy)*Rx/Ry - Ex)/Fx
+                      = [(Ey' - Dy)*Rx/Ry - (Ey - Dy)*Rx/Ry]/Fx
+                      = [Ey' - Ey]*(Rx/Ry)/Fx
+                      = [Ey' - Ey]*(Rx/Ry)/Fx
+                      = -PixelLengthY*(Rx/Ry)/Fx.
+   
+   Given c on a certain row, its pixel directly on its left (col), and the distance (offset) to it, can be found: 
+      col = floor(c)
+      offset = c - col
+  
+   The index of this pixel is
+      volumeIndex = row * colCount + col
+  
+   The projection kernel is defined by
+  
+                  LengthPerRow
+         /|\
+        / | \
+     __/  |  \__ 0
+       p0 p1 p2
+  
+   And thus
+      W_(rayIndex,volIndex) = (1 - offset) * lengthPerRow
+      W_(rayIndex,volIndex+1) = offset * lengthPerRow
+  
+  
+   Mainly horizontal rays (|Rx|<=|Ry|) are handled in a similar fashion:
+  
+      E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+      F = (0, -PixelLengthX)
+  
+      a = (Ex + b*Fx - Dx)/Rx = (Ex - Dx)/Rx
+      b = (Dy + a*Ry - Ey)/Fy = (Dy + (Ex - Dx)*Ry/Rx - Ey)/Fy
+      r = b
+      deltar = PixelLengthX*(Ry/Rx)/Fy.
+      row = floor(r+1/2)
+      offset = r - row
+      LengthPerCol = pixelLengthY * sqrt(Rx^2+Ry^2) / |Rx|
+  
+      W_(rayIndex,volIndex) = (1 - offset) * lengthPerCol
+      W_(rayIndex,volIndex+colcount) = offset * lengthPerCol
+*/
 template <typename Policy>
 void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
 {
-	// variables
-	float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, t;
-	float32 lengthPerRow, updatePerRow;
-	float32 lengthPerCol, updatePerCol;
-	bool switch_t;
-	int iAngle, iDetector, iVolumeIndex, iRayIndex;
-	int row, col, x1;
-	float32 P,x,x2;
+	// 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();
 
 	// loop angles
-	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-
-		// get theta
-		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
-		switch_t = false;
-		if (theta >= 7*PIdiv4) theta -= 2*PI;
-		if (theta >= 3*PIdiv4) {
-			theta -= PI;
-			switch_t = true;
+	for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+		// variables
+		float32 Dx, Dy, Ex, Ey, c, r, deltac, deltar, offset;
+		float32 RxOverRy, RyOverRx, lengthPerRow, lengthPerCol;
+		int iVolumeIndex, iRayIndex, row, col, iDetector;
+
+		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+		float32 detSize = sqrt(proj->fDetUX * proj->fDetUX + proj->fDetUY * proj->fDetUY);
+
+		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);
+			deltac = -pixelLengthY * RxOverRy * inv_pixelLengthX;
+		} else {
+			RyOverRx = proj->fRayY/proj->fRayX;
+			lengthPerCol = detSize * m_pVolumeGeometry->getPixelLengthY() * sqrt(proj->fRayY*proj->fRayY + proj->fRayX*proj->fRayX) / abs(proj->fRayX);
+			deltar = -pixelLengthX * RyOverRx * inv_pixelLengthY;
 		}
 
-		// precalculate sin, cos, 1/cos
-		sin_theta = sin(theta);
-		cos_theta = cos(theta);
-		inv_cos_theta = 1.0f / cos_theta; 
-		inv_sin_theta = 1.0f / sin_theta; 
-
-		// precalculate kernel limits
-		lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-		updatePerRow = sin_theta * inv_cos_theta;
-
-		// precalculate kernel limits
-		lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
-		updatePerCol = cos_theta * inv_sin_theta;
+		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
 		// loop detectors
 		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
@@ -95,79 +178,54 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
 			// POLICY: RAY PRIOR
 			if (!p.rayPrior(iRayIndex)) continue;
 	
-			// get t
-			t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
-			if (switch_t) {
-				t = -t;
-			}
+			Dx = proj->fDetSX + (iDetector+0.5f) * proj->fDetUX;
+			Dy = proj->fDetSY + (iDetector+0.5f) * proj->fDetUY;
+
+			bool isin = false;
 
 			// vertically
-			if (theta <= PIdiv4) {
-			
-				// calculate x for row 0
-				P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-				x = m_pVolumeGeometry->coordXToColF(P) - 0.5f;
+			if (vertical) {
+
+				// calculate c for row 0
+				c = (Dx + (Ey - Dy)*RxOverRy - Ex) * inv_pixelLengthX;
 
-				// for each row
-				for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
+				// loop rows
+				for (row = 0; row < rowCount; ++row, c += deltac) {
+
+					col = int(floor(c));
+					if (col < -1 || col >= colCount) { if (!isin) continue; else break; }
+					offset = c - float32(col);
+
+					iVolumeIndex = row * colCount + col;
+					if (col >= 0) { policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - offset) * lengthPerRow); }
 					
-					// get coords
-					x1 = int((x > 0.0f) ? x : x-1.0f);
-					x2 = x - x1; 
-					x += updatePerRow;
-
-					// add weights
-					if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1);
-						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-						if (p.pixelPrior(iVolumeIndex)) {
-							p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerRow);
-							p.pixelPosterior(iVolumeIndex);
-						}
-					}
-					if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridColCount()) {
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, x1+1);
-						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-						if (p.pixelPrior(iVolumeIndex)) {
-							p.addWeight(iRayIndex, iVolumeIndex, (x2) * lengthPerRow);
-							p.pixelPosterior(iVolumeIndex);
-						}
-					}
+					iVolumeIndex++;
+					if (col + 1 < colCount) { policy_weight(p, iRayIndex, iVolumeIndex, offset * lengthPerRow); }
+
+					isin = true;
 				}
 			}
 
 			// horizontally
-			else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
-
-				// calculate point P
-				P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-				x = m_pVolumeGeometry->coordYToRowF(P) - 0.5f;
-
-				// for each row
-				for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
-					// get coords
-					x1 = int((x > 0.0f) ? x : x-1.0f);
-					x2 = x - x1; 
-					x += updatePerCol;
-
-					// add weights
-					if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridRowCount()) {
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1, col);
-						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-						if (p.pixelPrior(iVolumeIndex)) {
-							p.addWeight(iRayIndex, iVolumeIndex, (1.0f - x2) * lengthPerCol);
-							p.pixelPosterior(iVolumeIndex);		
-						}
-					}
-					if (x1+1 >= 0 && x1+1 < m_pVolumeGeometry->getGridRowCount()) {
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(x1+1, col);
-						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
-						if (p.pixelPrior(iVolumeIndex)) {
-							p.addWeight(iRayIndex, iVolumeIndex, x2 * lengthPerCol);
-							p.pixelPosterior(iVolumeIndex);
-						}
-					}
+			else {
+
+				// calculate r for col 0
+				r = -(Dy + (Ex - Dx)*RyOverRx - Ey) * inv_pixelLengthY;
+
+				// loop columns
+				for (col = 0; col < colCount; ++col, r += deltar) {
+
+					row = int(floor(r));
+					if (row < -1 || row >= rowCount) { if (!isin) continue; else break; }
+					offset = r - float32(row);
+
+					iVolumeIndex = row * colCount + col;
+					if (row >= 0) { policy_weight(p, iRayIndex, iVolumeIndex, (1.0f - offset) * lengthPerCol); }
+
+					iVolumeIndex += colCount;
+					if (row + 1 < rowCount) { policy_weight(p, iRayIndex, iVolumeIndex, offset * lengthPerCol); }
+					
+					isin = true;
 				}
 			}
 	
@@ -177,5 +235,6 @@ void CParallelBeamLinearKernelProjector2D::projectBlock_internal(int _iProjFrom,
 		} // end loop detector
 	} // end loop angles
 
+	if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+		delete pVecProjectionGeometry;
 }
-
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h
index 373c7fa..908df1f 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.h
+++ b/include/astra/ParallelBeamStripKernelProjector2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #define _INC_ASTRA_PARALLELBEAMSTROKEKERNELPROJECTOR
 
 #include "ParallelProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
 #include "Float32Data2D.h"
 #include "Projector2D.h"
 
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index e8e3739..fdfcd90 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -29,7 +29,7 @@ template <typename Policy>
 void CParallelBeamStripKernelProjector2D::project(Policy& p)
 {
 	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
-	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+		                  0, m_pProjectionGeometry->getDetectorCount(), p);
 }
 
 template <typename Policy>
@@ -47,251 +47,244 @@ void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int
 }
 
 //----------------------------------------------------------------------------------------
-// PROJECT BLOCK
+/* PROJECT BLOCK
+   
+   Kernel limitations: isotropic pixels (PixelLengthX == PixelLengthY)
+  
+   For each angle/detector pair:
+   
+   Let DL=(DLx,DLy) denote the left of the detector (point) in volume coordinates, and
+   Let DR=(DRx,DRy) denote the right of the detector (point) in volume coordinates, and
+   let R=(Rx,Ry) denote the direction of the ray (vector).
+   
+   For mainly vertical rays (|Rx|<=|Ry|), 
+   let E=(Ex,Ey) denote the centre of the most upper left pixel:
+      E = (WindowMinX +  PixelLengthX/2, WindowMaxY - PixelLengthY/2),
+   and let F=(Fx,Fy) denote a vector to the next pixel 
+      F = (PixelLengthX, 0)
+   
+   The intersection of the left edge of the strip (DL+aR) with the centre line of the upper row of pixels (E+bF) is
+      { DLx + a*Rx = Ex + b*Fx
+      { DLy + a*Ry = Ey + b*Fy
+   Solving for (a,b) results in:
+      a = (Ey + b*Fy - DLy)/Ry
+        = (Ey - DLy)/Ry
+      b = (DLx + a*Rx - Ex)/Fx
+        = (DLx + (Ey - DLy)*Rx/Ry - Ex)/Fx
+  
+   Define cL as the x-value of the intersection of the left edge of the strip with the upper row in pixel coordinates. 
+      cL = b 
+  
+   cR, the x-value of the intersection of the right edge of the strip with the upper row in pixel coordinates can be found similarly.
+  
+   The intersection of the ray (DL+aR) with the left line of the second row of pixels (E'+bF) with
+      E'=(WindowMinX + PixelLengthX/2, WindowMaxY - 3*PixelLengthY/2)
+   expressed in x-value pixel coordinates is
+      cL' = (DLx + (Ey' - DLy)*Rx/Ry - Ex)/Fx.
+   And thus:
+      deltac = cL' - cL = (DLx + (Ey' - DLy)*Rx/Ry - Ex)/Fx - (DLx + (Ey - DLy)*Rx/Ry - Ex)/Fx
+                        = [(Ey' - DLy)*Rx/Ry - (Ey - DLy)*Rx/Ry]/Fx
+                        = [Ey' - Ey]*(Rx/Ry)/Fx
+                        = [Ey' - Ey]*(Rx/Ry)/Fx
+                        = -PixelLengthY*(Rx/Ry)/Fx.
+  
+   The projection weight for a certain pixel is defined by the area between two points of 
+  
+           _____       LengthPerRow
+         /|  |  |\
+        / |  |  | \
+     __/  |  |  |  \__ 0
+      -T -S  0  S  T
+   with S = 1/2 - 1/2*|Rx/Ry|, T = 1/2 + 1/2*|Rx/Ry|, and LengthPerRow = pixelLengthX * sqrt(Rx^2+Ry^2) / |Ry|
+   
+   For a certain row, all columns that are 'hit' by this kernel lie in the interval
+      (col_left, col_right) = (floor(cL-1/2+S), floor(cR+3/2-S))
+  
+   The offsets for both is
+      (offsetL, offsetR) = (cL - floor(col_left), cR - floor(col_left))
+  
+   The projection weight is found by the difference between the integrated values of the kernel
+           offset <= -T   Kernel = 0
+      -T < offset <= -S   Kernel = PixelArea/2*(T+offset)^2/(T-S)
+      -S < offset <=  S   Kernel = PixelArea/2 + offset
+       S < offset <=  T   Kernel = PixelArea - PixelArea/2*(T-offset)^2/(T-S)
+       T <= offset:       Kernel = PixelArea
+*/
 template <typename Policy>
 void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
 {
-	ASTRA_ASSERT(m_bIsInitialized);
-
-	// Some variables
-	float32 theta, t;
-	int row, col;
-	int iAngle;
-	int iDetector;
-	float32 res;
-	float32 PL, PLimitL, PLimitR;
-	float32 xL, xR, XLimitL, XLimitR;
-	int x1L,x1R;
-	float32 x2L, x2R, updateX;
-	int iVolumeIndex, iRayIndex;
-	
-	float32 sin_theta, cos_theta, inv_sin_theta, inv_cos_theta;
-	float32 fabs_sin_theta, fabs_cos_theta, fabs_inv_sin_theta, fabs_inv_cos_theta;
-	float32 PW, PH, DW, inv_PW, inv_PH;
-	float32 S, T, U, V, inv_4T;
+	// 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 pixelArea = pixelLengthX * pixelLengthY;
+	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();
+	const int detCount = pVecProjectionGeometry->getDetectorCount();
 
 	// loop angles
-	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
-		
-		// get values
-		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
-		bool switch_t = false;
-		if (theta >= 7*PIdiv4) theta -= 2*PI;
-		if (theta >= 3*PIdiv4) {
-			theta -= PI;
-			switch_t = true;
-		}
-
-		// Precalculate sin, cos, 1/cos
-		sin_theta = sin(theta);
-		cos_theta = cos(theta);
-		inv_cos_theta = 1.0f / cos_theta; 
-		inv_sin_theta = 1.0f / sin_theta;
-
-		fabs_sin_theta = (sin_theta < 0.0f) ? -sin_theta : sin_theta;
-		fabs_cos_theta = (cos_theta < 0.0f) ? -cos_theta : cos_theta;
-		fabs_inv_cos_theta = (inv_cos_theta < 0.0f) ? -inv_cos_theta : inv_cos_theta;
-		fabs_inv_sin_theta = (inv_sin_theta < 0.0f) ? -inv_sin_theta : inv_sin_theta;
-
-		// Other precalculations
-		PW = m_pVolumeGeometry->getPixelLengthX();
-		PH = m_pVolumeGeometry->getPixelLengthY();
-		DW = m_pProjectionGeometry->getDetectorWidth();
-		inv_PW = 1.0f / PW;
-		inv_PH = 1.0f / PH;
-
-		// [-45?,45?] and [135?,225?]
-		if (theta < PIdiv4) {
-
-			// Precalculate kernel limits
-			S = -0.5f * fabs_sin_theta * fabs_inv_cos_theta;
-			T = -S;
-			U = 1.0f + S;
-			V = 1.0f - S;
-			inv_4T = 0.25f / T;
-
-			updateX = sin_theta * inv_cos_theta;
-
-			// loop detectors
-			for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-			
-				iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
-				// POLICY: RAY PRIOR
-				if (!p.rayPrior(iRayIndex)) continue;
-
-				// get t
-				t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
-				if (switch_t) t = -t;
-			
-				// calculate left strip extremes (volume coordinates)
-				PL = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0) - DW*0.5f) * inv_cos_theta;
-				PLimitL = PL - 0.5f * fabs_sin_theta * fabs_inv_cos_theta * PH;
-				PLimitR = PLimitL + DW * inv_cos_theta + PH * fabs_sin_theta * fabs_inv_cos_theta; 
-
-				// calculate strip extremes (pixel coordinates)
-				XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-				XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-				xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-				xR = xL + (DW * inv_cos_theta) * inv_PW;
-	
-				// for each row
-				for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-				
-					// get strip extremes in column indices
-					x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f);
-					x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f);
-
-					// get coords w.r.t leftmost column hit by strip
-					x2L = xL - x1L; 
-					x2R = xR - x1L;
-					
-					// update strip extremes for the next row
-					XLimitL += updateX; 
-					XLimitR += updateX;
-					xL += updateX; 
-					xR += updateX; 
-
-					// for each affected col
-					for (col = x1L; col <= x1R; ++col) {
-
-						if (col < 0 || col >= m_pVolumeGeometry->getGridColCount()) { x2L -= 1.0f; x2R -= 1.0f;	continue; }
-
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
-						// POLICY: PIXEL PRIOR
-						if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-						
-						// right
-						if (x2R >= V)		res = 1.0f;
-						else if (x2R > U)	res = x2R - (x2R-U)*(x2R-U)*inv_4T;
-						else if (x2R >= T)	res = x2R;
-						else if (x2R > S)	res = (x2R-S)*(x2R-S) * inv_4T;
-						else				{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-								
-						// left
-						if (x2L <= S)		{} // - 0.0f
-						else if (x2L < T)	res -= (x2L-S)*(x2L-S) * inv_4T;
-						else if (x2L <= U)	res -= x2L;
-						else if (x2L < V)	res -= x2L - (x2L-U)*(x2L-U)*inv_4T;
-						else				{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-
-						// POLICY: ADD
-						p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
-
-						// POLICY: PIXEL POSTERIOR
-						p.pixelPosterior(iVolumeIndex);
-
-						x2L -= 1.0f;		
-						x2R -= 1.0f;
-
-					} // end col loop
-
-				} // end row loop
-
-				// POLICY: RAY POSTERIOR
-				p.rayPosterior(iRayIndex);
-
-			}	// end detector loop
-
-		// [45?,135?] and [225?,315?]
-		// horizontaly
+	for (int iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
+
+		// variables
+		float32 DLx, DLy, DRx, DRy, Ex, Ey, S, T, deltac, deltar, offsetL, offsetR, invTminS;
+		float32 res, RxOverRy, RyOverRx, cL, cR, rL, rR;
+		int iVolumeIndex, iRayIndex, iDetector;
+		int row, row_top, row_bottom, col, col_left, col_right;
+
+		const SParProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[iAngle];
+
+		bool vertical = fabs(proj->fRayX) < fabs(proj->fRayY);
+		if (vertical) {
+			RxOverRy = proj->fRayX/proj->fRayY;
+			deltac = -m_pVolumeGeometry->getPixelLengthY() * RxOverRy * inv_pixelLengthX;
+			S = 0.5f - 0.5f*fabs(RxOverRy);
+			T = 0.5f + 0.5f*fabs(RxOverRy);
+			invTminS = 1.0f / (T-S);
 		} else {
+			RyOverRx = proj->fRayY/proj->fRayX;
+			deltar = -m_pVolumeGeometry->getPixelLengthX() * RyOverRx * inv_pixelLengthY;
+			S = 0.5f - 0.5f*fabs(RyOverRx);
+			T = 0.5f + 0.5f*fabs(RyOverRx);
+			invTminS = 1.0f / (T-S);
+		}
 
-			// Precalculate kernel limits
-			S = -0.5f * fabs_cos_theta * fabs_inv_sin_theta;
-			T = -S;
-			U = 1.0f + S;
-			V = 1.0f - S;
-			inv_4T = 0.25f / T;
-
-			updateX = cos_theta * inv_sin_theta;
-
-			// loop detectors
-			for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
-			
-				iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
-				// POLICY: RAY PRIOR
-				if (!p.rayPrior(iRayIndex)) continue;
+		Ex = m_pVolumeGeometry->getWindowMinX() + pixelLengthX*0.5f;
+		Ey = m_pVolumeGeometry->getWindowMaxY() - pixelLengthY*0.5f;
 
-				// get t
-				t = m_pProjectionGeometry->indexToDetectorOffset(iDetector);
-				if (switch_t) t = -t;
+		// loop detectors
+		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
-				// calculate left strip extremes (volume coordinates)
-				PL = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0) + DW*0.5f) * inv_sin_theta;
-				PLimitL = PL + 0.5f * fabs_cos_theta * fabs_inv_sin_theta * PW;
-				PLimitR = PLimitL - DW * inv_sin_theta - PH * fabs_cos_theta * fabs_inv_sin_theta; 
-
-				// calculate strip extremes (pixel coordinates)
-				XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH;
-				XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH;
-				xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH;
-				xR = xL + (DW * fabs_inv_sin_theta) * inv_PH;
-
-				// for each col
-				for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
-					// get strip extremes in column indices
-					x1L = int((XLimitL > 0.0f) ? XLimitL : XLimitL-1.0f);
-					x1R = int((XLimitR > 0.0f) ? XLimitR : XLimitR-1.0f);
-
-					// get coords w.r.t leftmost column hit by strip
-					x2L = xL - x1L; 
-					x2R = xR - x1L;
-					
-					// update strip extremes for the next row
-					XLimitL += updateX; 
-					XLimitR += updateX;
-					xL += updateX; 
-					xR += updateX; 
+			iRayIndex = iAngle * detCount + iDetector;
 
-					// for each affected col
-					for (row = x1L; row <= x1R; ++row) {
-
-						if (row < 0 || row >= m_pVolumeGeometry->getGridRowCount()) { x2L -= 1.0f; x2R -= 1.0f;	continue; }
-
-						iVolumeIndex = m_pVolumeGeometry->pixelRowColToIndex(row, col);
-
-						// POLICY: PIXEL PRIOR
-						if (!p.pixelPrior(iVolumeIndex)) { x2L -= 1.0f; x2R -= 1.0f; continue; }
-
-						// right
-						if (x2R >= V)		res = 1.0f;
-						else if (x2R > U)	res = x2R - (x2R-U)*(x2R-U)*inv_4T;
-						else if (x2R >= T)	res = x2R;
-						else if (x2R > S)	res = (x2R-S)*(x2R-S) * inv_4T;
-						else				{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-								
-						// left
-						if (x2L <= S)		{} // - 0.0f
-						else if (x2L < T)	res -= (x2L-S)*(x2L-S) * inv_4T;
-						else if (x2L <= U)	res -= x2L;
-						else if (x2L < V)	res -= x2L - (x2L-U)*(x2L-U)*inv_4T;
-						else				{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-
-						// POLICY: ADD
-						p.addWeight(iRayIndex, iVolumeIndex, PW*PH * res);
-
-						// POLICY: PIXEL POSTERIOR
-						p.pixelPosterior(iVolumeIndex);
-
-						x2L -= 1.0f;		
-						x2R -= 1.0f;
-
-					} // end row loop
-
-				} // end col loop
-
-				// POLICY: RAY POSTERIOR
-				p.rayPosterior(iRayIndex);
-
-			} // end detector loop
-
-
-		} // end theta switch
+			// POLICY: RAY PRIOR
+			if (!p.rayPrior(iRayIndex)) continue;
+	
+			DLx = proj->fDetSX + iDetector * proj->fDetUX;
+			DLy = proj->fDetSY + iDetector * proj->fDetUY;
+			DRx = DLx + proj->fDetUX;
+			DRy = DLy + proj->fDetUY;
+
+			// vertically
+			if (vertical) {
+
+				// calculate cL and cR for row 0
+				cL = (DLx + (Ey - DLy)*RxOverRy - Ex) * inv_pixelLengthX;
+				cR = (DRx + (Ey - DRy)*RxOverRy - Ex) * inv_pixelLengthX;
+
+				if (cR < cL) {
+					float32 tmp = cL;
+					cL = cR;
+					cR = tmp;
+				}
+
+				// loop rows
+				for (row = 0; row < rowCount; ++row, cL += deltac, cR += deltac) {
+
+					col_left = int(cL-0.5f+S);
+					col_right = int(cR+1.5-S);
+
+					if (col_left < 0) col_left = 0; 
+					if (col_right > colCount-1) col_right = colCount-1; 
+
+					float32 tmp = float32(col_left);
+					offsetL = cL - tmp;
+					offsetR = cR - tmp;
+
+					// loop columns
+					for (col = col_left; col <= col_right; ++col, offsetL -= 1.0f, offsetR -= 1.0f) {
+
+						iVolumeIndex = row * colCount + col;
+						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
+						if (p.pixelPrior(iVolumeIndex)) {
+
+							// right ray edge
+							if (T <= offsetR)       res = 1.0f;
+							else if (S < offsetR)   res = 1.0f - 0.5f*(T-offsetR)*(T-offsetR)*invTminS;
+							else if (-S < offsetR)  res = 0.5f + offsetR;
+							else if (-T < offsetR)  res = 0.5f*(offsetR+T)*(offsetR+T)*invTminS;
+							else                    res = 0.0f;
+
+							// left ray edge
+							if (T <= offsetL)       res -= 1.0f;
+							else if (S < offsetL)   res -= 1.0f - 0.5f*(T-offsetL)*(T-offsetL)*invTminS;
+							else if (-S < offsetL)  res -= 0.5f + offsetL;
+							else if (-T < offsetL)  res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
+
+							p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+							p.pixelPosterior(iVolumeIndex);
+						}
+					}
+				}
+			}
+
+			// horizontally
+			else {
+
+				// calculate rL and rR for row 0
+				rL = -(DLy + (Ex - DLx)*RyOverRx - Ey) * inv_pixelLengthY;
+				rR = -(DRy + (Ex - DRx)*RyOverRx - Ey) * inv_pixelLengthY;
+
+				if (rR < rL) {
+					float32 tmp = rL;
+					rL = rR;
+					rR = tmp;
+				}
+
+				// loop columns
+				for (col = 0; col < colCount; ++col, rL += deltar, rR += deltar) {
+
+					row_top = int(rL-0.5f+S);
+					row_bottom = int(rR+1.5-S);
+
+					if (row_top < 0) row_top = 0; 
+					if (row_bottom > rowCount-1) row_bottom = rowCount-1; 
+
+					float32 tmp = float32(row_top);
+					offsetL = rL - tmp;
+					offsetR = rR - tmp;
+
+					// loop rows
+					for (row = row_top; row <= row_bottom; ++row, offsetL -= 1.0f, offsetR -= 1.0f) {
+
+						iVolumeIndex = row * colCount + col;
+						// POLICY: PIXEL PRIOR + ADD + POSTERIOR
+						if (p.pixelPrior(iVolumeIndex)) {
+
+							// right ray edge
+							if (T <= offsetR)       res = 1.0f;
+							else if (S < offsetR)   res = 1.0f - 0.5f*(T-offsetR)*(T-offsetR)*invTminS;
+							else if (-S < offsetR)  res = 0.5f + offsetR;
+							else if (-T < offsetR)  res = 0.5f*(offsetR+T)*(offsetR+T)*invTminS;
+							else                    res = 0.0f;
+
+							// left ray edge
+							if (T <= offsetL)       res -= 1.0f;
+							else if (S < offsetL)   res -= 1.0f - 0.5f*(T-offsetL)*(T-offsetL)*invTminS;
+							else if (-S < offsetL)  res -= 0.5f + offsetL;
+							else if (-T < offsetL)  res -= 0.5f*(offsetL+T)*(offsetL+T)*invTminS;
+
+							p.addWeight(iRayIndex, iVolumeIndex, pixelArea*res);
+							p.pixelPosterior(iVolumeIndex);
+						}
+					}
+				}
+			}
+	
+			// POLICY: RAY POSTERIOR
+			p.rayPosterior(iRayIndex);
+	
+		} // end loop detector
+	} // end loop angles
 
-	} // end angle loop
+	if (dynamic_cast<CParallelProjectionGeometry2D*>(m_pProjectionGeometry))
+		delete pVecProjectionGeometry;
 }
-
-
diff --git a/include/astra/ParallelProjectionGeometry2D.h b/include/astra/ParallelProjectionGeometry2D.h
index 963bb06..9625d35 100644
--- a/include/astra/ParallelProjectionGeometry2D.h
+++ b/include/astra/ParallelProjectionGeometry2D.h
@@ -29,6 +29,7 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 #define _INC_ASTRA_PARALLELPROJECTIONGEOMETRY2D
 
 #include "ProjectionGeometry2D.h"
+#include "ParallelVecProjectionGeometry2D.h"
 
 namespace astra
 {
@@ -82,8 +83,7 @@ public:
 	CParallelProjectionGeometry2D(int _iProjectionAngleCount, 
 								  int _iDetectorCount, 
 								  float32 _fDetectorWidth, 
-								  const float32* _pfProjectionAngles,
-								  const float32* _pfExtraDetectorOffsets = 0);
+								  const float32* _pfProjectionAngles);
 
 	/** Copy constructor. 
 	 */
@@ -115,8 +115,7 @@ public:
 	bool initialize(int _iProjectionAngleCount, 
 					int _iDetectorCount, 
 					float32 _fDetectorWidth, 
-					const float32* _pfProjectionAngles,
-					const float32* _pfExtraDetectorOffsets = 0);
+					const float32* _pfProjectionAngles);
 
 	/** Create a hard copy. 
 	*/
@@ -133,7 +132,7 @@ public:
 	 * @param _sType geometry type to compare to.
 	 * @return true if _sType == "parallel".
 	 */
-	 virtual bool isOfType(const std::string& _sType);
+	virtual bool isOfType(const std::string& _sType);
 
 	/** Get all settings in a Config object.
 	 *
@@ -150,7 +149,12 @@ public:
 	 *
 	 * @return a unit vector describing the direction
 	 */
-	 virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex = 0);
+	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex = 0);
+
+	/** Create a vector geom
+	*/
+	CParallelVecProjectionGeometry2D* toVectorGeometry();
+
 };
 
 } // namespace astra
diff --git a/include/astra/ParallelVecProjectionGeometry2D.h b/include/astra/ParallelVecProjectionGeometry2D.h
new file mode 100644
index 0000000..96f8a54
--- /dev/null
+++ b/include/astra/ParallelVecProjectionGeometry2D.h
@@ -0,0 +1,163 @@
+/*
+-----------------------------------------------------------------------
+Copyright: 2010-2015, iMinds-Vision Lab, University of Antwerp
+           2014-2015, CWI, Amsterdam
+
+Contact: astra@uantwerpen.be
+Website: http://sf.net/projects/astra-toolbox
+
+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/>.
+
+-----------------------------------------------------------------------
+$Id$
+*/
+
+#ifndef _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D
+#define _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D
+
+#include "ProjectionGeometry2D.h"
+#include "GeometryUtil2D.h"
+
+namespace astra
+{
+
+/**
+ * This class defines a 2D parallel beam projection geometry. 
+ *
+ * \par XML Configuration
+ * \astra_xml_item{DetectorCount, int, Number of detectors for each projection.}
+ *
+ * \par MATLAB example
+ * \astra_code{
+ *		proj_geom = astra_struct('parallel_vec');\n
+ *		proj_geom.DetectorCount = 512;\n
+ *		proj_geom.Vectors = V;\n
+ * }
+ *
+ * \par Vectors
+ * Vectors is a matrix containing the actual geometry. Each row corresponds
+ * to a single projection, and consists of:
+ * ( rayX, rayY, dX, dY, uX, uY)
+ *      ray: the ray direction
+ *      d  : the centre of the detector line
+ *      u  : the vector from detector pixel (0) to (1)
+ */
+class _AstraExport CParallelVecProjectionGeometry2D : public CProjectionGeometry2D
+{
+protected:
+
+	SParProjection *m_pProjectionAngles;
+
+public:
+
+	/** Default constructor. Sets all variables to zero. Note that this constructor leaves the object in an unusable state and must
+	 * be followed by a call to init(). 
+	 */
+	CParallelVecProjectionGeometry2D();
+
+	/** Constructor.
+	 *
+	 * @param _iProjectionAngleCount Number of projection angles.
+	 * @param _iDetectorCount Number of detectors, i.e., the number of detector measurements for each projection angle.
+	 * @param _pfProjectionAngles Pointer to an array of projection angles. The angles will be copied from this array. 
+	 */
+	CParallelVecProjectionGeometry2D(int _iProjectionAngleCount, 
+	                                int _iDetectorCount, 
+	                                const SParProjection* _pfProjectionAngles);
+
+	/** Copy constructor. 
+	 */
+	CParallelVecProjectionGeometry2D(const CParallelVecProjectionGeometry2D& _projGeom);
+
+	/** Assignment operator.
+	 */
+	CParallelVecProjectionGeometry2D& operator=(const CParallelVecProjectionGeometry2D& _other);
+
+	/** Destructor.
+	 */
+	virtual ~CParallelVecProjectionGeometry2D();
+
+	/** Initialize the geometry with a config object.
+	 *
+	 * @param _cfg Configuration Object
+	 * @return initialization successful?
+	 */
+	virtual bool initialize(const Config& _cfg);
+
+	/** Initialization. This function MUST be called after using the default constructor and MAY be called to 
+	 * reset a previously initialized object.
+	 *
+	 * @param _iProjectionAngleCount Number of projection angles.
+	 * @param _iDetectorCount Number of detectors, i.e., the number of detector measurements for each projection angle.
+	 * @param _pfProjectionAngles Pointer to an array of projection angles. The angles will be copied from this array.
+	 */
+	bool initialize(int _iProjectionAngleCount, 
+	                int _iDetectorCount, 
+	                const SParProjection* _pfProjectionAngles);
+
+	virtual bool _check();
+
+	/** Create a hard copy. 
+	*/
+	virtual CProjectionGeometry2D* clone();
+
+	/** Returns true if the type of geometry defined in this class is the one specified in _sType.
+	 *
+	 * @param _sType geometry type to compare to.
+	 * @return true if _sType == "fanflat_vec".
+	 */
+	 virtual bool isOfType(const std::string& _sType);
+
+    /** Return true if this geometry instance is the same as the one specified.
+	 *
+	 * @return true if this geometry instance is the same as the one specified.
+	 */
+	virtual bool isEqual(CProjectionGeometry2D*) const;
+
+	/** Get all settings in a Config object.
+	 *
+	 * @return Configuration Object.
+	 */
+	virtual Config* getConfiguration() const;
+
+
+	/** Get the value for t and theta, based upon the row and column index.
+	 *
+	 * @param _iRow		row index 
+	 * @param _iColumn	column index
+	 * @param _fT		output: value of t
+	 * @param _fTheta	output: value of theta, always lies within the [0,pi[ interval.
+	 */
+	virtual void getRayParams(int _iRow, int _iColumn, float32& _fT, float32& _fTheta) const;
+
+	/**
+	 * Returns a vector describing the direction of a ray belonging to a certain detector
+	 *
+	 * @param _iProjectionIndex index of projection
+	 * @param _iProjectionIndex index of detector
+	 *
+	 * @return a unit vector describing the direction
+	 */
+	virtual CVector3D getProjectionDirection(int _iProjectionIndex, int _iDetectorIndex);
+
+	const SParProjection* getProjectionVectors() const { return m_pProjectionAngles; }
+
+};
+
+} // namespace astra
+
+#endif /* _INC_ASTRA_PARALLELVECPROJECTIONGEOMETRY2D */
diff --git a/include/astra/ProjectionGeometry2D.h b/include/astra/ProjectionGeometry2D.h
index 393db77..504f588 100644
--- a/include/astra/ProjectionGeometry2D.h
+++ b/include/astra/ProjectionGeometry2D.h
@@ -64,10 +64,6 @@ protected:
 	 */
 	float32 m_fDetectorWidth;
 
-	/** An array of m_iProjectionAngleCount elements containing an extra detector offset for each projection.
-	 */
-	float32* m_pfExtraDetectorOffset;
-
 	/** Dynamically allocated array of projection angles. All angles are represented in radians and lie in 
 	 * the [0,2pi[ interval.
 	 */
@@ -93,8 +89,7 @@ protected:
 	CProjectionGeometry2D(int _iProjectionAngleCount, 
 						  int _iDetectorCount, 
 						  float32 _fDetectorWidth, 
-						  const float32* _pfProjectionAngles,
-						  const float32* _pfExtraDetectorOffsets = 0);
+						  const float32* _pfProjectionAngles);
 
 	/** Copy constructor. 
 	 */
@@ -120,8 +115,7 @@ protected:
 	bool _initialize(int _iProjectionAngleCount,
 					 int _iDetectorCount,
 					 float32 _fDetectorWidth,
-					 const float32* _pfProjectionAngles,
-					 const float32* _pfExtraDetectorOffsets = 0);
+					 const float32* _pfProjectionAngles);
 
 public:
 
@@ -201,9 +195,6 @@ public:
 	 */
 	float32 getProjectionAngleDegrees(int _iProjectionIndex) const;
 
-	float32 getExtraDetectorOffset(int iAngle) const;
-	const float32* getExtraDetectorOffset() const { return m_pfExtraDetectorOffset; }
-
 	/** Get the index coordinate of a point on a detector array.
 	 *
 	 * @param _fOffset	distance between the center of the detector array and a certain point
@@ -272,12 +263,6 @@ public:
 //----------------------------------------------------------------------------------------
 
 
-inline float32 CProjectionGeometry2D::getExtraDetectorOffset(int _iAngle) const
-{
-	return m_pfExtraDetectorOffset ? m_pfExtraDetectorOffset[_iAngle] : 0.0f;
-}
-
-
 // Get the initialization state.
 inline bool CProjectionGeometry2D::isInitialized() const
 {