diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h index d898061e3d..a9b779da18 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h @@ -1,101 +1,133 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ -/*=================================================================== - -This file is based heavily on a corresponding ITK filter. - -===================================================================*/ #ifndef __itkNonLocalMeansDenoisingFilter_h_ #define __itkNonLocalMeansDenoisingFilter_h_ #include "itkImageToImageFilter.h" #include "itkVectorImage.h" #include -#include +//#include namespace itk{ /** \class NonLocalMeansDenoisingFilter * \brief This class denoises a vectorimage according to the non local means procedure. * - * This Filter needs as an input the diffusion weigthed image and a related brainmask. - * Search- and comparisonradius need to be set! + * This Filter needs as an input a diffusion weigthed image. + * An input mask is optional to denoise only inside the mask range. All other voxels will be set to 0. */ template< class TPixelType > class NonLocalMeansDenoisingFilter : public ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > > { public: /** Typedefs */ typedef NonLocalMeansDenoisingFilter Self; typedef SmartPointer Pointer; typedef SmartPointer ConstPointer; typedef ImageToImageFilter< VectorImage < TPixelType, 3 >, VectorImage < TPixelType, 3 > > Superclass; typedef typename Superclass::InputImageType InputImageType; typedef typename Superclass::OutputImageType OutputImageType; typedef typename Superclass::OutputImageRegionType OutputImageRegionType; typedef Image MaskImageType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(NonLocalMeansDenoisingFilter, ImageToImageFilter) - /** Set/Get Macros */ + /** + * \brief Set flag to use joint information + */ itkSetMacro(UseJointInformation, bool) + /** + * \brief Set the searchradius + * + * The searchradius generates a neighborhood of size (2 * searchradius + 1)³. + * Default is 4. + */ itkSetMacro(SearchRadius, int) + /** + * \brief Set the comparisonradius + * + * The comparisonradius generates neighborhoods of size (2 * comparisonradius +1)³. + * Default is 1. + */ itkSetMacro(ComparisonRadius, int) + /** + * \brief Set the variance of the noise + * + * The variance of the noise needs to be estimated tu use this filter properly. + * Default is 1. + */ itkSetMacro(Variance, double) + /** + * \brief Set flag to use a rician adaption + * + * If this flag is true the filter uses a method which is optimized for Rician distributed noise. + */ itkSetMacro(UseRicianAdaption, bool) + /** + * \brief Get the amount of calculated Voxels + * + * Returns the number of calculated Voxels until yet, useful for the use of a progressbars. + */ itkGetMacro(CurrentVoxelCount, unsigned int) + /** \brief Set the input image. **/ void SetInputImage(const InputImageType* image); + /** + * \brief Set a denoising mask + * + * optional + * + * Set a mask to denoise only the masked area, all voxel outside this area will be set to 0. + */ void SetInputMask(MaskImageType* mask); protected: NonLocalMeansDenoisingFilter(); ~NonLocalMeansDenoisingFilter() {} - void PrintSelf(std::ostream& os, Indent indent) const; void BeforeThreadedGenerateData(); void ThreadedGenerateData( const OutputImageRegionType &outputRegionForThread, ThreadIdType); private: int m_SearchRadius; int m_ComparisonRadius; bool m_UseJointInformation; bool m_UseRicianAdaption; unsigned int m_CurrentVoxelCount; double m_Variance; typename MaskImageType::Pointer m_Mask; }; } #ifndef ITK_MANUAL_INSTANTIATION #include "itkNonLocalMeansDenoisingFilter.txx" #endif #endif //__itkNonLocalMeansDenoisingFilter_h_ diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx index b06ca89770..31269f5710 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx @@ -1,377 +1,370 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef __itkNonLocalMeansDenoisingFilter_txx #define __itkNonLocalMeansDenoisingFilter_txx #include #include #include #define _USE_MATH_DEFINES #include #include "itkImageRegionIterator.h" #include "itkNeighborhoodIterator.h" #include #include namespace itk { template< class TPixelType > NonLocalMeansDenoisingFilter< TPixelType > ::NonLocalMeansDenoisingFilter() : m_SearchRadius(4), m_ComparisonRadius(1), m_UseJointInformation(false), m_UseRicianAdaption(false), m_Variance(1), m_Mask(NULL) { this->SetNumberOfRequiredInputs( 1 ); } template< class TPixelType > void NonLocalMeansDenoisingFilter< TPixelType > ::BeforeThreadedGenerateData() { MITK_INFO << "SearchRadius: " << m_SearchRadius; MITK_INFO << "ComparisonRadius: " << m_ComparisonRadius; MITK_INFO << "Noisevariance: " << m_Variance; MITK_INFO << "Use Rician Adaption: " << std::boolalpha << m_UseRicianAdaption; MITK_INFO << "Use Joint Information: " << std::boolalpha << m_UseJointInformation; typename InputImageType::Pointer inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); if (m_Mask.IsNull()) { m_Mask = MaskImageType::New(); m_Mask->SetRegions(inputImagePointer->GetLargestPossibleRegion()); m_Mask->Allocate(); m_Mask->FillBuffer(1); } else { ImageRegionIterator< MaskImageType > mit(m_Mask, m_Mask->GetLargestPossibleRegion()); mit.GoToBegin(); typename MaskImageType::IndexType minIndex; typename MaskImageType::IndexType maxIndex; minIndex.Fill(10000); maxIndex.Fill(0); while (!mit.IsAtEnd()) { if (mit.Get()) { minIndex[0] = minIndex[0] < mit.GetIndex()[0] ? minIndex[0] : mit.GetIndex()[0]; minIndex[1] = minIndex[1] < mit.GetIndex()[1] ? minIndex[1] : mit.GetIndex()[1]; minIndex[2] = minIndex[2] < mit.GetIndex()[2] ? minIndex[2] : mit.GetIndex()[2]; maxIndex[0] = maxIndex[0] > mit.GetIndex()[0] ? maxIndex[0] : mit.GetIndex()[0]; maxIndex[1] = maxIndex[1] > mit.GetIndex()[1] ? maxIndex[1] : mit.GetIndex()[1]; maxIndex[2] = maxIndex[2] > mit.GetIndex()[2] ? maxIndex[2] : mit.GetIndex()[2]; } ++mit; } typename OutputImageType::SizeType size; size[0] = maxIndex[0] - minIndex[0]; size[1] = maxIndex[1] - minIndex[1]; size[2] = maxIndex[2] - minIndex[2]; typename OutputImageType::RegionType region (minIndex, size); typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); outputImage->SetRequestedRegion(region); } m_CurrentVoxelCount = 0; } template< class TPixelType > void NonLocalMeansDenoisingFilter< TPixelType > ::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, ThreadIdType ) { // initialize iterators typename OutputImageType::Pointer outputImage = static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); ImageRegionIterator< OutputImageType > oit(outputImage, outputRegionForThread); oit.GoToBegin(); // typename MaskImageType::Pointer maskImage = // static_cast< MaskImageType* >( this->ProcessObject::GetInput(1) ); ImageRegionIterator< MaskImageType > mit(m_Mask, outputRegionForThread); mit.GoToBegin(); typedef ImageRegionIteratorWithIndex InputIteratorType; typename InputImageType::Pointer inputImagePointer = NULL; inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); InputIteratorType git(inputImagePointer, outputRegionForThread ); git.GoToBegin(); // iterate over complete image region while( !git.IsAtEnd() ) { typename OutputImageType::PixelType outpix; outpix.SetSize (inputImagePointer->GetVectorLength()); if (mit.Get() != 0 && !this->GetAbortGenerateData()) { if(!m_UseJointInformation) { for (int i = 0; i < (int)inputImagePointer->GetVectorLength(); ++i) { double summw = 0; double sumj = 0; double w = 0; std::vector wj; std::vector p; typename InputIteratorType::IndexType index; index = git.GetIndex(); for (int x = index.GetElement(0) - m_SearchRadius; x <= index.GetElement(0) + m_SearchRadius; ++x) { for (int y = index.GetElement(1) - m_SearchRadius; y <= index.GetElement(1) + m_SearchRadius; ++y) { for (int z = index.GetElement(2) - m_SearchRadius; z <= index.GetElement(2) + m_SearchRadius; ++z) { typename InputIteratorType::IndexType indexV; indexV.SetElement(0, x); indexV.SetElement(1, y); indexV.SetElement(2, z); if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexV)) { TPixelType pixelJ = inputImagePointer->GetPixel(indexV)[i]; double sumk = 0; double size = 0; for (int xi = index.GetElement(0) - m_ComparisonRadius, xj = x - m_ComparisonRadius; xi <= index.GetElement(0) + m_ComparisonRadius; ++xi, ++xj) { for (int yi = index.GetElement(1) - m_ComparisonRadius, yj = y - m_ComparisonRadius; yi <= index.GetElement(1) + m_ComparisonRadius; ++yi, ++yj) { for (int zi = index.GetElement(2) - m_ComparisonRadius, zj = z - m_ComparisonRadius; zi <= index.GetElement(2) + m_ComparisonRadius; ++zi, ++zj) { typename InputIteratorType::IndexType indexI, indexJ; indexI.SetElement(0, xi); indexI.SetElement(1, yi); indexI.SetElement(2, zi); indexJ.SetElement(0, xj); indexJ.SetElement(1, yj); indexJ.SetElement(2, zj); // Compare neighborhoods ni & nj if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexI) && inputImagePointer->GetLargestPossibleRegion().IsInside(indexJ)) { int diff = inputImagePointer->GetPixel(indexI)[i] - inputImagePointer->GetPixel(indexJ)[i]; sumk += (double)(diff*diff); ++size; } } } } // weight all neighborhoods w = std::exp( - sumk / size / m_Variance); wj.push_back(w); if (m_UseRicianAdaption) { p.push_back((double)(pixelJ*pixelJ)); } else { p.push_back((double)(pixelJ)); } summw += w; } } } } for (unsigned int n = 0; n < wj.size(); ++n) { sumj += (wj[n]/summw) * p[n]; } if (m_UseRicianAdaption) { sumj -=2 * m_Variance; } if (sumj < 0) { sumj = 0; } TPixelType outval; if (m_UseRicianAdaption) { outval = std::floor(std::sqrt(sumj) + 0.5); } else { outval = std::floor(sumj + 0.5); } outpix.SetElement(i, outval); } } else { double Z = 0; itk::VariableLengthVector sumj; sumj.SetSize(inputImagePointer->GetVectorLength()); sumj.Fill(0.0); double w = 0; std::vector wj; std::vector > p; typename InputIteratorType::IndexType index; index = git.GetIndex(); for (int x = index.GetElement(0) - m_SearchRadius; x <= index.GetElement(0) + m_SearchRadius; ++x) { for (int y = index.GetElement(1) - m_SearchRadius; y <= index.GetElement(1) + m_SearchRadius; ++y) { for (int z = index.GetElement(2) - m_SearchRadius; z <= index.GetElement(2) + m_SearchRadius; ++z) { typename InputIteratorType::IndexType indexV; indexV.SetElement(0, x); indexV.SetElement(1, y); indexV.SetElement(2, z); if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexV)) { typename InputImageType::PixelType pixelJ = inputImagePointer->GetPixel(indexV); double sumk = 0; double size = 0; for (int xi = index.GetElement(0) - m_ComparisonRadius, xj = x - m_ComparisonRadius; xi <= index.GetElement(0) + m_ComparisonRadius; ++xi, ++xj) { for (int yi = index.GetElement(1) - m_ComparisonRadius, yj = y - m_ComparisonRadius; yi <= index.GetElement(1) + m_ComparisonRadius; ++yi, ++yj) { for (int zi = index.GetElement(2) - m_ComparisonRadius, zj = z - m_ComparisonRadius; zi <= index.GetElement(2) + m_ComparisonRadius; ++zi, ++zj) { typename InputIteratorType::IndexType indexI, indexJ; indexI.SetElement(0, xi); indexI.SetElement(1, yi); indexI.SetElement(2, zi); indexJ.SetElement(0, xj); indexJ.SetElement(1, yj); indexJ.SetElement(2, zj); // Compare neighborhoods ni & nj if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexI) && inputImagePointer->GetLargestPossibleRegion().IsInside(indexJ)) { typename InputImageType::PixelType diff = inputImagePointer->GetPixel(indexI) - inputImagePointer->GetPixel(indexJ); sumk += (double)(diff.GetSquaredNorm()); ++size; } } } } // weight all neighborhoods size *= inputImagePointer->GetVectorLength() + 1; w = std::exp( - (sumk / size) / m_Variance); wj.push_back(w); if (m_UseRicianAdaption) { itk::VariableLengthVector m; m.SetSize(inputImagePointer->GetVectorLength()); for (unsigned int i = 0; i < inputImagePointer->GetVectorLength(); ++i) { double sp = (double)(pixelJ.GetElement(i) * pixelJ.GetElement(i)); m.SetElement(i,sp); } p.push_back(m); } else ++size; { p.push_back(pixelJ); } Z += w; } } } } for (unsigned int n = 0; n < wj.size(); ++n) { sumj = sumj + ((wj[n]/Z) * p[n]); } if (m_UseRicianAdaption) { sumj = sumj - (2 * m_Variance); } for (unsigned int i = 0; i < inputImagePointer->GetVectorLength(); ++i) { double a = sumj.GetElement(i); if (a < 0) { a = 0; } TPixelType outval; if (m_UseRicianAdaption) { outval = std::floor(std::sqrt(a) + 0.5); } else { outval = std::floor(a + 0.5); } outpix.SetElement(i, outval); } } } else { outpix.Fill(0); } oit.Set(outpix); ++oit; ++m_CurrentVoxelCount; ++git; ++mit; } MITK_INFO << "One Thread finished calculation"; } -template< class TPixelType > -void NonLocalMeansDenoisingFilter< TPixelType >::PrintSelf(std::ostream& os, Indent indent) const -{ - Superclass::PrintSelf(os,indent); -} - template< class TPixelType > void NonLocalMeansDenoisingFilter< TPixelType >::SetInputImage(const InputImageType* image) { this->SetNthInput(0, const_cast< InputImageType* >(image)); } template< class TPixelType > void NonLocalMeansDenoisingFilter< TPixelType >::SetInputMask(MaskImageType* mask) { - //this->SetNthInput(1, const_cast< MaskImageType* >(mask)); m_Mask = mask; } } #endif // __itkNonLocalMeansDenoisingFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp index 5d83ff21ec..680b5c6295 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp @@ -1,124 +1,176 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkDiffusionImage.h" #include "mitkIOUtil.h" #include "mitkTestingMacros.h" #include "mitkTestFixture.h" #include "itkNonLocalMeansDenoisingFilter.h" class mitkNonLocalMeansDenoisingTestSuite : public mitk::TestFixture { CPPUNIT_TEST_SUITE(mitkNonLocalMeansDenoisingTestSuite); - MITK_TEST(Denoise_NLMr_1_1_shouldReturnTrue); - MITK_TEST(Denoise_NLMv_1_1_1_shouldReturnTrue); + MITK_TEST(Denoise_NLMg_shouldReturnTrue); + MITK_TEST(Denoise_NLMr_shouldReturnTrue); + MITK_TEST(Denoise_NLMv_shouldReturnTrue); + MITK_TEST(Denoise_NLMvr_shouldReturnTrue); CPPUNIT_TEST_SUITE_END(); private: /** Members used inside the different (sub-)tests. All members are initialized via setUp().*/ mitk::DiffusionImage::Pointer m_Image; mitk::DiffusionImage::Pointer m_ReferenceImage; mitk::DiffusionImage::Pointer m_DenoisedImage; itk::Image::Pointer m_ImageMask; itk::NonLocalMeansDenoisingFilter::Pointer m_DenoisingFilter; public: /** * @brief Setup Always call this method before each Test-case to ensure correct and new intialization of the used members for a new test case. (If the members are not used in a test, the method does not need to be called). */ void setUp() { //generate test images std::string imagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi.dwi"); - std::string maskPath = GetTestDataFilePath("DiffusionImaging/Denoising/denoising_mask.nrrd"); +// std::string maskPath = GetTestDataFilePath("DiffusionImaging/Denoising/denoising_mask.nrrd"); m_Image = static_cast*>( mitk::IOUtil::LoadImage(imagePath).GetPointer()); - mitk::Image::Pointer imageMask = static_cast( mitk::IOUtil::LoadImage(maskPath).GetPointer()); - mitk::CastToItkImage(imageMask, m_ImageMask); +// mitk::Image::Pointer imageMask = static_cast( mitk::IOUtil::LoadImage(maskPath).GetPointer()); +// mitk::CastToItkImage(imageMask, m_ImageMask); m_ReferenceImage = NULL; m_DenoisedImage = mitk::DiffusionImage::New(); //initialise Filter m_DenoisingFilter = itk::NonLocalMeansDenoisingFilter::New(); m_DenoisingFilter->SetInputImage(m_Image->GetVectorImage()); - m_DenoisingFilter->SetInputMask(m_ImageMask); +// m_DenoisingFilter->SetInputMask(m_ImageMask); m_DenoisingFilter->SetNumberOfThreads(1); m_DenoisingFilter->SetComparisonRadius(1); m_DenoisingFilter->SetSearchRadius(1); + m_DenoisingFilter->SetVariance(500); } void tearDown() { m_Image = NULL; m_ImageMask = NULL; m_ReferenceImage = NULL; m_DenoisingFilter = NULL; m_DenoisedImage = NULL; } - void Denoise_NLMr_1_1_shouldReturnTrue() + void Denoise_NLMg_shouldReturnTrue() { - std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMr_1-1.dwi"); + std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMg.dwi"); m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + m_DenoisingFilter->SetUseRicianAdaption(false); + m_DenoisingFilter->SetUseJointInformation(false); + try + { + m_DenoisingFilter->Update(); + } + catch(std::exception& e) + { + MITK_ERROR << e.what(); + } + + m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); + m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); + m_DenoisedImage->SetDirections(m_Image->GetDirections()); + m_DenoisedImage->InitializeFromVectorImage(); + + MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMg should always return the same result."); + } + + void Denoise_NLMr_shouldReturnTrue() + { + std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMr.dwi"); + m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + + m_DenoisingFilter->SetUseRicianAdaption(true); m_DenoisingFilter->SetUseJointInformation(false); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); m_DenoisedImage->SetDirections(m_Image->GetDirections()); m_DenoisedImage->InitializeFromVectorImage(); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMr should always return the same result."); } - void Denoise_NLMv_1_1_1_shouldReturnTrue() + void Denoise_NLMv_shouldReturnTrue() { - std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMv_1-1-1.dwi"); + std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMv.dwi"); m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); -// m_DenoisingFilter->SetChannelRadius(1); + m_DenoisingFilter->SetUseRicianAdaption(false); m_DenoisingFilter->SetUseJointInformation(true); try { m_DenoisingFilter->Update(); } catch(std::exception& e) { MITK_ERROR << e.what(); } m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); m_DenoisedImage->SetDirections(m_Image->GetDirections()); m_DenoisedImage->InitializeFromVectorImage(); MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMv should always return the same result."); } + void Denoise_NLMvr_shouldReturnTrue() + { + std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMvr.dwi"); + m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); + + m_DenoisingFilter->SetUseRicianAdaption(true); + m_DenoisingFilter->SetUseJointInformation(true); + try + { + m_DenoisingFilter->Update(); + } + catch(std::exception& e) + { + MITK_ERROR << e.what(); + } + + m_DenoisedImage->SetVectorImage(m_DenoisingFilter->GetOutput()); + m_DenoisedImage->SetReferenceBValue(m_Image->GetReferenceBValue()); + m_DenoisedImage->SetDirections(m_Image->GetDirections()); + m_DenoisedImage->InitializeFromVectorImage(); + + MITK_ASSERT_EQUAL( m_DenoisedImage, m_ReferenceImage, "NLMvr should always return the same result."); + } + }; MITK_TEST_SUITE_REGISTRATION(mitkNonLocalMeansDenoising) diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp index 49ce33de31..aab1587ca0 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp @@ -1,463 +1,479 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // Blueberry #include #include // Qmitk #include "QmitkDenoisingView.h" #include #include QmitkDenoisingWorker::QmitkDenoisingWorker(QmitkDenoisingView *view) : m_View(view) { } void QmitkDenoisingWorker::run() { if (m_View->m_ImageNode.IsNotNull()) { switch (m_View->m_SelectedFilter) { case QmitkDenoisingView::NOFILTERSELECTED: case QmitkDenoisingView::GAUSS: { break; } case QmitkDenoisingView::NLM: { try { m_View->m_CompletedCalculation = true; m_View->m_NonLocalMeansFilter->Update(); } catch (itk::ExceptionObject& e) { m_View->m_CompletedCalculation = false; MITK_ERROR << e.what(); } break; } } m_View->m_DenoisingThread.quit(); } } const std::string QmitkDenoisingView::VIEW_ID = "org.mitk.views.denoisingview"; QmitkDenoisingView::QmitkDenoisingView() : QmitkFunctionality() , m_Controls( 0 ) , m_ImageNode(NULL) , m_BrainMaskNode(NULL) , m_DenoisingWorker(this) , m_ThreadIsRunning(false) , m_NonLocalMeansFilter(NULL) , m_InputImage(NULL) , m_LastProgressCount(0) , m_MaxProgressCount(0) , m_SelectedFilter(NOFILTERSELECTED) { m_DenoisingWorker.moveToThread(&m_DenoisingThread); connect(&m_DenoisingThread, SIGNAL(started()), this, SLOT(BeforeThread())); connect(&m_DenoisingThread, SIGNAL(started()), &m_DenoisingWorker, SLOT(run())); connect(&m_DenoisingThread, SIGNAL(finished()), this, SLOT(AfterThread())); connect(&m_DenoisingThread, SIGNAL(terminated()), this, SLOT(AfterThread())); m_DenoisingTimer = new QTimer(this); } QmitkDenoisingView::~QmitkDenoisingView() { delete m_DenoisingTimer; } void QmitkDenoisingView::CreateQtPartControl( QWidget *parent ) { // build up qt view, unless already done if ( !m_Controls ) { // create GUI widgets from the Qt Designer's .ui file m_Controls = new Ui::QmitkDenoisingViewControls; m_Controls->setupUi( parent ); CreateConnections(); ResetParameterPanel(); } } void QmitkDenoisingView::CreateConnections() { if ( m_Controls ) { connect( (QObject*)(m_Controls->m_ApplyButton), SIGNAL(clicked()), this, SLOT(StartDenoising())); connect( (QObject*)(m_Controls->m_SelectFilterComboBox), SIGNAL(activated(int)), this, SLOT(SelectFilter(int))); connect( m_DenoisingTimer, SIGNAL(timeout()), this, SLOT(UpdateProgress())); } } void QmitkDenoisingView::Activated() { QmitkFunctionality::Activated(); m_Controls->m_SelectFilterComboBox->clear(); m_Controls->m_SelectFilterComboBox->insertItem(NOFILTERSELECTED, QString( QApplication::translate("QmitkDenoisingView", "Please select a filter", 0, QApplication::UnicodeUTF8) )); m_Controls->m_SelectFilterComboBox->insertItem(NLM, QString( QApplication::translate("QmitkDenoisingView", "Non-local means filter", 0, QApplication::UnicodeUTF8) )); m_Controls->m_SelectFilterComboBox->insertItem(GAUSS, QString( QApplication::translate("QmitkDenoisingView", "Discrete gaussian filter", 0, QApplication::UnicodeUTF8) )); } void QmitkDenoisingView::OnSelectionChanged( std::vector nodes ) { if (m_ThreadIsRunning) return; if (m_SelectedFilter != NOFILTERSELECTED) { m_Controls->m_InputImageLabel->setText("mandatory"); } else { m_Controls->m_InputImageLabel->setText("mandatory"); } m_Controls->m_InputBrainMaskLabel->setText("optional"); m_Controls->m_ApplyButton->setEnabled(false); m_ImageNode = NULL; m_BrainMaskNode = NULL; // iterate selection for( std::vector::iterator it = nodes.begin(); it != nodes.end(); ++it ) { mitk::DataNode::Pointer node = *it; if( node.IsNotNull() && dynamic_cast(node->GetData())) { m_Controls->m_InputImageLabel->setText(node->GetName().c_str()); m_ImageNode = node; } bool isBinary = false; node->GetBoolProperty("binary", isBinary); // look for a brainmask in selection if( node.IsNotNull() && static_cast(node->GetData()) && isBinary) { m_Controls->m_InputBrainMaskLabel->setText(node->GetName().c_str()); m_BrainMaskNode = node; } } // Preparation of GUI to start denoising if a filter is selected if (m_ImageNode.IsNotNull() && m_SelectedFilter != NOFILTERSELECTED) { m_Controls->m_ApplyButton->setEnabled(true); } } void QmitkDenoisingView::StartDenoising() { if (!m_ThreadIsRunning) { if (m_ImageNode.IsNotNull()) { m_LastProgressCount = 0; switch (m_SelectedFilter) { case NOFILTERSELECTED: { break; } case NLM: { // initialize NLM m_InputImage = dynamic_cast (m_ImageNode->GetData()); m_NonLocalMeansFilter = NonLocalMeansDenoisingFilterType::New(); if (m_BrainMaskNode.IsNotNull()) { // use brainmask if set m_ImageMask = dynamic_cast(m_BrainMaskNode->GetData()); itk::Image::Pointer itkMask = itk::Image::New(); mitk::CastToItkImage(m_ImageMask, itkMask); m_NonLocalMeansFilter->SetInputMask(itkMask); } m_NonLocalMeansFilter->SetNumberOfThreads(12); m_NonLocalMeansFilter->SetInputImage(m_InputImage->GetVectorImage()); m_NonLocalMeansFilter->SetUseRicianAdaption(m_Controls->m_RicianCheckbox->isChecked()); m_NonLocalMeansFilter->SetUseJointInformation(m_Controls->m_JointInformationCheckbox->isChecked()); m_NonLocalMeansFilter->SetSearchRadius(m_Controls->m_SpinBoxParameter1->value()); m_NonLocalMeansFilter->SetComparisonRadius(m_Controls->m_SpinBoxParameter2->value()); + m_NonLocalMeansFilter->SetVariance(m_Controls->m_DoubleSpinBoxParameter3->value()); // initialize the progressbar m_MaxProgressCount = m_InputImage->GetDimension(0) * m_InputImage->GetDimension(1) * m_InputImage->GetDimension(2); mitk::ProgressBar::GetInstance()->AddStepsToDo(m_MaxProgressCount); // start denoising in detached thread m_DenoisingThread.start(QThread::HighestPriority); break; } case GAUSS: { // initialize GAUSS and run m_InputImage = dynamic_cast (m_ImageNode->GetData()); ExtractFilterType::Pointer extractor = ExtractFilterType::New(); extractor->SetInput(m_InputImage->GetVectorImage()); ComposeFilterType::Pointer composer = ComposeFilterType::New(); for (unsigned int i = 0; i < m_InputImage->GetVectorImage()->GetVectorLength(); ++i) { extractor->SetIndex(i); extractor->Update(); m_GaussianFilter = GaussianFilterType::New(); m_GaussianFilter->SetVariance(m_Controls->m_SpinBoxParameter1->value()); if (m_BrainMaskNode.IsNotNull()) { m_ImageMask = dynamic_cast(m_BrainMaskNode->GetData()); itk::Image::Pointer itkMask = itk::Image::New(); mitk::CastToItkImage(m_ImageMask, itkMask); itk::MaskImageFilter , itk::Image >::Pointer maskImageFilter = itk::MaskImageFilter , itk::Image >::New(); maskImageFilter->SetInput(extractor->GetOutput()); maskImageFilter->SetMaskImage(itkMask); maskImageFilter->Update(); m_GaussianFilter->SetInput(maskImageFilter->GetOutput()); } else { m_GaussianFilter->SetInput(extractor->GetOutput()); } m_GaussianFilter->Update(); composer->SetInput(i, m_GaussianFilter->GetOutput()); } composer->Update(); DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage(composer->GetOutput()); image->SetReferenceBValue(m_InputImage->GetReferenceBValue()); image->SetDirections(m_InputImage->GetDirections()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_ImageNode->GetName().c_str(); imageNode->SetName((name+"_gauss_"+QString::number(m_Controls->m_SpinBoxParameter1->value())).toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); break; } } } } else { m_NonLocalMeansFilter->AbortGenerateDataOn(); m_CompletedCalculation = false; } } void QmitkDenoisingView::ResetParameterPanel() { m_Controls->m_DwiLabel->setEnabled(false); m_Controls->m_InputImageLabel->setEnabled(false); m_Controls->m_BrainMaskLabel->setEnabled(false); m_Controls->m_InputBrainMaskLabel->setEnabled(false); m_Controls->m_ParameterBox->hide(); m_Controls->m_LabelParameter_1->hide(); m_Controls->m_LabelParameter_2->hide(); m_Controls->m_LabelParameter_3->hide(); m_Controls->m_SpinBoxParameter1->hide(); m_Controls->m_SpinBoxParameter2->hide(); m_Controls->m_DoubleSpinBoxParameter3->hide(); m_Controls->m_RicianLabel->hide(); m_Controls->m_RicianCheckbox->hide(); m_Controls->m_JointInformationLabel->hide(); m_Controls->m_JointInformationCheckbox->hide(); m_Controls->m_ApplyButton->setEnabled(false); } void QmitkDenoisingView::SelectFilter(int filter) { if (m_ThreadIsRunning) return; //Prepare GUI this->ResetParameterPanel(); switch (filter) { case 0: { m_SelectedFilter = NOFILTERSELECTED; break; } case 1: { m_SelectedFilter = NLM; m_Controls->m_DwiLabel->setEnabled(true); m_Controls->m_InputImageLabel->setEnabled(true); m_Controls->m_BrainMaskLabel->setEnabled(true); m_Controls->m_InputBrainMaskLabel->setEnabled(true); m_Controls->m_ParameterBox->show(); m_Controls->m_LabelParameter_1->show(); m_Controls->m_LabelParameter_1->setText("Search Radius:"); m_Controls->m_LabelParameter_2->show(); m_Controls->m_LabelParameter_2->setText("Comparision Radius:"); m_Controls->m_LabelParameter_3->show(); m_Controls->m_LabelParameter_3->setText("Noise variance:"); m_Controls->m_SpinBoxParameter1->show(); m_Controls->m_SpinBoxParameter1->setValue(4); m_Controls->m_SpinBoxParameter2->show(); m_Controls->m_SpinBoxParameter2->setValue(1); m_Controls->m_DoubleSpinBoxParameter3->show(); m_Controls->m_DoubleSpinBoxParameter3->setValue(1.0); m_Controls->m_RicianLabel->show(); m_Controls->m_RicianCheckbox->show(); m_Controls->m_RicianCheckbox->setChecked(true); m_Controls->m_JointInformationLabel->show(); m_Controls->m_JointInformationCheckbox->show(); m_Controls->m_JointInformationCheckbox->setChecked(false); break; } case 2: { m_SelectedFilter = GAUSS; m_Controls->m_DwiLabel->setEnabled(true); m_Controls->m_InputImageLabel->setEnabled(true); m_Controls->m_BrainMaskLabel->setEnabled(true); m_Controls->m_InputBrainMaskLabel->setEnabled(true); m_Controls->m_ParameterBox->show(); m_Controls->m_LabelParameter_1->show(); m_Controls->m_LabelParameter_1->setText("Variance:"); m_Controls->m_SpinBoxParameter1->show(); m_Controls->m_SpinBoxParameter1->setValue(2); break; } } if (m_ImageNode.IsNull()) { if (m_SelectedFilter != NOFILTERSELECTED) m_Controls->m_InputImageLabel->setText("mandatory"); else m_Controls->m_InputImageLabel->setText("mandatory"); } if (m_ImageNode.IsNotNull()) { m_Controls->m_ApplyButton->setEnabled(false); switch(filter) { case NOFILTERSELECTED: { break; } case NLM: case GAUSS: { m_Controls->m_ApplyButton->setEnabled(true); break; } } } } void QmitkDenoisingView::BeforeThread() { m_ThreadIsRunning = true; // initialize timer to update the progressbar at each timestep m_DenoisingTimer->start(500); m_Controls->m_ParameterBox->setEnabled(false); m_Controls->m_ApplyButton->setText("Abort"); } void QmitkDenoisingView::AfterThread() { m_ThreadIsRunning = false; // stop timer to stop updates of progressbar m_DenoisingTimer->stop(); // make sure progressbar is finished mitk::ProgressBar::GetInstance()->Progress(m_MaxProgressCount); if (m_CompletedCalculation) { switch (m_SelectedFilter) { case NOFILTERSELECTED: case GAUSS: { break; } case NLM: { DiffusionImageType::Pointer image = DiffusionImageType::New(); image->SetVectorImage(m_NonLocalMeansFilter->GetOutput()); image->SetReferenceBValue(m_InputImage->GetReferenceBValue()); image->SetDirections(m_InputImage->GetDirections()); image->InitializeFromVectorImage(); mitk::DataNode::Pointer imageNode = mitk::DataNode::New(); imageNode->SetData( image ); QString name = m_ImageNode->GetName().c_str(); //TODO: Rician adaption & joint information in name - imageNode->SetName((name+"_NLM_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); + if (m_Controls->m_RicianCheckbox->isChecked() && !m_Controls->m_JointInformationCheckbox->isChecked()) + { + imageNode->SetName((name+"_NLMr_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); + } + else if(!m_Controls->m_RicianCheckbox->isChecked() && m_Controls->m_JointInformationCheckbox->isChecked()) + { + imageNode->SetName((name+"_NLMv_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); + } + else if(m_Controls->m_RicianCheckbox->isChecked() && m_Controls->m_JointInformationCheckbox->isChecked()) + { + imageNode->SetName((name+"_NLMvr_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); + } + else + { + imageNode->SetName((name+"_NLM_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); + } GetDefaultDataStorage()->Add(imageNode); break; } } } m_Controls->m_ParameterBox->setEnabled(true); m_Controls->m_ApplyButton->setText("Apply"); } void QmitkDenoisingView::UpdateProgress() { switch (m_SelectedFilter) { case NOFILTERSELECTED: case GAUSS: { break; } case NLM: { unsigned int currentProgressCount = m_NonLocalMeansFilter->GetCurrentVoxelCount(); mitk::ProgressBar::GetInstance()->Progress(currentProgressCount-m_LastProgressCount); m_LastProgressCount = currentProgressCount; break; } } } diff --git a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.h b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.h index 55bf451baf..ba63e1a8b3 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.h +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.h @@ -1,131 +1,133 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef _QMITKQmitkDenoisingView_H_INCLUDED #define _QMITKQmitkDenoisingView_H_INCLUDED #include #include #include "ui_QmitkDenoisingViewControls.h" #include #include #include #include #include #include #include #include #include #include +/** + * \class QmitkDenoisingView + * \brief View displaying details to denoise diffusionweighted images. + * + * \sa QmitkFunctionality + * \ingroup Functionalities + */ class QmitkDenoisingView; class QmitkDenoisingWorker : public QObject { Q_OBJECT public: QmitkDenoisingWorker(QmitkDenoisingView* view); public slots: void run(); private: QmitkDenoisingView* m_View; }; -/*! - \brief View displaying details to denoise diffusionweighted images. - \sa QmitkFunctionality - \ingroup Functionalities -*/ class QmitkDenoisingView : public QmitkFunctionality { // this is needed for all Qt objects that should have a Qt meta-object // (everything that derives from QObject and wants to have signal/slots) Q_OBJECT public: static const std::string VIEW_ID; QmitkDenoisingView(); virtual ~QmitkDenoisingView(); /** Typedefs */ typedef short DiffusionPixelType; typedef mitk::DiffusionImage< DiffusionPixelType > DiffusionImageType; typedef mitk::Image MaskImageType; typedef itk::NonLocalMeansDenoisingFilter< DiffusionPixelType > NonLocalMeansDenoisingFilterType; typedef itk::DiscreteGaussianImageFilter < itk::Image< DiffusionPixelType, 3>, itk::Image< DiffusionPixelType, 3> > GaussianFilterType; typedef itk::VectorImageToImageFilter < DiffusionPixelType > ExtractFilterType; typedef itk::ComposeImageFilter < itk::Image > ComposeFilterType; virtual void CreateQtPartControl(QWidget *parent); /// \brief Creation of the connections of main and control widget virtual void CreateConnections(); /// \brief Creation of the connections of the FilterComboBox virtual void Activated(); private slots: - void StartDenoising(); ///< prepares filter condition and starts thread for denoising - void SelectFilter(int filter); ///< updates which filter is selected - void BeforeThread(); ///< starts timer & disables all buttons while denoising - void AfterThread(); ///< stops timer & creates a new datanode of the denoised image - void UpdateProgress(); ///< updates the progressbar each timestep + void StartDenoising(); //< prepares filter condition and starts thread for denoising + void SelectFilter(int filter); //< updates which filter is selected + void BeforeThread(); //< starts timer & disables all buttons while denoising + void AfterThread(); //< stops timer & creates a new datanode of the denoised image + void UpdateProgress(); //< updates the progressbar each timestep private: /// \brief called by QmitkFunctionality when DataManager's selection has changed virtual void OnSelectionChanged( std::vector nodes ); void ResetParameterPanel(); Ui::QmitkDenoisingViewControls* m_Controls; mitk::DataNode::Pointer m_ImageNode; mitk::DataNode::Pointer m_BrainMaskNode; QmitkDenoisingWorker m_DenoisingWorker; QThread m_DenoisingThread; bool m_ThreadIsRunning; bool m_CompletedCalculation; NonLocalMeansDenoisingFilterType::Pointer m_NonLocalMeansFilter; GaussianFilterType::Pointer m_GaussianFilter; DiffusionImageType::Pointer m_InputImage; MaskImageType::Pointer m_ImageMask; QTimer* m_DenoisingTimer; unsigned int m_LastProgressCount; unsigned int m_MaxProgressCount; enum FilterType { NOFILTERSELECTED, NLM, GAUSS }m_SelectedFilter; friend class QmitkDenoisingWorker; }; #endif // _QmitkDenoisingView_H_INCLUDED