Add relaxation parameters to SIRT3D

6 files changed, 37 insertions, 2 deletions

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/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/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;
 
 	}