Merge pull request #105 from astamm/orientation-priors

Tool for computing orientation and weight priors for MCM anisotropic compartments

Anima/diffusion/mcm_tools/CMakeLists.txt

@@ -4,6 +4,7 @@ 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)

Anima/diffusion/mcm_tools/mcm_orientation_priors/CMakeLists.txt

@@ -0,0 +1,38 @@
+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,117 @@
+#include <animaMCMFileReader.h>
+#include <animaReadWriteFunctions.h>
+#include <animaMCMOrientationPriorsImageFilter.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 priors on the weights.", 
+        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<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));
+
+    return EXIT_SUCCESS;
+}

Anima/diffusion/mcm_tools/mcm_orientation_priors/animaMCMOrientationPriorsImageFilter.h

@@ -0,0 +1,83 @@
+#pragma once
+
+#include <animaMCMImage.h>
+#include <animaMultiCompartmentModel.h>
+
+#include <itkImageToImageFilter.h>
+#include <itkVectorImage.h>
+
+namespace anima
+{
+
+template <class TPixelType>
+class MCMOrientationPriorsImageFilter :
+public itk::ImageToImageFilter< anima::MCMImage <TPixelType, 3>, itk::VectorImage<TPixelType, 3> >
+{
+public:
+    /** Standard class type def */
+
+    using Self = MCMOrientationPriorsImageFilter;
+    using InputImageType = anima::MCMImage <TPixelType, 3>;
+    using OutputImageType = itk::VectorImage <TPixelType, 3>;
+    using Superclass = itk::ImageToImageFilter <InputImageType, OutputImageType >;
+    using Pointer = itk::SmartPointer<Self>;
+    using ConstPointer = itk::SmartPointer<const Self>;
+
+    /** 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 InputRegionType = typename InputImageType::RegionType;
+
+    using OutputImagePointer = typename OutputImageType::Pointer;
+    using OutputPixelType = typename OutputImageType::PixelType;
+
+    using MaskImageType = itk::Image <unsigned char, 3>;
+    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);}
+    unsigned int GetNumberOfAnisotropicCompartments() {return m_NumberOfAnisotropicCompartments;}
+
+protected:
+    MCMOrientationPriorsImageFilter()
+    {
+        m_ReferenceInputModels.clear();
+        m_MaskImages.clear();
+        m_NumberOfAnisotropicCompartments = 0;
+    }
+
+    virtual ~MCMOrientationPriorsImageFilter() {}
+
+    bool isZero(const itk::VariableLengthVector <double> &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 &region) ITK_OVERRIDE;
+
+private:
+    ITK_DISALLOW_COPY_AND_ASSIGN(MCMOrientationPriorsImageFilter);
+
+    std::vector <MCModelPointer> m_ReferenceInputModels;
+    std::vector <MaskImagePointer> m_MaskImages;
+    unsigned int m_NumberOfAnisotropicCompartments;
+};
+
+} // end namespace anima
+
+#include "animaMCMOrientationPriorsImageFilter.hxx"