Merge branch 'master' into aniso

15 files changed, 368 insertions, 451 deletions

diff --git a/src/ArtAlgorithm.cpp b/src/ArtAlgorithm.cpp
index b59bd93..526c263 100644
--- a/src/ArtAlgorithm.cpp
+++ b/src/ArtAlgorithm.cpp
@@ -156,8 +156,12 @@ bool CArtAlgorithm::initialize(const Config& _cfg)
 		return false;
 	}
 
+	// "Lambda" is replaced by the more descriptive "Relaxation"
 	m_fLambda = _cfg.self.getOptionNumerical("Lambda", 1.0f);
-	CC.markOptionParsed("Lambda");
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation", m_fLambda);
+	if (!_cfg.self.hasOption("Relaxation"))
+		CC.markOptionParsed("Lambda");
+	CC.markOptionParsed("Relaxation");
 
 	// success
 	m_bIsInitialized = _check();
@@ -232,7 +236,7 @@ map<string,boost::any> CArtAlgorithm::getInformation()
 {
 	map<string, boost::any> res;
 	res["RayOrder"] = getInformation("RayOrder");
-	res["Lambda"] = getInformation("Lambda");
+	res["Relaxation"] = getInformation("Relaxation");
 	return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res);
 };
 
@@ -240,7 +244,7 @@ map<string,boost::any> CArtAlgorithm::getInformation()
 // Information - Specific
 boost::any CArtAlgorithm::getInformation(std::string _sIdentifier) 
 {
-	if (_sIdentifier == "Lambda")	{ return m_fLambda; }
+	if (_sIdentifier == "Relaxation")	{ return m_fLambda; }
 	if (_sIdentifier == "RayOrder") {
 		vector<float32> res;
 		for (int i = 0; i < m_iRayCount; i++) {
diff --git a/src/AstraObjectManager.cpp b/src/AstraObjectManager.cpp
index c49f273..46eae4b 100644
--- a/src/AstraObjectManager.cpp
+++ b/src/AstraObjectManager.cpp
@@ -31,13 +31,13 @@ $Id$
 
 namespace astra {
 
-int CAstraIndexManager::m_iPreviousIndex = 0;
-
-DEFINE_SINGLETON(CAstraObjectManager<CProjector2D>);
-DEFINE_SINGLETON(CAstraObjectManager<CProjector3D>);
-DEFINE_SINGLETON(CAstraObjectManager<CFloat32Data2D>);
-DEFINE_SINGLETON(CAstraObjectManager<CFloat32Data3D>);
-DEFINE_SINGLETON(CAstraObjectManager<CAlgorithm>);
-DEFINE_SINGLETON(CAstraObjectManager<CSparseMatrix>);
+DEFINE_SINGLETON(CProjector2DManager);
+DEFINE_SINGLETON(CProjector3DManager);
+DEFINE_SINGLETON(CData2DManager);
+DEFINE_SINGLETON(CData3DManager);
+DEFINE_SINGLETON(CAlgorithmManager);
+DEFINE_SINGLETON(CMatrixManager);
+
+DEFINE_SINGLETON(CAstraIndexManager);
 
 } // end namespace
diff --git a/src/CompositeGeometryManager.cpp b/src/CompositeGeometryManager.cpp
index 96b28e9..084ba8c 100644
--- a/src/CompositeGeometryManager.cpp
+++ b/src/CompositeGeometryManager.cpp
@@ -45,14 +45,15 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 
 #include <cstring>
 #include <sstream>
+#include <climits>
 
 #ifndef USE_PTHREADS
 #include <boost/thread/mutex.hpp>
 #include <boost/thread.hpp>
 #endif
 
-namespace astra {
 
+namespace astra {
 
 SGPUParams* CCompositeGeometryManager::s_params = 0;
 
@@ -98,6 +99,9 @@ CCompositeGeometryManager::CCompositeGeometryManager()
 
 bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div, TJobSet &split)
 {
+	int maxBlockDim = astraCUDA3d::maxBlockDimension();
+	ASTRA_DEBUG("Found max block dim %d", maxBlockDim);
+
 	split.clear();
 
 	for (TJobSet::const_iterator i = jobs.begin(); i != jobs.end(); ++i)
@@ -111,7 +115,22 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
 		//    b. split input part
 		//    c. create jobs for new (input,output) subparts
 
-		TPartList splitOutput = pOutput->split(maxSize/3, div);
+		TPartList splitOutput;
+		pOutput->splitZ(splitOutput, maxSize/3, UINT_MAX, div);
+#if 0
+		TPartList splitOutput2;
+		for (TPartList::iterator i_out = splitOutput.begin(); i_out != splitOutput.end(); ++i_out) {
+			boost::shared_ptr<CPart> outputPart = *i_out;
+			outputPart.get()->splitX(splitOutput2, UINT_MAX, UINT_MAX, 1);
+		}
+		splitOutput.clear();
+		for (TPartList::iterator i_out = splitOutput2.begin(); i_out != splitOutput2.end(); ++i_out) {
+			boost::shared_ptr<CPart> outputPart = *i_out;
+					outputPart.get()->splitY(splitOutput, UINT_MAX, UINT_MAX, 1);
+		}
+		splitOutput2.clear();
+#endif
+
 
 		for (TJobList::const_iterator j = L.begin(); j != L.end(); ++j)
 		{
@@ -139,8 +158,21 @@ bool CCompositeGeometryManager::splitJobs(TJobSet &jobs, size_t maxSize, int div
 
 				size_t remainingSize = ( maxSize - outputPart->getSize() ) / 2;
 
-				TPartList splitInput = input->split(remainingSize, 1);
+				TPartList splitInput;
+				input->splitZ(splitInput, remainingSize, maxBlockDim, 1);
 				delete input;
+				TPartList splitInput2;
+				for (TPartList::iterator i_in = splitInput.begin(); i_in != splitInput.end(); ++i_in) {
+					boost::shared_ptr<CPart> inputPart = *i_in;
+					inputPart.get()->splitX(splitInput2, UINT_MAX, maxBlockDim, 1);
+				}
+				splitInput.clear();
+				for (TPartList::iterator i_in = splitInput2.begin(); i_in != splitInput2.end(); ++i_in) {
+					boost::shared_ptr<CPart> inputPart = *i_in;
+					inputPart.get()->splitY(splitInput, UINT_MAX, maxBlockDim, 1);
+				}
+				splitInput2.clear();
+
 				ASTRA_DEBUG("Input split into %d parts", splitInput.size());
 
 				for (TPartList::iterator i_in = splitInput.begin();
@@ -327,10 +359,14 @@ static size_t ceildiv(size_t a, size_t b) {
 	return (a + b - 1) / b;
 }
 
-static size_t computeVerticalSplit(size_t maxBlock, int div, size_t sliceCount)
+static size_t computeLinearSplit(size_t maxBlock, int div, size_t sliceCount)
 {
 	size_t blockSize = maxBlock;
-	size_t blockCount = ceildiv(sliceCount, blockSize);
+	size_t blockCount;
+	if (sliceCount <= blockSize)
+		blockCount = 1;
+	else
+		blockCount = ceildiv(sliceCount, blockSize);
 
 	// Increase number of blocks to be divisible by div
 	size_t divCount = div * ceildiv(blockCount, div);
@@ -410,7 +446,17 @@ SPar3DProjection* getProjectionVectors(const CParallelVecProjectionGeometry3D* p
 
 
 template<class V>
-static void translateProjectionVectors(V* pProjs, int count, double dv)
+static void translateProjectionVectorsU(V* pProjs, int count, double du)
+{
+	for (int i = 0; i < count; ++i) {
+		pProjs[i].fDetSX += du * pProjs[i].fDetUX;
+		pProjs[i].fDetSY += du * pProjs[i].fDetUY;
+		pProjs[i].fDetSZ += du * pProjs[i].fDetUZ;
+	}
+}
+
+template<class V>
+static void translateProjectionVectorsV(V* pProjs, int count, double dv)
 {
 	for (int i = 0; i < count; ++i) {
 		pProjs[i].fDetSX += dv * pProjs[i].fDetVX;
@@ -420,8 +466,58 @@ static void translateProjectionVectors(V* pProjs, int count, double dv)
 }
 
 
+static CProjectionGeometry3D* getSubProjectionGeometryU(const CProjectionGeometry3D* pProjGeom, int u, int size)
+{
+	// First convert to vectors, then translate, then convert into new object
 
-static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry3D* pProjGeom, int v, int size)
+	const CConeProjectionGeometry3D* conegeom = dynamic_cast<const CConeProjectionGeometry3D*>(pProjGeom);
+	const CParallelProjectionGeometry3D* par3dgeom = dynamic_cast<const CParallelProjectionGeometry3D*>(pProjGeom);
+	const CParallelVecProjectionGeometry3D* parvec3dgeom = dynamic_cast<const CParallelVecProjectionGeometry3D*>(pProjGeom);
+	const CConeVecProjectionGeometry3D* conevec3dgeom = dynamic_cast<const CConeVecProjectionGeometry3D*>(pProjGeom);
+
+	if (conegeom || conevec3dgeom) {
+		SConeProjection* pConeProjs;
+		if (conegeom) {
+			pConeProjs = getProjectionVectors<SConeProjection>(conegeom);
+		} else {
+			pConeProjs = getProjectionVectors<SConeProjection>(conevec3dgeom);
+		}
+
+		translateProjectionVectorsU(pConeProjs, pProjGeom->getProjectionCount(), u);
+
+		CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
+		                                                              pProjGeom->getDetectorRowCount(),
+		                                                              size,
+		                                                              pConeProjs);
+
+
+		delete[] pConeProjs;
+		return ret;
+	} else {
+		assert(par3dgeom || parvec3dgeom);
+		SPar3DProjection* pParProjs;
+		if (par3dgeom) {
+			pParProjs = getProjectionVectors<SPar3DProjection>(par3dgeom);
+		} else {
+			pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom);
+		}
+
+		translateProjectionVectorsU(pParProjs, pProjGeom->getProjectionCount(), u);
+
+		CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
+		                                                                  pProjGeom->getDetectorRowCount(),
+		                                                                  size,
+		                                                                  pParProjs);
+
+		delete[] pParProjs;
+		return ret;
+	}
+
+}
+
+
+
+static CProjectionGeometry3D* getSubProjectionGeometryV(const CProjectionGeometry3D* pProjGeom, int v, int size)
 {
 	// First convert to vectors, then translate, then convert into new object
 
@@ -438,7 +534,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry
 			pConeProjs = getProjectionVectors<SConeProjection>(conevec3dgeom);
 		}
 
-		translateProjectionVectors(pConeProjs, pProjGeom->getProjectionCount(), v);
+		translateProjectionVectorsV(pConeProjs, pProjGeom->getProjectionCount(), v);
 
 		CProjectionGeometry3D* ret = new CConeVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
 		                                                              size,
@@ -457,7 +553,7 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry
 			pParProjs = getProjectionVectors<SPar3DProjection>(parvec3dgeom);
 		}
 
-		translateProjectionVectors(pParProjs, pProjGeom->getProjectionCount(), v);
+		translateProjectionVectorsV(pParProjs, pProjGeom->getProjectionCount(), v);
 
 		CProjectionGeometry3D* ret = new CParallelVecProjectionGeometry3D(pProjGeom->getProjectionCount(),
 		                                                                  size,
@@ -476,17 +572,110 @@ static CProjectionGeometry3D* getSubProjectionGeometry(const CProjectionGeometry
 // - each no bigger than maxSize
 // - number of sub-parts is divisible by div
 // - maybe all approximately the same size?
-CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::split(size_t maxSize, int div)
+void CCompositeGeometryManager::CVolumePart::splitX(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
+{
+	if (true) {
+		// Split in vertical direction only at first, until we figure out
+		// a model for splitting in other directions
+
+		size_t sliceSize = ((size_t) pGeom->getGridSliceCount()) * pGeom->getGridRowCount();
+		int sliceCount = pGeom->getGridColCount();
+		size_t m = std::min(maxSize / sliceSize, maxDim);
+		size_t blockSize = computeLinearSplit(m, div, sliceCount);
+
+		int rem = sliceCount % blockSize;
+
+		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
+
+		for (int x = -(rem / 2); x < sliceCount; x += blockSize) {
+			int newsubX = x;
+			if (newsubX < 0) newsubX = 0;
+			int endX = x + blockSize;
+			if (endX > sliceCount) endX = sliceCount;
+			int size = endX - newsubX;
+
+			CVolumePart *sub = new CVolumePart();
+			sub->subX = this->subX + newsubX;
+			sub->subY = this->subY;
+			sub->subZ = this->subZ;
+
+			ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub);
+
+			double shift = pGeom->getPixelLengthX() * newsubX;
+
+			sub->pData = pData;
+			sub->pGeom = new CVolumeGeometry3D(size,
+			                                   pGeom->getGridRowCount(),
+			                                   pGeom->getGridSliceCount(),
+			                                   pGeom->getWindowMinX() + shift,
+			                                   pGeom->getWindowMinY(),
+			                                   pGeom->getWindowMinZ(),
+			                                   pGeom->getWindowMinX() + shift + size * pGeom->getPixelLengthX(),
+			                                   pGeom->getWindowMaxY(),
+			                                   pGeom->getWindowMaxZ());
+
+			out.push_back(boost::shared_ptr<CPart>(sub));
+		}
+	}
+}
+
+void CCompositeGeometryManager::CVolumePart::splitY(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
 {
-	TPartList ret;
+	if (true) {
+		// Split in vertical direction only at first, until we figure out
+		// a model for splitting in other directions
+
+		size_t sliceSize = ((size_t) pGeom->getGridColCount()) * pGeom->getGridSliceCount();
+		int sliceCount = pGeom->getGridRowCount();
+		size_t m = std::min(maxSize / sliceSize, maxDim);
+		size_t blockSize = computeLinearSplit(m, div, sliceCount);
+
+		int rem = sliceCount % blockSize;
+
+		ASTRA_DEBUG("From %d to %d step %d", -(rem / 2), sliceCount, blockSize);
+
+		for (int y = -(rem / 2); y < sliceCount; y += blockSize) {
+			int newsubY = y;
+			if (newsubY < 0) newsubY = 0;
+			int endY = y + blockSize;
+			if (endY > sliceCount) endY = sliceCount;
+			int size = endY - newsubY;
 
+			CVolumePart *sub = new CVolumePart();
+			sub->subX = this->subX;
+			sub->subY = this->subY + newsubY;
+			sub->subZ = this->subZ;
+
+			ASTRA_DEBUG("VolumePart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub);
+
+			double shift = pGeom->getPixelLengthY() * newsubY;
+
+			sub->pData = pData;
+			sub->pGeom = new CVolumeGeometry3D(pGeom->getGridColCount(),
+			                                   size,
+			                                   pGeom->getGridSliceCount(),
+			                                   pGeom->getWindowMinX(),
+			                                   pGeom->getWindowMinY() + shift,
+			                                   pGeom->getWindowMinZ(),
+			                                   pGeom->getWindowMaxX(),
+			                                   pGeom->getWindowMinY() + shift + size * pGeom->getPixelLengthY(),
+			                                   pGeom->getWindowMaxZ());
+
+			out.push_back(boost::shared_ptr<CPart>(sub));
+		}
+	}
+}
+
+void CCompositeGeometryManager::CVolumePart::splitZ(CCompositeGeometryManager::TPartList& out, size_t maxSize, size_t maxDim, int div)
+{
 	if (true) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
 		size_t sliceSize = ((size_t) pGeom->getGridColCount()) * pGeom->getGridRowCount();
 		int sliceCount = pGeom->getGridSliceCount();
-		size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount);
+		size_t m = std::min(maxSize / sliceSize, maxDim);
+		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = sliceCount % blockSize;
 
@@ -519,11 +708,9 @@ CCompositeGeometryManager::TPartList CCompositeGeometryManager::CVolumePart::spl
 			                                   pGeom->getWindowMaxY(),
 			                                   pGeom->getWindowMinZ() + shift + size * pGeom->getPixelLengthZ());
 
-			ret.push_back(boost::shared_ptr<CPart>(sub));
+			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
 	}
-
-	return ret;
 }
 
 CCompositeGeometryManager::CVolumePart* CCompositeGeometryManager::CVolumePart::clone() const
@@ -630,7 +817,7 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re
 	if (_vmin == _vmax) {
 		sub->pGeom = 0;
 	} else {
-		sub->pGeom = getSubProjectionGeometry(pGeom, _vmin, _vmax - _vmin);
+		sub->pGeom = getSubProjectionGeometryV(pGeom, _vmin, _vmax - _vmin);
 	}
 
 	ASTRA_DEBUG("Reduce projection from %d - %d to %d - %d", this->subZ, this->subZ + pGeom->getDetectorRowCount(), this->subZ + _vmin, this->subZ + _vmax);
@@ -639,17 +826,58 @@ CCompositeGeometryManager::CPart* CCompositeGeometryManager::CProjectionPart::re
 }
 
 
-CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart::split(size_t maxSize, int div)
+void CCompositeGeometryManager::CProjectionPart::splitX(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div)
+{
+	if (true) {
+		// Split in vertical direction only at first, until we figure out
+		// a model for splitting in other directions
+
+		size_t sliceSize = ((size_t) pGeom->getDetectorRowCount()) * pGeom->getProjectionCount();
+		int sliceCount = pGeom->getDetectorColCount();
+		size_t m = std::min(maxSize / sliceSize, maxDim);
+		size_t blockSize = computeLinearSplit(m, div, sliceCount);
+
+		int rem = sliceCount % blockSize;
+
+		for (int x = -(rem / 2); x < sliceCount; x += blockSize) {
+			int newsubX = x;
+			if (newsubX < 0) newsubX = 0;
+			int endX = x + blockSize;
+			if (endX > sliceCount) endX = sliceCount;
+			int size = endX - newsubX;
+
+			CProjectionPart *sub = new CProjectionPart();
+			sub->subX = this->subX + newsubX;
+			sub->subY = this->subY;
+			sub->subZ = this->subZ;
+
+			ASTRA_DEBUG("ProjectionPart split %d %d %d -> %p", sub->subX, sub->subY, sub->subZ, (void*)sub);
+
+			sub->pData = pData;
+
+			sub->pGeom = getSubProjectionGeometryU(pGeom, newsubX, size);
+
+			out.push_back(boost::shared_ptr<CPart>(sub));
+		}
+	}
+}
+
+void CCompositeGeometryManager::CProjectionPart::splitY(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div)
 {
-	TPartList ret;
+	// TODO
+	out.push_back(boost::shared_ptr<CPart>(clone()));
+}
 
+void CCompositeGeometryManager::CProjectionPart::splitZ(CCompositeGeometryManager::TPartList &out, size_t maxSize, size_t maxDim, int div)
+{
 	if (true) {
 		// Split in vertical direction only at first, until we figure out
 		// a model for splitting in other directions
 
 		size_t sliceSize = ((size_t) pGeom->getDetectorColCount()) * pGeom->getProjectionCount();
 		int sliceCount = pGeom->getDetectorRowCount();
-		size_t blockSize = computeVerticalSplit(maxSize / sliceSize, div, sliceCount);
+		size_t m = std::min(maxSize / sliceSize, maxDim);
+		size_t blockSize = computeLinearSplit(m, div, sliceCount);
 
 		int rem = sliceCount % blockSize;
 
@@ -669,14 +897,12 @@ CCompositeGeometryManager::TPartList CCompositeGeometryManager::CProjectionPart:
 
 			sub->pData = pData;
 
-			sub->pGeom = getSubProjectionGeometry(pGeom, newsubZ, size);
+			sub->pGeom = getSubProjectionGeometryV(pGeom, newsubZ, size);
 
-			ret.push_back(boost::shared_ptr<CPart>(sub));
+			out.push_back(boost::shared_ptr<CPart>(sub));
 		}
 	}
 
-	return ret;
-
 }
 
 CCompositeGeometryManager::CProjectionPart* CCompositeGeometryManager::CProjectionPart::clone() const
diff --git a/src/CudaDataOperationAlgorithm.cpp b/src/CudaDataOperationAlgorithm.cpp
index 15886a4..82b676b 100644
--- a/src/CudaDataOperationAlgorithm.cpp
+++ b/src/CudaDataOperationAlgorithm.cpp
@@ -76,7 +76,7 @@ bool CCudaDataOperationAlgorithm::initialize(const Config& _cfg)
 	node = _cfg.self.getSingleNode("DataId");
 	ASTRA_CONFIG_CHECK(node, "CCudaDataOperationAlgorithm", "No DataId tag specified.");
 	vector<string> data = node.getContentArray();
-	for (vector<string>::iterator it = data.begin(); it != data.end(); it++){
+	for (vector<string>::iterator it = data.begin(); it != data.end(); ++it){
 		int id = StringUtil::stringToInt(*it);
 		m_pData.push_back(dynamic_cast<CFloat32Data2D*>(CData2DManager::getSingleton().get(id)));
 	}
diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp
index 5a1910c..2798434 100644
--- a/src/CudaReconstructionAlgorithm2D.cpp
+++ b/src/CudaReconstructionAlgorithm2D.cpp
@@ -329,6 +329,20 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry()
 }
 
 //----------------------------------------------------------------------------------------
+
+void CCudaReconstructionAlgorithm2D::initCUDAAlgorithm()
+{
+	bool ok;
+
+	ok = setupGeometry();
+	ASTRA_ASSERT(ok);
+
+	ok = m_pAlgo->allocateBuffers();
+	ASTRA_ASSERT(ok);
+}
+
+
+//----------------------------------------------------------------------------------------
 // Iterate
 void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 {
@@ -339,13 +353,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations)
 	const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
 
 	if (!m_bAlgoInit) {
-
-		ok = setupGeometry();
-		ASTRA_ASSERT(ok);
-
-		ok = m_pAlgo->allocateBuffers();
-		ASTRA_ASSERT(ok);
-
+		initCUDAAlgorithm();
 		m_bAlgoInit = true;
 	}
 
diff --git a/src/CudaSartAlgorithm.cpp b/src/CudaSartAlgorithm.cpp
index d202847..bf97224 100644
--- a/src/CudaSartAlgorithm.cpp
+++ b/src/CudaSartAlgorithm.cpp
@@ -107,7 +107,8 @@ bool CCudaSartAlgorithm::initialize(const Config& _cfg)
 		CC.markOptionParsed("ProjectionOrderList");
 	}
 
-
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+	CC.markOptionParsed("Relaxation");
 
 	return true;
 }
@@ -123,12 +124,26 @@ bool CCudaSartAlgorithm::initialize(CProjector2D* _pProjector,
 	if (!m_bIsInitialized)
 		return false;
 
+	m_fLambda = 1.0f;
+
 	m_pAlgo = new astraCUDA::SART();
 	m_bAlgoInit = false;
 
 	return true;
 }
 
+//----------------------------------------------------------------------------------------
+
+void CCudaSartAlgorithm::initCUDAAlgorithm()
+{
+	CCudaReconstructionAlgorithm2D::initCUDAAlgorithm();
+
+	astraCUDA::SART* pSart = dynamic_cast<astraCUDA::SART*>(m_pAlgo);
+
+	pSart->setRelaxation(m_fLambda);
+}
+
+
 
 } // namespace astra
 
diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp
index 33e381a..c8dc677 100644
--- a/src/CudaSirtAlgorithm.cpp
+++ b/src/CudaSirtAlgorithm.cpp
@@ -50,6 +50,8 @@ CCudaSirtAlgorithm::CCudaSirtAlgorithm()
 
 	m_pMinMask = 0;
 	m_pMaxMask = 0;
+
+	m_fLambda = 1.0f;
 }
 
 //----------------------------------------------------------------------------------------
@@ -86,6 +88,8 @@ bool CCudaSirtAlgorithm::initialize(const Config& _cfg)
 	}
 	CC.markOptionParsed("MaxMaskId");
 
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+	CC.markOptionParsed("Relaxation");
 
 	m_pAlgo = new astraCUDA::SIRT();
 	m_bAlgoInit = false;
@@ -108,41 +112,30 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector,
 
 	m_pAlgo = new astraCUDA::SIRT();
 	m_bAlgoInit = false;
+	m_fLambda = 1.0f;
 
 	return true;
 }
 
 //----------------------------------------------------------------------------------------
-// Iterate
-void CCudaSirtAlgorithm::run(int _iNrIterations)
+
+void CCudaSirtAlgorithm::initCUDAAlgorithm()
 {
-	// check initialized
-	ASTRA_ASSERT(m_bIsInitialized);
+	CCudaReconstructionAlgorithm2D::initCUDAAlgorithm();
 
-	if (!m_bAlgoInit) {
-		// We only override the initialisation step to copy the min/max masks
+	astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo);
 
-		bool ok = setupGeometry();
+	if (m_pMinMask || m_pMaxMask) {
+		const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
+		const float *pfMinMaskData = 0;
+		const float *pfMaxMaskData = 0;
+		if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst();
+		if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst();
+		bool ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount());
 		ASTRA_ASSERT(ok);
-
-		ok = m_pAlgo->allocateBuffers();
-		ASTRA_ASSERT(ok);
-
-		if (m_pMinMask || m_pMaxMask) {
-			const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry();
-			astraCUDA::SIRT* pSirt = dynamic_cast<astraCUDA::SIRT*>(m_pAlgo);
-			const float *pfMinMaskData = 0;
-			const float *pfMaxMaskData = 0;
-			if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst();
-			if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst();
-			ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount());
-			ASTRA_ASSERT(ok);
-		}
-
-		m_bAlgoInit = true;
 	}
 
-	CCudaReconstructionAlgorithm2D::run(_iNrIterations);
+	pSirt->setRelaxation(m_fLambda);
 }
 
 
diff --git a/src/CudaSirtAlgorithm3D.cpp b/src/CudaSirtAlgorithm3D.cpp
index 605c470..c819f8e 100644
--- a/src/CudaSirtAlgorithm3D.cpp
+++ b/src/CudaSirtAlgorithm3D.cpp
@@ -56,6 +56,7 @@ CCudaSirtAlgorithm3D::CCudaSirtAlgorithm3D()
 	m_iGPUIndex = -1;
 	m_iVoxelSuperSampling = 1;
 	m_iDetectorSuperSampling = 1;
+	m_fLambda = 1.0f;
 }
 
 //----------------------------------------------------------------------------------------
@@ -128,6 +129,8 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
 		return false;
 	}
 
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation");
+
 	initializeFromProjector();
 
 	// Deprecated options
@@ -135,6 +138,7 @@ bool CCudaSirtAlgorithm3D::initialize(const Config& _cfg)
 	m_iDetectorSuperSampling = (int)_cfg.self.getOptionNumerical("DetectorSuperSampling", m_iDetectorSuperSampling);
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUindex", m_iGPUIndex);
 	m_iGPUIndex = (int)_cfg.self.getOptionNumerical("GPUIndex", m_iGPUIndex);
+
 	CC.markOptionParsed("VoxelSuperSampling");
 	CC.markOptionParsed("DetectorSuperSampling");
 	CC.markOptionParsed("GPUIndex");
@@ -164,6 +168,8 @@ bool CCudaSirtAlgorithm3D::initialize(CProjector3D* _pProjector,
 		clear();
 	}
 
+	m_fLambda = 1.0f;
+
 	// required classes
 	m_pProjector = _pProjector;
 	m_pSinogram = _pSinogram;
@@ -224,6 +230,8 @@ void CCudaSirtAlgorithm3D::run(int _iNrIterations)
 
 		ASTRA_ASSERT(ok);
 
+		m_pSirt->setRelaxation(m_fLambda);
+
 		m_bAstraSIRTInit = true;
 
 	}
diff --git a/src/FilteredBackProjectionAlgorithm.cpp b/src/FilteredBackProjectionAlgorithm.cpp
index c195578..ccbfec6 100644
--- a/src/FilteredBackProjectionAlgorithm.cpp
+++ b/src/FilteredBackProjectionAlgorithm.cpp
@@ -117,12 +117,10 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 		int angleCount = projectionIndex.size();
 		int detectorCount = m_pProjector->getProjectionGeometry()->getDetectorCount();
 
+		// TODO: There is no need to allocate this. Better just
+		// create the CFloat32ProjectionData2D object directly, and use its
+		// memory.
 		float32 * sinogramData2D = new float32[angleCount* detectorCount];
-		float32 ** sinogramData = new float32*[angleCount];
-		for (int i = 0; i < angleCount; i++)
-		{
-			sinogramData[i] = &(sinogramData2D[i * detectorCount]);
-		}
 
 		float32 * projectionAngles = new float32[angleCount];
 		float32 detectorWidth = m_pProjector->getProjectionGeometry()->getDetectorWidth();
@@ -130,6 +128,8 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 		for (int i = 0; i < angleCount; i ++) {
 			if (projectionIndex[i] > m_pProjector->getProjectionGeometry()->getProjectionAngleCount() -1 )
 			{
+				delete[] sinogramData2D;
+				delete[] projectionAngles;
 				ASTRA_ERROR("Invalid Projection Index");
 				return false;
 			} else {
@@ -139,7 +139,6 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 				{
 					sinogramData2D[i*detectorCount+ iDetector] = m_pSinogram->getData2D()[orgIndex][iDetector];
 				}
-//				sinogramData[i] = m_pSinogram->getSingleProjectionData(projectionIndex[i]);
 				projectionAngles[i] = m_pProjector->getProjectionGeometry()->getProjectionAngle((int)projectionIndex[i] );
 
 			}
@@ -148,6 +147,9 @@ bool CFilteredBackProjectionAlgorithm::initialize(const Config& _cfg)
 		CParallelProjectionGeometry2D * pg = new CParallelProjectionGeometry2D(angleCount, detectorCount,detectorWidth,projectionAngles);
 		m_pProjector = new CParallelBeamLineKernelProjector2D(pg,m_pReconstruction->getGeometry());
 		m_pSinogram = new CFloat32ProjectionData2D(pg, sinogramData2D);
+
+		delete[] sinogramData2D;
+		delete[] projectionAngles;
 	}
 
 	// TODO: check that the angles are linearly spaced between 0 and pi
diff --git a/src/Float32VolumeData3DMemory.cpp b/src/Float32VolumeData3DMemory.cpp
index af45cb9..14adb1a 100644
--- a/src/Float32VolumeData3DMemory.cpp
+++ b/src/Float32VolumeData3DMemory.cpp
@@ -136,7 +136,6 @@ CFloat32VolumeData2D * CFloat32VolumeData3DMemory::fetchSliceZ(int _iSliceIndex)
 	CFloat32VolumeData2D* res = new CFloat32VolumeData2D(&volGeom);
 
 	// copy data
-	int iSliceCount = m_pGeometry->getGridSliceCount();
 	float * pfTargetData = res->getData();
 	for(int iRowIndex = 0; iRowIndex < iRowCount; iRowIndex++)
 	{
diff --git a/src/PluginAlgorithm.cpp b/src/PluginAlgorithm.cpp
index 9fc511a..1bcfbdb 100644
--- a/src/PluginAlgorithm.cpp
+++ b/src/PluginAlgorithm.cpp
@@ -26,376 +26,11 @@ along with the ASTRA Toolbox. If not, see <http://www.gnu.org/licenses/>.
 $Id$
 */
 
-#ifdef ASTRA_PYTHON
-
 #include "astra/PluginAlgorithm.h"
-#include "astra/Logging.h"
-#include "astra/Utilities.h"
-#include <boost/algorithm/string.hpp>
-#include <boost/algorithm/string/split.hpp>
-#include <iostream>
-#include <fstream>
-#include <string>
-
-#include <Python.h>
-#include "bytesobject.h"
 
 namespace astra {
 
+CPluginAlgorithmFactory *CPluginAlgorithmFactory::m_factory = 0;
 
-
-void logPythonError(){
-    if(PyErr_Occurred()){
-        PyObject *ptype, *pvalue, *ptraceback;
-        PyErr_Fetch(&ptype, &pvalue, &ptraceback);
-        PyErr_NormalizeException(&ptype, &pvalue, &ptraceback);
-        PyObject *traceback = PyImport_ImportModule("traceback");
-        if(traceback!=NULL){
-            PyObject *exc;
-            if(ptraceback==NULL){
-                exc = PyObject_CallMethod(traceback,"format_exception_only","OO",ptype, pvalue);
-            }else{
-                exc = PyObject_CallMethod(traceback,"format_exception","OOO",ptype, pvalue, ptraceback);
-            }
-            if(exc!=NULL){
-                PyObject *six = PyImport_ImportModule("six");
-                if(six!=NULL){
-                    PyObject *iter = PyObject_GetIter(exc);
-                    if(iter!=NULL){
-                        PyObject *line;
-                        std::string errStr = "";
-                        while(line = PyIter_Next(iter)){
-                            PyObject *retb = PyObject_CallMethod(six,"b","O",line);
-                            if(retb!=NULL){
-                                errStr += std::string(PyBytes_AsString(retb));
-                                Py_DECREF(retb);
-                            }
-                            Py_DECREF(line);
-                        }
-                        ASTRA_ERROR("%s",errStr.c_str());
-                        Py_DECREF(iter);
-                    }
-                    Py_DECREF(six);
-                }
-                Py_DECREF(exc);
-            }
-            Py_DECREF(traceback);
-        }
-        if(ptype!=NULL) Py_DECREF(ptype);
-        if(pvalue!=NULL) Py_DECREF(pvalue);
-        if(ptraceback!=NULL) Py_DECREF(ptraceback);
-    }
-}
-
-
-CPluginAlgorithm::CPluginAlgorithm(PyObject* pyclass){
-    instance = PyObject_CallObject(pyclass, NULL);
-    if(instance==NULL) logPythonError();
-}
-
-CPluginAlgorithm::~CPluginAlgorithm(){
-    if(instance!=NULL){
-        Py_DECREF(instance);
-        instance = NULL;
-    }
-}
-
-bool CPluginAlgorithm::initialize(const Config& _cfg){
-    if(instance==NULL) return false;
-    PyObject *cfgDict = XMLNode2dict(_cfg.self);
-    PyObject *retVal = PyObject_CallMethod(instance, "astra_init", "O",cfgDict);
-    Py_DECREF(cfgDict);
-    if(retVal==NULL){
-        logPythonError();
-        return false;
-    }
-    m_bIsInitialized = true;
-    Py_DECREF(retVal);
-    return m_bIsInitialized;
-}
-
-void CPluginAlgorithm::run(int _iNrIterations){
-    if(instance==NULL) return;
-    PyGILState_STATE state = PyGILState_Ensure();
-    PyObject *retVal = PyObject_CallMethod(instance, "run", "i",_iNrIterations);
-    if(retVal==NULL){
-        logPythonError();
-    }else{
-        Py_DECREF(retVal);
-    }
-    PyGILState_Release(state);
 }
 
-void fixLapackLoading(){
-    // When running in Matlab, we need to force numpy
-    // to use its internal lapack library instead of
-    // Matlab's MKL library to avoid errors. To do this,
-    // we set Python's dlopen flags to RTLD_NOW|RTLD_DEEPBIND
-    // and import 'numpy.linalg.lapack_lite' here. We reset
-    // Python's dlopen flags afterwards.
-    PyObject *sys = PyImport_ImportModule("sys");
-    if(sys!=NULL){
-        PyObject *curFlags = PyObject_CallMethod(sys,"getdlopenflags",NULL);
-        if(curFlags!=NULL){
-            PyObject *retVal = PyObject_CallMethod(sys, "setdlopenflags", "i",10);
-            if(retVal!=NULL){
-                PyObject *lapack = PyImport_ImportModule("numpy.linalg.lapack_lite");
-                if(lapack!=NULL){
-                    Py_DECREF(lapack);
-                }
-                PyObject_CallMethod(sys, "setdlopenflags", "O",curFlags);
-                Py_DECREF(retVal);
-            }
-            Py_DECREF(curFlags);
-        }
-        Py_DECREF(sys);
-    }
-}
-
-CPluginAlgorithmFactory::CPluginAlgorithmFactory(){
-    if(!Py_IsInitialized()){
-        Py_Initialize();
-        PyEval_InitThreads();
-    }
-#ifndef _MSC_VER
-    if(astra::running_in_matlab) fixLapackLoading();
-#endif
-    pluginDict = PyDict_New();
-    inspect = PyImport_ImportModule("inspect");
-    six = PyImport_ImportModule("six");
-}
-
-CPluginAlgorithmFactory::~CPluginAlgorithmFactory(){
-    if(pluginDict!=NULL){
-        Py_DECREF(pluginDict);
-    }
-    if(inspect!=NULL) Py_DECREF(inspect);
-    if(six!=NULL) Py_DECREF(six);
-}
-
-PyObject * getClassFromString(std::string str){
-    std::vector<std::string> items;
-    boost::split(items, str, boost::is_any_of("."));
-    PyObject *pyclass = PyImport_ImportModule(items[0].c_str());
-    if(pyclass==NULL){
-        logPythonError();
-        return NULL;
-    }
-    PyObject *submod = pyclass;
-    for(unsigned int i=1;i<items.size();i++){
-        submod = PyObject_GetAttrString(submod,items[i].c_str());
-        Py_DECREF(pyclass);
-        pyclass = submod;
-        if(pyclass==NULL){
-            logPythonError();
-            return NULL;
-        }
-    }
-    return pyclass;
-}
-
-bool CPluginAlgorithmFactory::registerPlugin(std::string name, std::string className){
-    PyObject *str = PyBytes_FromString(className.c_str());
-    PyDict_SetItemString(pluginDict, name.c_str(), str);
-    Py_DECREF(str);
-    return true;
-}
-
-bool CPluginAlgorithmFactory::registerPlugin(std::string className){
-    PyObject *pyclass = getClassFromString(className);
-    if(pyclass==NULL) return false;
-    bool ret = registerPluginClass(pyclass);
-    Py_DECREF(pyclass);
-    return ret;
-}
-
-bool CPluginAlgorithmFactory::registerPluginClass(std::string name, PyObject * className){
-    PyDict_SetItemString(pluginDict, name.c_str(), className);
-    return true;
-}
-
-bool CPluginAlgorithmFactory::registerPluginClass(PyObject * className){
-    PyObject *astra_name = PyObject_GetAttrString(className,"astra_name");
-    if(astra_name==NULL){
-        logPythonError();
-        return false;
-    }
-    PyObject *retb = PyObject_CallMethod(six,"b","O",astra_name);
-    if(retb!=NULL){
-        PyDict_SetItemString(pluginDict,PyBytes_AsString(retb),className);
-        Py_DECREF(retb);
-    }else{
-        logPythonError();
-    }
-    Py_DECREF(astra_name);
-    return true;
-}
-
-CPluginAlgorithm * CPluginAlgorithmFactory::getPlugin(std::string name){
-    PyObject *className = PyDict_GetItemString(pluginDict, name.c_str());
-    if(className==NULL) return NULL;
-    CPluginAlgorithm *alg = NULL;
-    if(PyBytes_Check(className)){
-        std::string str = std::string(PyBytes_AsString(className));
-    	PyObject *pyclass = getClassFromString(str);
-        if(pyclass!=NULL){
-            alg = new CPluginAlgorithm(pyclass);
-            Py_DECREF(pyclass);
-        }
-    }else{
-        alg = new CPluginAlgorithm(className);
-    }
-    return alg;
-}
-
-PyObject * CPluginAlgorithmFactory::getRegistered(){
-    Py_INCREF(pluginDict);
-    return pluginDict;
-}
-
-std::map<std::string, std::string> CPluginAlgorithmFactory::getRegisteredMap(){
-    std::map<std::string, std::string> ret;
-    PyObject *key, *value;
-    Py_ssize_t pos = 0;
-    while (PyDict_Next(pluginDict, &pos, &key, &value)) {
-        PyObject *keystr = PyObject_Str(key);
-        if(keystr!=NULL){
-            PyObject *valstr = PyObject_Str(value);
-            if(valstr!=NULL){
-                PyObject * keyb = PyObject_CallMethod(six,"b","O",keystr);
-                if(keyb!=NULL){
-                    PyObject * valb = PyObject_CallMethod(six,"b","O",valstr);
-                    if(valb!=NULL){
-                        ret[PyBytes_AsString(keyb)] = PyBytes_AsString(valb);
-                        Py_DECREF(valb);
-                    }
-                    Py_DECREF(keyb);
-                }
-                Py_DECREF(valstr);
-            }
-            Py_DECREF(keystr);
-        }
-        logPythonError();
-    }
-    return ret;
-}
-
-std::string CPluginAlgorithmFactory::getHelp(std::string name){
-    PyObject *className = PyDict_GetItemString(pluginDict, name.c_str());
-    if(className==NULL){
-        ASTRA_ERROR("Plugin %s not found!",name.c_str());
-        PyErr_Clear();
-        return "";
-    }
-    std::string ret = "";
-    PyObject *pyclass;
-    if(PyBytes_Check(className)){
-        std::string str = std::string(PyBytes_AsString(className));
-        pyclass = getClassFromString(str);
-    }else{
-        pyclass = className;
-    }
-    if(pyclass==NULL) return "";
-    if(inspect!=NULL && six!=NULL){
-        PyObject *retVal = PyObject_CallMethod(inspect,"getdoc","O",pyclass);
-        if(retVal!=NULL){
-            if(retVal!=Py_None){
-                PyObject *retb = PyObject_CallMethod(six,"b","O",retVal);
-                if(retb!=NULL){
-                    ret = std::string(PyBytes_AsString(retb));
-                    Py_DECREF(retb);
-                }
-            }
-            Py_DECREF(retVal);
-        }else{
-            logPythonError();
-        }
-    }
-    if(PyBytes_Check(className)){
-        Py_DECREF(pyclass);
-    }
-    return ret;
-}
-
-DEFINE_SINGLETON(CPluginAlgorithmFactory);
-
-#if PY_MAJOR_VERSION >= 3
-PyObject * pyStringFromString(std::string str){
-    return PyUnicode_FromString(str.c_str());
-}
-#else
-PyObject * pyStringFromString(std::string str){
-    return PyBytes_FromString(str.c_str());
-}
-#endif
-
-PyObject* stringToPythonValue(std::string str){
-    if(str.find(";")!=std::string::npos){
-        std::vector<std::string> rows, row;
-        boost::split(rows, str, boost::is_any_of(";"));
-        PyObject *mat = PyList_New(rows.size());
-        for(unsigned int i=0; i<rows.size(); i++){
-            boost::split(row, rows[i], boost::is_any_of(","));
-            PyObject *rowlist = PyList_New(row.size());
-            for(unsigned int j=0;j<row.size();j++){
-                PyList_SetItem(rowlist, j, PyFloat_FromDouble(StringUtil::stringToDouble(row[j])));
-            }
-            PyList_SetItem(mat, i, rowlist);
-        }
-        return mat;
-    }
-    if(str.find(",")!=std::string::npos){
-        std::vector<std::string> vec;
-        boost::split(vec, str, boost::is_any_of(","));
-        PyObject *veclist = PyList_New(vec.size());
-        for(unsigned int i=0;i<vec.size();i++){
-            PyList_SetItem(veclist, i, PyFloat_FromDouble(StringUtil::stringToDouble(vec[i])));
-        }
-        return veclist;
-    }
-    try{
-        return PyLong_FromLong(StringUtil::stringToInt(str));
-    }catch(const StringUtil::bad_cast &){
-        try{
-            return PyFloat_FromDouble(StringUtil::stringToDouble(str));
-        }catch(const StringUtil::bad_cast &){
-            return pyStringFromString(str);
-        }
-    }
-}
-
-PyObject* XMLNode2dict(XMLNode node){
-    PyObject *dct = PyDict_New();
-    PyObject *opts = PyDict_New();
-    if(node.hasAttribute("type")){
-        PyObject *obj = pyStringFromString(node.getAttribute("type").c_str());
-        PyDict_SetItemString(dct, "type", obj);
-        Py_DECREF(obj);
-    }
-    std::list<XMLNode> nodes = node.getNodes();
-    std::list<XMLNode>::iterator it = nodes.begin();
-    while(it!=nodes.end()){
-        XMLNode subnode = *it;
-        if(subnode.getName()=="Option"){
-            PyObject *obj;
-            if(subnode.hasAttribute("value")){
-                obj = stringToPythonValue(subnode.getAttribute("value"));
-            }else{
-                obj = stringToPythonValue(subnode.getContent());
-            }
-            PyDict_SetItemString(opts, subnode.getAttribute("key").c_str(), obj);
-            Py_DECREF(obj);
-        }else{
-            PyObject *obj = stringToPythonValue(subnode.getContent());
-            PyDict_SetItemString(dct, subnode.getName().c_str(), obj);
-            Py_DECREF(obj);
-        }
-        ++it;
-    }
-    PyDict_SetItemString(dct, "options", opts);
-    Py_DECREF(opts);
-    return dct;
-}
-
-}
-#endif
diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp
index 9346160..8df3342 100644
--- a/src/SartAlgorithm.cpp
+++ b/src/SartAlgorithm.cpp
@@ -151,6 +151,9 @@ bool CSartAlgorithm::initialize(const Config& _cfg)
 		CC.markOptionParsed("ProjectionOrderList");
 	}
 
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+	CC.markOptionParsed("Relaxation");
+
 	// create data objects
 	m_pTotalRayLength = new CFloat32ProjectionData2D(m_pProjector->getProjectionGeometry());
 	m_pTotalPixelWeight = new CFloat32VolumeData2D(m_pProjector->getVolumeGeometry());
@@ -246,6 +249,7 @@ map<string,boost::any> CSartAlgorithm::getInformation()
 {
 	map<string, boost::any> res;
 	res["ProjectionOrder"] = getInformation("ProjectionOrder");
+	res["Relaxation"] = getInformation("Relaxation");
 	return mergeMap<string,boost::any>(CReconstructionAlgorithm2D::getInformation(), res);
 };
 
@@ -253,6 +257,8 @@ map<string,boost::any> CSartAlgorithm::getInformation()
 // Information - Specific
 boost::any CSartAlgorithm::getInformation(std::string _sIdentifier) 
 {
+	if (_sIdentifier == "Relaxation")
+		return m_fLambda;
 	if (_sIdentifier == "ProjectionOrder") {
 		vector<float32> res;
 		for (int i = 0; i < m_iProjectionCount; i++) {
@@ -272,9 +278,8 @@ void CSartAlgorithm::run(int _iNrIterations)
 
 	m_bShouldAbort = false;
 
-	int iIteration = 0;
-
 	// data projectors
+	CDataProjectorInterface* pFirstForwardProjector;
 	CDataProjectorInterface* pForwardProjector;
 	CDataProjectorInterface* pBackProjector;
 
@@ -286,13 +291,13 @@ void CSartAlgorithm::run(int _iNrIterations)
 			m_pProjector, 
 			SinogramMaskPolicy(m_pSinogramMask),														// sinogram mask
 			ReconstructionMaskPolicy(m_pReconstructionMask),											// reconstruction mask
-			SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength),	// SIRT backprojection
+			SIRTBPPolicy(m_pReconstruction, m_pDiffSinogram, m_pTotalPixelWeight, m_pTotalRayLength, m_fLambda),	// SIRT backprojection
 			m_bUseSinogramMask, m_bUseReconstructionMask, true // options on/off
 		); 
 
 	// first time forward projection data projector,
 	// also computes total pixel weight and total ray length
-	pForwardProjector = dispatchDataProjector(
+	pFirstForwardProjector = dispatchDataProjector(
 			m_pProjector, 
 			SinogramMaskPolicy(m_pSinogramMask),														// sinogram mask
 			ReconstructionMaskPolicy(m_pReconstructionMask),											// reconstruction mask
@@ -303,16 +308,30 @@ void CSartAlgorithm::run(int _iNrIterations)
 			m_bUseSinogramMask, m_bUseReconstructionMask, true											 // options on/off
 		);
 
+	// forward projection data projector
+	pForwardProjector = dispatchDataProjector(
+			m_pProjector,
+			SinogramMaskPolicy(m_pSinogramMask),														// sinogram mask
+			ReconstructionMaskPolicy(m_pReconstructionMask),											// reconstruction mask
+			CombinePolicy<DiffFPPolicy, TotalPixelWeightPolicy>(					// 2 basic operations
+				DiffFPPolicy(m_pReconstruction, m_pDiffSinogram, m_pSinogram),								// forward projection with difference calculation
+				TotalPixelWeightPolicy(m_pTotalPixelWeight)),												// calculate the total pixel weights
+			m_bUseSinogramMask, m_bUseReconstructionMask, true											 // options on/off
+		);
+
 
 
 	// iteration loop
-	for (; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) {
+	for (int iIteration = 0; iIteration < _iNrIterations && !m_bShouldAbort; ++iIteration) {
 
 		int iProjection = m_piProjectionOrder[m_iIterationCount % m_iProjectionCount];
 	
 		// forward projection and difference calculation
 		m_pTotalPixelWeight->setData(0.0f);
-		pForwardProjector->projectSingleProjection(iProjection);
+		if (iIteration < m_iProjectionCount)
+			pFirstForwardProjector->projectSingleProjection(iProjection);
+		else
+			pForwardProjector->projectSingleProjection(iProjection);
 		// backprojection
 		pBackProjector->projectSingleProjection(iProjection);
 		// update iteration count
@@ -325,6 +344,7 @@ void CSartAlgorithm::run(int _iNrIterations)
 	}
 
 
+	ASTRA_DELETE(pFirstForwardProjector);
 	ASTRA_DELETE(pForwardProjector);
 	ASTRA_DELETE(pBackProjector);
 
diff --git a/src/SirtAlgorithm.cpp b/src/SirtAlgorithm.cpp
index d9f3a65..ff25648 100644
--- a/src/SirtAlgorithm.cpp
+++ b/src/SirtAlgorithm.cpp
@@ -76,6 +76,7 @@ void CSirtAlgorithm::_clear()
 	m_pDiffSinogram = NULL;
 	m_pTmpVolume = NULL;
 
+	m_fLambda = 1.0f;
 	m_iIterationCount = 0;
 }
 
@@ -91,6 +92,7 @@ void CSirtAlgorithm::clear()
 	ASTRA_DELETE(m_pDiffSinogram);
 	ASTRA_DELETE(m_pTmpVolume);
 
+	m_fLambda = 1.0f;
 	m_iIterationCount = 0;
 }
 
@@ -128,6 +130,9 @@ bool CSirtAlgorithm::initialize(const Config& _cfg)
 		return false;
 	}
 
+	m_fLambda = _cfg.self.getOptionNumerical("Relaxation", 1.0f);
+	CC.markOptionParsed("Relaxation");
+
 	// init data objects and data projectors
 	_init();
 
@@ -152,6 +157,8 @@ bool CSirtAlgorithm::initialize(CProjector2D* _pProjector,
 	m_pSinogram = _pSinogram;
 	m_pReconstruction = _pReconstruction;
 
+	m_fLambda = 1.0f;
+
 	// init data objects and data projectors
 	_init();
 
@@ -248,7 +255,7 @@ void CSirtAlgorithm::run(int _iNrIterations)
 			x = 1.0f / x;
 		else
 			x = 0.0f;
-		pfT[i] = x;
+		pfT[i] = m_fLambda * x;
 	}
 	pfT = m_pTotalRayLength->getData();
 	for (int i = 0; i < m_pTotalRayLength->getSize(); ++i) {
@@ -296,7 +303,7 @@ void CSirtAlgorithm::run(int _iNrIterations)
 		m_pTmpVolume->setData(0.0f);
 		pBackProjector->project();
 
-		// divide by pixel weights
+		// multiply with relaxation factor divided by pixel weights
 		(*m_pTmpVolume) *= (*m_pTotalPixelWeight);
 		(*m_pReconstruction) += (*m_pTmpVolume);
 
diff --git a/src/Utilities.cpp b/src/Utilities.cpp
index 4b80503..c9740bf 100644
--- a/src/Utilities.cpp
+++ b/src/Utilities.cpp
@@ -42,7 +42,7 @@ namespace StringUtil {
 
 int stringToInt(const std::string& s)
 {
-	double i;
+	int i;
 	std::istringstream iss(s);
 	iss.imbue(std::locale::classic());
 	iss >> i;
diff --git a/src/XMLNode.cpp b/src/XMLNode.cpp
index 40a9b22..cf268c2 100644
--- a/src/XMLNode.cpp
+++ b/src/XMLNode.cpp
@@ -158,7 +158,7 @@ vector<string> XMLNode::getContentArray() const
 	vector<string> res(iSize);
 	// loop all list item nodes
 	list<XMLNode> nodes = getNodes("ListItem");
-	for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) {
+	for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) {
 		int iIndex = it->getAttributeNumerical("index");
 		string sValue = it->getAttribute("value");
 		ASTRA_ASSERT(iIndex < iSize);
@@ -290,7 +290,7 @@ vector<float32> XMLNode::getOptionNumericalArray(string _sKey) const
 	if (!hasOption(_sKey)) return vector<float32>();
 
 	list<XMLNode> nodes = getNodes("Option");
-	for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); it++) {
+	for (list<XMLNode>::iterator it = nodes.begin(); it != nodes.end(); ++it) {
 		if (it->getAttribute("key") == _sKey) {
 			vector<float32> vals = it->getContentNumericalArray();
 			return vals;