Add relaxation parameters to SIRT, SART

6 files changed, 37 insertions, 6 deletions

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 91cc206..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
@@ -113,6 +115,10 @@ protected:
 	CFloat32VolumeData2D* m_pMinMask;
 	CFloat32VolumeData2D* m_pMaxMask;
 
+	/** Relaxation factor
+	 */
+	float m_fLambda;
+
 	virtual void initCUDAAlgorithm();
 };
 
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