From 5edb35edc2c721b458334a65512b534912c2c542 Mon Sep 17 00:00:00 2001
From: Willem Jan Palenstijn <Willem.Jan.Palenstijn@cwi.nl>
Date: Wed, 23 Mar 2016 15:30:56 +0100
Subject: Add relaxation parameters to SIRT, SART
---
cuda/2d/sart.cu | 7 +++++--
cuda/2d/sart.h | 4 ++++
cuda/2d/sirt.cu | 8 ++++++++
cuda/2d/sirt.h | 4 ++++
include/astra/CudaSartAlgorithm.h | 10 ++++++++++
include/astra/CudaSirtAlgorithm.h | 6 ++++++
include/astra/DataProjectorPolicies.h | 4 +++-
include/astra/DataProjectorPolicies.inl | 6 ++++--
include/astra/SartAlgorithm.h | 8 ++++++--
include/astra/SirtAlgorithm.h | 9 ++++++++-
src/CudaSartAlgorithm.cpp | 17 ++++++++++++++++-
src/CudaSirtAlgorithm.cpp | 6 ++++++
src/SartAlgorithm.cpp | 8 +++++++-
src/SirtAlgorithm.cpp | 11 +++++++++--
14 files changed, 96 insertions(+), 12 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/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 {
* <Option key="UseMinConstraint" value="yes"/>\n
* <Option key="UseMaxConstraint" value="yes"/>\n
* <Option key="MaxConstraintValue" value="1024"/>\n
+ * <Option key="Relaxation" value="1"/>\n
* </Algorithm>
* }
*
@@ -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/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 7beb30e..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,6 +112,7 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector,
m_pAlgo = new astraCUDA::SIRT();
m_bAlgoInit = false;
+ m_fLambda = 1.0f;
return true;
}
@@ -130,6 +135,7 @@ void CCudaSirtAlgorithm::initCUDAAlgorithm()
ASTRA_ASSERT(ok);
}
+ pSirt->setRelaxation(m_fLambda);
}
diff --git a/src/SartAlgorithm.cpp b/src/SartAlgorithm.cpp
index 9346160..403f851 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++) {
@@ -286,7 +292,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);
--
cgit v1.2.3