From 16430239d04ff738a21146c410918c285552543f Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Wed, 23 Mar 2016 15:50:24 +0100
Subject: Add relaxation parameters to SIRT3D
---
cuda/3d/astra3d.cu | 13 +++++++++++++
cuda/3d/astra3d.h | 2 ++
cuda/3d/sirt3d.cu | 8 +++++++-
cuda/3d/sirt3d.h | 5 +++++
include/astra/CudaSirtAlgorithm3D.h | 3 ++-
src/CudaSirtAlgorithm3D.cpp | 8 ++++++++
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;
}
--
cgit v1.2.3