From 495903529d473a9968c1333d5a515e3b94732f0b Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:29:43 +0100 Subject: Move CUDA algorithm initialization to its own function --- include/astra/CudaReconstructionAlgorithm2D.h | 3 +++ include/astra/CudaSirtAlgorithm.h | 4 +-- src/CudaReconstructionAlgorithm2D.cpp | 22 +++++++++++------ src/CudaSirtAlgorithm.cpp | 35 +++++++++------------------ 4 files changed, 31 insertions(+), 33 deletions(-) diff --git a/include/astra/CudaReconstructionAlgorithm2D.h b/include/astra/CudaReconstructionAlgorithm2D.h index dc93a1a..bb5f2a7 100644 --- a/include/astra/CudaReconstructionAlgorithm2D.h +++ b/include/astra/CudaReconstructionAlgorithm2D.h @@ -141,6 +141,9 @@ protected: */ bool setupGeometry(); + /** Initialize CUDA algorithm. For internal use only. + */ + virtual void initCUDAAlgorithm(); /** The internally used CUDA algorithm object */ diff --git a/include/astra/CudaSirtAlgorithm.h b/include/astra/CudaSirtAlgorithm.h index 929ac30..91cc206 100644 --- a/include/astra/CudaSirtAlgorithm.h +++ b/include/astra/CudaSirtAlgorithm.h @@ -93,8 +93,6 @@ public: */ virtual bool initialize(const Config& _cfg); - virtual void run(int _iNrIterations); - /** Initialize class. * * @param _pProjector Projector Object. (Optional) @@ -114,6 +112,8 @@ public: protected: CFloat32VolumeData2D* m_pMinMask; CFloat32VolumeData2D* m_pMaxMask; + + virtual void initCUDAAlgorithm(); }; // inline functions diff --git a/src/CudaReconstructionAlgorithm2D.cpp b/src/CudaReconstructionAlgorithm2D.cpp index 5a1910c..2798434 100644 --- a/src/CudaReconstructionAlgorithm2D.cpp +++ b/src/CudaReconstructionAlgorithm2D.cpp @@ -328,6 +328,20 @@ bool CCudaReconstructionAlgorithm2D::setupGeometry() return true; } +//---------------------------------------------------------------------------------------- + +void CCudaReconstructionAlgorithm2D::initCUDAAlgorithm() +{ + bool ok; + + ok = setupGeometry(); + ASTRA_ASSERT(ok); + + ok = m_pAlgo->allocateBuffers(); + ASTRA_ASSERT(ok); +} + + //---------------------------------------------------------------------------------------- // Iterate void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) @@ -339,13 +353,7 @@ void CCudaReconstructionAlgorithm2D::run(int _iNrIterations) const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); if (!m_bAlgoInit) { - - ok = setupGeometry(); - ASTRA_ASSERT(ok); - - ok = m_pAlgo->allocateBuffers(); - ASTRA_ASSERT(ok); - + initCUDAAlgorithm(); m_bAlgoInit = true; } diff --git a/src/CudaSirtAlgorithm.cpp b/src/CudaSirtAlgorithm.cpp index 33e381a..7beb30e 100644 --- a/src/CudaSirtAlgorithm.cpp +++ b/src/CudaSirtAlgorithm.cpp @@ -113,36 +113,23 @@ bool CCudaSirtAlgorithm::initialize(CProjector2D* _pProjector, } //---------------------------------------------------------------------------------------- -// Iterate -void CCudaSirtAlgorithm::run(int _iNrIterations) -{ - // check initialized - ASTRA_ASSERT(m_bIsInitialized); - if (!m_bAlgoInit) { - // We only override the initialisation step to copy the min/max masks +void CCudaSirtAlgorithm::initCUDAAlgorithm() +{ + CCudaReconstructionAlgorithm2D::initCUDAAlgorithm(); - bool ok = setupGeometry(); - ASTRA_ASSERT(ok); + astraCUDA::SIRT* pSirt = dynamic_cast(m_pAlgo); - ok = m_pAlgo->allocateBuffers(); + if (m_pMinMask || m_pMaxMask) { + const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); + const float *pfMinMaskData = 0; + const float *pfMaxMaskData = 0; + if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); + if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); + bool ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); ASTRA_ASSERT(ok); - - if (m_pMinMask || m_pMaxMask) { - const CVolumeGeometry2D& volgeom = *m_pReconstruction->getGeometry(); - astraCUDA::SIRT* pSirt = dynamic_cast(m_pAlgo); - const float *pfMinMaskData = 0; - const float *pfMaxMaskData = 0; - if (m_pMinMask) pfMinMaskData = m_pMinMask->getDataConst(); - if (m_pMaxMask) pfMaxMaskData = m_pMaxMask->getDataConst(); - ok = pSirt->uploadMinMaxMasks(pfMinMaskData, pfMaxMaskData, volgeom.getGridColCount()); - ASTRA_ASSERT(ok); - } - - m_bAlgoInit = true; } - CCudaReconstructionAlgorithm2D::run(_iNrIterations); } -- cgit v1.2.3 From f03ceb16d2dbde0c43e8c90683c5feafe01e5356 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Wed, 23 Mar 2016 15:30:47 +0100 Subject: Rename ART lambda option to Relaxation --- include/astra/ArtAlgorithm.h | 4 ++-- src/ArtAlgorithm.cpp | 10 +++++++--- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/include/astra/ArtAlgorithm.h b/include/astra/ArtAlgorithm.h index d232b95..1ad9f3f 100644 --- a/include/astra/ArtAlgorithm.h +++ b/include/astra/ArtAlgorithm.h @@ -59,7 +59,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{Lamda, float, 1, The relaxation factor.} + * \astra_xml_item_option{Relaxation, float, 1, The relaxation factor.} * \astra_xml_item_option{RayOrder, string, "sequential", the order in which the rays are updated. 'sequential' or 'custom'} * \astra_xml_item_option{RayOrderList, n by 2 vector of float, not used, if RayOrder='custom': use this ray order. Each row consist of a projection id and detector id.} * @@ -73,7 +73,7 @@ namespace astra { * cfg.option.UseMinConstraint = 'yes';\n * cfg.option.UseMaxConstraint = 'yes';\n * cfg.option.MaxConstraintValue = 1024;\n - * cfg.option.Lamda = 0.7;\n + * cfg.option.Relaxation = 0.7;\n * cfg.option.RayOrder = 'custom';\n * cfg.option.RayOrderList = [0\,0; 0\,2; 1\,0];\n * alg_id = astra_mex_algorithm('create'\, cfg);\n diff --git a/src/ArtAlgorithm.cpp b/src/ArtAlgorithm.cpp index b59bd93..526c263 100644 --- a/src/ArtAlgorithm.cpp +++ b/src/ArtAlgorithm.cpp @@ -156,8 +156,12 @@ bool CArtAlgorithm::initialize(const Config& _cfg) return false; } + // "Lambda" is replaced by the more descriptive "Relaxation" m_fLambda = _cfg.self.getOptionNumerical("Lambda", 1.0f); - CC.markOptionParsed("Lambda"); + m_fLambda = _cfg.self.getOptionNumerical("Relaxation", m_fLambda); + if (!_cfg.self.hasOption("Relaxation")) + CC.markOptionParsed("Lambda"); + CC.markOptionParsed("Relaxation"); // success m_bIsInitialized = _check(); @@ -232,7 +236,7 @@ map CArtAlgorithm::getInformation() { map res; res["RayOrder"] = getInformation("RayOrder"); - res["Lambda"] = getInformation("Lambda"); + res["Relaxation"] = getInformation("Relaxation"); return mergeMap(CReconstructionAlgorithm2D::getInformation(), res); }; @@ -240,7 +244,7 @@ map CArtAlgorithm::getInformation() // Information - Specific boost::any CArtAlgorithm::getInformation(std::string _sIdentifier) { - if (_sIdentifier == "Lambda") { return m_fLambda; } + if (_sIdentifier == "Relaxation") { return m_fLambda; } if (_sIdentifier == "RayOrder") { vector res; for (int i = 0; i < m_iRayCount; i++) { -- cgit v1.2.3 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 From 16430239d04ff738a21146c410918c285552543f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn 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(D_pixelWeight, D_maskData, dims); } + processVol3D(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(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