From 5edb35edc2c721b458334a65512b534912c2c542 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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(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(D_pixelWeight, D_maskData, pixelPitch, dims); } + // Also fold the relaxation factor into pixel weights + processVol(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(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(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 CSartAlgorithm::getInformation() { map res; res["ProjectionOrder"] = getInformation("ProjectionOrder"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -253,6 +257,8 @@ map CSartAlgorithm::getInformation() // Information - Specific boost::any CSartAlgorithm::getInformation(std::string _sIdentifier) { + if (_sIdentifier == "Relaxation") + return m_fLambda; if (_sIdentifier == "ProjectionOrder") { vector 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