Merge branch 'Inria-Empenn:master' into fusion

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,7 +2,9 @@ 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_orientation_priors)
 add_subdirectory(mcm_scalar_maps)
 add_subdirectory(mt_estimation_validation)
-add_subdirectory(test_averaging)
+add_subdirectory(test_averaging)

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,113 @@
+#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",
+        "A string specifying the name of the file that stores the input MCM image.",
+        true, "", "input mcm image", cmd);
+    TCLAP::ValueArg<std::string> outArg(
+        "o", "outputfile",
+        "A string specifying the name of the file that stores the output image with combined anisotropic weights.",
+        true, "", "output combined 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/mcm_tools/mcm_orientation_priors/CMakeLists.txt

@@ -0,0 +1,30 @@
+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 AnimaStatisticalDistributions)
+
+  # ############################################################################
+  # install
+  # ############################################################################
+
+  set_exe_install_rules(${PROJECT_NAME})
+
+endif()

Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriors.cxx

@@ -0,0 +1,133 @@
+#include <animaMCMFileReader.h>
+#include <animaMCMOrientationPriorsImageFilter.h>
+#include <animaReadWriteFunctions.h>
+
+#include <itkTimeProbe.h>
+
+#include <tclap/CmdLine.h>
+
+int main(int argc, char **argv)
+{
+    TCLAP::CmdLine cmd("INRIA / IRISA - VisAGeS Team", ' ', ANIMA_VERSION);
+
+    TCLAP::ValueArg<std::string> inArg(
+        "i", "input-mcms",
+        "A text file specifying a list of MCM images in a common geometry.",
+        true, "", "input MCM images", cmd);
+    TCLAP::ValueArg<std::string> 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<std::string> outWeightsArg(
+        "w", "output-weights",
+        "A string specifying the filename for the output vector image that will store fractions prior.",
+        true, "", "output weights priors", cmd);
+
+    TCLAP::ValueArg<std::string> 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<std::string> meanOrientationsArg(
+        "", "mean-orientations",
+        "A string specifying the basename for the output vector images that will store mean orientations.",
+        false, "", "output mean orientation images", cmd);
+    TCLAP::ValueArg<std::string> meanFractionsArg(
+        "", "mean-fractions",
+        "A string specifying the filename for the output vector image that will store fractions mean.",
+        false, "", "output mean fractions image", cmd);
+    TCLAP::ValueArg<unsigned int> 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<double>;
+    MainFilterType::Pointer mainFilter = MainFilterType::New();
+
+    using MCMReaderType = anima::MCMFileReader<double, 3>;
+    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 + 1 << ": " << 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<MaskImageType>(tmpStr));
+        }
+    }
+
+    mainFilter->SetNumberOfWorkUnits(nbThreadsArg.getValue());
+    mainFilter->Update();
+
+    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<OutputImageType>(fileName, mainFilter->GetOutput(i));
+    }
+
+    anima::writeImage<OutputImageType>(outWeightsArg.getValue(), mainFilter->GetOutput(numberOfAnisotropicCompartments));
+
+    if (meanOrientationsArg.getValue() != "")
+    {
+        for (unsigned int i = 0; i < numberOfAnisotropicCompartments; ++i)
+        {
+            std::string fileName = meanOrientationsArg.getValue() + "_" + std::to_string(i) + ".nrrd";
+            anima::writeImage<OutputImageType>(fileName, mainFilter->GetMeanOrientationImage(i));
+        }
+    }
+
+    if (meanFractionsArg.getValue() != "")
+        anima::writeImage<OutputImageType>(meanFractionsArg.getValue(), mainFilter->GetMeanFractionsImage());
+
+    return EXIT_SUCCESS;
+}