From 2a68dc6fc6e5c1d2232937139282b06011003988 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Sat, 14 Jul 2018 20:21:30 +0200 Subject: [PATCH] Cleaning, renaming and making Rician to Gaussian noise tool compile only when Arb is linked with Anima. Otherwise, unstable behaviour. --- Anima/filtering/CMakeLists.txt | 2 +- .../rician_noise_to_gaussian/CMakeLists.txt | 37 --- .../animaRiceToGaussianImageFilter.h | 104 -------- .../animaRiceToGaussianImageFilter.hxx | 247 ----------------- .../animaRicianToGaussianNoise.cxx | 250 ------------------ 5 files changed, 1 insertion(+), 639 deletions(-) delete mode 100644 Anima/filtering/rician_noise_to_gaussian/CMakeLists.txt delete mode 100644 Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.h delete mode 100644 Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.hxx delete mode 100644 Anima/filtering/rician_noise_to_gaussian/animaRicianToGaussianNoise.cxx diff --git a/Anima/filtering/CMakeLists.txt b/Anima/filtering/CMakeLists.txt index 75c36a1d5..20d208af2 100644 --- a/Anima/filtering/CMakeLists.txt +++ b/Anima/filtering/CMakeLists.txt @@ -7,7 +7,7 @@ project(ANIMA-FILTERING) add_subdirectory(bias_correction) add_subdirectory(denoising) add_subdirectory(dti_tools) -add_subdirectory(rician_noise_to_gaussian) +add_subdirectory(rician_to_gaussian_noise) add_subdirectory(noise_generator) add_subdirectory(nyul_standardization) add_subdirectory(regularization) diff --git a/Anima/filtering/rician_noise_to_gaussian/CMakeLists.txt b/Anima/filtering/rician_noise_to_gaussian/CMakeLists.txt deleted file mode 100644 index cf7eac2a5..000000000 --- a/Anima/filtering/rician_noise_to_gaussian/CMakeLists.txt +++ /dev/null @@ -1,37 +0,0 @@ -if(BUILD_TOOLS) - -project(animaRicianToGaussianNoise) - -## ############################################################################# -## List Sources -## ############################################################################# - -list_source_files(${PROJECT_NAME} - ${CMAKE_CURRENT_SOURCE_DIR} - ) - -## ############################################################################# -## add executable -## ############################################################################# - -add_executable(${PROJECT_NAME} - ${${PROJECT_NAME}_CFILES} - ) - - -## ############################################################################# -## Link -## ############################################################################# - -target_link_libraries(${PROJECT_NAME} - ${ITKIO_LIBRARIES} - AnimaSpecialFunctions - ) - -## ############################################################################# -## install -## ############################################################################# - -set_exe_install_rules(${PROJECT_NAME}) - -endif() diff --git a/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.h b/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.h deleted file mode 100644 index 24c5f4d7b..000000000 --- a/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.h +++ /dev/null @@ -1,104 +0,0 @@ -#pragma once - -#include - -#include -#include - -namespace anima -{ -template -class RiceToGaussianImageFilter : -public anima::MaskedImageToImageFilter,itk::Image > -{ -public: - /** Standard class typedefs. */ - typedef RiceToGaussianImageFilter Self; - - typedef itk::Image InputImageType; - typedef typename InputImageType::PixelType InputPixelType; - - typedef itk::Image OutputImageType; - typedef typename OutputImageType::PixelType OutputPixelType; - - typedef anima::MaskedImageToImageFilter Superclass; - - typedef itk::SmartPointer Pointer; - typedef itk::SmartPointer ConstPointer; - - /** Method for creation through the object factory. */ - itkNewMacro(Self) - - /** Run-time type information (and related methods) */ - itkTypeMacro(RiceToGaussianImageFilter, MaskedImageToImageFilter) - - /** Superclass typedefs. */ - typedef typename Superclass::InputImageRegionType InputImageRegionType; - typedef typename Superclass::OutputImageRegionType OutputImageRegionType; - typedef typename Superclass::MaskImageType MaskImageType; - typedef typename MaskImageType::SizeType SizeType; - - typedef typename InputImageType::Pointer InputPointerType; - typedef typename OutputImageType::Pointer OutputPointerType; - - itkSetMacro(MaximumNumberOfIterations, unsigned int) - itkSetMacro(Epsilon, double) - itkSetMacro(Sigma, double) - itkSetMacro(MeanImage, InputPointerType) - itkSetMacro(VarianceImage, InputPointerType) - - itkSetMacro(Scale, double) - itkGetConstMacro(Scale, double) - - itkSetMacro(Alpha, double) - itkGetConstMacro(Alpha, double) - - OutputPointerType GetLocationImage() {return this->GetOutput(0);} - OutputPointerType GetScaleImage() {return this->GetOutput(1);} - OutputPointerType GetGaussianImage() {return this->GetOutput(2);} - -protected: - RiceToGaussianImageFilter() - : Superclass() - { - m_MaximumNumberOfIterations = 100; - m_Epsilon = 1.0e-8; - m_Sigma = 1.0; - m_Scale = 0.0; - m_Alpha = 0.005; - m_ThreadScaleSamples.clear(); - m_NeighborWeights.clear(); - m_Radius.Fill(0); - m_MeanImage = NULL; - m_VarianceImage = NULL; - - this->SetNumberOfRequiredOutputs(3); - this->SetNthOutput(0, this->MakeOutput(0)); - this->SetNthOutput(1, this->MakeOutput(1)); - this->SetNthOutput(2, this->MakeOutput(2)); - } - - virtual ~RiceToGaussianImageFilter() {} - - void BeforeThreadedGenerateData(void) ITK_OVERRIDE; - void ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, - itk::ThreadIdType threadId) ITK_OVERRIDE; - void AfterThreadedGenerateData(void) ITK_OVERRIDE; - -private: - ITK_DISALLOW_COPY_AND_ASSIGN(RiceToGaussianImageFilter); - - unsigned int m_MaximumNumberOfIterations; - double m_Epsilon, m_Sigma, m_Scale, m_Alpha; - std::vector > m_ThreadScaleSamples; - SizeType m_Radius; - std::vector m_NeighborWeights; - InputPointerType m_MeanImage, m_VarianceImage; - - boost::math::normal_distribution<> m_NormalDistribution; - boost::math::rayleigh_distribution<> m_RayleighDistribution; -}; - -} // end of namespace anima - -#include "animaRiceToGaussianImageFilter.hxx" diff --git a/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.hxx b/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.hxx deleted file mode 100644 index af55812a4..000000000 --- a/Anima/filtering/rician_noise_to_gaussian/animaRiceToGaussianImageFilter.hxx +++ /dev/null @@ -1,247 +0,0 @@ -#pragma once -#include "animaRiceToGaussianImageFilter.h" -#include -#include - -#include -#include -#include -#include - -namespace anima -{ - -template -void -RiceToGaussianImageFilter -::BeforeThreadedGenerateData(void) -{ - // Compute spatial weights of neighbors beforehand - typename InputImageType::SpacingType spacing = this->GetInput()->GetSpacing(); - typedef itk::GaussianOperator GaussianOperatorType; - std::vector gaussianKernels(ImageDimension); - - for (unsigned int i = 0;i < ImageDimension;++i) - { - unsigned int reverse_i = ImageDimension - i - 1; - double stddev = m_Sigma / spacing[i]; - gaussianKernels[reverse_i].SetDirection(i); - gaussianKernels[reverse_i].SetVariance(stddev * stddev); - gaussianKernels[reverse_i].SetMaximumError(1e-3); - gaussianKernels[reverse_i].CreateDirectional(); - gaussianKernels[reverse_i].ScaleCoefficients(1.0e4); - m_Radius[i] = gaussianKernels[reverse_i].GetRadius(i); - } - - m_NeighborWeights.clear(); - for (unsigned int i = 0;i < gaussianKernels[0].Size();++i) - { - for (unsigned int j = 0;j < gaussianKernels[1].Size();++j) - { - for (unsigned int k = 0;k < gaussianKernels[2].Size();++k) - { - double weight = gaussianKernels[0][i] * gaussianKernels[1][j] * gaussianKernels[2][k]; - m_NeighborWeights.push_back(weight); - } - } - } - - // Initialize thread containers for global scale estimation - unsigned int numThreads = this->GetNumberOfThreads(); - m_ThreadScaleSamples.resize(numThreads); - this->Superclass::BeforeThreadedGenerateData(); -} - -template -void -RiceToGaussianImageFilter -::ThreadedGenerateData(const OutputImageRegionType& outputRegionForThread, - itk::ThreadIdType threadId) -{ - typedef itk::ConstNeighborhoodIterator InputIteratorType; - typedef itk::ConstNeighborhoodIterator MaskIteratorType; - typedef itk::ImageRegionConstIterator InputSimpleIteratorType; - typedef itk::ImageRegionIterator OutputIteratorType; - - InputIteratorType inputItr(m_Radius, this->GetInput(), outputRegionForThread); - unsigned int neighborhoodSize = inputItr.Size(); - - MaskIteratorType maskItr(m_Radius, this->GetComputationMask(), outputRegionForThread); - - InputSimpleIteratorType meanItr, varItr; - if (m_MeanImage) - meanItr = InputSimpleIteratorType(m_MeanImage, outputRegionForThread); - if (m_VarianceImage) - varItr = InputSimpleIteratorType(m_VarianceImage, outputRegionForThread); - - OutputIteratorType locationItr(this->GetOutput(0), outputRegionForThread); - OutputIteratorType scaleItr(this->GetOutput(1), outputRegionForThread); - OutputIteratorType signalItr(this->GetOutput(2), outputRegionForThread); - - std::vector samples, weights; - bool isInBounds; - typename InputImageType::IndexType currentIndex, neighborIndex; - typename InputImageType::PointType currentPoint, neighborPoint; - - while (!maskItr.IsAtEnd()) - { - // Discard voxels outside of the brain - if (maskItr.GetCenterPixel() == 0) - { - locationItr.Set(0); - scaleItr.Set(0); - signalItr.Set(0); - - ++inputItr; - ++maskItr; - if (m_MeanImage) - ++meanItr; - if (m_VarianceImage) - ++varItr; - ++locationItr; - ++scaleItr; - ++signalItr; - - continue; - } - - // Get signal at current central voxel - double inputSignal = inputItr.GetCenterPixel(); - - // Rice-corrupted signals should all be positive - if (inputSignal <= 0) - { - locationItr.Set(0); - scaleItr.Set(0); - signalItr.Set(0); - - ++inputItr; - ++maskItr; - if (m_MeanImage) - ++meanItr; - if (m_VarianceImage) - ++varItr; - ++locationItr; - ++scaleItr; - ++signalItr; - - continue; - } - - // Estimation of location and scale - double location = 0; - double scale = m_Scale; - - if (m_MeanImage && m_VarianceImage) - { - // If mean and variance images are available, use them instead of neighborhood - double sigmaValue = std::sqrt(varItr.Get()); - double rValue = meanItr.Get() / sigmaValue; - double k1Value = 0; - double thetaValue = anima::FixedPointFinder(rValue, 1, k1Value); - k1Value = anima::KummerFunction(-thetaValue * thetaValue / 2.0, -0.5, 1); - - scale = sigmaValue / std::sqrt(anima::XiFunction(thetaValue * thetaValue, 1, k1Value, M_PI / 2.0)); - location = thetaValue * scale; - } - else - { - // Use neighbors to create samples - samples.clear(); - weights.clear(); - - for (unsigned int i = 0; i < neighborhoodSize; ++i) - { - double tmpVal = static_cast(inputItr.GetPixel(i, isInBounds)); - - if (isInBounds && !std::isnan(tmpVal) && std::isfinite(tmpVal)) - { - if (maskItr.GetPixel(i) != maskItr.GetCenterPixel()) - continue; - - double weight = m_NeighborWeights[i]; - - if (weight < m_Epsilon) - continue; - - samples.push_back(tmpVal); - weights.push_back(weight); - } - } - - if (samples.size() == 1) - location = inputSignal; - else - anima::GetRiceParameters(samples, weights, location, scale); - } - - // Transform Rice signal in Gaussian signal - double outputSignal = inputSignal; - double snrValue = location / scale; - - if (!std::isfinite(snrValue) || std::isnan(snrValue)) - { - std::cout << snrValue << " " << location << " " << scale << std::endl; - itkExceptionMacro("Estimated SNR is invalid"); - } - - if (snrValue <= m_Epsilon) - outputSignal = 0.0; - else if (snrValue <= 0.1) - { - double unifSignal = boost::math::cdf(m_RayleighDistribution, inputSignal / scale); - - if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha) - unifSignal = boost::math::cdf(m_RayleighDistribution, snrValue); - - outputSignal = scale * boost::math::quantile(m_NormalDistribution, unifSignal); - } - else if (snrValue <= 600) // if SNR if > 600 keep signal as is, else... - { - double unifSignal = anima::GetRiceCDF(inputSignal, location, scale); - - if (unifSignal >= 1.0 - m_Alpha || unifSignal <= m_Alpha) - unifSignal = anima::GetRiceCDF(location, location, scale); - - outputSignal = location + scale * boost::math::quantile(m_NormalDistribution, unifSignal); - } - - m_ThreadScaleSamples[threadId].push_back(scale); - - locationItr.Set(static_cast(location)); - scaleItr.Set(static_cast(scale)); - signalItr.Set(static_cast(outputSignal)); - - ++inputItr; - ++maskItr; - if (m_MeanImage) - ++meanItr; - if (m_VarianceImage) - ++varItr; - ++locationItr; - ++scaleItr; - ++signalItr; - } -} - -template -void -RiceToGaussianImageFilter -::AfterThreadedGenerateData(void) -{ - if (m_Scale == 0) - { - std::vector scaleSamples; - for (unsigned int i = 0;i < this->GetNumberOfThreads();++i) - scaleSamples.insert(scaleSamples.end(), m_ThreadScaleSamples[i].begin(), m_ThreadScaleSamples[i].end()); - - m_Scale = anima::GetMedian(scaleSamples); - } - - m_ThreadScaleSamples.clear(); - m_NeighborWeights.clear(); - - this->Superclass::AfterThreadedGenerateData(); -} - -} // end of namespace anima diff --git a/Anima/filtering/rician_noise_to_gaussian/animaRicianToGaussianNoise.cxx b/Anima/filtering/rician_noise_to_gaussian/animaRicianToGaussianNoise.cxx deleted file mode 100644 index abf107393..000000000 --- a/Anima/filtering/rician_noise_to_gaussian/animaRicianToGaussianNoise.cxx +++ /dev/null @@ -1,250 +0,0 @@ -#include -#include -#include - -#include - -#include -#include -#include - -//Update progression of the process -void eventCallback (itk::Object* caller, const itk::EventObject& event, void* clientData) -{ - itk::ProcessObject * processObject = (itk::ProcessObject*) caller; - std::cout<<"\033[K\rProgression: "<<(int)(processObject->GetProgress() * 100)<<"%"< inputArg("i", "input", "Input Rice-corrupted image", true, "", "Input Rice-corrupted image", cmd); - TCLAP::ValueArg outputArg("o", "output", "Output Gaussian-corrupted image", true, "", "Output Gaussian-corrupted image", cmd); - TCLAP::ValueArg bvecArg("g","grad","Input gradients",true,"","Input gradients",cmd); - TCLAP::ValueArg bvalArg("b","bval","Input b-values",true,"","Input b-values",cmd); - - TCLAP::ValueArg scaleArg("s", "scale", "Output scale image", false, "", "Output scale image", cmd); - TCLAP::ValueArg maskArg("m", "mask", "Segmentation mask", false, "", "Segmentation mask", cmd); - - TCLAP::ValueArg epsArg("E", "epsilon", "Minimal absolute value difference betweem fixed point iterations (default: 1e-8)", false, 1.0e-8, "Minimal absolute value difference betweem fixed point iterations", cmd); - TCLAP::ValueArg sigmaArg("S", "sigma", "Gaussian standard deviation for defining neighbor weights (default: 1.0)", false, 1.0, "Gaussian standard deviation for defining neighbor weights", cmd); - - TCLAP::ValueArg maxiterArg("I", "maxiter", "Maximum number of iterations (default: 100)", false, 100, "Maximum number of iterations", cmd); - TCLAP::ValueArg nbpArg("T", "nbp", "Number of threads to run on -> default : automatically determine", false, itk::MultiThreader::GetGlobalDefaultNumberOfThreads(), "Number of threads", cmd); - - try - { - cmd.parse(ac,av); - } - catch (TCLAP::ArgException& e) - { - std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl; - return EXIT_FAILURE; - } - - const unsigned int ImageDimension = 3; - - typedef anima::GradientFileReader ,double> GFReaderType; - GFReaderType gfReader; - gfReader.SetGradientFileName(bvecArg.getValue()); - gfReader.SetBValueBaseString(bvalArg.getValue()); - gfReader.Update(); - - GFReaderType::GradientVectorType directions = gfReader.GetGradients(); - GFReaderType::BValueVectorType bvalues = gfReader.GetBValues(); - - typedef itk::Image Image4DType; - typedef anima::RiceToGaussianImageFilter RiceToGaussianImageFilterType; - typedef RiceToGaussianImageFilterType::InputImageType InputImageType; - typedef RiceToGaussianImageFilterType::OutputImageType OutputImageType; - typedef RiceToGaussianImageFilterType::MaskImageType MaskImageType; - typedef itk::SubtractImageFilter SubtractFilterType; - typedef itk::MultiplyImageFilter MultiplyImageFilterType; - typedef itk::NaryAddImageFilter NaryAddImageFilterType; - - Image4DType::Pointer input = anima::readImage(inputArg.getValue()); - std::vector inputData; - inputData = anima::getImagesFromHigherDimensionImage(input); - - unsigned int numberOfImages = inputData.size(); - unsigned int numberOfB0Images = 0; - unsigned int indexOfFirstB0Image = 0; - - for (unsigned int i = 0;i < numberOfImages;++i) - { - if (bvalues[i] <= 10) - { - if (numberOfB0Images == 0) - indexOfFirstB0Image = i; - ++numberOfB0Images; - } - } - - std::cout << "Number of B0 images: " << numberOfB0Images << std::endl; - - itk::CStyleCommand::Pointer callback; - RiceToGaussianImageFilterType::Pointer mainFilter; - - // Estimate scale parameter from B0 image(s) - if (numberOfB0Images == 1) - { - mainFilter = RiceToGaussianImageFilterType::New(); - mainFilter->SetInput(inputData[indexOfFirstB0Image]); - mainFilter->SetMaximumNumberOfIterations(maxiterArg.getValue()); - mainFilter->SetEpsilon(epsArg.getValue()); - mainFilter->SetSigma(sigmaArg.getValue()); - mainFilter->SetNumberOfThreads(nbpArg.getValue()); - - if (maskArg.getValue() != "") - mainFilter->SetComputationMask(anima::readImage(maskArg.getValue())); - - callback = itk::CStyleCommand::New(); - callback->SetCallback(eventCallback); - mainFilter->AddObserver(itk::ProgressEvent(), callback); - - std::cout << "First run of Gaussianizer on the B0 image for estimating the scale parameter..." << std::endl; - } - else - { - NaryAddImageFilterType::Pointer addFilter = NaryAddImageFilterType::New(); - unsigned int pos = 0; - for (unsigned int i = 0;i < numberOfImages;++i) - { - if (bvalues[i] > 10) - continue; - - addFilter->SetInput(pos, inputData[i]); - ++pos; - } - - addFilter->SetNumberOfThreads(nbpArg.getValue()); - addFilter->Update(); - - MultiplyImageFilterType::Pointer mulFilter = MultiplyImageFilterType::New(); - mulFilter->SetInput1(addFilter->GetOutput()); - mulFilter->SetConstant2(1.0 / (double)numberOfB0Images); - mulFilter->SetNumberOfThreads(nbpArg.getValue()); - mulFilter->Update(); - - InputImageType::Pointer meanImage = mulFilter->GetOutput(); - - addFilter = NaryAddImageFilterType::New(); - pos = 0; - for (unsigned int i = 0;i < numberOfImages;++i) - { - if (bvalues[i] > 10) - continue; - - SubtractFilterType::Pointer subFilter = SubtractFilterType::New(); - subFilter->SetInput1(inputData[i]); - subFilter->SetInput2(meanImage); - subFilter->SetNumberOfThreads(nbpArg.getValue()); - subFilter->Update(); - - mulFilter = MultiplyImageFilterType::New(); - mulFilter->SetInput1(subFilter->GetOutput()); - mulFilter->SetInput2(subFilter->GetOutput()); - mulFilter->SetNumberOfThreads(nbpArg.getValue()); - mulFilter->Update(); - - addFilter->SetInput(pos, mulFilter->GetOutput()); - ++pos; - } - - addFilter->SetNumberOfThreads(nbpArg.getValue()); - addFilter->Update(); - - mulFilter = MultiplyImageFilterType::New(); - mulFilter->SetInput1(addFilter->GetOutput()); - mulFilter->SetConstant2(1.0 / (numberOfB0Images - 1.0)); - mulFilter->SetNumberOfThreads(nbpArg.getValue()); - mulFilter->Update(); - - InputImageType::Pointer varianceImage = mulFilter->GetOutput(); - - mainFilter = RiceToGaussianImageFilterType::New(); - mainFilter->SetInput(inputData[indexOfFirstB0Image]); - mainFilter->SetMaximumNumberOfIterations(maxiterArg.getValue()); - mainFilter->SetEpsilon(epsArg.getValue()); - mainFilter->SetSigma(sigmaArg.getValue()); - mainFilter->SetMeanImage(meanImage); - mainFilter->SetVarianceImage(varianceImage); - mainFilter->SetNumberOfThreads(nbpArg.getValue()); - - if (maskArg.getValue() != "") - mainFilter->SetComputationMask(anima::readImage(maskArg.getValue())); - - callback = itk::CStyleCommand::New(); - callback->SetCallback(eventCallback); - mainFilter->AddObserver(itk::ProgressEvent(), callback); - - std::cout << "First run of Gaussianizer using the mean and variance of the B0 images for estimating the scale parameter..." << std::endl; - } - - try - { - mainFilter->Update(); - } - catch (itk::ExceptionObject & err) - { - std::cerr << err << std::endl; - return EXIT_FAILURE; - } - - double scale = mainFilter->GetScale(); - - if (scaleArg.getValue() != "") - anima::writeImage(scaleArg.getValue(), mainFilter->GetScaleImage()); - - std::cout << "\nDone. Estimated scale parameter: " << scale << std::endl; - - // Transform Rice to Gaussian data - - for (unsigned int i = 0;i < numberOfImages;++i) - { - std::cout << "\nRun Gaussianizer for Image "<< i+1 << " with fixed scale parameter..." << std::endl; - - mainFilter = RiceToGaussianImageFilterType::New(); - mainFilter->SetInput(inputData[i]); - mainFilter->SetMaximumNumberOfIterations(maxiterArg.getValue()); - mainFilter->SetEpsilon(epsArg.getValue()); - mainFilter->SetSigma(sigmaArg.getValue()); - mainFilter->SetScale(scale); - mainFilter->SetNumberOfThreads(nbpArg.getValue()); - - if (maskArg.getValue() != "") - mainFilter->SetComputationMask(anima::readImage(maskArg.getValue())); - - callback = itk::CStyleCommand::New(); - callback->SetCallback(eventCallback); - mainFilter->AddObserver(itk::ProgressEvent(), callback); - - try - { - mainFilter->Update(); - } - catch (itk::ExceptionObject & err) - { - std::cerr << err << std::endl; - return EXIT_FAILURE; - } - - anima::writeImage(outputArg.getValue() + to_string(i + 1, 3) + ".nrrd", mainFilter->GetGaussianImage()); - } - - return EXIT_SUCCESS; -}