From 17aaad4ee9db9885cede8fca5dfb3ef1fdf60042 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Fri, 27 Oct 2023 14:33:43 +0200 Subject: [PATCH 1/4] First draft of MCM orientation priors tool. --- Anima/diffusion/mcm_tools/CMakeLists.txt | 1 + .../mcm_orientation_priors/CMakeLists.txt | 37 +++ .../animaMCMOrientationPriors.cxx | 111 +++++++ .../animaMCMOrientationPriorsImageFilter.h | 83 +++++ .../animaMCMOrientationPriorsImageFilter.hxx | 308 ++++++++++++++++++ 5 files changed, 540 insertions(+) create mode 100644 Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt create mode 100644 Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx create mode 100644 Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h create mode 100644 Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx diff --git a/Anima/diffusion/mcm_tools/CMakeLists.txt b/Anima/diffusion/mcm_tools/CMakeLists.txt index 911579ba5..f565d95b4 100644 --- a/Anima/diffusion/mcm_tools/CMakeLists.txt +++ b/Anima/diffusion/mcm_tools/CMakeLists.txt @@ -3,6 +3,7 @@ add_subdirectory(generate_isoradius_ddi_surface) add_subdirectory(get_scalar_map_from_ddi) add_subdirectory(mcm_average_images) add_subdirectory(mcm_merge_block_images) +add_subdirectory(mcm_orientation_priors) add_subdirectory(mcm_scalar_maps) add_subdirectory(mt_estimation_validation) add_subdirectory(test_averaging) \ No newline at end of file diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt b/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt new file mode 100644 index 000000000..9bbb3ca82 --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt @@ -0,0 +1,37 @@ +if(BUILD_TOOLS) + +project(animaMCMOrientationPriors) + +## ############################################################################# +## 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} + ${TinyXML2_LIBRARY} + AnimaMCM + ) + +## ############################################################################# +## install +## ############################################################################# + +set_exe_install_rules(${PROJECT_NAME}) + +endif() diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx new file mode 100644 index 000000000..f44b7b2fe --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx @@ -0,0 +1,111 @@ +#include +#include +#include + +#include +#include + +int main(int argc, char **argv) +{ + TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS Team", ' ',ANIMA_VERSION); + + TCLAP::ValueArg inArg( + "i", "input-mcms", + "A text file specifying a list of MCM images in a common geometry.", + true, "", "input MCM images", cmd + ); + TCLAP::ValueArg outOrientationArg( + "o", "output-orientation", + "A string specifying the basename for the output vector images that will store priors on the orientations.", + true, "", "output orientation priors",cmd + ); + TCLAP::ValueArg outWeightsArg( + "w", "output-weights", + "A string specifying the filename for the output vector image that will store priors on the weights.", + true, "", "output weights priors", cmd + ); + + TCLAP::ValueArg maskArg( + "m", "input-masks", + "A text file specifying a list of mask images in the same common geometry as the input MCM images (default: none).", + false, "", "input mask images", cmd + ); + TCLAP::ValueArg nbThreadsArg( + "T", "nthreads", + "An integer value specifying the number of threads to run on (default: all cores).", + false, itk::MultiThreaderBase::GetGlobalDefaultNumberOfThreads(), "number of threads", 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 MainFilterType = anima::MCMOrientationPriorsImageFilter ; + MainFilterType::Pointer mainFilter = MainFilterType::New(); + + using MCMReaderType = anima::MCMFileReader ; + using MaskImageType = MainFilterType::MaskImageType; + using OutputImageType = MainFilterType::OutputImageType; + + // Load MCM images + std::ifstream inputFile(inArg.getValue().c_str()); + + if (!inputFile.is_open()) + { + std::cerr << "Please provide usable file with input MCMs" << std::endl; + return EXIT_FAILURE; + } + + unsigned int nbOfImages = 0; + + while (!inputFile.eof()) + { + char tmpStr[2048]; + inputFile.getline(tmpStr,2048); + + if (strcmp(tmpStr,"") == 0) + continue; + + std::cout << "Loading image " << nbOfImages << " : " << tmpStr << std::endl; + + MCMReaderType mcmReader; + mcmReader.SetFileName(tmpStr); + mcmReader.Update(); + + mainFilter->SetInput(nbOfImages, mcmReader.GetModelVectorImage()); + + nbOfImages++; + } + + std::ifstream masksIn; + if (maskArg.getValue() != "") + masksIn.open(maskArg.getValue()); + + if (masksIn.is_open()) + { + char tmpStr[2048]; + while (!masksIn.eof()) + { + masksIn.getline(tmpStr,2048); + + if (strcmp(tmpStr,"") == 0) + continue; + + mainFilter->AddMaskImage(anima::readImage (tmpStr)); + } + } + + mainFilter->SetNumberOfWorkUnits(nbThreadsArg.getValue()); + mainFilter->Update(); + + anima::writeImage (outOrientationArg.getValue(),mainFilter->GetOutput(0)); + anima::writeImage (outWeightsArg.getValue(),mainFilter->GetOutput(1)); + + return EXIT_SUCCESS; +} diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h new file mode 100644 index 000000000..313af24ac --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h @@ -0,0 +1,83 @@ +#pragma once + +#include +#include + +#include +#include + +namespace anima +{ + +template +class MCMOrientationPriorsImageFilter : +public itk::ImageToImageFilter< anima::MCMImage , itk::VectorImage > +{ +public: + /** Standard class type def */ + + using Self = MCMOrientationPriorsImageFilter; + using InputImageType = anima::MCMImage ; + using OutputImageType = itk::VectorImage ; + using Superclass = itk::ImageToImageFilter ; + using Pointer = itk::SmartPointer; + using ConstPointer = itk::SmartPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self) + + /** Run-time type information (and related methods) */ + itkTypeMacro(MCMOrientationPriorsImageFilter, ImageToImageFilter) + + using InputImagePointer = typename InputImageType::ConstPointer; + using OutputImagePointer = typename OutputImageType::Pointer; + + using InputRegionType = typename InputImageType::RegionType; + using InputIndexType = typename InputImageType::IndexType; + using PixelType = typename InputImageType::PixelType; + + using MaskImageType = itk::Image ; + using MaskImagePointer = MaskImageType::Pointer; + + // Multi-compartment models typedefs + using MCModelType = anima::MultiCompartmentModel; + using MCModelPointer = MCModelType::Pointer; + + void AddMaskImage(MaskImageType *maskImage) {m_MaskImages.push_back(maskImage);} + +protected: + MCMOrientationPriorsImageFilter() + { + m_ReferenceInputModels.clear(); + m_MaskImages.clear(); + m_NumberOfAnisotropicCompartments = 0; + } + + virtual ~MCMOrientationPriorsImageFilter() {} + + bool isZero(const itk::VariableLengthVector &value) const + { + for (unsigned int i = 0;i < value.GetNumberOfElements();++i) + { + if (value[i] != 0) + return false; + } + + return true; + } + + void GenerateOutputInformation() ITK_OVERRIDE; + void BeforeThreadedGenerateData() ITK_OVERRIDE; + void DynamicThreadedGenerateData(const InputRegionType ®ion) ITK_OVERRIDE; + +private: + ITK_DISALLOW_COPY_AND_ASSIGN(MCMOrientationPriorsImageFilter); + + std::vector m_ReferenceInputModels; + std::vector m_MaskImages; + unsigned int m_NumberOfAnisotropicCompartments; +}; + +} // end namespace anima + +#include "animaMCMOrientationPriorsImageFilter.hxx" diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx new file mode 100644 index 000000000..3a13ac97e --- /dev/null +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx @@ -0,0 +1,308 @@ +#pragma once +#include "animaMCMOrientationPriorsImageFilter.h" +#include + +#include +#include + +namespace anima +{ + +template +void +MCMOrientationPriorsImageFilter +::GenerateOutputInformation() +{ + //------------------------------- + // Creates the additional outputs + //------------------------------- + unsigned int prevNum = this->GetNumberOfOutputs(); + + InputImageType *input = const_cast (this->GetInput(0)); + MCModelPointer tmpMCM = input->GetDescriptionModel(); + unsigned int numIsoCompartments = tmpMCM->GetNumberOfIsotropicCompartments(); + unsigned int numCompartments = tmpMCM->GetNumberOfCompartments(); + m_NumberOfAnisotropicCompartments = numCompartments - numIsoCompartments; + unsigned int numOutputs = m_NumberOfAnisotropicCompartments + 1; + + this->SetNumberOfIndexedOutputs(numOutputs); + + for (unsigned int i = prevNum;i < numOutputs;++i) + this->SetNthOutput(i,this->MakeOutput(i).GetPointer()); + + //--------------------------------------------------- + // Call the superclass' implementation of this method + //--------------------------------------------------- + Superclass::GenerateOutputInformation(); + + //-------------------------------------------------- + // Set length of vectors in the output vector images + //-------------------------------------------------- + for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + { + OutputImageType *output = this->GetOutput(i); + output->SetVectorLength(4); // AST: now fixed for Watson + } + + OutputImageType *output = this->GetOutput(m_NumberOfAnisotropicCompartments); + output->SetVectorLength(m_NumberOfAnisotropicCompartments); +} + +template +void +MCMOrientationPriorsImageFilter +::BeforeThreadedGenerateData() +{ + this->Superclass::BeforeThreadedGenerateData(); + + m_ReferenceInputModels.resize(this->GetNumberOfIndexedInputs()); + + InputImageType *input = const_cast (this->GetInput(0)); + MCModelPointer tmpMCM = input->GetDescriptionModel(); + unsigned int numIsoCompartments = tmpMCM->GetNumberOfIsotropicCompartments(); + unsigned int numCompartments = tmpMCM->GetNumberOfCompartments(); + + for (unsigned int i = 0;i < this->GetNumberOfIndexedInputs();++i) + { + InputImageType *input = const_cast (this->GetInput(i)); + MCModelPointer model = input->GetDescriptionModel(); + + if (model->GetNumberOfIsotropicCompartments() != numIsoCompartments || model->GetNumberOfCompartments() != numCompartments) + itkExceptionMacro("All input MCM images should have the same compartment structure."); + + m_ReferenceInputModels[i] = model; + } +} + +template +void +MCMOrientationPriorsImageFilter +::DynamicThreadedGenerateData(const InputRegionType ®ion) +{ + unsigned int numInputs = this->GetNumberOfIndexedInputs(); + unsigned int numOutputs = this->GetNumberOfIndexedOutputs(); + + using InputIteratorType = itk::ImageRegionConstIterator ; + std::vector inputIterators(numInputs); + for (unsigned int i = 0; i < numInputs; ++i) + inputIterators[i] = InputIteratorType(this->GetInput(i), region); + + using OutputIteratorType = itk::ImageRegionIterator ; + std::vector outputIterators(numOutputs); + for (unsigned int i = 0; i < numInputs; ++i) + outputIterators[i] = OutputIteratorType(this->GetOutput(i), region); + + using MaskIteratorType = itk::ImageRegionConstIterator ; + std::vector maskIterators; + if (m_MaskImages.size() == numInputs) + { + for (unsigned int i = 0; i < numInputs; ++i) + maskIterators.push_back(MaskIteratorType(m_MaskImages[i], region)); + } + unsigned int numMasks = maskIterators.size(); + + std::vector workInputModels(numInputs); + for (unsigned int i = 0;i < numInputs;++i) + workInputModels[i] = m_ReferenceInputModels[i]->Clone(); + + using ClusterFilterType = anima::SpectralClusteringFilter; + using OrientationVectorType = itk::Vector; + ClusterFilterType clusterFilter; + ClusterFilterType::MatrixType squaredDistanceMatrix; + std::vector workAllOrientations(numInputs * m_NumberOfAnisotropicCompartments); + std::vector workInputOrientations; + std::vector workAllWeights(numInputs * m_NumberOfAnisotropicCompartments); + vnl_matrix workInputWeights; + std::vector usefulInputsForWeights(numInputs, false); + std::vector workAllMemberships(numInputs * m_NumberOfAnisotropicCompartments); + OrientationVectorType workSphericalOrientation, workCartesianOrientation; + workSphericalOrientation[2] = 1.0; + std::vector clusterMembers; + // anima::WatsonDistribution watsonDistribution; + // anima::DirichletDistribution dirichletDistribution; + + std::vector< itk::VariableLengthVector > outputOrientationValues(m_NumberOfAnisotropicCompartments); + for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + outputOrientationValues[i].SetSize(4); + itk::VariableLengthVector outputWeightValues; + outputWeightValues.SetSize(m_NumberOfAnisotropicCompartments); + + while (!inputIterators[0].IsAtEnd()) + { + std::fill(workInputWeights.begin(),workInputWeights.end(),0.0); + + unsigned int numUsefulInputs = 0; + for (unsigned int i = 0;i < numInputs;++i) + { + if (isZero(inputIterators[i].Get())) + continue; + + if (numMasks == numInputs) + { + if (maskIterators[i].Get() == 0) + continue; + } + + workInputModels[i]->SetModelVector(inputIterators[i].Get()); + + ++numUsefulInputs; + } + + if (numUsefulInputs == 0) + { + for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + { + outputOrientationValues[i].Fill(0.0); + outputIterators[i].Set(outputOrientationValues[i]); + ++outputIterators[i]; + } + + outputWeightValues.Fill(0.0); + outputIterators[m_NumberOfAnisotropicCompartments].Set(outputWeightValues); + ++outputIterators[m_NumberOfAnisotropicCompartments]; + + for (unsigned int i = 0;i < numInputs;++i) + ++inputIterators[i]; + + if (numMasks == numInputs) + { + for (unsigned int i = 0;i < numMasks;++i) + ++maskIterators[i]; + } + + continue; + } + + // Actual work comes here + + // Grab list of all orientations, memberships and weights from all MCM models + for (unsigned int i = 0;i < numInputs;++i) + { + unsigned int numIsoCompartments = workInputModels[i]->GetNumberOfIsotropicCompartments(); + + double sumOfAnisotropicWeights = 0.0; + for (unsigned int j = 0;j < m_NumberOfAnisotropicCompartments;++j) + { + workSphericalOrientation[0] = workInputModels[i]->GetCompartment(j + numIsoCompartments)->GetOrientationTheta(); + workSphericalOrientation[1] = workInputModels[i]->GetCompartment(j + numIsoCompartments)->GetOrientationPhi(); + anima::TransformSphericalToCartesianCoordinates(workSphericalOrientation, workCartesianOrientation); + workAllOrientations[i * m_NumberOfAnisotropicCompartments + j] = workCartesianOrientation; + workAllMemberships[i * m_NumberOfAnisotropicCompartments + j] = i; + sumOfAnisotropicWeights += workInputModels[i]->GetCompartmentWeight(j + numIsoCompartments); + } + + for (unsigned int j = 0;j < m_NumberOfAnisotropicCompartments;++j) + workAllWeights[i * m_NumberOfAnisotropicCompartments + j] = workInputModels[i]->GetCompartmentWeight(j + numIsoCompartments) / sumOfAnisotropicWeights; + } + + // Compute squared distance matrix between orientations + unsigned int numOrientations = workAllOrientations.size(); + squaredDistanceMatrix.set_size(numOrientations, numOrientations); + squaredDistanceMatrix.fill_diagonal(0.0); + for (unsigned int i = 0;i < numOrientations - 1;++i) + { + workCartesianOrientation = workAllOrientations[i]; + for (unsigned int j = i + 1;j < numOrientations;++j) + { + double cosValue = anima::ComputeScalarProduct(workAllOrientations[j], workCartesianOrientation); + cosValue = std::abs(cosValue); + if (cosValue > 1.0) + cosValue = 1.0; + double distanceValue = std::acos(cosValue) / M_PI; + squaredDistanceMatrix.put(i, j, distanceValue * distanceValue); + if (j != i) + squaredDistanceMatrix.put(i, j, distanceValue * distanceValue); + } + } + + // Perform clustering based on orientations + clusterFilter.SetInputData(squaredDistanceMatrix); + clusterFilter.SetNbClass(m_NumberOfAnisotropicCompartments); + clusterFilter.Update(); + + // Characterize each cluster by + // - A Watson distribution for orientations, + // - A Beta distribution for the weight. + for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + { + clusterMembers = clusterFilter.GetClassMembers(i); + unsigned int clusterSize = clusterMembers.size(); + workInputOrientations.resize(clusterSize); + for (unsigned int j = 0;j < clusterSize;++j) + workInputOrientations[j] = workAllOrientations[clusterMembers[j]]; + + // watsonDistribution.Fit(workInputOrientations, "mle"); + // workCartesianOrientation = watsonDistribution.GetMeanAxis(); + // for (unsigned int j = 0;j < 3;++j) + // outputOrientationValues[i][j] = workCartesianOrientation[j]; + // outputOrientationValues[i][3] = watsonDistribution.GetConcentrationParameter(); + outputIterators[i].Set(outputOrientationValues[i]); + ++outputIterators[i]; + } + + std::fill(usefulInputsForWeights.begin(), usefulInputsForWeights.end(), false); + unsigned int nbUsefulInputsForWeights = 0; + for (unsigned int i = 0;i < numInputs;++i) + { + unsigned int numCompartments = 0; + for (unsigned int j = 0;j < m_NumberOfAnisotropicCompartments;++j) + { + clusterMembers = clusterFilter.GetClassMembers(j); + for (unsigned int k = 0;k < clusterMembers.size();++k) + { + if (workAllMemberships[clusterMembers[k]] == i) + { + numCompartments++; + break; + } + } + } + + if (numCompartments == m_NumberOfAnisotropicCompartments) + { + usefulInputsForWeights[i] = true; + nbUsefulInputsForWeights++; + } + } + + workInputWeights.set_size(nbUsefulInputsForWeights, m_NumberOfAnisotropicCompartments); + unsigned int pos = 0; + for (unsigned int i = 0;i < numInputs;++i) + { + unsigned int numCompartments = 0; + for (unsigned int j = 0;j < m_NumberOfAnisotropicCompartments;++j) + { + clusterMembers = clusterFilter.GetClassMembers(j); + for (unsigned int k = 0;k < clusterMembers.size();++k) + { + if (workAllMemberships[clusterMembers[k]] == i) + { + workInputWeights.put(pos, j, workAllWeights[clusterMembers[k]]); + break; + } + } + } + + if (usefulInputsForWeights[i]) + pos++; + } + + // dirichletDistribution.Fit(workInputWeights, "mle"); + // workAllWeights = dirichletDistribution.GetConcentrationParameters(); + // for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + // outputWeightValues[i] = workAllWeights[i]; + outputIterators[m_NumberOfAnisotropicCompartments].Set(outputWeightValues); + ++outputIterators[m_NumberOfAnisotropicCompartments]; + + for (unsigned int i = 0;i < numInputs;++i) + ++inputIterators[i]; + + if (numMasks == numInputs) + { + for (unsigned int i = 0;i < numMasks;++i) + ++maskIterators[i]; + } + } +} + +} // end namespace anima From dc08eac4e7f181714c2a70a4a2a38f2e890292da Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Fri, 27 Oct 2023 16:31:58 +0200 Subject: [PATCH 2/4] MCM orientation priors tool seem to work. Needs to link with stat distribution after PR#104 is merged. --- .../animaMCMOrientationPriors.cxx | 10 ++++- .../animaMCMOrientationPriorsImageFilter.h | 1 + .../animaMCMOrientationPriorsImageFilter.hxx | 44 ++++++++++++++----- 3 files changed, 43 insertions(+), 12 deletions(-) diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx index f44b7b2fe..8b04eb382 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx @@ -104,8 +104,14 @@ int main(int argc, char **argv) mainFilter->SetNumberOfWorkUnits(nbThreadsArg.getValue()); mainFilter->Update(); - anima::writeImage (outOrientationArg.getValue(),mainFilter->GetOutput(0)); - anima::writeImage (outWeightsArg.getValue(),mainFilter->GetOutput(1)); + unsigned int numberOfAnisotropicCompartments = mainFilter->GetNumberOfAnisotropicCompartments(); + for (unsigned int i = 0;i < numberOfAnisotropicCompartments;++i) + { + std::string fileName = outOrientationArg.getValue() + "_" + std::to_string(i) + ".nrrd"; + anima::writeImage (fileName, mainFilter->GetOutput(i)); + } + + anima::writeImage (outWeightsArg.getValue(), mainFilter->GetOutput(numberOfAnisotropicCompartments)); return EXIT_SUCCESS; } diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h index 313af24ac..53643b39f 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h @@ -44,6 +44,7 @@ public itk::ImageToImageFilter< anima::MCMImage , itk::VectorImag using MCModelPointer = MCModelType::Pointer; void AddMaskImage(MaskImageType *maskImage) {m_MaskImages.push_back(maskImage);} + unsigned int GetNumberOfAnisotropicCompartments() {return m_NumberOfAnisotropicCompartments;} protected: MCMOrientationPriorsImageFilter() diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx index 3a13ac97e..f531a2f8c 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx @@ -89,7 +89,7 @@ MCMOrientationPriorsImageFilter using OutputIteratorType = itk::ImageRegionIterator ; std::vector outputIterators(numOutputs); - for (unsigned int i = 0; i < numInputs; ++i) + for (unsigned int i = 0; i < numOutputs; ++i) outputIterators[i] = OutputIteratorType(this->GetOutput(i), region); using MaskIteratorType = itk::ImageRegionConstIterator ; @@ -173,8 +173,6 @@ MCMOrientationPriorsImageFilter continue; } - // Actual work comes here - // Grab list of all orientations, memberships and weights from all MCM models for (unsigned int i = 0;i < numInputs;++i) { @@ -218,11 +216,10 @@ MCMOrientationPriorsImageFilter // Perform clustering based on orientations clusterFilter.SetInputData(squaredDistanceMatrix); clusterFilter.SetNbClass(m_NumberOfAnisotropicCompartments); + clusterFilter.SetVerbose(false); clusterFilter.Update(); - // Characterize each cluster by - // - A Watson distribution for orientations, - // - A Beta distribution for the weight. + // Characterize orientations in each cluster by a Watson distribution for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) { clusterMembers = clusterFilter.GetClassMembers(i); @@ -236,10 +233,15 @@ MCMOrientationPriorsImageFilter // for (unsigned int j = 0;j < 3;++j) // outputOrientationValues[i][j] = workCartesianOrientation[j]; // outputOrientationValues[i][3] = watsonDistribution.GetConcentrationParameter(); + for (unsigned int j = 0;j < 3;++j) + outputOrientationValues[i][j] = workInputOrientations[0][j]; outputIterators[i].Set(outputOrientationValues[i]); ++outputIterators[i]; } + // Characterize anisotropic weights in each cluster by a Dirichlet distribution + + // First, find MCMs whose anisotropic compartments did spread in all clusters std::fill(usefulInputsForWeights.begin(), usefulInputsForWeights.end(), false); unsigned int nbUsefulInputsForWeights = 0; for (unsigned int i = 0;i < numInputs;++i) @@ -265,11 +267,32 @@ MCMOrientationPriorsImageFilter } } + // Next, compute Dirichlet priors in each cluster + if (nbUsefulInputsForWeights < 2) + { + outputWeightValues.Fill(1.0); + outputIterators[m_NumberOfAnisotropicCompartments].Set(outputWeightValues); + ++outputIterators[m_NumberOfAnisotropicCompartments]; + + for (unsigned int i = 0;i < numInputs;++i) + ++inputIterators[i]; + + if (numMasks == numInputs) + { + for (unsigned int i = 0;i < numMasks;++i) + ++maskIterators[i]; + } + + continue; + } + workInputWeights.set_size(nbUsefulInputsForWeights, m_NumberOfAnisotropicCompartments); unsigned int pos = 0; for (unsigned int i = 0;i < numInputs;++i) { - unsigned int numCompartments = 0; + if (!usefulInputsForWeights[i]) + continue; + for (unsigned int j = 0;j < m_NumberOfAnisotropicCompartments;++j) { clusterMembers = clusterFilter.GetClassMembers(j); @@ -282,15 +305,16 @@ MCMOrientationPriorsImageFilter } } } - - if (usefulInputsForWeights[i]) - pos++; + + pos++; } // dirichletDistribution.Fit(workInputWeights, "mle"); // workAllWeights = dirichletDistribution.GetConcentrationParameters(); // for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) // outputWeightValues[i] = workAllWeights[i]; + for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) + outputWeightValues[i] = workInputWeights.get(0, i); outputIterators[m_NumberOfAnisotropicCompartments].Set(outputWeightValues); ++outputIterators[m_NumberOfAnisotropicCompartments]; From 9ac75c49982c3f87ebe1b42c3840a8f01eec4fc3 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Fri, 27 Oct 2023 16:50:18 +0200 Subject: [PATCH 3/4] Some cleaning and bug fix in distance matrix computation. --- .../animaMCMOrientationPriorsImageFilter.h | 7 +++---- .../animaMCMOrientationPriorsImageFilter.hxx | 7 ++++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h index 53643b39f..168982e07 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h @@ -30,11 +30,10 @@ public itk::ImageToImageFilter< anima::MCMImage , itk::VectorImag itkTypeMacro(MCMOrientationPriorsImageFilter, ImageToImageFilter) using InputImagePointer = typename InputImageType::ConstPointer; - using OutputImagePointer = typename OutputImageType::Pointer; - using InputRegionType = typename InputImageType::RegionType; - using InputIndexType = typename InputImageType::IndexType; - using PixelType = typename InputImageType::PixelType; + + using OutputImagePointer = typename OutputImageType::Pointer; + using OutputPixelType = typename OutputImageType::PixelType; using MaskImageType = itk::Image ; using MaskImagePointer = MaskImageType::Pointer; diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx index f531a2f8c..e95f7b341 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx @@ -107,6 +107,7 @@ MCMOrientationPriorsImageFilter using ClusterFilterType = anima::SpectralClusteringFilter; using OrientationVectorType = itk::Vector; + ClusterFilterType clusterFilter; ClusterFilterType::MatrixType squaredDistanceMatrix; std::vector workAllOrientations(numInputs * m_NumberOfAnisotropicCompartments); @@ -121,10 +122,10 @@ MCMOrientationPriorsImageFilter // anima::WatsonDistribution watsonDistribution; // anima::DirichletDistribution dirichletDistribution; - std::vector< itk::VariableLengthVector > outputOrientationValues(m_NumberOfAnisotropicCompartments); + std::vector outputOrientationValues(m_NumberOfAnisotropicCompartments); for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) outputOrientationValues[i].SetSize(4); - itk::VariableLengthVector outputWeightValues; + OutputPixelType outputWeightValues; outputWeightValues.SetSize(m_NumberOfAnisotropicCompartments); while (!inputIterators[0].IsAtEnd()) @@ -209,7 +210,7 @@ MCMOrientationPriorsImageFilter double distanceValue = std::acos(cosValue) / M_PI; squaredDistanceMatrix.put(i, j, distanceValue * distanceValue); if (j != i) - squaredDistanceMatrix.put(i, j, distanceValue * distanceValue); + squaredDistanceMatrix.put(j, i, distanceValue * distanceValue); } } From 3b3def66a22291ae2d66df7a4480adab5e1de66a Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Fri, 27 Oct 2023 18:55:37 +0200 Subject: [PATCH 4/4] Working tool for MCM orientation priors. --- .../mcm_orientation_priors/CMakeLists.txt | 1 + .../animaMCMOrientationPriors.cxx | 2 +- .../animaMCMOrientationPriorsImageFilter.hxx | 27 +++++++++---------- 3 files changed, 14 insertions(+), 16 deletions(-) diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt b/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt index 9bbb3ca82..8db408adc 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt @@ -26,6 +26,7 @@ target_link_libraries(${PROJECT_NAME} ${ITKIO_LIBRARIES} ${TinyXML2_LIBRARY} AnimaMCM + AnimaStatisticalDistributions ) ## ############################################################################# diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx index 8b04eb382..8f31243f8 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx @@ -72,7 +72,7 @@ int main(int argc, char **argv) if (strcmp(tmpStr,"") == 0) continue; - std::cout << "Loading image " << nbOfImages << " : " << tmpStr << std::endl; + std::cout << "Loading image " << nbOfImages + 1 << ": " << tmpStr << std::endl; MCMReaderType mcmReader; mcmReader.SetFileName(tmpStr); diff --git a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx index e95f7b341..ad94a0de0 100644 --- a/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx +++ b/Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.hxx @@ -1,6 +1,8 @@ #pragma once #include "animaMCMOrientationPriorsImageFilter.h" #include +#include +#include #include #include @@ -114,13 +116,14 @@ MCMOrientationPriorsImageFilter std::vector workInputOrientations; std::vector workAllWeights(numInputs * m_NumberOfAnisotropicCompartments); vnl_matrix workInputWeights; + std::vector workConcentrationParameters; std::vector usefulInputsForWeights(numInputs, false); std::vector workAllMemberships(numInputs * m_NumberOfAnisotropicCompartments); OrientationVectorType workSphericalOrientation, workCartesianOrientation; workSphericalOrientation[2] = 1.0; std::vector clusterMembers; - // anima::WatsonDistribution watsonDistribution; - // anima::DirichletDistribution dirichletDistribution; + anima::WatsonDistribution watsonDistribution; + anima::DirichletDistribution dirichletDistribution; std::vector outputOrientationValues(m_NumberOfAnisotropicCompartments); for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) @@ -130,8 +133,6 @@ MCMOrientationPriorsImageFilter while (!inputIterators[0].IsAtEnd()) { - std::fill(workInputWeights.begin(),workInputWeights.end(),0.0); - unsigned int numUsefulInputs = 0; for (unsigned int i = 0;i < numInputs;++i) { @@ -229,13 +230,11 @@ MCMOrientationPriorsImageFilter for (unsigned int j = 0;j < clusterSize;++j) workInputOrientations[j] = workAllOrientations[clusterMembers[j]]; - // watsonDistribution.Fit(workInputOrientations, "mle"); - // workCartesianOrientation = watsonDistribution.GetMeanAxis(); - // for (unsigned int j = 0;j < 3;++j) - // outputOrientationValues[i][j] = workCartesianOrientation[j]; - // outputOrientationValues[i][3] = watsonDistribution.GetConcentrationParameter(); + watsonDistribution.Fit(workInputOrientations, "mle"); + workCartesianOrientation = watsonDistribution.GetMeanAxis(); for (unsigned int j = 0;j < 3;++j) - outputOrientationValues[i][j] = workInputOrientations[0][j]; + outputOrientationValues[i][j] = workCartesianOrientation[j]; + outputOrientationValues[i][3] = watsonDistribution.GetConcentrationParameter(); outputIterators[i].Set(outputOrientationValues[i]); ++outputIterators[i]; } @@ -310,12 +309,10 @@ MCMOrientationPriorsImageFilter pos++; } - // dirichletDistribution.Fit(workInputWeights, "mle"); - // workAllWeights = dirichletDistribution.GetConcentrationParameters(); - // for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) - // outputWeightValues[i] = workAllWeights[i]; + dirichletDistribution.Fit(workInputWeights, "mle"); + workConcentrationParameters = dirichletDistribution.GetConcentrationParameters(); for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i) - outputWeightValues[i] = workInputWeights.get(0, i); + outputWeightValues[i] = workConcentrationParameters[i]; outputIterators[m_NumberOfAnisotropicCompartments].Set(outputWeightValues); ++outputIterators[m_NumberOfAnisotropicCompartments];