A dedicated library for statistical distributions

Anima/diffusion/mcm/CMakeLists.txt

@@ -24,6 +24,7 @@ add_library(${PROJECT_NAME}
 target_link_libraries(${PROJECT_NAME}
   AnimaMCMBase
   AnimaSpecialFunctions
+  AnimaStatisticalDistributions
   ITKOptimizers
   ITKCommon
   ${TinyXML2_LIBRARY}

Anima/diffusion/mcm/animaNODDICompartment.cxx

@@ -494,7 +494,8 @@ void NODDICompartment::UpdateKappaValues()
     double dawsonValue = anima::EvaluateDawsonIntegral(std::sqrt(kappa), true);
     m_Tau1 = (1.0 / dawsonValue - 1.0) / (2.0 * kappa);
     m_Tau1Deriv = (1.0 - (1.0 - dawsonValue * (2.0 * kappa - 1.0)) / (2.0 * dawsonValue * dawsonValue)) / (2.0 * kappa * kappa);
-    anima::GetStandardWatsonSHCoefficients(kappa,m_WatsonSHCoefficients,m_WatsonSHCoefficientDerivatives);
+    m_WatsonDistribution.SetConcentrationParameter(kappa);
+    m_WatsonDistribution.GetStandardWatsonSHCoefficients(m_WatsonSHCoefficients, m_WatsonSHCoefficientDerivatives);
 
     m_ModifiedConcentration = false;
 }

Anima/diffusion/mcm/animaNODDICompartment.h

@@ -2,6 +2,7 @@
 
 #include <animaBaseCompartment.h>
 #include <AnimaMCMExport.h>
+#include <animaWatsonDistribution.h>
 
 namespace anima
 {
@@ -84,6 +85,12 @@ class ANIMAMCM_EXPORT NODDICompartment : public BaseCompartment
 
         m_CurrentBValue = -1.0;
         m_CurrentGradient.fill(0.0);
+
+        itk::Vector<double,3> meanAxis;
+        meanAxis[0] = 0.0;
+        meanAxis[1] = 0.0;
+        meanAxis[2] = 1.0;
+        m_WatsonDistribution.SetMeanAxis(meanAxis);
     }
 
     virtual ~NODDICompartment() {}
@@ -112,6 +119,8 @@ class ANIMAMCM_EXPORT NODDICompartment : public BaseCompartment
     double m_Tau1, m_Tau1Deriv;
     double m_ExtraAxonalSignal, m_IntraAxonalSignal;
     double m_IntraAngleDerivative, m_IntraKappaDerivative, m_IntraAxialDerivative;
+
+    anima::WatsonDistribution m_WatsonDistribution;
 };
 
 } //end namespace anima

Anima/diffusion/mcm_tools/CMakeLists.txt

@@ -2,6 +2,7 @@ add_subdirectory(dwi_simulation_from_mcm)
 add_subdirectory(generate_isoradius_ddi_surface)
 add_subdirectory(get_scalar_map_from_ddi)
 add_subdirectory(mcm_average_images)
+add_subdirectory(mcm_merge_anisotropic_weights)
 add_subdirectory(mcm_merge_block_images)
 add_subdirectory(mcm_scalar_maps)
 add_subdirectory(mt_estimation_validation)

Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/CMakeLists.txt

@@ -0,0 +1,37 @@
+if(BUILD_TOOLS)
+
+project(animaMCMMergeAnisotropicWeights)
+
+## #############################################################################
+## 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}
+  AnimaMCM
+  )
+
+## #############################################################################
+## install
+## #############################################################################
+
+set_exe_install_rules(${PROJECT_NAME})
+
+endif()

Anima/diffusion/mcm_tools/mcm_merge_anisotropic_weights/animaMCMMergeAnisotropicWeights.cxx

@@ -0,0 +1,106 @@
+#include <animaMCMFileReader.h>
+#include <animaMCMImage.h>
+#include <animaMultiCompartmentModel.h>
+#include <animaReadWriteFunctions.h>
+#include <itkImageRegionConstIterator.h>
+#include <itkImageRegionIterator.h>
+
+#include <tclap/CmdLine.h>
+
+int main(int argc, char **argv)
+{
+    TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS/Empenn Team", ' ', ANIMA_VERSION);
+
+    TCLAP::ValueArg<std::string>   inArg("i", "inputfile", "Input MCM image", true, "", "input mcm image", cmd);
+    TCLAP::ValueArg<std::string>  outArg("o", "outputfile", "Output image with combined anisotropic weights", true, "", "output conbined weight image", cmd);
+
+    try
+    {
+        cmd.parse(argc,argv);
+    }
+    catch (TCLAP::ArgException& e)
+    {
+        std::cerr << "Error: " << e.error() << "for argument " << e.argId() << std::endl;
+        return EXIT_FAILURE;
+    }
+
+    using InputImageType = anima::MCMImage <double, 3>;
+    using OutputImageType = itk::VectorImage <double, 3>;
+    using InputPixelType = InputImageType::PixelType;
+    using OutputPixelType = OutputImageType::PixelType;
+    using MCModelType = anima::MultiCompartmentModel;
+    using MCModelPointer = MCModelType::Pointer;
+
+    // Read input sample
+    anima::MCMFileReader <double, 3> mcmReader;
+    mcmReader.SetFileName(inArg.getValue());
+    mcmReader.Update();
+    InputImageType::Pointer mcmImage = mcmReader.GetModelVectorImage();
+    MCModelPointer mcmPtr = mcmImage->GetDescriptionModel()->Clone();
+    InputPixelType mcmValue;
+    unsigned int numCompartments = mcmPtr->GetNumberOfCompartments();
+    unsigned int numIsoCompartments = mcmPtr->GetNumberOfIsotropicCompartments();
+
+    using InputImageIteratorType = itk::ImageRegionConstIterator <InputImageType>;
+    using OutputImageIteratorType = itk::ImageRegionIterator <OutputImageType>;
+    InputImageIteratorType inItr(mcmImage, mcmImage->GetLargestPossibleRegion());
+
+    // Initialize output image
+    unsigned int nbOutputComponents = numIsoCompartments + 1;
+    OutputImageType::Pointer outputImage = OutputImageType::New();
+    outputImage->SetRegions(mcmImage->GetLargestPossibleRegion());
+    outputImage->CopyInformation(mcmImage);
+    outputImage->SetVectorLength(nbOutputComponents);
+    outputImage->Allocate();
+
+    OutputPixelType outputValue(nbOutputComponents);
+    outputValue.Fill(0.0);
+    outputImage->FillBuffer(outputValue);
+
+    std::cout << "- Number of compartments in input image: " << numCompartments << std::endl;
+    std::cout << "- Number of compartments in output image: " << nbOutputComponents << std::endl;
+
+    OutputImageIteratorType outItr(outputImage, outputImage->GetLargestPossibleRegion());
+
+	while (!outItr.IsAtEnd())
+	{
+        mcmValue = inItr.Get();
+
+        bool backgroundVoxel = true;
+        for (unsigned int i = 0;i < mcmValue.GetNumberOfElements();++i)
+        {
+            if (mcmValue[i] != 0.0)
+            {
+                backgroundVoxel = false;
+                break;
+            }
+        }
+
+        if (backgroundVoxel)
+        {
+            ++inItr;
+            ++outItr;
+            continue;
+        }
+
+        mcmPtr->SetModelVector(mcmValue);
+
+        for (unsigned int i = 0;i < numIsoCompartments;++i)
+            outputValue[i] = mcmPtr->GetCompartmentWeight(i);
+
+        double anisoWeight = 0.0;
+        for (unsigned int i = numIsoCompartments;i < numCompartments;++i)
+            anisoWeight += mcmPtr->GetCompartmentWeight(i);
+
+        outputValue[numIsoCompartments] = anisoWeight;
+
+        outItr.Set(outputValue);
+
+        ++inItr;
+        ++outItr;
+	}
+
+    anima::writeImage <OutputImageType> (outArg.getValue(), outputImage);
+
+	return EXIT_SUCCESS;
+}

Anima/diffusion/tractography/CMakeLists.txt

@@ -29,6 +29,7 @@ target_link_libraries(${PROJECT_NAME}
   AnimaDataIO
   AnimaOptimizers
   AnimaSHTools
+  AnimaStatisticalDistributions
   AnimaMCM
   AnimaMCMBase
   )

Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.h

@@ -1,5 +1,7 @@
 #pragma once
 
+#include <animaWatsonDistribution.h>
+
 #include <itkVectorImage.h>
 #include <itkImage.h>
 
@@ -221,6 +223,12 @@ class BaseProbabilisticTractographyImageFilter : public itk::ProcessObject
     //! Computes additional scalar maps that are model dependent to add to the output
     virtual void ComputeAdditionalScalarMaps() {}
 
+    //! Holds an object of class WatsonDistribution for sampling new directions
+    anima::WatsonDistribution m_WatsonDistribution;
+
+    //! Holds a sample of size 1 of directions
+    DirectionVectorType m_SampleOfDirections;
+
 private:
     ITK_DISALLOW_COPY_AND_ASSIGN(BaseProbabilisticTractographyImageFilter);
 

Anima/diffusion/tractography/animaBaseProbabilisticTractographyImageFilter.hxx

@@ -58,6 +58,8 @@ BaseProbabilisticTractographyImageFilter <TInputModelImageType>
 
     m_HighestProcessedSeed = 0;
     m_ProgressReport = 0;
+
+    m_SampleOfDirections.resize(1);
 }
 
 template <class TInputModelImageType>

Anima/diffusion/tractography/animaDTIProbabilisticTractographyImageFilter.cxx

@@ -4,7 +4,6 @@
 #include <animaVectorOperations.h>
 #include <animaDistributionSampling.h>
 #include <animaVMFDistribution.h>
-#include <animaWatsonDistribution.h>
 #include <animaLogarithmFunctions.h>
 #include <animaBaseTensorTools.h>
 
@@ -64,7 +63,10 @@ DTIProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi
 //    else
 //        anima::SampleFromVMFDistribution(concentrationParameter,sampling_direction,resVec,random_generator);
 
-    anima::SampleFromWatsonDistribution(concentrationParameter,sampling_direction,resVec,3,random_generator);
+    m_WatsonDistribution.SetMeanAxis(sampling_direction);
+    m_WatsonDistribution.SetConcentrationParameter(concentrationParameter);
+    m_WatsonDistribution.Random(m_SampleOfDirections, random_generator);
+    resVec = m_SampleOfDirections[0];
 
     if (is2d)
     {
@@ -75,10 +77,14 @@ DTIProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi
     if (LC > m_ThresholdForProlateTensor)
     {
 //        log_prior = anima::safe_log(anima::ComputeVMFPdf(resVec, oldDirection, this->GetKappaOfPriorDistribution()));
-        log_prior = anima::safe_log(anima::EvaluateWatsonPDF(resVec, oldDirection, this->GetKappaOfPriorDistribution()));
+        m_WatsonDistribution.SetMeanAxis(oldDirection);
+        m_WatsonDistribution.SetConcentrationParameter(this->GetKappaOfPriorDistribution());
+        log_prior = m_WatsonDistribution.GetLogDensity(resVec);
 
 //        log_proposal = anima::safe_log(anima::ComputeVMFPdf(resVec, sampling_direction, concentrationParameter));
-        log_proposal = anima::safe_log(anima::EvaluateWatsonPDF(resVec, sampling_direction, concentrationParameter));
+        m_WatsonDistribution.SetMeanAxis(sampling_direction);
+        m_WatsonDistribution.SetConcentrationParameter(concentrationParameter);
+        log_proposal = m_WatsonDistribution.GetLogDensity(resVec);
     }
 
     if (anima::ComputeScalarProduct(oldDirection, resVec) < 0)
@@ -199,19 +205,23 @@ double DTIProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0
     if (noiseValue > 0)
         concentrationParameter = b0Value / std::sqrt(noiseValue);
 
+    m_WatsonDistribution.SetMeanAxis(newDirection);
+
     if (LC > m_ThresholdForProlateTensor)
     {
         Vector3DType dtiPrincipalDirection(0.0);
 
         this->GetDTIPrincipalDirection(modelValue, dtiPrincipalDirection, is2d);
-        logLikelihood = std::log(anima::EvaluateWatsonPDF(dtiPrincipalDirection, newDirection, concentrationParameter));
+        m_WatsonDistribution.SetConcentrationParameter(concentrationParameter);
+        logLikelihood = m_WatsonDistribution.GetLogDensity(dtiPrincipalDirection);
     }
     else
     {
         Vector3DType dtiMinorDirection(0.0);
 
         this->GetDTIMinorDirection(modelValue, dtiMinorDirection);
-        logLikelihood = std::log(anima::EvaluateWatsonPDF(dtiMinorDirection, newDirection, - concentrationParameter));
+        m_WatsonDistribution.SetConcentrationParameter(-concentrationParameter);
+        logLikelihood = m_WatsonDistribution.GetLogDensity(dtiMinorDirection);
     }
 
     double resVal = log_prior + logLikelihood - log_proposal;

Anima/diffusion/tractography/animaMCMProbabilisticTractographyImageFilter.cxx

@@ -122,7 +122,10 @@ MCMProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi
     //    else
     //        anima::SampleFromVMFDistribution(chosenKappa,sampling_direction,resVec,random_generator);
 
-    anima::SampleFromWatsonDistribution(chosenKappa,sampling_direction,resVec,3,random_generator);
+    m_WatsonDistribution.SetMeanAxis(sampling_direction);
+    m_WatsonDistribution.SetConcentrationParameter(chosenKappa);
+    m_WatsonDistribution.Random(m_SampleOfDirections, random_generator);
+    resVec = m_SampleOfDirections[0];
 
     if (is2d)
     {
@@ -133,14 +136,18 @@ MCMProbabilisticTractographyImageFilter::ProposeNewDirection(Vector3DType &oldDi
     if (effectiveNumDirs > 0)
     {
         //        log_prior = anima::safe_log(anima::ComputeVMFPdf(resVec, oldDirection, this->GetKappaOfPriorDistribution()));
-        log_prior = anima::safe_log(anima::EvaluateWatsonPDF(resVec, oldDirection, this->GetKappaOfPriorDistribution()));
+        m_WatsonDistribution.SetMeanAxis(oldDirection);
+        m_WatsonDistribution.SetConcentrationParameter(this->GetKappaOfPriorDistribution());
+        log_prior = m_WatsonDistribution.GetLogDensity(resVec);
 
         log_proposal = 0;
         double sumWeights = 0;
         for (unsigned int i = 0;i < effectiveNumDirs;++i)
         {
             //            log_proposal += mixtureWeights[i] * anima::ComputeVMFPdf(resVec, maximaMCM[i], kappaValues[i]);
-            log_proposal += mixtureWeights[i] * anima::EvaluateWatsonPDF(resVec, maximaMCM[i], kappaValues[i]);
+            m_WatsonDistribution.SetMeanAxis(maximaMCM[i]);
+            m_WatsonDistribution.SetConcentrationParameter(kappaValues[i]);
+            log_proposal += mixtureWeights[i] * m_WatsonDistribution.GetDensity(resVec);
             sumWeights += mixtureWeights[i];
         }
 
@@ -290,6 +297,9 @@ double MCMProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0
 
     double concentrationParameter = b0Value / std::sqrt(noiseValue);
 
+    m_WatsonDistribution.SetMeanAxis(newDirection);
+    m_WatsonDistribution.SetConcentrationParameter(concentrationParameter);
+
     bool oneTested = false;
     for (unsigned int i = workModel->GetNumberOfIsotropicCompartments();i < workModel->GetNumberOfCompartments();++i)
     {
@@ -306,7 +316,7 @@ double MCMProbabilisticTractographyImageFilter::ComputeLogWeightUpdate(double b0
             anima::Normalize(tmpVec,tmpVec);
         }
 
-        double tmpVal = std::log(anima::EvaluateWatsonPDF(tmpVec, newDirection, concentrationParameter));
+        double tmpVal = m_WatsonDistribution.GetLogDensity(tmpVec);
         if ((tmpVal > logLikelihood)||(!oneTested))
         {
             logLikelihood = tmpVal;