diff --git a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h index 4b05697352..1e6aa1c007 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.h @@ -1,110 +1,114 @@ /*=================================================================== 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 #include #include #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! */ 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; typedef VectorImageToImageFilter < TPixelType > ImageExtractorType; typedef ChangeInformationImageFilter < MaskImageType > ChangeInformationType; typedef ExtractImageFilter < MaskImageType, MaskImageType > ExtractImageFilterType; typedef LabelStatisticsImageFilter < MaskImageType, MaskImageType > LabelStatisticsFilterType; typedef InvertIntensityImageFilter < MaskImageType, MaskImageType > InvertImageFilterType; typedef StatisticsImageFilter < MaskImageType > StatisticsFilterType; /** Method for creation through the object factory. */ itkFactorylessNewMacro(Self) itkCloneMacro(Self) /** Runtime information support. */ itkTypeMacro(NonLocalMeansDenoisingFilter, ImageToImageFilter) /** Set/Get Macros */ itkSetMacro(UseJointInformation, bool) - itkSetMacro(ChannelRadius, int) itkSetMacro(SearchRadius, int) itkSetMacro(ComparisonRadius, int) + itkSetMacro(Variance, double) + itkSetMacro(UseRicianAdaption, bool) itkGetMacro(CurrentVoxelCount, unsigned int) + void SetInputImage(const InputImageType* image); void SetInputMask(const 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; - int m_ChannelRadius; - VariableLengthVector< double > m_Deviations; +// int m_ChannelRadius; +// VariableLengthVector< double > m_Deviations; bool m_UseJointInformation; + bool m_UseRicianAdaption; unsigned int m_CurrentVoxelCount; + double m_Variance; }; } #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 39147e4b89..2bc2f7206e 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx +++ b/Modules/DiffusionImaging/DiffusionCore/Algorithms/itkNonLocalMeansDenoisingFilter.txx @@ -1,289 +1,392 @@ /*=================================================================== 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(5), + m_ComparisonRadius(1), + m_UseJointInformation(true), + m_UseRicianAdaption(false), + m_Variance(536.87) { this->SetNumberOfRequiredInputs( 2 ); } template< class TPixelType > void NonLocalMeansDenoisingFilter< TPixelType > ::BeforeThreadedGenerateData() { - typename OutputImageType::Pointer outputImage = - static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); - typename OutputImageType::PixelType px; - px.SetSize(1); - px.SetElement(0,0); - outputImage->FillBuffer(px); - - typename InputImageType::Pointer inImage = static_cast< InputImageType* >(this->ProcessObject::GetInput(0)); - typename MaskImageType::Pointer mask = static_cast< MaskImageType* >(this->ProcessObject::GetInput(1)); - int size = inImage->GetVectorLength(); - m_Deviations.SetSize(size); - typename ImageExtractorType::Pointer extractor = ImageExtractorType::New(); - extractor->SetInput(inImage); - - // calculate max value of mask, for correct inversion - typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); - statisticsFilter->SetInput(mask); - statisticsFilter->Update(); - - // invert mask, to mask the backround - typename InvertImageFilterType::Pointer inverter = InvertImageFilterType::New(); - inverter->SetInput(mask); - inverter->SetMaximum(statisticsFilter->GetMaximum()); - inverter->Update(); - - // make sure inverted mask has same origin is the brainmask - typename ChangeInformationType::Pointer changeMaskFilter = ChangeInformationType::New(); - changeMaskFilter->ChangeOriginOn(); - changeMaskFilter->SetInput(inverter->GetOutput()); - changeMaskFilter->SetOutputOrigin(mask->GetOrigin()); - changeMaskFilter->Update(); - typename MaskImageType::Pointer invertedMask = changeMaskFilter->GetOutput(); - typename MaskImageType::PointType imageOrigin = inImage->GetOrigin(); - typename MaskImageType::PointType maskOrigin = invertedMask->GetOrigin(); - long offset[3]; - - typedef itk::ContinuousIndex ContinousIndexType; - ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; - - inImage->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); - inImage->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); - - // make sure there is no misalignment between mask and image - for ( unsigned int i = 0; i < 3; ++i ) - { - double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); - if ( fabs( misalignment ) > mitk::eps ) - { - itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); - } - - double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; - offset[i] = (int) indexCoordDistance + inImage->GetBufferedRegion().GetIndex()[i]; - } - - // calculate for each channel the stddev - for ( int i = 0; i < size; ++i) - { - /// extract channel i of the input - extractor->SetIndex(i); - extractor->Update(); - - // adapt mask to the image - typename ChangeInformationType::Pointer adaptMaskFilter; - adaptMaskFilter = ChangeInformationType::New(); - adaptMaskFilter->ChangeOriginOn(); - adaptMaskFilter->ChangeRegionOn(); - adaptMaskFilter->SetInput( invertedMask ); - adaptMaskFilter->SetOutputOrigin( extractor->GetOutput()->GetOrigin() /*image->GetOrigin()*/ ); - adaptMaskFilter->SetOutputOffset( offset ); - adaptMaskFilter->Update(); - - // extract backround as the ROI - typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); - typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); - extractImageFilter->SetInput( extractor->GetOutput() ); - extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); - extractImageFilter->Update(); - - // calculate statistics of ROI - typename MaskImageType::Pointer adaptedImage = extractImageFilter->GetOutput(); - typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); - labelStatisticsFilter->SetInput(adaptedImage); - labelStatisticsFilter->SetLabelInput(adaptedMaskImage); - labelStatisticsFilter->UseHistogramsOff(); - labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); - labelStatisticsFilter->Update(); - - // save the stddev of each channel - m_Deviations.SetElement(i, labelStatisticsFilter->GetSigma(1)); - } +// typename OutputImageType::Pointer outputImage = +// static_cast< OutputImageType * >(this->ProcessObject::GetOutput(0)); +// typename OutputImageType::PixelType px; +// px.SetSize(1); +// px.SetElement(0,0); +// outputImage->FillBuffer(px); + +// typename InputImageType::Pointer inImage = static_cast< InputImageType* >(this->ProcessObject::GetInput(0)); +// typename MaskImageType::Pointer mask = static_cast< MaskImageType* >(this->ProcessObject::GetInput(1)); +// int size = inImage->GetVectorLength(); +// m_Deviations.SetSize(size); +// typename ImageExtractorType::Pointer extractor = ImageExtractorType::New(); +// extractor->SetInput(inImage); + +// // calculate max value of mask, for correct inversion +// typename StatisticsFilterType::Pointer statisticsFilter = StatisticsFilterType::New(); +// statisticsFilter->SetInput(mask); +// statisticsFilter->Update(); + +// // invert mask, to mask the backround +// typename InvertImageFilterType::Pointer inverter = InvertImageFilterType::New(); +// inverter->SetInput(mask); +// inverter->SetMaximum(statisticsFilter->GetMaximum()); +// inverter->Update(); + +// // make sure inverted mask has same origin is the brainmask +// typename ChangeInformationType::Pointer changeMaskFilter = ChangeInformationType::New(); +// changeMaskFilter->ChangeOriginOn(); +// changeMaskFilter->SetInput(inverter->GetOutput()); +// changeMaskFilter->SetOutputOrigin(mask->GetOrigin()); +// changeMaskFilter->Update(); +// typename MaskImageType::Pointer invertedMask = changeMaskFilter->GetOutput(); +// typename MaskImageType::PointType imageOrigin = inImage->GetOrigin(); +// typename MaskImageType::PointType maskOrigin = invertedMask->GetOrigin(); +// long offset[3]; + +// typedef itk::ContinuousIndex ContinousIndexType; +// ContinousIndexType maskOriginContinousIndex, imageOriginContinousIndex; + +// inImage->TransformPhysicalPointToContinuousIndex(maskOrigin, maskOriginContinousIndex); +// inImage->TransformPhysicalPointToContinuousIndex(imageOrigin, imageOriginContinousIndex); + +// // make sure there is no misalignment between mask and image +// for ( unsigned int i = 0; i < 3; ++i ) +// { +// double misalignment = maskOriginContinousIndex[i] - floor( maskOriginContinousIndex[i] + 0.5 ); +// if ( fabs( misalignment ) > mitk::eps ) +// { +// itkExceptionMacro( << "Pixels/voxels of mask and image are not sufficiently aligned! (Misalignment: " << misalignment << ")" ); +// } + +// double indexCoordDistance = maskOriginContinousIndex[i] - imageOriginContinousIndex[i]; +// offset[i] = (int) indexCoordDistance + inImage->GetBufferedRegion().GetIndex()[i]; +// } + +// // calculate for each channel the stddev +// for ( int i = 0; i < size; ++i) +// { +// /// extract channel i of the input +// extractor->SetIndex(i); +// extractor->Update(); + +// // adapt mask to the image +// typename ChangeInformationType::Pointer adaptMaskFilter; +// adaptMaskFilter = ChangeInformationType::New(); +// adaptMaskFilter->ChangeOriginOn(); +// adaptMaskFilter->ChangeRegionOn(); +// adaptMaskFilter->SetInput( invertedMask ); +// adaptMaskFilter->SetOutputOrigin( extractor->GetOutput()->GetOrigin() /*image->GetOrigin()*/ ); +// adaptMaskFilter->SetOutputOffset( offset ); +// adaptMaskFilter->Update(); + +// // extract backround as the ROI +// typename MaskImageType::Pointer adaptedMaskImage = adaptMaskFilter->GetOutput(); +// typename ExtractImageFilterType::Pointer extractImageFilter = ExtractImageFilterType::New(); +// extractImageFilter->SetInput( extractor->GetOutput() ); +// extractImageFilter->SetExtractionRegion( adaptedMaskImage->GetBufferedRegion() ); +// extractImageFilter->Update(); + +// // calculate statistics of ROI +// typename MaskImageType::Pointer adaptedImage = extractImageFilter->GetOutput(); +// typename LabelStatisticsFilterType::Pointer labelStatisticsFilter = LabelStatisticsFilterType::New(); +// labelStatisticsFilter->SetInput(adaptedImage); +// labelStatisticsFilter->SetLabelInput(adaptedMaskImage); +// labelStatisticsFilter->UseHistogramsOff(); +// labelStatisticsFilter->GetOutput()->SetRequestedRegion( adaptedMaskImage->GetLargestPossibleRegion() ); +// labelStatisticsFilter->Update(); + +// // save the stddev of each channel +// m_Deviations.SetElement(i, labelStatisticsFilter->GetSigma(1)); +// } 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(); typedef ImageRegionIteratorWithIndex InputIteratorType; typename InputImageType::Pointer inputImagePointer = NULL; inputImagePointer = static_cast< InputImageType * >( this->ProcessObject::GetInput(0) ); InputIteratorType git(inputImagePointer, outputRegionForThread ); - InputIteratorType njit(inputImagePointer, outputRegionForThread ); - InputIteratorType niit(inputImagePointer, outputRegionForThread ); - InputIteratorType hit(inputImagePointer, outputRegionForThread); git.GoToBegin(); + m_Variance = 536.87; + // iterate over complete image region while( !git.IsAtEnd() ) { typename OutputImageType::PixelType outpix; outpix.SetSize (inputImagePointer->GetVectorLength()); - for (int i = 0; i < (int)inputImagePointer->GetVectorLength(); ++i) + if(!m_UseJointInformation) + { + for (int i = 0; i < (int)inputImagePointer->GetVectorLength(); ++i) + { + double Z = 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)); + } + Z += w; + } + } + } + } + for (unsigned int n = 0; n < wj.size(); ++n) + { + sumj += (wj[n]/Z) * p[n]; + } + double df; + if (m_UseRicianAdaption) + { + df = sumj - (2 * m_Variance); + } + else + { + df = sumj; + } + if (df < 0) + { + df = 0; + } + + TPixelType outval = std::floor(std::sqrt(df) + 0.5); + outpix.SetElement(i, outval); + } + } + + else { double Z = 0; - double sumj = 0; + itk::VariableLengthVector sumj; + sumj.SetSize(inputImagePointer->GetVectorLength()); + sumj.Fill(0.0); double w = 0; - double variance = m_Deviations.GetElement(i) * m_Deviations.GetElement(i); std::vector wj; - std::vector p; + std::vector > p; + typename InputIteratorType::IndexType index; + index = git.GetIndex(); - for (int x = git.GetIndex().GetElement(0) - m_SearchRadius; x <= git.GetIndex().GetElement(0) + m_SearchRadius; ++x) + for (int x = index.GetElement(0) - m_SearchRadius; x <= index.GetElement(0) + m_SearchRadius; ++x) { - for (int y = git.GetIndex().GetElement(1) - m_SearchRadius; y <= git.GetIndex().GetElement(1) + m_SearchRadius; ++y) + for (int y = index.GetElement(1) - m_SearchRadius; y <= index.GetElement(1) + m_SearchRadius; ++y) { - for (int z = git.GetIndex().GetElement(2) - m_SearchRadius; z <= git.GetIndex().GetElement(2) + m_SearchRadius; ++z) + for (int z = index.GetElement(2) - m_SearchRadius; z <= index.GetElement(2) + m_SearchRadius; ++z) { - typename InputIteratorType::IndexType idx; - idx.SetElement(0, x); - idx.SetElement(1, y); - idx.SetElement(2, z); - if (inputImagePointer->GetLargestPossibleRegion().IsInside(idx)) + typename InputIteratorType::IndexType indexV; + indexV.SetElement(0, x); + indexV.SetElement(1, y); + indexV.SetElement(2, z); + if (inputImagePointer->GetLargestPossibleRegion().IsInside(indexV)) { - hit.SetIndex(idx); - TPixelType pixelJ = hit.Get()[i]; + typename InputImageType::PixelType pixelJ = inputImagePointer->GetPixel(indexV); double sumk = 0; double size = 0; - for (int xi = git.GetIndex().GetElement(0) - m_ComparisonRadius, xj = hit.GetIndex().GetElement(0) - m_ComparisonRadius; xi <= git.GetIndex().GetElement(0) + m_ComparisonRadius; ++xi, ++xj) + for (int xi = index.GetElement(0) - m_ComparisonRadius, xj = x - m_ComparisonRadius; xi <= index.GetElement(0) + m_ComparisonRadius; ++xi, ++xj) { - for (int yi = git.GetIndex().GetElement(1) - m_ComparisonRadius, yj = hit.GetIndex().GetElement(1) - m_ComparisonRadius; yi <= git.GetIndex().GetElement(1) + m_ComparisonRadius; ++yi, ++yj) + for (int yi = index.GetElement(1) - m_ComparisonRadius, yj = y - m_ComparisonRadius; yi <= index.GetElement(1) + m_ComparisonRadius; ++yi, ++yj) { - for (int zi = git.GetIndex().GetElement(2) - m_ComparisonRadius, zj = hit.GetIndex().GetElement(2) - m_ComparisonRadius; zi <= git.GetIndex().GetElement(2) + m_ComparisonRadius; ++zi, ++zj) + 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)) { - niit.SetIndex(indexI); - njit.SetIndex(indexJ); - if (m_UseJointInformation) - { - // if filter should use joint information. A 4d Neighborhood including surrounding channels is used. - for (int d = i - m_ChannelRadius; d <= i + m_ChannelRadius; ++d) - { - if (d >= 0 && d < (int)inputImagePointer->GetVectorLength()) - { - int diff = niit.Get()[d] - njit.Get()[d]; - sumk += (double)(diff*diff); - ++size; - } - } - } - else - { - int diff = niit.Get()[i] - njit.Get()[i]; - sumk += (double)(diff*diff); - ++size; - } + typename InputImageType::PixelType diff = inputImagePointer->GetPixel(indexI) - inputImagePointer->GetPixel(indexJ); + sumk += (double)(diff.GetSquaredNorm()); + ++size; } } } } // weight all neighborhoods - w = std::exp( - (sumk / size) / variance); + w = std::exp( - (sumk / size) / m_Variance); wj.push_back(w); - p.push_back((double)(pixelJ*pixelJ)); + 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 + { + p.push_back(pixelJ); + } Z += w; } } } } for (unsigned int n = 0; n < wj.size(); ++n) { sumj += (wj[n]/Z) * p[n]; } - double df = sumj - (2 * variance); - if (df < 0) + typename InputImageType::PixelType df; + if (m_UseRicianAdaption) { - df = 0; + df = sumj - (2 * m_Variance); + } + else + { + df = sumj; + } + for (unsigned int i = 0; i < inputImagePointer->GetVectorLength(); ++i) + { + double a = df.GetElement(i); + if (a < 0) + { + a = 0; + } + TPixelType outval = std::floor(std::sqrt(a) + 0.5); + outpix.SetElement(i, outval); } - - TPixelType outval = std::floor(std::sqrt(df) + 0.5); - outpix.SetElement(i, outval); } + oit.Set(outpix); ++oit; ++m_CurrentVoxelCount; ++git; } 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(const MaskImageType* mask) { this->SetNthInput(1, const_cast< MaskImageType* >(mask)); } } #endif // __itkNonLocalMeansDenoisingFilter_txx diff --git a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp index e3c70f3b87..5d83ff21ec 100644 --- a/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp +++ b/Modules/DiffusionImaging/DiffusionCore/Testing/mitkNonLocalMeansDenoisingTest.cpp @@ -1,124 +1,124 @@ /*=================================================================== 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); 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"); 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); 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->SetNumberOfThreads(1); m_DenoisingFilter->SetComparisonRadius(1); m_DenoisingFilter->SetSearchRadius(1); } void tearDown() { m_Image = NULL; m_ImageMask = NULL; m_ReferenceImage = NULL; m_DenoisingFilter = NULL; m_DenoisedImage = NULL; } void Denoise_NLMr_1_1_shouldReturnTrue() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMr_1-1.dwi"); m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); 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() { std::string referenceImagePath = GetTestDataFilePath("DiffusionImaging/Denoising/test_multi_NLMv_1-1-1.dwi"); m_ReferenceImage = static_cast*>( mitk::IOUtil::LoadImage(referenceImagePath).GetPointer()); - m_DenoisingFilter->SetChannelRadius(1); +// m_DenoisingFilter->SetChannelRadius(1); 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."); } }; 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 e91810749f..4128697d60 100644 --- a/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp +++ b/Plugins/org.mitk.gui.qt.diffusionimaging/src/internal/QmitkDenoisingView.cpp @@ -1,500 +1,500 @@ /*=================================================================== 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() && m_View->m_BrainMaskNode.IsNotNull()) { switch (m_View->m_SelectedFilter) { case 0: case 3: { break; } case 1: case 2: { try { m_View->m_NoExceptionThrown = true; m_View->m_NonLocalMeansFilter->Update(); } catch (itk::ExceptionObject& e) { m_View->m_NoExceptionThrown = 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_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(NLMR, QString( QApplication::translate("QmitkDenoisingView", "Non local means filter", 0, QApplication::UnicodeUTF8) )); m_Controls->m_SelectFilterComboBox->insertItem(NLMV, QString( QApplication::translate("QmitkDenoisingView", "Non local means filter with joint information", 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"); m_Controls->m_InputBrainMaskLabel->setText("mandatory"); } else { m_Controls->m_InputImageLabel->setText("mandatory"); } 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_BrainMaskNode.IsNotNull()) { if(m_SelectedFilter != NOFILTERSELECTED) { m_Controls->m_ApplyButton->setEnabled(true); } } else if (m_ImageNode.IsNotNull() && m_SelectedFilter == GAUSS) { m_Controls->m_ApplyButton->setEnabled(true); } } void QmitkDenoisingView::StartDenoising() { if (m_ImageNode.IsNotNull() && m_BrainMaskNode.IsNotNull() && m_SelectedFilter != GAUSS) { m_LastProgressCount = 0; switch (m_SelectedFilter) { case NOFILTERSELECTED: case GAUSS: { break; } case NLMR: { // initialize NLMr m_InputImage = dynamic_cast (m_ImageNode->GetData()); m_ImageMask = dynamic_cast(m_BrainMaskNode->GetData()); itk::Image::Pointer itkMask = itk::Image::New(); mitk::CastToItkImage(m_ImageMask, itkMask); m_NonLocalMeansFilter = NonLocalMeansDenoisingFilterType::New(); m_NonLocalMeansFilter->SetNumberOfThreads(12); m_NonLocalMeansFilter->SetInputImage(m_InputImage->GetVectorImage()); m_NonLocalMeansFilter->SetInputMask(itkMask); m_NonLocalMeansFilter->SetUseJointInformation(false); m_NonLocalMeansFilter->SetSearchRadius(m_Controls->m_SpinBoxParameter1->value()); m_NonLocalMeansFilter->SetComparisonRadius(m_Controls->m_SpinBoxParameter2->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 NLMV: { // initialize NLMv m_InputImage = dynamic_cast (m_ImageNode->GetData()); m_ImageMask = dynamic_cast(m_BrainMaskNode->GetData()); itk::Image::Pointer itkMask = itk::Image::New(); mitk::CastToItkImage(m_ImageMask, itkMask); m_NonLocalMeansFilter = NonLocalMeansDenoisingFilterType::New(); - m_NonLocalMeansFilter->SetNumberOfThreads(12); + m_NonLocalMeansFilter->SetNumberOfThreads(1); m_NonLocalMeansFilter->SetInputImage(m_InputImage->GetVectorImage()); m_NonLocalMeansFilter->SetInputMask(itkMask); m_NonLocalMeansFilter->SetUseJointInformation(true); m_NonLocalMeansFilter->SetSearchRadius(m_Controls->m_SpinBoxParameter1->value()); m_NonLocalMeansFilter->SetComparisonRadius(m_Controls->m_SpinBoxParameter2->value()); - m_NonLocalMeansFilter->SetChannelRadius(m_Controls->m_SpinBoxParameter3->value()); +// m_NonLocalMeansFilter->SetChannelRadius(m_Controls->m_SpinBoxParameter3->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; } } } else if(m_SelectedFilter == GAUSS && m_ImageNode.IsNotNull()) { // initialize GAUSS 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->SetInput(extractor->GetOutput()); m_GaussianFilter->SetVariance(m_Controls->m_SpinBoxParameter1->value()); 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); } } void QmitkDenoisingView::ResetParameterPanel() { m_Controls->m_DwiLabel->setEnabled(false); m_Controls->m_InputImageLabel->setEnabled(false); m_Controls->m_BrainMaskLabel->hide(); m_Controls->m_InputBrainMaskLabel->hide(); 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_SpinBoxParameter3->hide(); } 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 = NLMR; m_Controls->m_DwiLabel->setEnabled(true); m_Controls->m_InputImageLabel->setEnabled(true); m_Controls->m_BrainMaskLabel->show(); m_Controls->m_InputBrainMaskLabel->show(); 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_SpinBoxParameter1->show(); m_Controls->m_SpinBoxParameter1->setValue(1); m_Controls->m_SpinBoxParameter2->show(); m_Controls->m_SpinBoxParameter2->setValue(1); break; } case 2: { m_SelectedFilter = NLMV; m_Controls->m_DwiLabel->setEnabled(true); m_Controls->m_InputImageLabel->setEnabled(true); m_Controls->m_BrainMaskLabel->show(); m_Controls->m_InputBrainMaskLabel->show(); 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("Number of neighboring gradients:"); m_Controls->m_SpinBoxParameter1->show(); m_Controls->m_SpinBoxParameter1->setValue(1); m_Controls->m_SpinBoxParameter2->show(); m_Controls->m_SpinBoxParameter2->setValue(1); m_Controls->m_SpinBoxParameter3->show(); m_Controls->m_SpinBoxParameter3->setValue(1); break; } case 3: { m_SelectedFilter = GAUSS; m_Controls->m_DwiLabel->setEnabled(true); m_Controls->m_InputImageLabel->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); m_Controls->m_LabelParameter_2->hide(); m_Controls->m_SpinBoxParameter2->hide(); } } 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 NLMR: case NLMV: { if (m_BrainMaskNode.IsNotNull()) m_Controls->m_ApplyButton->setEnabled(true); break; } 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_LabelParameter_1->setEnabled(false); m_Controls->m_LabelParameter_2->setEnabled(false); m_Controls->m_LabelParameter_3->setEnabled(false); m_Controls->m_SpinBoxParameter1->setEnabled(false); m_Controls->m_SpinBoxParameter2->setEnabled(false); m_Controls->m_SpinBoxParameter3->setEnabled(false); m_Controls->m_SelectFilterComboBox->setEnabled(false); m_Controls->m_ApplyButton->setEnabled(false); } 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_NoExceptionThrown) { switch (m_SelectedFilter) { case NOFILTERSELECTED: case GAUSS: { break; } case NLMR: { 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(); imageNode->SetName((name+"_NLMr_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())).toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); break; } case NLMV: { 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(); imageNode->SetName((name+"_NLMv_"+QString::number(m_Controls->m_SpinBoxParameter1->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter2->value())+"-"+QString::number(m_Controls->m_SpinBoxParameter3->value())).toStdString().c_str()); GetDefaultDataStorage()->Add(imageNode); m_Controls->m_LabelParameter_3->setEnabled(true); m_Controls->m_SpinBoxParameter3->setEnabled(true); break; } } } m_Controls->m_LabelParameter_1->setEnabled(true); m_Controls->m_LabelParameter_2->setEnabled(true); m_Controls->m_SpinBoxParameter1->setEnabled(true); m_Controls->m_SpinBoxParameter2->setEnabled(true); m_Controls->m_SelectFilterComboBox->setEnabled(true); m_Controls->m_ApplyButton->setEnabled(true); } void QmitkDenoisingView::UpdateProgress() { switch (m_SelectedFilter) { case NOFILTERSELECTED: case GAUSS: { break; } case NLMR: case NLMV: { unsigned int currentProgressCount = m_NonLocalMeansFilter->GetCurrentVoxelCount(); mitk::ProgressBar::GetInstance()->Progress(currentProgressCount-m_LastProgressCount); m_LastProgressCount = currentProgressCount; break; } } }