Merge pull request #35 from wjp/relaxation

Add relaxation factor option to SIRT, SART

24 files changed, 173 insertions, 52 deletions

diff --git a/cuda/2d/sart.cu b/cuda/2d/sart.cu
index e5cb5bb..c8608a3 100644
--- a/cuda/2d/sart.cu
+++ b/cuda/2d/sart.cu
@@ -71,6 +71,8 @@ SART::SART() : ReconAlgo()
 	projectionCount = 0;
 	iteration = 0;
 	customOrder = false;
+
+	fRelaxation = 1.0f;
 }
 
 
@@ -98,6 +100,7 @@ void SART::reset()
 	projectionCount = 0;
 	iteration = 0;
 	customOrder = false;
+	fRelaxation = 1.0f;
 
 	ReconAlgo::reset();
 }
@@ -200,10 +203,10 @@ bool SART::iterate(unsigned int iterations)
 			// BP, mask, and add back
 			// TODO: Try putting the masking directly in the BP
 			zeroVolumeData(D_tmpData, tmpPitch, dims);
-			callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, 1.0f);
+			callBP_SART(D_tmpData, tmpPitch, D_projData, projPitch, angle, fRelaxation);
 			processVol<opAddMul>(D_volumeData, D_maskData, D_tmpData, volumePitch, dims);
 		} else {
-			callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, 1.0f);
+			callBP_SART(D_volumeData, volumePitch, D_projData, projPitch, angle, fRelaxation);
 		}
 
 		if (useMinConstraint)
diff --git a/cuda/2d/sart.h b/cuda/2d/sart.h
index 7dcd641..eff9ecf 100644
--- a/cuda/2d/sart.h
+++ b/cuda/2d/sart.h
@@ -50,6 +50,8 @@ public:
 
 	virtual float computeDiffNorm();
 
+	void setRelaxation(float r) { fRelaxation = r; }
+
 protected:
 	void reset();
 	bool precomputeWeights();
@@ -78,6 +80,8 @@ protected:
 	// Geometry-specific precomputed data
 	float* D_lineWeight;
 	unsigned int linePitch;
+
+	float fRelaxation;
 };
 
 }
diff --git a/cuda/2d/sirt.cu b/cuda/2d/sirt.cu
index 162ee77..4baaccb 100644
--- a/cuda/2d/sirt.cu
+++ b/cuda/2d/sirt.cu
@@ -50,6 +50,8 @@ SIRT::SIRT() : ReconAlgo()
 	D_minMaskData = 0;
 	D_maxMaskData = 0;
 
+	fRelaxation = 1.0f;
+
 	freeMinMaxMasks = false;
 }
 
@@ -83,6 +85,8 @@ void SIRT::reset()
 	useVolumeMask = false;
 	useSinogramMask = false;
 
+	fRelaxation = 1.0f;
+
 	ReconAlgo::reset();
 }
 
@@ -139,6 +143,9 @@ bool SIRT::precomputeWeights()
 		processVol<opMul>(D_pixelWeight, D_maskData, pixelPitch, dims);
 	}
 
+	// Also fold the relaxation factor into pixel weights
+	processVol<opMul>(D_pixelWeight, fRelaxation, pixelPitch, dims);
+
 	return true;
 }
 
@@ -253,6 +260,7 @@ bool SIRT::iterate(unsigned int iterations)
 
 		callBP(D_tmpData, tmpPitch, D_projData, projPitch, 1.0f);
 
+		// pixel weights also contain the volume mask and relaxation factor
 		processVol<opAddMul>(D_volumeData, D_pixelWeight, D_tmpData, volumePitch, dims);
 
 		if (useMinConstraint)
diff --git a/cuda/2d/sirt.h b/cuda/2d/sirt.h
index 21094a1..bc913f4 100644
--- a/cuda/2d/sirt.h
+++ b/cuda/2d/sirt.h
@@ -53,6 +53,8 @@ public:
 	bool uploadMinMaxMasks(const float* minMaskData, const float* maxMaskData,
 	                       unsigned int pitch);
 
+	void setRelaxation(float r) { fRelaxation = r; }
+
 	virtual bool iterate(unsigned int iterations);
 
 	virtual float computeDiffNorm();
@@ -81,6 +83,8 @@ protected:
 	unsigned int minMaskPitch;
 	float* D_maxMaskData;
 	unsigned int maxMaskPitch;
+
+	float fRelaxation;
 };
 
 bool doSIRT(float* D_volumeData, unsigned int volumePitch,
diff --git a/cuda/3d/astra3d.cu b/cuda/3d/astra3d.cu
index 8328229..5670873 100644
--- a/cuda/3d/astra3d.cu
+++ b/cuda/3d/astra3d.cu
@@ -267,6 +267,7 @@ public:
 	float fOriginDetectorDistance;
 	float fSourceZ;
 	float fDetSize;
+	float fRelaxation;
 
 	SConeProjection* projs;
 	SPar3DProjection* parprojs;
@@ -311,6 +312,8 @@ AstraSIRT3d::AstraSIRT3d()
 	pData->parprojs = 0;
 	pData->fOutputScale = 1.0f;
 
+	pData->fRelaxation = 1.0f;
+
 	pData->initialized = false;
 	pData->setStartReconstruction = false;
 
@@ -389,6 +392,14 @@ bool AstraSIRT3d::enableSuperSampling(unsigned int iVoxelSuperSampling,
 	return true;
 }
 
+void AstraSIRT3d::setRelaxation(float r)
+{
+	if (pData->initialized)
+		return;
+
+	pData->fRelaxation = r;
+}
+
 bool AstraSIRT3d::enableVolumeMask()
 {
 	if (pData->initialized)
@@ -448,6 +459,8 @@ bool AstraSIRT3d::init()
 	if (!ok)
 		return false;
 
+	pData->sirt.setRelaxation(pData->fRelaxation);
+
 	pData->D_volumeData = allocateVolumeData(pData->dims);
 	ok = pData->D_volumeData.ptr;
 	if (!ok)
diff --git a/cuda/3d/astra3d.h b/cuda/3d/astra3d.h
index 2782994..2137587 100644
--- a/cuda/3d/astra3d.h
+++ b/cuda/3d/astra3d.h
@@ -68,6 +68,8 @@ public:
 	bool enableSuperSampling(unsigned int iVoxelSuperSampling,
 	                         unsigned int iDetectorSuperSampling);
 
+	void setRelaxation(float r);
+
 	// Enable volume/sinogram masks
 	//
 	// This may optionally be called before init().
diff --git a/cuda/3d/sirt3d.cu b/cuda/3d/sirt3d.cu
index 484521e..713944b 100644
--- a/cuda/3d/sirt3d.cu
+++ b/cuda/3d/sirt3d.cu
@@ -59,6 +59,8 @@ SIRT::SIRT() : ReconAlgo3D()
 
 	useMinConstraint = false;
 	useMaxConstraint = false;
+
+	fRelaxation = 1.0f;
 }
 
 
@@ -89,6 +91,8 @@ void SIRT::reset()
 	useVolumeMask = false;
 	useSinogramMask = false;
 
+	fRelaxation = 1.0f;
+
 	ReconAlgo3D::reset();
 }
 
@@ -196,6 +200,8 @@ bool SIRT::precomputeWeights()
 		// scale pixel weights with mask to zero out masked pixels
 		processVol3D<opMul>(D_pixelWeight, D_maskData, dims);
 	}
+	processVol3D<opMul>(D_pixelWeight, fRelaxation, dims);
+
 
 	return true;
 }
@@ -307,7 +313,7 @@ bool SIRT::iterate(unsigned int iterations)
 	}
 #endif
 
-
+		// pixel weights also contain the volume mask and relaxation factor
 		processVol3D<opAddMul>(D_volumeData, D_tmpData, D_pixelWeight, dims);
 
 		if (useMinConstraint)
diff --git a/cuda/3d/sirt3d.h b/cuda/3d/sirt3d.h
index bb3864a..5e93deb 100644
--- a/cuda/3d/sirt3d.h
+++ b/cuda/3d/sirt3d.h
@@ -48,6 +48,9 @@ public:
 	// init should be called after setting all geometry
 	bool init();
 
+	// Set relaxation factor. This may be called after init and before iterate.
+	void setRelaxation(float r) { fRelaxation = r; }
+
 	// setVolumeMask should be called after init and before iterate,
 	// but only if enableVolumeMask was called before init.
 	// It may be called again after iterate.
@@ -91,6 +94,8 @@ protected:
 	float fMinConstraint;
 	float fMaxConstraint;
 
+	float fRelaxation;
+
 	cudaPitchedPtr D_maskData;
 	cudaPitchedPtr D_smaskData;
 
diff --git a/include/astra/ArtAlgorithm.h b/include/astra/ArtAlgorithm.h
index d232b95..1ad9f3f 100644
--- a/include/astra/ArtAlgorithm.h
+++ b/include/astra/ArtAlgorithm.h
@@ -59,7 +59,7 @@ namespace astra {
  * \astra_xml_item_option{MinConstraintValue, float, 0, Minimum constraint value.}
  * \astra_xml_item_option{UseMaxConstraint, bool, false, Use maximum value constraint.}
  * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.}
- * \astra_xml_item_option{Lamda, float, 1, The relaxation factor.}
+ * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.}
  * \astra_xml_item_option{RayOrder, string, "sequential", the order in which the rays are updated. 'sequential' or 'custom'}
  * \astra_xml_item_option{RayOrderList, n by 2 vector of float, not used, if RayOrder='custom': use this ray order.  Each row consist of a projection id and detector id.}
  * 
@@ -73,7 +73,7 @@ namespace astra {
  *		cfg.option.UseMinConstraint = 'yes';\n 
  *		cfg.option.UseMaxConstraint = 'yes';\n
  *		cfg.option.MaxConstraintValue = 1024;\n
- *		cfg.option.Lamda = 0.7;\n
+ *		cfg.option.Relaxation = 0.7;\n
  *		cfg.option.RayOrder = 'custom';\n
  *		cfg.option.RayOrderList = [0\,0; 0\,2; 1\,0];\n
  *		alg_id = astra_mex_algorithm('create'\, cfg);\n
diff --git a/include/astra/CudaReconstructionAlgorithm2D.h b/include/astra/CudaReconstructionAlgorithm2D.h
index dc93a1a..bb5f2a7 100644
--- a/include/astra/CudaReconstructionAlgorithm2D.h
+++ b/include/astra/CudaReconstructionAlgorithm2D.h
@@ -141,6 +141,9 @@ protected:
 	 */
 	bool setupGeometry();
 
+	/** Initialize CUDA algorithm. For internal use only.
+	 */
+	virtual void initCUDAAlgorithm();
 
 	/** The internally used CUDA algorithm object
 	 */ 
diff --git a/include/astra/CudaSartAlgorithm.h b/include/astra/CudaSartAlgorithm.h
index c22dc4f..e5e18f0 100644
--- a/include/astra/CudaSartAlgorithm.h
+++ b/include/astra/CudaSartAlgorithm.h
@@ -46,6 +46,7 @@ namespace astra {
  * \astra_xml_item{ProjectionDataId, integer, Identifier of a projection data object as it is stored in the DataManager.}
  * \astra_xml_item{ReconstructionDataId, integer, Identifier of a volume data object as it is stored in the DataManager.}
  * \astra_xml_item_option{ReconstructionMaskId, integer, not used, Identifier of a volume data object that acts as a reconstruction mask. 0 = reconstruct on this pixel. 1 = don't reconstruct on this pixel.}
+ * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.}
  *
  * \par MATLAB example
  * \astra_code{
@@ -53,6 +54,7 @@ namespace astra {
  *		cfg.ProjectionDataId = sino_id;\n
  *		cfg.ReconstructionDataId = recon_id;\n
  *		cfg.option.ReconstructionMaskId = mask_id;\n
+ *		cfg.option.Relaxation = 1.0;\n
  *		alg_id = astra_mex_algorithm('create'\, cfg);\n
  *		astra_mex_algorithm('iterate'\, alg_id\, 10);\n
  *		astra_mex_algorithm('delete'\, alg_id);\n
@@ -97,6 +99,14 @@ public:
 	 * @return description string
 	 */
 	virtual std::string description() const;
+
+protected:
+
+	/** Relaxation factor
+	 */
+	float m_fLambda;
+
+	virtual void initCUDAAlgorithm();
 };
 
 // inline functions
diff --git a/include/astra/CudaSirtAlgorithm.h b/include/astra/CudaSirtAlgorithm.h
index 929ac30..3cd8acc 100644
--- a/include/astra/CudaSirtAlgorithm.h
+++ b/include/astra/CudaSirtAlgorithm.h
@@ -53,6 +53,7 @@ namespace astra {
  * \astra_xml_item{ProjectionDataId, integer, Identifier of a projection data object as it is stored in the DataManager.}
  * \astra_xml_item{ReconstructionDataId, integer, Identifier of a volume data object as it is stored in the DataManager.}
  * \astra_xml_item_option{ReconstructionMaskId, integer, not used, Identifier of a volume data object that acts as a reconstruction mask. 0 = reconstruct on this pixel. 1 = don't reconstruct on this pixel.}
+ * \astra_xml_item_option{Relaxation, float, 1, Relaxation parameter.}
  *
  * \par MATLAB example
  * \astra_code{
@@ -62,6 +63,7 @@ namespace astra {
  *		cfg.ProjectionDataId = sino_id;\n
  *		cfg.ReconstructionDataId = recon_id;\n
  *		cfg.option.ReconstructionMaskId = mask_id;\n
+ *		cfg.option.Relaxation = 1.0;\n
  *		alg_id = astra_mex_algorithm('create'\, cfg);\n
  *		astra_mex_algorithm('iterate'\, alg_id\, 10);\n
  *		astra_mex_algorithm('delete'\, alg_id);\n
@@ -93,8 +95,6 @@ public:
 	 */
 	virtual bool initialize(const Config& _cfg);
 
-	virtual void run(int _iNrIterations);
-
 	/** Initialize class.
 	 *
 	 * @param _pProjector		Projector Object. (Optional)
@@ -114,6 +114,12 @@ public:
 protected:
 	CFloat32VolumeData2D* m_pMinMask;
 	CFloat32VolumeData2D* m_pMaxMask;
+
+	/** Relaxation factor
+	 */
+	float m_fLambda;
+
+	virtual void initCUDAAlgorithm();
 };
 
 // inline functions
diff --git a/include/astra/CudaSirtAlgorithm3D.h b/include/astra/CudaSirtAlgorithm3D.h
index 379720e..60191cd 100644
--- a/include/astra/CudaSirtAlgorithm3D.h
+++ b/include/astra/CudaSirtAlgorithm3D.h
@@ -50,7 +50,7 @@ class AstraSIRT3d;
  *
  * The update step of pixel \f$v_j\f$ for iteration \f$k\f$ is given by:
  * \f[
- *	v_j^{(k+1)} = v_j^{(k)} + \alpha \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}}
+ *	v_j^{(k+1)} = v_j^{(k)} + \lambda \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}}
  * \f]
  *
  * \par XML Configuration
@@ -175,6 +175,7 @@ protected:
 	bool m_bAstraSIRTInit;
 	int m_iDetectorSuperSampling;
 	int m_iVoxelSuperSampling;
+	float m_fLambda;
 
 	void initializeFromProjector();
 };
diff --git a/include/astra/DataProjectorPolicies.h b/include/astra/DataProjectorPolicies.h
index c258f5c..acfb36f 100644
--- a/include/astra/DataProjectorPolicies.h
+++ b/include/astra/DataProjectorPolicies.h
@@ -319,10 +319,12 @@ class SIRTBPPolicy {
 	CFloat32ProjectionData2D* m_pTotalRayLength;
 	CFloat32VolumeData2D* m_pTotalPixelWeight;
 
+	float m_fRelaxation;
+
 public:
 
 	FORCEINLINE SIRTBPPolicy();
-	FORCEINLINE SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pTotalPixelWeight, CFloat32ProjectionData2D* _pTotalRayLength); 
+	FORCEINLINE SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, CFloat32ProjectionData2D* _pSinogram, CFloat32VolumeData2D* _pTotalPixelWeight, CFloat32ProjectionData2D* _pTotalRayLength, float _fRelaxation);
 	FORCEINLINE ~SIRTBPPolicy();
 
 	FORCEINLINE bool rayPrior(int _iRayIndex);
diff --git a/include/astra/DataProjectorPolicies.inl b/include/astra/DataProjectorPolicies.inl
index 0c0eddd..f300761 100644
--- a/include/astra/DataProjectorPolicies.inl
+++ b/include/astra/DataProjectorPolicies.inl
@@ -712,12 +712,14 @@ SIRTBPPolicy::SIRTBPPolicy()
 SIRTBPPolicy::SIRTBPPolicy(CFloat32VolumeData2D* _pReconstruction, 
 						   CFloat32ProjectionData2D* _pSinogram, 
 						   CFloat32VolumeData2D* _pTotalPixelWeight, 
-						   CFloat32ProjectionData2D* _pTotalRayLength) 
+						   CFloat32ProjectionData2D* _pTotalRayLength,
+						   float _fRelaxation)
 {
 	m_pReconstruction = _pReconstruction;
 	m_pSinogram = _pSinogram;
 	m_pTotalPixelWeight = _pTotalPixelWeight;
 	m_pTotalRayLength = _pTotalRayLength;
+	m_fRelaxation = _fRelaxation;
 }
 //----------------------------------------------------------------------------------------	
 SIRTBPPolicy::~SIRTBPPolicy() 
@@ -739,7 +741,7 @@ void SIRTBPPolicy::addWeight(int _iRayIndex, int _iVolumeIndex, float32 _fWeight
 {
 	float32 fGammaBeta = m_pTotalPixelWeight->getData()[_iVolumeIndex] * m_pTotalRayLength->getData()[_iRayIndex];
 	if ((fGammaBeta > 0.001f) || (fGammaBeta < -0.001f)) {
-		m_pReconstruction->getData()[_iVolumeIndex] += _fWeight * m_pSinogram->getData()[_iRayIndex] / fGammaBeta;
+		m_pReconstruction->getData()[_iVolumeIndex] += _fWeight * m_fRelaxation * m_pSinogram->getData()[_iRayIndex] / fGammaBeta;
 	}
 }
 //----------------------------------------------------------------------------------------
diff --git a/include/astra/SartAlgorithm.h b/include/astra/SartAlgorithm.h
index eb4c61e..cdae029 100644
--- a/include/astra/SartAlgorithm.h
+++ b/include/astra/SartAlgorithm.h
@@ -49,7 +49,7 @@ namespace astra {
  *
  * The update step of pixel \f$v_j\f$ for projection \f$phi\f$ and iteration \f$k\f$ is given by:
  * \f[
- *	v_j^{(k+1)} = v_j^{(k)} + \frac{\sum_{p_i \in P_\phi} \left(  \lambda \frac{p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}} {\sum_{r=1}^{N}w_{ir} }    \right)} {\sum_{p_i \in P_\phi}w_{ij}}
+ *	v_j^{(k+1)} = v_j^{(k)} + \lambda \frac{\sum_{p_i \in P_\phi} \left(  \frac{p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}} {\sum_{r=1}^{N}w_{ir} }    \right)} {\sum_{p_i \in P_\phi}w_{ij}}
  * \f]
  *
  * \par XML Configuration
@@ -64,6 +64,7 @@ namespace astra {
  * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.}
  * \astra_xml_item_option{ProjectionOrder, string, "sequential", the order in which the projections are updated. 'sequential', 'random' or 'custom'}
  * \astra_xml_item_option{ProjectionOrderList, vector of float, not used, if ProjectionOrder='custom': use this order.}
+ * \astra_xml_item_option{Relaxation, float, 1, The relaxation parameter.}
  *
  * \par MATLAB example
  * \astra_code{
@@ -76,7 +77,8 @@ namespace astra {
  *		cfg.option.UseMaxConstraint = 'yes';\n
  *		cfg.option.MaxConstraintValue = 1024;\n
  *		cfg.option.ProjectionOrder = 'custom';\n
-*		cfg.option.ProjectionOrderList = randperm(100);\n
+ *		cfg.option.ProjectionOrderList = randperm(100);\n
+ *		cfg.option.Relaxation = 1.0;\n
  *		alg_id = astra_mex_algorithm('create'\, cfg);\n
  *		astra_mex_algorithm('iterate'\, alg_id\, 10);\n
  *		astra_mex_algorithm('delete'\, alg_id);\n
@@ -215,6 +217,8 @@ protected:
 	//< Current index in the projection order array.
 	int m_iCurrentProjection;
 
+	//< Relaxation parameter
+	float m_fLambda;
 };
 
 // inline functions
diff --git a/include/astra/SirtAlgorithm.h b/include/astra/SirtAlgorithm.h
index 05b3fa9..8044d09 100644
--- a/include/astra/SirtAlgorithm.h
+++ b/include/astra/SirtAlgorithm.h
@@ -49,7 +49,7 @@ namespace astra {
  *
  * The update step of pixel \f$v_j\f$ for iteration \f$k\f$ is given by:
  * \f[
- *	v_j^{(k+1)} = v_j^{(k)} + \alpha \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}}
+ *	v_j^{(k+1)} = v_j^{(k)} + \lambda \sum_{i=1}^{M} \left( \frac{w_{ij}\left( p_i - \sum_{r=1}^{N} w_{ir}v_r^{(k)}\right)}{\sum_{k=1}^{N} w_{ik}} \right) \frac{1}{\sum_{l=1}^{M}w_{lj}}
  * \f]
  *
  * \par XML Configuration
@@ -62,6 +62,7 @@ namespace astra {
  * \astra_xml_item_option{MinConstraintValue, float, 0, Minimum constraint value.}
  * \astra_xml_item_option{UseMaxConstraint, bool, false, Use maximum value constraint.}
  * \astra_xml_item_option{MaxConstraintValue, float, 255, Maximum constraint value.}
+ * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.}
  *
  * \par XML Example
  * \astra_code{
@@ -74,6 +75,7 @@ namespace astra {
  *		&lt;Option key="UseMinConstraint" value="yes"/&gt;\n
  *		&lt;Option key="UseMaxConstraint" value="yes"/&gt;\n
  *		&lt;Option key="MaxConstraintValue" value="1024"/&gt;\n
+ *		&lt;Option key="Relaxation" value="1"/&gt;\n
  *		&lt;/Algorithm&gt;
  * }
  *
@@ -88,6 +90,7 @@ namespace astra {
  *		cfg.option.UseMinConstraint = 'yes';\n 
  *		cfg.option.UseMaxConstraint = 'yes';\n
  *		cfg.option.MaxConstraintValue = 1024;\n
+ *		cfg.option.Relaxation = 1.0;\n
  *		alg_id = astra_mex_algorithm('create'\, cfg);\n
  *		astra_mex_algorithm('iterate'\, alg_id\, 10);\n
  *		astra_mex_algorithm('delete'\, alg_id);\n
@@ -136,6 +139,10 @@ protected:
 	 */
 	int m_iIterationCount;
 
+	/** Relaxation parameter
+	 */
+	float m_fLambda;
+
 public:
 	
 	// type of the algorithm, needed to register with CAlgorithmFactory
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/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/SartAlgorithm.cpp b/src/SartAlgorithm.cpp
index f80df61..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++) {
@@ -285,7 +291,7 @@ 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
 		); 
 
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);