From 23d910316bf423a4d5b611f45ecaa64842bef844 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Wed, 7 Mar 2018 15:50:49 +0100 Subject: [PATCH] Add filter for converting Rice-corrupted images into Gaussian-corrupted images --- Anima/cmake/ANIMAConfig.cmake.in | 20 ++ Anima/cmake/ANIMAExports.cmake.in | 1 + Anima/cmake/ANIMAUse.cmake.in | 36 +++ Anima/cmake/AnimaDependencies.cmake | 40 +++ Anima/cmake/AnimaModulesSetup.cmake | 4 + Anima/cmake/FindARB.cmake | 29 ++ Anima/cmake/FindFLINT.cmake | 29 ++ Anima/cmake/FindGMP.cmake | 29 ++ Anima/cmake/FindMPFR.cmake | 29 ++ .../animaDTIEstimationImageFilter.hxx | 4 +- .../animaMCMEstimatorImageFilter.hxx | 6 +- Anima/filtering/CMakeLists.txt | 5 + Anima/filtering/gaussianizer/CMakeLists.txt | 37 +++ .../gaussianizer/animaGaussianizer.cxx | 250 ++++++++++++++++++ .../animaRiceToGaussianImageFilter.h | 104 ++++++++ .../animaRiceToGaussianImageFilter.hxx | 247 +++++++++++++++++ Anima/math-tools/CMakeLists.txt | 4 + .../arb_special_functions/CMakeLists.txt | 43 +++ .../animaARBBesselFunctions.cxx | 57 ++++ .../animaARBBesselFunctions.h | 33 +++ .../animaARBKummerFunctions.cxx | 37 +++ .../animaARBKummerFunctions.h | 11 + .../common_tools/convert_shape/CMakeLists.txt | 2 +- .../convert_shape/animaConvertShape.cxx | 9 +- .../data_io/animaGradientFileReader.hxx | 6 +- .../matrix_operations/animaVectorOperations.h | 16 ++ .../animaVectorOperations.hxx | 60 +++++ .../animaBesselFunctions.cxx | 11 +- .../special_functions/animaBesselFunctions.h | 3 + .../animaKummerFunctions.cxx | 1 + .../animaChiDistribution.h | 28 ++ .../animaChiDistribution.hxx | 183 +++++++++++++ CMakeLists.txt | 29 ++ 33 files changed, 1388 insertions(+), 15 deletions(-) create mode 100644 Anima/cmake/FindARB.cmake create mode 100644 Anima/cmake/FindFLINT.cmake create mode 100644 Anima/cmake/FindGMP.cmake create mode 100644 Anima/cmake/FindMPFR.cmake create mode 100644 Anima/filtering/gaussianizer/CMakeLists.txt create mode 100644 Anima/filtering/gaussianizer/animaGaussianizer.cxx create mode 100644 Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.h create mode 100644 Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.hxx create mode 100644 Anima/math-tools/arb_special_functions/CMakeLists.txt create mode 100644 Anima/math-tools/arb_special_functions/animaARBBesselFunctions.cxx create mode 100644 Anima/math-tools/arb_special_functions/animaARBBesselFunctions.h create mode 100644 Anima/math-tools/arb_special_functions/animaARBKummerFunctions.cxx create mode 100644 Anima/math-tools/arb_special_functions/animaARBKummerFunctions.h create mode 100644 Anima/math-tools/statistical_distributions/animaChiDistribution.h create mode 100644 Anima/math-tools/statistical_distributions/animaChiDistribution.hxx diff --git a/Anima/cmake/ANIMAConfig.cmake.in b/Anima/cmake/ANIMAConfig.cmake.in index d0811cceb..84177572c 100644 --- a/Anima/cmake/ANIMAConfig.cmake.in +++ b/Anima/cmake/ANIMAConfig.cmake.in @@ -24,6 +24,26 @@ set(NLOPT_INCLUDE_DIR "@NLOPT_INCLUDE_DIR@") set(NLOPT_LIBRARY_DIR "@NLOPT_LIBRARY_DIR@") set(NLOPT_LIBRARY "@NLOPT_LIBRARY@") +set(USE_GMP "@USE_GMP@") +set(GMP_INCLUDE_DIR "@GMP_INCLUDE_DIR@") +set(GMP_LIBRARY_DIR "@GMP_LIBRARY_DIR@") +set(GMP_LIBRARY "@GMP_LIBRARY@") + +set(USE_MPFR "@USE_MPFR@") +set(MPFR_INCLUDE_DIR "@MPFR_INCLUDE_DIR@") +set(MPFR_LIBRARY_DIR "@MPFR_LIBRARY_DIR@") +set(MPFR_LIBRARY "@MPFR_LIBRARY@") + +set(USE_FLINT "@USE_FLINT@") +set(FLINT_INCLUDE_DIR "@FLINT_INCLUDE_DIR@") +set(FLINT_LIBRARY_DIR "@FLINT_LIBRARY_DIR@") +set(FLINT_LIBRARY "@FLINT_LIBRARY@") + +set(USE_ARB "@USE_ARB@") +set(ARB_INCLUDE_DIR "@ARB_INCLUDE_DIR@") +set(ARB_LIBRARY_DIR "@ARB_LIBRARY_DIR@") +set(ARB_LIBRARY "@ARB_LIBRARY@") + set(TinyXML2_INCLUDE_DIR "@TinyXML2_INCLUDE_DIR@") set(TinyXML2_LIBRARY_DIR "@TinyXML2_LIBRARY_DIR@") set(TinyXML2_LIBRARY "@TinyXML2_LIBRARY@") diff --git a/Anima/cmake/ANIMAExports.cmake.in b/Anima/cmake/ANIMAExports.cmake.in index 019a1f0b8..080a275cd 100644 --- a/Anima/cmake/ANIMAExports.cmake.in +++ b/Anima/cmake/ANIMAExports.cmake.in @@ -1,4 +1,5 @@ include_directories(@ANIMA_BINARY_DIR@) +include_directories(@ANIMA_BINARY_DIR@/math-tools/arb_special_functions) include_directories(@ANIMA_BINARY_DIR@/math-tools/data_io) include_directories(@ANIMA_BINARY_DIR@/math-tools/optimizers) include_directories(@ANIMA_BINARY_DIR@/math-tools/multi_compartment_base) diff --git a/Anima/cmake/ANIMAUse.cmake.in b/Anima/cmake/ANIMAUse.cmake.in index bcbc1485f..2b31cfee8 100644 --- a/Anima/cmake/ANIMAUse.cmake.in +++ b/Anima/cmake/ANIMAUse.cmake.in @@ -37,6 +37,42 @@ if(NOT ANIMA_USE_FILE_INCLUDED) set(NLOPT_FOUND 1) endif(USE_NLOPT) + set(USE_GMP ${USE_GMP}) + if (USE_GMP) + set(GMP_INCLUDE_DIR ${GMP_INCLUDE_DIR}) + include_directories(${GMP_INCLUDE_DIR}) + set(GMP_LIBRARY_DIR ${GMP_LIBRARY_DIR}) + set(GMP_LIBRARY ${GMP_LIBRARY}) + set(GMP_FOUND 1) + endif(USE_GMP) + + set(USE_MPFR ${USE_MPFR}) + if (USE_MPFR) + set(MPFR_INCLUDE_DIR ${MPFR_INCLUDE_DIR}) + include_directories(${MPFR_INCLUDE_DIR}) + set(MPFR_LIBRARY_DIR ${MPFR_LIBRARY_DIR}) + set(MPFR_LIBRARY ${MPFR_LIBRARY}) + set(MPFR_FOUND 1) + endif(USE_MPFR) + + set(USE_FLINT ${USE_FLINT}) + if (USE_FLINT) + set(FLINT_INCLUDE_DIR ${FLINT_INCLUDE_DIR}) + include_directories(${FLINT_INCLUDE_DIR}) + set(FLINT_LIBRARY_DIR ${FLINT_LIBRARY_DIR}) + set(FLINT_LIBRARY ${FLINT_LIBRARY}) + set(FLINT_FOUND 1) + endif(USE_FLINT) + + set(USE_ARB ${USE_ARB}) + if (USE_ARB) + set(ARB_INCLUDE_DIR ${ARB_INCLUDE_DIR}) + include_directories(${ARB_INCLUDE_DIR}) + set(ARB_LIBRARY_DIR ${ARB_LIBRARY_DIR}) + set(ARB_LIBRARY ${ARB_LIBRARY}) + set(ARB_FOUND 1) + endif(USE_ARB) + set(TinyXML2_INCLUDE_DIR ${TinyXML2_INCLUDE_DIR}) include_directories(${TinyXML2_INCLUDE_DIR}) set(TinyXML2_LIBRARY_DIR ${TinyXML2_LIBRARY_DIR}) diff --git a/Anima/cmake/AnimaDependencies.cmake b/Anima/cmake/AnimaDependencies.cmake index 4bfe32386..6853abf9a 100644 --- a/Anima/cmake/AnimaDependencies.cmake +++ b/Anima/cmake/AnimaDependencies.cmake @@ -30,6 +30,46 @@ if(USE_NLOPT) find_package(NLOPT REQUIRED) endif() +# GMP +option(USE_GMP + "Use GMP external library (necessary for some special functions)" + OFF + ) + +if(USE_GMP) + find_package(GMP REQUIRED) +endif() + +# MPFR +option(USE_MPFR + "Use MPFR external library (necessary for some special functions)" + OFF + ) + +if(USE_MPFR) + find_package(MPFR REQUIRED) +endif() + +# FLINT +option(USE_FLINT + "Use FLINT external library (necessary for some special functions)" + OFF + ) + +if(USE_FLINT) + find_package(FLINT REQUIRED) +endif() + +# ARB +option(USE_ARB + "Use ARB external library (necessary for some special functions)" + OFF + ) + +if(USE_ARB) + find_package(ARB REQUIRED) +endif() + # TinyXML2 if (BUILD_MODULE_REGISTRATION OR BUILD_MODULE_DIFFUSION) find_package(TinyXML2 REQUIRED) diff --git a/Anima/cmake/AnimaModulesSetup.cmake b/Anima/cmake/AnimaModulesSetup.cmake index 6b22c029e..5e0f924ff 100644 --- a/Anima/cmake/AnimaModulesSetup.cmake +++ b/Anima/cmake/AnimaModulesSetup.cmake @@ -30,6 +30,10 @@ set(${PROJECT_NAME}_INCLUDE_DIRS "${Boost_INCLUDE_DIR}" "${TCLAP_INCLUDE_DIR}" "${NLOPT_INCLUDE_DIR}" + "${GMP_INCLUDE_DIR}" + "${MPFR_INCLUDE_DIR}" + "${FLINT_INCLUDE_DIR}" + "${ARB_INCLUDE_DIR}" "${TinyXML2_INCLUDE_DIR}" ) diff --git a/Anima/cmake/FindARB.cmake b/Anima/cmake/FindARB.cmake new file mode 100644 index 000000000..bdfb10b5e --- /dev/null +++ b/Anima/cmake/FindARB.cmake @@ -0,0 +1,29 @@ +# Try to find the arb library +# Once done this will define +# +# ARB_FOUND - system has arb and it can be used +# ARB_INCLUDE_DIR - directory where the header file can be found +# ARB_LIBRARY - Path where arb library file can be found +# + +SET (ARB_FOUND FALSE) + +FIND_PATH(ARB_INCLUDE_DIR arb.h + /usr/include + /usr/local/include +) + +set(ARB_LIBRARY_DIR "" CACHE PATH "ARB library folder") + +FIND_LIBRARY(ARB_LIBRARY NAMES arb libarb + HINTS "${ARB_LIBRARY_DIR}" +) + +mark_as_advanced(ARB_LIBRARY) + +IF( EXISTS "${ARB_INCLUDE_DIR}" AND EXISTS "${ARB_LIBRARY}" ) + set(ARB_FOUND TRUE) +else() + message(FATAL_ERROR "ARB not found") +endif() + diff --git a/Anima/cmake/FindFLINT.cmake b/Anima/cmake/FindFLINT.cmake new file mode 100644 index 000000000..1a2bdfdd1 --- /dev/null +++ b/Anima/cmake/FindFLINT.cmake @@ -0,0 +1,29 @@ +# Try to find the flint library +# Once done this will define +# +# FLINT_FOUND - system has flint and it can be used +# FLINT_INCLUDE_DIR - directory where the header file can be found +# FLINT_LIBRARY - Path where flint library file can be found +# + +SET (FLINT_FOUND FALSE) + +FIND_PATH(FLINT_INCLUDE_DIR flint/flint.h + /usr/include + /usr/local/include +) + +set(FLINT_LIBRARY_DIR "" CACHE PATH "FLINT library folder") + +FIND_LIBRARY(FLINT_LIBRARY NAMES flint libflint + HINTS "${FLINT_LIBRARY_DIR}" +) + +mark_as_advanced(FLINT_LIBRARY) + +IF( EXISTS "${FLINT_INCLUDE_DIR}" AND EXISTS "${FLINT_LIBRARY}" ) + set(FLINT_FOUND TRUE) +else() + message(FATAL_ERROR "FLINT not found") +endif() + diff --git a/Anima/cmake/FindGMP.cmake b/Anima/cmake/FindGMP.cmake new file mode 100644 index 000000000..9a72747df --- /dev/null +++ b/Anima/cmake/FindGMP.cmake @@ -0,0 +1,29 @@ +# Try to find the gmp library +# Once done this will define +# +# GMP_FOUND - system has gmp and it can be used +# GMP_INCLUDE_DIR - directory where the header file can be found +# GMP_LIBRARY - Path where gmp library file can be found +# + +SET (GMP_FOUND FALSE) + +FIND_PATH(GMP_INCLUDE_DIR gmp.h + /usr/include + /usr/local/include +) + +set(GMP_LIBRARY_DIR "" CACHE PATH "GMP library folder") + +FIND_LIBRARY(GMP_LIBRARY NAMES gmp libgmp + HINTS "${GMP_LIBRARY_DIR}" +) + +mark_as_advanced(GMP_LIBRARY) + +IF( EXISTS "${GMP_INCLUDE_DIR}" AND EXISTS "${GMP_LIBRARY}" ) + set(GMP_FOUND TRUE) +else() + message(FATAL_ERROR "GMP not found") +endif() + diff --git a/Anima/cmake/FindMPFR.cmake b/Anima/cmake/FindMPFR.cmake new file mode 100644 index 000000000..e34b85cd1 --- /dev/null +++ b/Anima/cmake/FindMPFR.cmake @@ -0,0 +1,29 @@ +# Try to find the mpfr library +# Once done this will define +# +# MPFR_FOUND - system has mpfr and it can be used +# MPFR_INCLUDE_DIR - directory where the header file can be found +# MPFR_LIBRARY - Path where mpfr library file can be found +# + +SET (MPFR_FOUND FALSE) + +FIND_PATH(MPFR_INCLUDE_DIR mpfr.h + /usr/include + /usr/local/include +) + +set(MPFR_LIBRARY_DIR "" CACHE PATH "MPFR library folder") + +FIND_LIBRARY(MPFR_LIBRARY NAMES mpfr libmpfr + HINTS "${MPFR_LIBRARY_DIR}" +) + +mark_as_advanced(MPFR_LIBRARY) + +IF( EXISTS "${MPFR_INCLUDE_DIR}" AND EXISTS "${MPFR_LIBRARY}" ) + set(MPFR_FOUND TRUE) +else() + message(FATAL_ERROR "MPFR not found") +endif() + diff --git a/Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.hxx b/Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.hxx index 60739bb43..ab44a2a7c 100644 --- a/Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.hxx +++ b/Anima/diffusion/dti_estimator/animaDTIEstimationImageFilter.hxx @@ -214,14 +214,12 @@ DTIEstimationImageFilter for (unsigned int i = 0;i < numInputs;++i) { dwi[i] = inIterators[i].Get(); - lnDwi[i] = std::log(std::max(1.0e-6,dwi[i])); + lnDwi[i] = (dwi[i] <= 0) ? 0 : std::log(dwi[i]); } for (unsigned int i = 0;i < m_NumberOfComponents;++i) - { for (unsigned int j = 0;j < numInputs;++j) resVec[i] += m_InitialMatrixSolver(i+1,j) * lnDwi[j]; - } anima::GetTensorFromVectorRepresentation(resVec,data.workTensor); diff --git a/Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.hxx b/Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.hxx index 7174fe8cb..332cf9fec 100644 --- a/Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.hxx +++ b/Anima/diffusion/mcm_estimator/animaMCMEstimatorImageFilter.hxx @@ -479,7 +479,7 @@ MCMEstimatorImageFilter // Load DWI for (unsigned int i = 0;i < m_NumberOfImages;++i) observedSignals[i] = inIterators[i].Get(); - + int moseValue = -1; bool estimateNonIsoCompartments = false; if (m_ExternalMoseVolume) @@ -818,7 +818,7 @@ MCMEstimatorImageFilter workVec = mcmUpdateValue->GetParametersAsVector(); for (unsigned int i = 0;i < dimension;++i) p[i] = workVec[i]; - + double costValue = this->PerformSingleOptimization(p,cost,lowerBounds,upperBounds); this->GetProfiledInformation(cost,mcmUpdateValue,b0Value,sigmaSqValue); @@ -832,7 +832,7 @@ MCMEstimatorImageFilter if (m_CompartmentType == Stick) return; - + // We're done with ball and stick, next up is ball and zeppelin // - First create model mcmCreator->SetCompartmentType(Zeppelin); diff --git a/Anima/filtering/CMakeLists.txt b/Anima/filtering/CMakeLists.txt index 9586ff701..c24c6e5ad 100644 --- a/Anima/filtering/CMakeLists.txt +++ b/Anima/filtering/CMakeLists.txt @@ -7,6 +7,11 @@ project(ANIMA-FILTERING) add_subdirectory(bias_correction) add_subdirectory(denoising) add_subdirectory(dti_tools) + +if (USE_GMP AND GMP_FOUND AND USE_MPFR AND MPFR_FOUND AND USE_FLINT AND FLINT_FOUND AND USE_ARB AND ARB_FOUND) + add_subdirectory(gaussianizer) +endif() + add_subdirectory(noise_generator) add_subdirectory(nyul_standardization) add_subdirectory(regularization) diff --git a/Anima/filtering/gaussianizer/CMakeLists.txt b/Anima/filtering/gaussianizer/CMakeLists.txt new file mode 100644 index 000000000..dbb2cfbde --- /dev/null +++ b/Anima/filtering/gaussianizer/CMakeLists.txt @@ -0,0 +1,37 @@ +if(BUILD_TOOLS) + +project(animaGaussianizer) + +## ############################################################################# +## 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} + AnimaArbSpecialFunctions + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/filtering/gaussianizer/animaGaussianizer.cxx b/Anima/filtering/gaussianizer/animaGaussianizer.cxx new file mode 100644 index 000000000..abf107393 --- /dev/null +++ b/Anima/filtering/gaussianizer/animaGaussianizer.cxx @@ -0,0 +1,250 @@ +#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; +} diff --git a/Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.h b/Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.h new file mode 100644 index 000000000..24c5f4d7b --- /dev/null +++ b/Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.h @@ -0,0 +1,104 @@ +#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/gaussianizer/animaRiceToGaussianImageFilter.hxx b/Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.hxx new file mode 100644 index 000000000..bd673b186 --- /dev/null +++ b/Anima/filtering/gaussianizer/animaRiceToGaussianImageFilter.hxx @@ -0,0 +1,247 @@ +#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::GetKummerM(-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/math-tools/CMakeLists.txt b/Anima/math-tools/CMakeLists.txt index d7327b8a9..d33dfb981 100644 --- a/Anima/math-tools/CMakeLists.txt +++ b/Anima/math-tools/CMakeLists.txt @@ -4,6 +4,10 @@ project(ANIMA-MATHS) # Here go the add_subdirectories, no code should be at the root of the project ################################################################################ +if (USE_GMP AND GMP_FOUND AND USE_MPFR AND MPFR_FOUND AND USE_FLINT AND FLINT_FOUND AND USE_ARB AND ARB_FOUND) + add_subdirectory(arb_special_functions) +endif() + add_subdirectory(arithmetic) add_subdirectory(common_tools) diff --git a/Anima/math-tools/arb_special_functions/CMakeLists.txt b/Anima/math-tools/arb_special_functions/CMakeLists.txt new file mode 100644 index 000000000..354d60905 --- /dev/null +++ b/Anima/math-tools/arb_special_functions/CMakeLists.txt @@ -0,0 +1,43 @@ +project(AnimaArbSpecialFunctions) + +## ############################################################################# +## List Sources +## ############################################################################# + +list_source_files(${PROJECT_NAME} + ${CMAKE_CURRENT_SOURCE_DIR} + ) + + +## ############################################################################# +## add lib +## ############################################################################# + +add_library(${PROJECT_NAME} + ${${PROJECT_NAME}_CFILES} + ) + +## ############################################################################# +## Link +## ############################################################################# + +target_link_libraries(${PROJECT_NAME} + ${ARB_LIBRARY} + ) + + +################################################################################ +# Auto generate the export file for the libs +################################################################################ + +generate_export_header(${PROJECT_NAME} + STATIC_DEFINE ${PROJECT_NAME}_BUILT_AS_STATIC + EXPORT_FILE_NAME "${PROJECT_NAME}Export.h" + ) + + +## ############################################################################# +## install +## ############################################################################# + +set_lib_install_rules(${PROJECT_NAME}) diff --git a/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.cxx b/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.cxx new file mode 100644 index 000000000..5de843614 --- /dev/null +++ b/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.cxx @@ -0,0 +1,57 @@ +#include "animaARBBesselFunctions.h" +#include +#include +#include + +namespace anima +{ + +double GetScaledBesselI(const unsigned int N, const double x) +{ + arb_t inputBall, outputBall, dofBall; + arb_init(inputBall); + arb_init(outputBall); + arb_init(dofBall); + arb_set_d(inputBall, x); + arb_set_si(dofBall, N); + + unsigned int precision = 64; + arb_hypgeom_bessel_i_scaled(outputBall, dofBall, inputBall, precision); + + while (arb_can_round_arf(outputBall, 53, MPFR_RNDN) == 0) + { + precision *= 2; + arb_hypgeom_bessel_i_scaled(outputBall, dofBall, inputBall, precision); + } + + double resVal = arf_get_d(arb_midref(outputBall), MPFR_RNDN); + + arb_clear(inputBall); + arb_clear(outputBall); + arb_clear(dofBall); + + return resVal; +} + +double GetMarcumQ(const unsigned int M, const double a, const double b) +{ + if (M < 1) + { + std::cerr << "The order M should be greater or equal to 1." << std::endl; + exit(-1); + } + + MarcumQIntegrand integrand; + integrand.SetAValue(a); + integrand.SetBValue(b); + integrand.SetMValue(M); + + double resVal = b * b * boost::math::quadrature::gauss::integrate(integrand, 0.0, 1.0); + + if (M > 1) + resVal *= std::pow(b / a, M - 1.0); + + return 1.0 - resVal; +} + +} // end namespace anima diff --git a/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.h b/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.h new file mode 100644 index 000000000..106f074d5 --- /dev/null +++ b/Anima/math-tools/arb_special_functions/animaARBBesselFunctions.h @@ -0,0 +1,33 @@ +#pragma once + +#include + +#include "AnimaArbSpecialFunctionsExport.h" + +namespace anima +{ + +//! Computes exp(-x) I_N(x) using the ARB library +ANIMAARBSPECIALFUNCTIONS_EXPORT double GetScaledBesselI(const unsigned int N, const double x); + +class MarcumQIntegrand +{ +public: + void SetAValue(double val) {m_AValue = val;} + void SetBValue(double val) {m_BValue = val;} + void SetMValue(double val) {m_MValue = val;} + + double operator() (const double t) + { + double tmpVal = (m_MValue == 1) ? t : std::pow(t, (double)m_MValue); + return tmpVal * std::exp(-(m_BValue * t - m_AValue) * (m_BValue * t - m_AValue) / 2.0) * GetScaledBesselI(m_MValue - 1, m_AValue * m_BValue * t); + } + +private: + double m_AValue, m_BValue; + unsigned int m_MValue; +}; + +ANIMAARBSPECIALFUNCTIONS_EXPORT double GetMarcumQ(const unsigned int M, const double a, const double b); + +} // end of namespace anima diff --git a/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.cxx b/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.cxx new file mode 100644 index 000000000..13e8e9720 --- /dev/null +++ b/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.cxx @@ -0,0 +1,37 @@ +#include "animaARBKummerFunctions.h" +#include + +namespace anima +{ + +double GetKummerM(const double x, const double a, const double b) +{ + arb_t inputBall, outputBall, aBall, bBall; + arb_init(inputBall); + arb_init(outputBall); + arb_init(aBall); + arb_init(bBall); + arb_set_d(inputBall, x); + arb_set_d(aBall, a); + arb_set_d(bBall, b); + + unsigned int precision = 64; + arb_hypgeom_m(outputBall, aBall, bBall, inputBall, 0, precision); + + while (arb_can_round_arf(outputBall, 53, MPFR_RNDN) == 0) + { + precision *= 2; + arb_hypgeom_m(outputBall, aBall, bBall, inputBall, 0, precision); + } + + double resVal = arf_get_d(arb_midref(outputBall), MPFR_RNDN); + + arb_clear(inputBall); + arb_clear(outputBall); + arb_clear(aBall); + arb_clear(bBall); + + return resVal; +} + +} // end of namespace anima diff --git a/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.h b/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.h new file mode 100644 index 000000000..5aaf4e7e5 --- /dev/null +++ b/Anima/math-tools/arb_special_functions/animaARBKummerFunctions.h @@ -0,0 +1,11 @@ +#pragma once + +#include "AnimaArbSpecialFunctionsExport.h" + +namespace anima +{ + +//! Computes Kummer's M functions (i.e. confluent hypergeometric function 1F1) using the ARB library +ANIMAARBSPECIALFUNCTIONS_EXPORT double GetKummerM(const double x, const double a, const double b); + +} // end namespace anima diff --git a/Anima/math-tools/common_tools/convert_shape/CMakeLists.txt b/Anima/math-tools/common_tools/convert_shape/CMakeLists.txt index f4f1c133c..160ac3916 100644 --- a/Anima/math-tools/common_tools/convert_shape/CMakeLists.txt +++ b/Anima/math-tools/common_tools/convert_shape/CMakeLists.txt @@ -1,6 +1,6 @@ if(BUILD_TOOLS AND USE_VTK AND VTK_FOUND) -project(animaConvertShapes) +project(animaConvertShape) ## ############################################################################# ## List Sources diff --git a/Anima/math-tools/common_tools/convert_shape/animaConvertShape.cxx b/Anima/math-tools/common_tools/convert_shape/animaConvertShape.cxx index 30237039f..c772ef1d4 100644 --- a/Anima/math-tools/common_tools/convert_shape/animaConvertShape.cxx +++ b/Anima/math-tools/common_tools/convert_shape/animaConvertShape.cxx @@ -5,9 +5,12 @@ int main(int ac, const char** av) { - TCLAP::CmdLine cmd("animaConverthape can be used to convert a VTK, VTP, CSV or FDS (medInria fiber format) file into another format. It does not check before writing fds " - "that the file is actually made only of fibers.\n" - "INRIA / IRISA - VisAGeS Team",' ',ANIMA_VERSION); + TCLAP::CmdLine cmd("animaConvertShape can be used to convert shape data " + "back and forth between the following file formats: " + "VTK, VTP, CSV or FDS (medInria fiber format). It " + "does not check before writing fds that the file is " + "actually made only of fibers.\n" + "INRIA / IRISA - VisAGeS Team (and MOX - Politecnico di Milano for CSV support)",' ',ANIMA_VERSION); TCLAP::ValueArg inputArg("i","in","input data filename",true,"","input data",cmd); TCLAP::ValueArg outputArg("o","out","output data filename",true,"","output data",cmd); diff --git a/Anima/math-tools/data_io/animaGradientFileReader.hxx b/Anima/math-tools/data_io/animaGradientFileReader.hxx index d217dca5a..cd637e664 100644 --- a/Anima/math-tools/data_io/animaGradientFileReader.hxx +++ b/Anima/math-tools/data_io/animaGradientFileReader.hxx @@ -91,6 +91,7 @@ Update() startIndex = 0; std::string testExt = m_GradientFileName.substr(startIndex); + if (testExt.find("bvec") != std::string::npos) this->ReadGradientsFromBVecFile(); else @@ -209,7 +210,7 @@ ReadGradientsFromTextFile() gradTmp[0] = a0; gradTmp[1] = a1; gradTmp[2] = a2; - + m_Gradients.push_back(gradTmp); } @@ -240,7 +241,8 @@ ReadBValues() startIndex = 0; std::string testExt = m_BValueBaseString.substr(startIndex); - if (testExt.find("bval") != testExt.size()) + + if (testExt.find("bval") != std::string::npos) { std::string workStr; do { diff --git a/Anima/math-tools/matrix_operations/animaVectorOperations.h b/Anima/math-tools/matrix_operations/animaVectorOperations.h index cb85f5c76..6341ecb88 100644 --- a/Anima/math-tools/matrix_operations/animaVectorOperations.h +++ b/Anima/math-tools/matrix_operations/animaVectorOperations.h @@ -10,6 +10,22 @@ namespace anima { +//***** GetMedian() *****// +// Main +template double GetMedian(const VectorType &data, const unsigned int NDimension); +// For itkVector +template double GetMedian(const itk::Vector &data); +// For itkVariableLengthVector +template double GetMedian(const itk::VariableLengthVector &data); +// For itkPoint +template double GetMedian(const itk::Point &data); +// For vnl_vector +template double GetMedian(const vnl_vector &data); +// For vnl_vector_fixed +template double GetMedian(const vnl_vector_fixed &data); +// For std::vector +template double GetMedian(const std::vector &data); + /******* Main function ComputeEuclideanDistance *******/ // Main template double ComputeEuclideanDistance(const VectorType &x1, const VectorType &x2, const unsigned int NDimension); diff --git a/Anima/math-tools/matrix_operations/animaVectorOperations.hxx b/Anima/math-tools/matrix_operations/animaVectorOperations.hxx index f95b2a5a1..de2802492 100644 --- a/Anima/math-tools/matrix_operations/animaVectorOperations.hxx +++ b/Anima/math-tools/matrix_operations/animaVectorOperations.hxx @@ -9,6 +9,66 @@ namespace anima { +/******* Main function GetMedian *******/ +// Main +template double GetMedian(const VectorType &data, const unsigned int NDimension) +{ + std::vector array(NDimension, 0); + + for (unsigned int i = 0;i < NDimension;++i) + array[i] = data[i]; + + std::nth_element(array.begin(), array.begin() + NDimension / 2, array.end()); + + double median = array[NDimension / 2]; + + if (NDimension % 2 == 0) + { + median += *std::max_element(array.begin(), array.begin() + NDimension / 2); + median /= 2.0; + } + + return median; +} +// For itkVector +template double GetMedian(const itk::Vector &data) +{ + return GetMedian(data, NDimension); +} +// For itkVariableLengthVector +template double GetMedian(const itk::VariableLengthVector &data) +{ + unsigned int NDimension = data.GetSize(); + return GetMedian(data, NDimension); +} +// For itkPoint +template double GetMedian(const itk::Point &data) +{ + return GetMedian(data, NDimension); +} +// For vnl_vector +template double GetMedian(const vnl_vector &data) +{ + unsigned int NDimension = data.size(); + return GetMedian(data, NDimension); +} +// For vnl_vector_fixed +template +double GetMedian(const vnl_vector_fixed &data) +{ + return GetMedian(data, NDimension); +} +// For std::vector +template double GetMedian(const std::vector &data) +{ + unsigned int NDimension = data.size(); + return GetMedian(data, NDimension); +} +/******************************************************/ + + + + /******* Main function ComputeEuclideanDistance *******/ // Main template double ComputeEuclideanDistance(const VectorType &x1, const VectorType &x2, const unsigned int NDimension) diff --git a/Anima/math-tools/special_functions/animaBesselFunctions.cxx b/Anima/math-tools/special_functions/animaBesselFunctions.cxx index 44cef8625..86421fdd6 100644 --- a/Anima/math-tools/special_functions/animaBesselFunctions.cxx +++ b/Anima/math-tools/special_functions/animaBesselFunctions.cxx @@ -9,19 +9,24 @@ namespace anima { +double scaled_bessel_i(unsigned int N, double x) +{ + double logValue = log_bessel_i(N, x) - x; + return std::exp(logValue); +} + double log_bessel_i(unsigned int N, double x) { if (x < std::sqrt((double)N+1) / 100) { - // tgamma(N) = factorial(N-1) if (N == 0) return -std::log(std::tgamma(N+1)); else return -std::log(std::tgamma(N+1)) + N * std::log(x / 2.0); } - if (x <= std::max(9.23 * N + 15.934, 11.26 * N - 236.21)) // before was 600; now garantees an absolute error of maximum 1e-4 between true function and approximation - return std::log(boost::math::cyl_bessel_i(N,x)); + if (x <= 600) + return std::log(boost::math::cyl_bessel_i(N, x)); // // Too big value, switching to approximation // // Works less and less when orders go up but above that barrier we get non valid values diff --git a/Anima/math-tools/special_functions/animaBesselFunctions.h b/Anima/math-tools/special_functions/animaBesselFunctions.h index 93a5d604b..deb87953b 100644 --- a/Anima/math-tools/special_functions/animaBesselFunctions.h +++ b/Anima/math-tools/special_functions/animaBesselFunctions.h @@ -5,6 +5,9 @@ namespace anima { +//! Computes exp(-x) I_N(x) +ANIMASPECIALFUNCTIONS_EXPORT double scaled_bessel_i(unsigned int N, double x); + //! Computes a lower bound of the modified Bessel function of the first kind: I_{N} (N >= 0) ANIMASPECIALFUNCTIONS_EXPORT double bessel_i_lower_bound(unsigned int N, double x); diff --git a/Anima/math-tools/special_functions/animaKummerFunctions.cxx b/Anima/math-tools/special_functions/animaKummerFunctions.cxx index 28492cd61..bb5ffefa4 100644 --- a/Anima/math-tools/special_functions/animaKummerFunctions.cxx +++ b/Anima/math-tools/special_functions/animaKummerFunctions.cxx @@ -1,5 +1,6 @@ #include #include +#include namespace anima { diff --git a/Anima/math-tools/statistical_distributions/animaChiDistribution.h b/Anima/math-tools/statistical_distributions/animaChiDistribution.h new file mode 100644 index 000000000..37a4e53d5 --- /dev/null +++ b/Anima/math-tools/statistical_distributions/animaChiDistribution.h @@ -0,0 +1,28 @@ +#pragma once + +#include + +namespace anima +{ + +//! Implementation of Koay and Basser, Analytically exact correction scheme for signal extraction from noisy magnitude MR signals, Journal of Magnetic Resonance (2006). +double LowerBound(const unsigned int N, const double k0ValueSq, const double betaNSq); +double XiFunction(const double thetaSq, const unsigned int N, const double k1Value, const double betaNSq); +double GFunction(const double thetaSq, const double rSq, const unsigned int N, const double k1Value, const double betaNSq); +double KFunction(const double theta, const double r, const unsigned int N, const double betaNSq, double &k1Value); +double FixedPointFinder(const double r, const unsigned int N, double &k1Value, const unsigned int maximumNumberOfIterations = 500, const double epsilon = 1.0e-9); + +//! Implementation of Koay, Ozarslan and Basser, A signal transformational framework for breaking the noise floor and its applications in MRI, Journal of Magnetic Resonance (2009). +double GFunction2(const double etaSq, const double mSq, const double sigmaSq, const unsigned int N, const double k1Value, const double betaNSq); +double KFunction2(const double eta, const double m, const double sigma, const unsigned int N, const double betaNSq, double &k1Value); +double FixedPointFinder2(const double m, const double sigma, const unsigned int N, double &k1Value, const unsigned int maximumNumberOfIterations = 500, const double epsilon = 1.0e-9); + +//! Implementation of both Koay technique for Rice parameter estimation +void GetRiceParameters(const std::vector &samples, const std::vector &weights, double &location, double &scale); + +//! In-house function for evaluating the cumulative distribution function of a Rice distribution based on Gaussian quadrature integral approximation of the Marcum Q function. +double GetRiceCDF(const double x, const double location, const double scale); + +} // end of namespace + +#include "animaChiDistribution.hxx" diff --git a/Anima/math-tools/statistical_distributions/animaChiDistribution.hxx b/Anima/math-tools/statistical_distributions/animaChiDistribution.hxx new file mode 100644 index 000000000..194b658f0 --- /dev/null +++ b/Anima/math-tools/statistical_distributions/animaChiDistribution.hxx @@ -0,0 +1,183 @@ +#pragma once + +#include "animaChiDistribution.h" +#include +#include + +namespace anima +{ + +inline double factorial(unsigned int x) +{ + return (x <= 1) ? 1 : x * factorial(x - 1); +} + +inline double double_factorial(unsigned int x) +{ + return (x <= 1) ? 1 : x * double_factorial(x - 2); +} + +double LowerBound(const unsigned int N, const double k0ValueSq, const double betaNSq) +{ + double xiValue = 2.0 * N - betaNSq * k0ValueSq; + double insideValue = 2.0 * N / xiValue - 1.0; + return std::sqrt(std::max(insideValue, 0.0)); +} + +double XiFunction(const double thetaSq, const unsigned int N, const double k1Value, const double betaNSq) +{ + double resVal = 2.0 * N + thetaSq - betaNSq * k1Value * k1Value; + + if (resVal > 1.0) + { + std::cerr << "The xi function cannot take values greater than 1." << std::endl; + std::cout << thetaSq << " " << k1Value << " " << anima::GetKummerM(-thetaSq / 2.0, -0.5, N) << std::endl; + exit(-1); + } + + return resVal; +} + +double GFunction(const double thetaSq, const double rSq, const unsigned int N, const double k1Value, const double betaNSq) +{ + double insideValue = anima::XiFunction(thetaSq, N, k1Value, betaNSq) * (1.0 + rSq) - 2 * N; + return std::sqrt(std::max(insideValue, 0.0)); +} + +double KFunction(const double theta, const double r, const unsigned int N, const double betaNSq, double &k1Value) +{ + double thetaSq = theta * theta; + double rSq = r * r; + k1Value = anima::GetKummerM(-thetaSq / 2.0, -0.5, N); + double k2Value = anima::GetKummerM(-thetaSq / 2.0, 0.5, N + 1); + double gValue = anima::GFunction(thetaSq, rSq, N, k1Value, betaNSq); + double num = gValue * (gValue - theta); + double denom = theta * (1.0 + rSq) * (1.0 - betaNSq / (2.0 * N) * k1Value * k2Value) - gValue; + return theta - num / denom; +} + +double FixedPointFinder(const double r, const unsigned int N, double &k1Value, const unsigned int maximumNumberOfIterations, const double epsilon) +{ + double k0Value = anima::GetKummerM(0.0, -0.5, N); + double k0ValueSq = k0Value * k0Value; + double betaN = std::sqrt(M_PI / 2.0) * anima::double_factorial(2 * N - 1) / (std::pow(2, N - 1) * anima::factorial(N - 1)); + double betaNSq = betaN * betaN; + + double lb = anima::LowerBound(N, k0ValueSq, betaNSq); + if (r <= lb) + return 0.0; + + unsigned int counter = maximumNumberOfIterations; + + double t0 = r - lb; + double t1 = anima::KFunction(t0, r, N, betaNSq, k1Value); + + while (std::abs(t0 - t1) > epsilon) + { + t0 = t1; + t1 = anima::KFunction(t0, r, N, betaNSq, k1Value); + --counter; + + if (counter == 0) + break; + } + + return t1; +} + +double GFunction2(const double etaSq, const double mSq, const double sigmaSq, const unsigned int N, const double k1Value, const double betaNSq) +{ + double insideValue = mSq + (anima::XiFunction(etaSq / sigmaSq, N, k1Value, betaNSq) - 2 * N) * sigmaSq; + return std::sqrt(std::max(insideValue, 0.0)); +} + +double KFunction2(const double eta, const double m, const double sigma, const unsigned int N, const double betaNSq, double &k1Value) +{ + double mSq = m * m; + double etaSq = eta * eta; + double sigmaSq = sigma * sigma; + double thetaSq = etaSq / sigmaSq; + k1Value = anima::GetKummerM(-thetaSq / 2.0, -0.5, N); + double k2Value = anima::GetKummerM(-thetaSq / 2.0, 0.5, N + 1); + double gValue = anima::GFunction2(etaSq, mSq, sigmaSq, N, k1Value, betaNSq); + double num = gValue * (gValue - eta); + double denom = eta * (1.0 - betaNSq / (2.0 * N) * k1Value * k2Value) - gValue; + return eta - num / denom; +} + +double FixedPointFinder2(const double m, const double sigma, const unsigned int N, double &k1Value, const unsigned int maximumNumberOfIterations, const double epsilon) +{ + unsigned int counter = maximumNumberOfIterations; + double betaN = std::sqrt(M_PI / 2.0) * anima::double_factorial(2 * N - 1) / (std::pow(2, N - 1) * anima::factorial(N - 1)); + double betaNSq = betaN * betaN; + double delta = betaN * sigma - m; + + if (delta == 0) + return 0; + + double mCorrected = (delta > 0) ? betaN * sigma + delta : m; + + double t0 = mCorrected; + double t1 = anima::KFunction2(t0, mCorrected, sigma, N, betaNSq, k1Value); + + while (std::abs(t0 - t1) > epsilon) + { + t0 = t1; + t1 = anima::KFunction2(t0, mCorrected, sigma, N, betaNSq, k1Value); + --counter; + + if (counter == 0) + break; + } + + return (delta > 0) ? -t1 : t1; +} + +void GetRiceParameters(const std::vector &samples, const std::vector &weights, double &location, double &scale) +{ + unsigned int sampleSize = samples.size(); + double k1Value = 0; + double meanValue = 0; + double sumWeights = 0, sumSqWeights = 0; + + for (unsigned int i = 0;i < sampleSize;++i) + { + double weightValue = weights[i]; + meanValue += weightValue * samples[i]; + sumWeights += weightValue; + sumSqWeights += weightValue * weightValue; + } + + meanValue /= sumWeights; + + if (scale == 0) + { + double sigmaValue = 0; + + for (unsigned int i = 0;i < sampleSize;++i) + sigmaValue += weights[i] * (samples[i] - meanValue) * (samples[i] - meanValue); + + sigmaValue /= (sumWeights * sumWeights - sumSqWeights); + sigmaValue *= sumWeights; + sigmaValue = std::sqrt(sigmaValue); + + double rValue = meanValue / sigmaValue; + double thetaValue = FixedPointFinder(rValue, 1, k1Value); + k1Value = anima::GetKummerM(-thetaValue * thetaValue / 2.0, -0.5, 1); + + scale = sigmaValue / std::sqrt(anima::XiFunction(thetaValue * thetaValue, 1, k1Value, M_PI / 2.0)); + location = thetaValue * scale; + return; + } + + location = FixedPointFinder2(meanValue, scale, 1, k1Value); +} + +double GetRiceCDF(const double x, const double location, const double scale) +{ + double a = location / scale; + double b = x / scale; + return 1.0 - anima::GetMarcumQ(1, a, b); +} + +} // end of namespace anima diff --git a/CMakeLists.txt b/CMakeLists.txt index bed194bd0..d98d5a467 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -19,6 +19,7 @@ option(BUILD_SHARED_LIBS "Build shared libraries" ON) option(USE_VTK "Build VTK dependencies" ON) option(USE_RPI "Build RPI dependencies" ON) option(USE_NLOPT "Build NLOPT dependencies" ON) +option(USE_ARB "Build ARB dependencies" ON) option(BUILD_ANIMA_TOOLS "Build ANIMA tools" ON) option(BUILD_ANIMA_TESTING "Build ANIMA testing executables" OFF) option(BUILD_ANIMA_DOCUMENTATION "Build ANIMA doxygen" OFF) @@ -95,6 +96,34 @@ if (USE_NLOPT) endif() endif() +# GMP +if (USE_GMP) + find_package(GMP REQUIRED) + set(GMP_SRC_DIR ${GMP_INCLUDE_DIR}) + set(GMP_BUILD_DIR ${GMP_LIBRARY_DIR}) +endif() + +# MPFR +if (USE_MPFR) + find_package(MPFR REQUIRED) + set(MPFR_SRC_DIR ${MPFR_INCLUDE_DIR}) + set(MPFR_BUILD_DIR ${MPFR_LIBRARY_DIR}) +endif() + +# FLINT +if (USE_FLINT) + find_package(FLINT REQUIRED) + set(FLINT_SRC_DIR ${FLINT_INCLUDE_DIR}) + set(FLINT_BUILD_DIR ${FLINT_LIBRARY_DIR}) +endif() + +# ARB +if (USE_ARB) + find_package(ARB REQUIRED) + set(ARB_SRC_DIR ${ARB_INCLUDE_DIR}) + set(ARB_BUILD_DIR ${ARB_LIBRARY_DIR}) +endif() + # VTK if (USE_VTK) option(USE_SYSTEM_VTK "Use system installed VTK" OFF)