From 3b3def66a22291ae2d66df7a4480adab5e1de66a Mon Sep 17 00:00:00 2001
From: Aymeric Stamm <aymeric.stamm@cnrs.fr>
Date: Fri, 27 Oct 2023 18:55:37 +0200
Subject: [PATCH] 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 <animaSpectralClusteringFilter.h>
+#include <animaDirichletDistribution.h>
+#include <animaWatsonDistribution.h>
 
 #include <itkImageRegionConstIterator.h>
 #include <itkImageRegionIterator.h>
@@ -114,13 +116,14 @@ MCMOrientationPriorsImageFilter <TPixelType>
     std::vector<OrientationVectorType> workInputOrientations;
     std::vector<double> workAllWeights(numInputs * m_NumberOfAnisotropicCompartments);
     vnl_matrix<double> workInputWeights;
+    std::vector<double> workConcentrationParameters;
     std::vector<bool> usefulInputsForWeights(numInputs, false);
     std::vector<unsigned int> workAllMemberships(numInputs * m_NumberOfAnisotropicCompartments);
     OrientationVectorType workSphericalOrientation, workCartesianOrientation;
     workSphericalOrientation[2] = 1.0;
     std::vector<unsigned int> clusterMembers;
-    // anima::WatsonDistribution watsonDistribution;
-    // anima::DirichletDistribution dirichletDistribution;
+    anima::WatsonDistribution watsonDistribution;
+    anima::DirichletDistribution dirichletDistribution;
 
     std::vector<OutputPixelType> outputOrientationValues(m_NumberOfAnisotropicCompartments);
     for (unsigned int i = 0;i < m_NumberOfAnisotropicCompartments;++i)
@@ -130,8 +133,6 @@ MCMOrientationPriorsImageFilter <TPixelType>
 
     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 <TPixelType>
             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 <TPixelType>
             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];