Remove code duplication in templated projectors

10 files changed, 158 insertions, 2274 deletions

diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.h b/include/astra/FanFlatBeamLineKernelProjector2D.h
index 154504a..8ead007 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.h
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.h
@@ -178,6 +178,13 @@ public:
 
 	float32 angleBetweenVectors(float32 _fAX, float32 _fAY, float32 _fBX, float32 _fBY);
 
+protected:
+	/** Internal policy-based projection of a range of angles and range.
+ 	 * (_i*From is inclusive, _i*To exclusive) */
+	template <typename Policy>
+	void projectBlock_internal(int _iProjFrom, int _iProjTo,
+	                           int _iDetFrom, int _iDetTo, Policy& _policy);
+
 };
 
 //----------------------------------------------------------------------------------------
diff --git a/include/astra/FanFlatBeamLineKernelProjector2D.inl b/include/astra/FanFlatBeamLineKernelProjector2D.inl
index e686532..7ea2749 100644
--- a/include/astra/FanFlatBeamLineKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamLineKernelProjector2D.inl
@@ -29,11 +29,32 @@ $Id$
 
 using namespace astra;
 
-//----------------------------------------------------------------------------------------
-// PROJECT ALL 
 template <typename Policy>
 void CFanFlatBeamLineKernelProjector2D::project(Policy& p)
 {
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CFanFlatBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      _iDetector, _iDetector + 1, p);
+}
+
+//----------------------------------------------------------------------------------------
+// PROJECT BLOCK
+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;
@@ -47,7 +68,7 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p)
 	const SFanProjection * proj = 0;
 
 	// loop angles
-	for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
+	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 
 		// get theta
 		if (pProjectionGeometry) {
@@ -69,7 +90,7 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p)
 		}
 
 		// loop detectors
-		for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
+		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
 			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
@@ -273,468 +294,3 @@ void CFanFlatBeamLineKernelProjector2D::project(Policy& p)
 		} // end loop detector
 	} // end loop angles
 }
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
-template <typename Policy>
-void CFanFlatBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, 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, iDetector, x1;
-	bool switch_t;
-
-	const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
-	const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry);
-
-	float32 old_theta, theta, alpha;
-	const SFanProjection * proj = 0;
-
-	//get theta
-	if (pProjectionGeometry) {
-		old_theta = pProjectionGeometry->getProjectionAngle(_iProjection);
-	}
-	else if (pVecProjectionGeometry) {
-		proj = &pVecProjectionGeometry->getProjectionVectors()[_iProjection];
-		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;
-	}
-	
-	// loop detectors
-	for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		
-		iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
-		// POLICY: RAY PRIOR
-		if (!p.rayPrior(iRayIndex)) continue;
-
-		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();
-
-		// precalculate S and T
-		S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
-		T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
-
-		// 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;
-
-			// 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;
-
-				if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
-
-				// 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);
-						}
-					}
-				}
-
-				// 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 (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);
-						}
-					}
-				}
-			}
-		}
-
-		// horizontally
-		else {
-
-			// calculate point P
-			P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-			x = (m_pVolumeGeometry->getWindowMaxY() - P) * 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;
-
-				if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
-
-				// 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);
-						}
-					}
-				}
-
-				// 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);
-						}
-					}					
-				}
-
-				// 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);
-						}
-					}
-				}
-			}
-		} // end loop col
-
-		// POLICY: RAY POSTERIOR
-		p.rayPosterior(iRayIndex);
-
-	} // end loop detector
-
-}
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
-template <typename Policy>
-void CFanFlatBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, 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, x1;
-	bool switch_t;
-
-	const CFanFlatProjectionGeometry2D* pProjectionGeometry = dynamic_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
-	const CFanFlatVecProjectionGeometry2D* pVecProjectionGeometry = dynamic_cast<CFanFlatVecProjectionGeometry2D*>(m_pProjectionGeometry);
-
-	iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-	// POLICY: RAY PRIOR
-	if (!p.rayPrior(iRayIndex)) return;
-
-	float32 old_theta, theta, alpha;
-	if (pProjectionGeometry) {
-		old_theta = pProjectionGeometry->getProjectionAngle(_iProjection);
-		t = -pProjectionGeometry->indexToDetectorOffset(_iDetector);
-		alpha = atan(t / pProjectionGeometry->getSourceDetectorDistance());
-		t = sin(alpha) * pProjectionGeometry->getOriginSourceDistance();
-	}
-	else if (pVecProjectionGeometry) {
-		const SFanProjection * proj = &pVecProjectionGeometry->getProjectionVectors()[_iProjection];
-		old_theta = atan2(-proj->fSrcX, proj->fSrcY);
-		if (old_theta < 0) old_theta += 2*PI;
-
-		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);
-	}
-
-	switch_t = false;
-	if (old_theta >= 7*PIdiv4) old_theta -= 2*PI;
-	if (old_theta >= 3*PIdiv4) {
-		old_theta -= PI;
-		switch_t = true;
-	}
-	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();
-
-	// precalculate S and T
-	S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
-	T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
-
-	// 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;
-
-		// 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;
-
-			if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
-
-			// 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);
-					}
-				}
-			}
-
-			// 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 (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);
-					}
-				}
-			}
-		}
-	}
-
-	// horizontally
-	else {
-
-		// calculate point P
-		P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-		x = (m_pVolumeGeometry->getWindowMaxY() - P) * 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;
-
-			if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
-
-			// 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);
-					}
-				}
-			}
-
-			// 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);
-					}
-				}					
-			}
-
-			// 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);
-					}
-				}
-			}
-		}
-	} // end loop col
-
-	// POLICY: RAY POSTERIOR
-	p.rayPosterior(iRayIndex);
-}
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.h b/include/astra/FanFlatBeamStripKernelProjector2D.h
index 7fed4c8..6b95427 100644
--- a/include/astra/FanFlatBeamStripKernelProjector2D.h
+++ b/include/astra/FanFlatBeamStripKernelProjector2D.h
@@ -174,6 +174,12 @@ public:
 	 */
 	virtual std::string getType();
 
+protected:
+	/** Internal policy-based projection of a range of angles and range.
+ 	 * (_i*From is inclusive, _i*To exclusive) */
+	template <typename Policy>
+	void projectBlock_internal(int _iProjFrom, int _iProjTo,
+	                           int _iDetFrom, int _iDetTo, Policy& _policy);
 };
 
 //----------------------------------------------------------------------------------------
diff --git a/include/astra/FanFlatBeamStripKernelProjector2D.inl b/include/astra/FanFlatBeamStripKernelProjector2D.inl
index e3f8b29..f95a065 100644
--- a/include/astra/FanFlatBeamStripKernelProjector2D.inl
+++ b/include/astra/FanFlatBeamStripKernelProjector2D.inl
@@ -29,11 +29,32 @@ $Id$
 
 using namespace astra;
 
-//----------------------------------------------------------------------------------------
-// PROJECT ALL 
 template <typename Policy>
 void CFanFlatBeamStripKernelProjector2D::project(Policy& p)
 {
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CFanFlatBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CFanFlatBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      _iDetector, _iDetector + 1, p);
+}
+
+//----------------------------------------------------------------------------------------
+// PROJECT BLOCK
+template <typename Policy>
+void CFanFlatBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
+{
 	ASTRA_ASSERT(m_bIsInitialized);
 
 	// Some variables
@@ -65,7 +86,7 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p)
 	}
 
 	// loop angles
-	for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
+	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 		
 		// get values
 		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
@@ -84,7 +105,7 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p)
 		if (theta < PIdiv4) {
 
 			// loop detectors
-			for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
+			for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 
 				iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
@@ -345,617 +366,3 @@ void CFanFlatBeamStripKernelProjector2D::project(Policy& p)
 	delete[] sin_alpha;
 }
 
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
-template <typename Policy>
-void CFanFlatBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
-{
-	ASTRA_ASSERT(m_bIsInitialized);
-
-	// Some variables
-	float32 theta;
-	int row, col;
-	int iDetector;
-	float32 res;
-	int x1L, x1R;
-	float32 x2L, x2R;
-	int iVolumeIndex, iRayIndex;
-	
-	CFanFlatProjectionGeometry2D* projgeom = static_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
-
-	// Other precalculations
-	float32 PW = m_pVolumeGeometry->getPixelLengthX();
-	float32 PH = m_pVolumeGeometry->getPixelLengthY();
-	float32 DW = m_pProjectionGeometry->getDetectorWidth();
-	float32 inv_PW = 1.0f / PW;
-	float32 inv_PH = 1.0f / PH;
-
-	// calculate alpha's
-	float32 alpha;
-	float32* cos_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1];
-	float32* sin_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1];
-	for (int i = 0; i < m_pProjectionGeometry->getDetectorCount() + 1; ++i) {
-		alpha = -atan((i - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW / projgeom->getSourceDetectorDistance());
-		cos_alpha[i] = cos(alpha);
-		sin_alpha[i] = sin(alpha);
-	}
-
-		
-	// get values
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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
-	float32 sin_theta = sin(theta);
-	float32 cos_theta = cos(theta);
-
-	// [-45?,45?] and [135?,225?]
-	if (theta < PIdiv4) {
-
-		// loop detectors
-		for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-
-			iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
-			// POLICY: RAY PRIOR
-			if (!p.rayPrior(iRayIndex)) continue;
-
-			float32 sin_theta_left, cos_theta_left;
-			float32 sin_theta_right, cos_theta_right;
-
-			// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
-			float32 t_l, t_r;
-			if (!switch_t) {
-				sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1];
-				sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector];
-
-				cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1];
-				cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector];
-
-				t_l = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance();
-				t_r = sin_alpha[iDetector] * projgeom->getOriginSourceDistance();
-
-			} else {
-				sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector];
-				sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1];
-
-				cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector];
-				cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1];
-
-				t_l = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance();
-				t_r = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance();	
-			}
-
-			float32 inv_cos_theta_left = 1.0f / cos_theta_left; 
-			float32 inv_cos_theta_right = 1.0f / cos_theta_right; 
-
-			float32 updateX_left = sin_theta_left * inv_cos_theta_left;
-			float32 updateX_right = sin_theta_right * inv_cos_theta_right;
-
-			// Precalculate kernel limits
-			float32 S_l = -0.5f * updateX_left;
-			if (S_l > 0) {S_l = -S_l;}
-			float32 T_l = -S_l;
-			float32 U_l = 1.0f + S_l;
-			float32 V_l = 1.0f - S_l;
-			float32 inv_4T_l = 0.25f / T_l;
-
-			float32 S_r = -0.5f * updateX_right;
-			if (S_r > 0) {S_r = -S_r;}
-			float32 T_r = -S_r;
-			float32 U_r = 1.0f + S_r;
-			float32 V_r = 1.0f - S_r;
-			float32 inv_4T_r = 0.25f / T_r;
-
-			// calculate strip extremes (volume coordinates)
-			float32 PL = (t_l - sin_theta_left * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_left;
-			float32 PR = (t_r - sin_theta_right * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_right;
-			float32 PLimitL = PL + S_l * PH;
-			float32 PLimitR = PR - S_r * PH;
-			
-			// calculate strip extremes (pixel coordinates)
-			float32 XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-			float32 XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-			float32 xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-			float32 xR = (PR - m_pVolumeGeometry->getWindowMinX()) * 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_left; 
-				XLimitR += updateX_right;
-				xL += updateX_left; 
-				xR += updateX_right; 
-
-				// 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_r)			res = 1.0f;
-					else if (x2R > U_r)		res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
-					else if (x2R >= T_r)	res = x2R;
-					else if (x2R > S_r)		res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
-					else					{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-							
-					// left
-					if (x2L <= S_l)			{}
-					else if (x2L < T_l)		res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
-					else if (x2L <= U_l)	res -= x2L;
-					else if (x2L < V_l)		res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
-					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
-	} else {
-
-		// loop detectors
-		for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-
-			iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + iDetector;
-
-			// POLICY: RAY PRIOR
-			if (!p.rayPrior(iRayIndex)) continue;
-
-			// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
-			float32 sin_theta_left, cos_theta_left;
-			float32 sin_theta_right, cos_theta_right;
-			float32 t_l, t_r;
-			if (!switch_t) {
-				sin_theta_left = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector];
-				sin_theta_right = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1];
-
-				cos_theta_left = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector];
-				cos_theta_right = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1];
-
-				t_l = sin_alpha[iDetector] * projgeom->getOriginSourceDistance();
-				t_r = sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance();
-
-			} else {
-				sin_theta_left = sin_theta * cos_alpha[iDetector+1] + cos_theta * sin_alpha[iDetector+1];
-				sin_theta_right = sin_theta * cos_alpha[iDetector] + cos_theta * sin_alpha[iDetector];
-
-				cos_theta_left = cos_theta * cos_alpha[iDetector+1] - sin_theta * sin_alpha[iDetector+1];
-				cos_theta_right = cos_theta * cos_alpha[iDetector] - sin_theta * sin_alpha[iDetector];
-
-				t_l = -sin_alpha[iDetector+1] * projgeom->getOriginSourceDistance();
-				t_r = -sin_alpha[iDetector] * projgeom->getOriginSourceDistance();	
-			}
-
-			float32 inv_sin_theta_left = 1.0f / sin_theta_left;
-			float32 inv_sin_theta_right = 1.0f / sin_theta_right;
-
-			float32 updateX_left = cos_theta_left * inv_sin_theta_left;
-			float32 updateX_right = cos_theta_right * inv_sin_theta_right;
-
-			// Precalculate kernel limits
-			float32 S_l = -0.5f * updateX_left;
-			if (S_l > 0) { S_l = -S_l; }
-			float32 T_l = -S_l;
-			float32 U_l = 1.0f + S_l;
-			float32 V_l = 1.0f - S_l;
-			float32 inv_4T_l = 0.25f / T_l;
-
-			float32 S_r = -0.5f * updateX_right;
-			if (S_r > 0) { S_r = -S_r; }
-			float32 T_r = -S_r;
-			float32 U_r = 1.0f + S_r;
-			float32 V_r = 1.0f - S_r;
-			float32 inv_4T_r = 0.25f / T_r;
-
-			// calculate strip extremes (volume coordinates)
-			float32 PL = (t_l - cos_theta_left * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_left;
-			float32 PR = (t_r - cos_theta_right * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_right;
-			float32 PLimitL = PL - S_l * PW;
-			float32 PLimitR = PR + S_r * PW;
-			
-			// calculate strip extremes (pixel coordinates)
-			float32 XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH;
-			float32 XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH;
-			float32 xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH;
-			float32 xR = (m_pVolumeGeometry->getWindowMaxY() - PR) * 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_left; 
-				XLimitR += updateX_right;
-				xL += updateX_left; 
-				xR += updateX_right; 
-
-				// for each affected row
-				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_r)			res = 1.0f;
-					else if (x2R > U_r)		res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
-					else if (x2R >= T_r)	res = x2R;
-					else if (x2R > S_r)		res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
-					else					{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-							
-					// left
-					if (x2L <= S_l)			{}
-					else if (x2L < T_l)		res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
-					else if (x2L <= U_l)	res -= x2L;
-					else if (x2L < V_l)		res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
-					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
-
-	} // end theta switch
-
-	delete[] cos_alpha;
-	delete[] sin_alpha;
-}
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
-template <typename Policy>
-void CFanFlatBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
-{
-	ASTRA_ASSERT(m_bIsInitialized);
-
-	// Some variables
-	float32 theta;
-	int row, col;
-	float32 res;
-	int x1L, x1R;
-	float32 x2L, x2R;
-	int iVolumeIndex, iRayIndex;
-	
-	CFanFlatProjectionGeometry2D* projgeom = static_cast<CFanFlatProjectionGeometry2D*>(m_pProjectionGeometry);
-
-	// Other precalculations
-	float32 PW = m_pVolumeGeometry->getPixelLengthX();
-	float32 PH = m_pVolumeGeometry->getPixelLengthY();
-	float32 DW = m_pProjectionGeometry->getDetectorWidth();
-	float32 inv_PW = 1.0f / PW;
-	float32 inv_PH = 1.0f / PH;
-
-	// calculate alpha's
-	float32 alpha;
-	float32* cos_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1];
-	float32* sin_alpha = new float32[m_pProjectionGeometry->getDetectorCount() + 1];
-	for (int i = 0; i < m_pProjectionGeometry->getDetectorCount() + 1; ++i) {
-		alpha = -atan((i - m_pProjectionGeometry->getDetectorCount()*0.5f) * DW / projgeom->getSourceDetectorDistance());
-		cos_alpha[i] = cos(alpha);
-		sin_alpha[i] = sin(alpha);
-	}
-
-		
-	// get values
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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
-	float32 sin_theta = sin(theta);
-	float32 cos_theta = cos(theta);
-
-	// [-45?,45?] and [135?,225?]
-	if (theta < PIdiv4) {
-
-		iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-		// POLICY: RAY PRIOR
-		if (!p.rayPrior(iRayIndex)) { delete[] cos_alpha; delete[] sin_alpha; return; }
-
-		float32 sin_theta_left, cos_theta_left;
-		float32 sin_theta_right, cos_theta_right;
-
-		// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
-		float32 t_l, t_r;
-		if (!switch_t) {
-			sin_theta_left = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1];
-			sin_theta_right = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector];
-
-			cos_theta_left = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1];
-			cos_theta_right = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector];
-
-			t_l = sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance();
-			t_r = sin_alpha[_iDetector] * projgeom->getOriginSourceDistance();
-
-		} else {
-			sin_theta_left = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector];
-			sin_theta_right = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1];
-
-			cos_theta_left = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector];
-			cos_theta_right = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1];
-
-			t_l = -sin_alpha[_iDetector] * projgeom->getOriginSourceDistance();
-			t_r = -sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance();	
-		}
-
-		float32 inv_cos_theta_left = 1.0f / cos_theta_left; 
-		float32 inv_cos_theta_right = 1.0f / cos_theta_right; 
-
-		float32 updateX_left = sin_theta_left * inv_cos_theta_left;
-		float32 updateX_right = sin_theta_right * inv_cos_theta_right;
-
-		// Precalculate kernel limits
-		float32 S_l = -0.5f * updateX_left;
-		if (S_l > 0) {S_l = -S_l;}
-		float32 T_l = -S_l;
-		float32 U_l = 1.0f + S_l;
-		float32 V_l = 1.0f - S_l;
-		float32 inv_4T_l = 0.25f / T_l;
-
-		float32 S_r = -0.5f * updateX_right;
-		if (S_r > 0) {S_r = -S_r;}
-		float32 T_r = -S_r;
-		float32 U_r = 1.0f + S_r;
-		float32 V_r = 1.0f - S_r;
-		float32 inv_4T_r = 0.25f / T_r;
-
-		// calculate strip extremes (volume coordinates)
-		float32 PL = (t_l - sin_theta_left * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_left;
-		float32 PR = (t_r - sin_theta_right * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta_right;
-		float32 PLimitL = PL + S_l * PH;
-		float32 PLimitR = PR - S_r * PH;
-		
-		// calculate strip extremes (pixel coordinates)
-		float32 XLimitL = (PLimitL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-		float32 XLimitR = (PLimitR - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-		float32 xL = (PL - m_pVolumeGeometry->getWindowMinX()) * inv_PW;
-		float32 xR = (PR - m_pVolumeGeometry->getWindowMinX()) * 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_left; 
-			XLimitR += updateX_right;
-			xL += updateX_left; 
-			xR += updateX_right; 
-
-			// 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_r)			res = 1.0f;
-				else if (x2R > U_r)		res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
-				else if (x2R >= T_r)	res = x2R;
-				else if (x2R > S_r)		res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
-				else					{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-						
-				// left
-				if (x2L <= S_l)			{}
-				else if (x2L < T_l)		res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
-				else if (x2L <= U_l)	res -= x2L;
-				else if (x2L < V_l)		res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
-				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);
-
-	// [45?,135?] and [225?,315?]
-	// horizontaly
-	} else {
-
-		iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-		// POLICY: RAY PRIOR
-		if (!p.rayPrior(iRayIndex)) { delete[] cos_alpha; delete[] sin_alpha; return; }
-
-		// get theta_l = alpha_left + theta and theta_r = alpha_right + theta
-		float32 sin_theta_left, cos_theta_left;
-		float32 sin_theta_right, cos_theta_right;
-		float32 t_l, t_r;
-		if (!switch_t) {
-			sin_theta_left = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector];
-			sin_theta_right = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1];
-
-			cos_theta_left = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector];
-			cos_theta_right = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1];
-
-			t_l = sin_alpha[_iDetector] * projgeom->getOriginSourceDistance();
-			t_r = sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance();
-
-		} else {
-			sin_theta_left = sin_theta * cos_alpha[_iDetector+1] + cos_theta * sin_alpha[_iDetector+1];
-			sin_theta_right = sin_theta * cos_alpha[_iDetector] + cos_theta * sin_alpha[_iDetector];
-
-			cos_theta_left = cos_theta * cos_alpha[_iDetector+1] - sin_theta * sin_alpha[_iDetector+1];
-			cos_theta_right = cos_theta * cos_alpha[_iDetector] - sin_theta * sin_alpha[_iDetector];
-
-			t_l = -sin_alpha[_iDetector+1] * projgeom->getOriginSourceDistance();
-			t_r = -sin_alpha[_iDetector] * projgeom->getOriginSourceDistance();	
-		}
-
-		float32 inv_sin_theta_left = 1.0f / sin_theta_left;
-		float32 inv_sin_theta_right = 1.0f / sin_theta_right;
-
-		float32 updateX_left = cos_theta_left * inv_sin_theta_left;
-		float32 updateX_right = cos_theta_right * inv_sin_theta_right;
-
-		// Precalculate kernel limits
-		float32 S_l = -0.5f * updateX_left;
-		if (S_l > 0) { S_l = -S_l; }
-		float32 T_l = -S_l;
-		float32 U_l = 1.0f + S_l;
-		float32 V_l = 1.0f - S_l;
-		float32 inv_4T_l = 0.25f / T_l;
-
-		float32 S_r = -0.5f * updateX_right;
-		if (S_r > 0) { S_r = -S_r; }
-		float32 T_r = -S_r;
-		float32 U_r = 1.0f + S_r;
-		float32 V_r = 1.0f - S_r;
-		float32 inv_4T_r = 0.25f / T_r;
-
-		// calculate strip extremes (volume coordinates)
-		float32 PL = (t_l - cos_theta_left * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_left;
-		float32 PR = (t_r - cos_theta_right * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta_right;
-		float32 PLimitL = PL - S_l * PW;
-		float32 PLimitR = PR + S_r * PW;
-		
-		// calculate strip extremes (pixel coordinates)
-		float32 XLimitL = (m_pVolumeGeometry->getWindowMaxY() - PLimitL) * inv_PH;
-		float32 XLimitR = (m_pVolumeGeometry->getWindowMaxY() - PLimitR) * inv_PH;
-		float32 xL = (m_pVolumeGeometry->getWindowMaxY() - PL) * inv_PH;
-		float32 xR = (m_pVolumeGeometry->getWindowMaxY() - PR) * 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_left; 
-			XLimitR += updateX_right;
-			xL += updateX_left; 
-			xR += updateX_right; 
-
-			// for each affected row
-			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_r)			res = 1.0f;
-				else if (x2R > U_r)		res = x2R - (x2R-U_r)*(x2R-U_r)*inv_4T_r;
-				else if (x2R >= T_r)	res = x2R;
-				else if (x2R > S_r)		res = (x2R-S_r)*(x2R-S_r) * inv_4T_r;
-				else					{ x2L -= 1.0f; x2R -= 1.0f;	continue; }
-						
-				// left
-				if (x2L <= S_l)			{}
-				else if (x2L < T_l)		res -= (x2L-S_l)*(x2L-S_l) * inv_4T_l;
-				else if (x2L <= U_l)	res -= x2L;
-				else if (x2L < V_l)		res -= x2L - (x2L-U_l)*(x2L-U_l)*inv_4T_l;
-				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 theta switch
-
-	delete[] cos_alpha;
-	delete[] sin_alpha;
-}
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.h b/include/astra/ParallelBeamLineKernelProjector2D.h
index 9e2cc5a..64b5ef3 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.h
+++ b/include/astra/ParallelBeamLineKernelProjector2D.h
@@ -173,6 +173,13 @@ public:
 	 */
 	virtual std::string getType();
 
+
+protected:
+	/** Internal policy-based projection of a range of angles and range.
+ 	 * (_i*From is inclusive, _i*To exclusive) */
+	template <typename Policy>
+	void projectBlock_internal(int _iProjFrom, int _iProjTo,
+	                           int _iDetFrom, int _iDetTo, Policy& _policy);
 };
 
 inline std::string CParallelBeamLineKernelProjector2D::getType() 
diff --git a/include/astra/ParallelBeamLineKernelProjector2D.inl b/include/astra/ParallelBeamLineKernelProjector2D.inl
index 08bbe5f..149ac46 100644
--- a/include/astra/ParallelBeamLineKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLineKernelProjector2D.inl
@@ -26,12 +26,32 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 $Id$
 */
 
+template <typename Policy>
+void CParallelBeamLineKernelProjector2D::project(Policy& p)
+{
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      _iDetector, _iDetector + 1, p);
+}
 
 
 //----------------------------------------------------------------------------------------
-// PROJECT ALL 
+// PROJECT BLOCK
 template <typename Policy>
-void CParallelBeamLineKernelProjector2D::project(Policy& p)
+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;
@@ -40,7 +60,7 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p)
 	bool switch_t;
 
 	// loop angles
-	for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
+	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 
 		// get theta
 		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
@@ -72,7 +92,7 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p)
 		T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
 
 		// loop detectors
-		for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
+		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
 			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
@@ -267,465 +287,3 @@ void CParallelBeamLineKernelProjector2D::project(Policy& p)
 
 }
 
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
-template <typename Policy>
-void CParallelBeamLineKernelProjector2D::projectSingleProjection(int _iProjection, 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, iDetector, x1;
-	bool switch_t;
-
-	// get theta
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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_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);
-
-	// loop detectors
-	for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		
-		iRayIndex = _iProjection * 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) {
-		
-			// calculate x for row 0
-			P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-			x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
-
-			// get coords
-			int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-			float nextx2 = x - nextx1;
-
-			// for each row
-			for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-
-				x1 = nextx1;
-				x2 = nextx2;
-
-				nextx2 += updatePerRow;
-				while (nextx2 >= 1.0f) {
-					nextx2 -= 1.0f;
-					nextx1++;
-				}
-				while (nextx2 < 0.0f) {
-					nextx2 += 1.0f;
-					nextx1--;
-				}
-
-				if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
-
-				// 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);
-						}
-					}
-				}
-
-				// 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 {
-					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);
-						}
-					}
-				}
-			}
-		}
-
-		// 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->getWindowMaxY() - P) * inv_pixelLengthY;
-
-			// get coords
-			int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-			float nextx2 = x - nextx1;
-
-			// for each col
-			for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-			
-				x1 = nextx1;
-				x2 = nextx2;
-
-				nextx2 += updatePerCol;
-				while (nextx2 >= 1.0f) {
-					nextx2 -= 1.0f;
-					nextx1++;
-				}
-				while (nextx2 < 0.0f) {
-					nextx2 += 1.0f;
-					nextx1--;
-				}
-
-				if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
-
-				// 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);
-						}
-					}
-				}
-
-				// 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);
-						}
-					}					
-				}
-
-				// down
-				else {
-					I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-					if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-						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->getGridColCount()) {
-						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);
-						}
-					}
-				}
-			}
-		} // end loop col
-
-		// POLICY: RAY POSTERIOR
-		p.rayPosterior(iRayIndex);
-
-	} // end loop detector
-}
-
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
-template <typename Policy>
-void CParallelBeamLineKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
-{
-	int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-	// POLICY: RAY PRIOR
-	if (!p.rayPrior(iRayIndex)) return;
-
-	// variables
-	float32 t, I, P, x, x2;
-	int iVolumeIndex, row, col, x1;
-
-	// get theta
-	float32 theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	bool switch_t = false;
-	if (theta >= 7*PIdiv4) theta -= 2*PI;
-	if (theta >= 3*PIdiv4) {
-		theta -= PI;
-		switch_t = true;
-	}
-
-	// get t
-	t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector);
-	if (switch_t) t = -t;
-
-	// vertically
-	if (theta <= PIdiv4) {
-	
-		float32 sin_theta = sin(theta);
-		float32 inv_cos_theta = 1.0f / cos(theta); 
-
-		// precalculate kernel limits
-		float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-		float32 updatePerRow = sin_theta * inv_cos_theta;
-		float32 inv_pixelLengthX = 1.0f / m_pVolumeGeometry->getPixelLengthX();
-		float32 S = 0.5f - 0.5f * ((updatePerRow < 0) ? -updatePerRow : updatePerRow);
-
-		// calculate x for row 0
-		P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-		x = (P - m_pVolumeGeometry->getWindowMinX()) * inv_pixelLengthX;
-
-		// get coords
-		int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-		float nextx2 = x - nextx1;
-
-		// for each row
-		for (row = 0; row < m_pVolumeGeometry->getGridRowCount(); ++row) {
-
-			x1 = nextx1;
-			x2 = nextx2;
-
-			nextx2 += updatePerRow;
-			while (nextx2 >= 1.0f) {
-				nextx2 -= 1.0f;
-				nextx1++;
-			}
-			while (nextx2 < 0.0f) {
-				nextx2 += 1.0f;
-				nextx1--;
-			}
-
-			if (x1 < -1 || x1 > m_pVolumeGeometry->getGridColCount()) continue;
-
-			// 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);
-					}
-				}
-			}
-
-			// 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 {
-				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);
-					}
-				}
-			}
-		}
-	}
-
-	// horizontally
-	else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
-
-		float32 cos_theta = cos(theta);
-		float32 inv_sin_theta = 1.0f / sin(theta); 
-
-		// precalculate kernel limits
-		float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
-		float32 updatePerCol = cos_theta * inv_sin_theta;
-		float32 inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
-		float32 T = 0.5f - 0.5f * ((updatePerCol < 0) ? -updatePerCol : updatePerCol);
-
-		// calculate point P
-		P = (t - cos_theta * m_pVolumeGeometry->pixelColToCenterX(0)) * inv_sin_theta;
-		x = (m_pVolumeGeometry->getWindowMaxY() - P) * inv_pixelLengthY;
-
-		// get coords
-		int nextx1 = int((x > 0.0f) ? x : x-1.0f);
-		float nextx2 = x - nextx1;
-
-		// for each col
-		for (col = 0; col < m_pVolumeGeometry->getGridColCount(); ++col) {
-
-			x1 = nextx1;
-			x2 = nextx2;
-
-			nextx2 += updatePerCol;
-			while (nextx2 >= 1.0f) {
-				nextx2 -= 1.0f;
-				nextx1++;
-			}
-			while (nextx2 < 0.0f) {
-				nextx2 += 1.0f;
-				nextx1--;
-			}
-
-			if (x1 < -1 || x1 > m_pVolumeGeometry->getGridRowCount()) continue;
-
-			// 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);
-					}
-				}
-			}
-
-			// 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);
-					}
-				}					
-			}
-
-			// down
-			else {
-				I = (1.5f - T - x2) / (1.0f - 2.0f*T) * lengthPerCol;
-
-				if (x1 >= 0 && x1 < m_pVolumeGeometry->getGridColCount()) {
-					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->getGridColCount()) {
-					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);
-					}
-				}
-			}
-		}
-	} // end loop col
-
-	// POLICY: RAY POSTERIOR
-	p.rayPosterior(iRayIndex);
-
-}
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.h b/include/astra/ParallelBeamLinearKernelProjector2D.h
index ac5899e..0a011f1 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.h
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.h
@@ -177,6 +177,15 @@ public:
 	 */
 	virtual std::string getType();
 
+
+protected:
+	/** Internal policy-based projection of a range of angles and range.
+ 	 * (_i*From is inclusive, _i*To exclusive) */
+	template <typename Policy>
+	void projectBlock_internal(int _iProjFrom, int _iProjTo,
+	                           int _iDetFrom, int _iDetTo, Policy& _policy);
+
+
 };
 
 //----------------------------------------------------------------------------------------
diff --git a/include/astra/ParallelBeamLinearKernelProjector2D.inl b/include/astra/ParallelBeamLinearKernelProjector2D.inl
index 66b1cb2..d545f30 100644
--- a/include/astra/ParallelBeamLinearKernelProjector2D.inl
+++ b/include/astra/ParallelBeamLinearKernelProjector2D.inl
@@ -27,11 +27,31 @@ $Id$
 */
 
 
+template <typename Policy>
+void CParallelBeamLinearKernelProjector2D::project(Policy& p)
+{
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamLinearKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      _iDetector, _iDetector + 1, p);
+}
 
 //----------------------------------------------------------------------------------------
-// PROJECT ALL 
+// PROJECT BLOCK
 template <typename Policy>
-void CParallelBeamLinearKernelProjector2D::project(Policy& p)
+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;
@@ -43,7 +63,7 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p)
 	float32 P,x,x2;
 
 	// loop angles
-	for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
+	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 
 		// get theta
 		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
@@ -71,7 +91,7 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p)
 		inv_pixelLengthY = 1.0f / m_pVolumeGeometry->getPixelLengthY();
 
 		// loop detectors
-		for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
+		for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
 			iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
@@ -162,255 +182,3 @@ void CParallelBeamLinearKernelProjector2D::project(Policy& p)
 
 }
 
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
-template <typename Policy>
-void CParallelBeamLinearKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
-{
-	// variables
-	float32 theta, sin_theta, cos_theta, inv_sin_theta, inv_cos_theta, t;
-	float32 lengthPerRow, updatePerRow, inv_pixelLengthX;
-	float32 lengthPerCol, updatePerCol, inv_pixelLengthY;
-	bool switch_t;
-	int iDetector, iVolumeIndex, iRayIndex;
-	int row, col, x1;
-	float32 P,x,x2;
-
-	// get theta
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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; 
-
-	// 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();
-
-	// loop detectors
-	for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		
-		iRayIndex = _iProjection * 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) {
-		
-			// calculate x for row 0
-			P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-			x = m_pVolumeGeometry->coordXToColF(P) - 0.5f;
-
-			// 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;
-
-				// 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);
-					}
-				}
-			}
-		}
-
-		// 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);
-					}
-				}
-			}
-		}
-
-		// POLICY: RAY POSTERIOR
-		p.rayPosterior(iRayIndex);
-
-	} // end loop detector
-}
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
-template <typename Policy>
-void CParallelBeamLinearKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
-{
-	int iVolumeIndex, iRayIndex;
-
-	iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-	// POLICY: RAY PRIOR
-	if (!p.rayPrior(iRayIndex)) return;
-
-	// get theta
-	float32 theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	bool switch_t = false;
-	if (theta >= 7*PIdiv4) theta -= 2*PI;
-	if (theta >= 3*PIdiv4) {
-		theta -= PI;
-		switch_t = true;
-	}
-
-	// get t
-	float32 t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector);
-	if (switch_t) {
-		t = -t;
-	}
-
-	// vertically
-	if (theta <= PIdiv4) {
-	
-		// precalculate sin, 1/cos
-		float32 sin_theta = sin(theta);
-		float32 inv_cos_theta = 1.0f / cos(theta); 
-
-		// precalculate kernel limits
-		float32 lengthPerRow = m_pVolumeGeometry->getPixelLengthY() * inv_cos_theta;
-		float32 updatePerRow = sin_theta * inv_cos_theta;
-
-		int row, x1;
-		float32 P,x,x2;
-
-		// calculate x for row 0
-		P = (t - sin_theta * m_pVolumeGeometry->pixelRowToCenterY(0)) * inv_cos_theta;
-		x = m_pVolumeGeometry->coordXToColF(P) - 0.5f;
-
-		// 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;
-
-			// 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);
-				}
-			}
-		}
-	}
-
-	// horizontally
-	else if (PIdiv4 <= theta && theta <= 3*PIdiv4) {
-
-		// precalculate cos 1/sin
-		float32 cos_theta = cos(theta);
-		float32 inv_sin_theta = 1.0f / sin(theta); 
-
-		// precalculate kernel limits
-		float32 lengthPerCol = m_pVolumeGeometry->getPixelLengthX() * inv_sin_theta;
-		float32 updatePerCol = cos_theta * inv_sin_theta;
-
-		int col, x1;
-		float32 P,x,x2;
-
-		// 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);
-				}
-			}
-		}
-	}
-
-	// POLICY: RAY POSTERIOR
-	p.rayPosterior(iRayIndex);
-}
diff --git a/include/astra/ParallelBeamStripKernelProjector2D.h b/include/astra/ParallelBeamStripKernelProjector2D.h
index de056de..e3808ca 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.h
+++ b/include/astra/ParallelBeamStripKernelProjector2D.h
@@ -176,6 +176,13 @@ protected:
 	 */
 	virtual std::string getType();
 
+
+	/** 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/ParallelBeamStripKernelProjector2D.inl b/include/astra/ParallelBeamStripKernelProjector2D.inl
index c55fa8e..a9168ee 100644
--- a/include/astra/ParallelBeamStripKernelProjector2D.inl
+++ b/include/astra/ParallelBeamStripKernelProjector2D.inl
@@ -26,11 +26,31 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 $Id$
 */
 
+template <typename Policy>
+void CParallelBeamStripKernelProjector2D::project(Policy& p)
+{
+	projectBlock_internal(0, m_pProjectionGeometry->getProjectionAngleCount(),
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      0, m_pProjectionGeometry->getDetectorCount(), p);
+}
+
+template <typename Policy>
+void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
+{
+	projectBlock_internal(_iProjection, _iProjection + 1,
+	                      _iDetector, _iDetector + 1, p);
+}
 
 //----------------------------------------------------------------------------------------
-// PROJECT ALL 
+// PROJECT BLOCK
 template <typename Policy>
-void CParallelBeamStripKernelProjector2D::project(Policy& p)
+void CParallelBeamStripKernelProjector2D::projectBlock_internal(int _iProjFrom, int _iProjTo, int _iDetFrom, int _iDetTo, Policy& p)
 {
 	ASTRA_ASSERT(m_bIsInitialized);
 
@@ -52,7 +72,7 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p)
 	float32 S, T, U, V, inv_4T;
 
 	// loop angles
-	for (iAngle = 0; iAngle < m_pProjectionGeometry->getProjectionAngleCount(); ++iAngle) {
+	for (iAngle = _iProjFrom; iAngle < _iProjTo; ++iAngle) {
 		
 		// get values
 		theta = m_pProjectionGeometry->getProjectionAngle(iAngle);
@@ -94,7 +114,7 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p)
 			updateX = sin_theta * inv_cos_theta;
 
 			// loop detectors
-			for (iDetector = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
+			for (iDetector = _iDetFrom; iDetector < _iDetTo; ++iDetector) {
 			
 				iRayIndex = iAngle * m_pProjectionGeometry->getDetectorCount() + iDetector;
 
@@ -276,464 +296,3 @@ void CParallelBeamStripKernelProjector2D::project(Policy& p)
 }
 
 
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE PROJECTION
-template <typename Policy>
-void CParallelBeamStripKernelProjector2D::projectSingleProjection(int _iProjection, Policy& p)
-{
-	// Some variables
-	float32 theta, t;
-	int row, col, 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 values
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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 = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		
-			iRayIndex = _iProjection * 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
-	} else {
-
-		// 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 = 0; iDetector < m_pProjectionGeometry->getDetectorCount(); ++iDetector) {
-		
-			iRayIndex = _iProjection * 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 - 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; 
-
-				// 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
-}
-
-
-//----------------------------------------------------------------------------------------
-// PROJECT SINGLE RAY
-template <typename Policy>
-void CParallelBeamStripKernelProjector2D::projectSingleRay(int _iProjection, int _iDetector, Policy& p)
-{
-	int iRayIndex = _iProjection * m_pProjectionGeometry->getDetectorCount() + _iDetector;
-
-	// POLICY: RAY PRIOR
-	if (!p.rayPrior(iRayIndex)) return;
-
-	// Some variables
-	float32 theta, t;
-	int row, col;
-	float32 res;
-	float32 PL, PLimitL, PLimitR;
-	float32 xL, xR, XLimitL, XLimitR;
-	int x1L,x1R;
-	float32 x2L, x2R, updateX;
-	int iVolumeIndex;
-	
-	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 values
-	theta = m_pProjectionGeometry->getProjectionAngle(_iProjection);
-	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;
-		
-		// 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
-
-	// [45?,135?] and [225?,315?]
-	// horizontaly
-	} else {
-
-		// 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;
-
-		// get t
-		t = m_pProjectionGeometry->indexToDetectorOffset(_iDetector);
-		if (switch_t) t = -t;
-	
-		// 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; 
-
-			// 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
-
-	} // end theta switch
-
-	// POLICY: RAY POSTERIOR
-	p.rayPosterior(iRayIndex);
-
-}