Skip to content

Commit

Permalink
Merge pull request #104 from MariePoirierGit/dirichlet_distribution
Browse files Browse the repository at this point in the history
A dedicated library for statistical distributions
  • Loading branch information
astamm authored Oct 27, 2023
2 parents 9a055f1 + 9c134ba commit 815d8b8
Show file tree
Hide file tree
Showing 41 changed files with 1,623 additions and 486 deletions.
1 change: 1 addition & 0 deletions Anima/diffusion/mcm/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ add_library(${PROJECT_NAME}
target_link_libraries(${PROJECT_NAME}
AnimaMCMBase
AnimaSpecialFunctions
AnimaStatisticalDistributions
ITKOptimizers
ITKCommon
${TinyXML2_LIBRARY}
Expand Down
3 changes: 2 additions & 1 deletion Anima/diffusion/mcm/animaNODDICompartment.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
Expand Down
9 changes: 9 additions & 0 deletions Anima/diffusion/mcm/animaNODDICompartment.h
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include <animaBaseCompartment.h>
#include <AnimaMCMExport.h>
#include <animaWatsonDistribution.h>

namespace anima
{
Expand Down Expand Up @@ -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() {}
Expand Down Expand Up @@ -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
1 change: 1 addition & 0 deletions Anima/diffusion/mcm_tools/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Original file line number Diff line number Diff line change
@@ -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()
Original file line number Diff line number Diff line change
@@ -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;
}
1 change: 1 addition & 0 deletions Anima/diffusion/tractography/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,7 @@ target_link_libraries(${PROJECT_NAME}
AnimaDataIO
AnimaOptimizers
AnimaSHTools
AnimaStatisticalDistributions
AnimaMCM
AnimaMCMBase
)
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#pragma once

#include <animaWatsonDistribution.h>

#include <itkVectorImage.h>
#include <itkImage.h>

Expand Down Expand Up @@ -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);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -58,6 +58,8 @@ BaseProbabilisticTractographyImageFilter <TInputModelImageType>

m_HighestProcessedSeed = 0;
m_ProgressReport = 0;

m_SampleOfDirections.resize(1);
}

template <class TInputModelImageType>
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,6 @@
#include <animaVectorOperations.h>
#include <animaDistributionSampling.h>
#include <animaVMFDistribution.h>
#include <animaWatsonDistribution.h>
#include <animaLogarithmFunctions.h>
#include <animaBaseTensorTools.h>

Expand Down Expand Up @@ -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)
{
Expand All @@ -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)
Expand Down Expand Up @@ -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;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand All @@ -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];
}

Expand Down Expand Up @@ -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)
{
Expand All @@ -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;
Expand Down
Loading

0 comments on commit 815d8b8

Please sign in to comment.