Repository: connectomicslab/connectomemapper3 (language: python)

Issue:
[enh] Merge last update with debug parcellation (unassigned voxels in insula) and integration msmt_csd to master branch



Discussion / hints:

Fix patch:
diff --git a/.readthedocs.yml b/.readthedocs.yml
index 094871b69..3d486a9e8 100644
--- a/.readthedocs.yml
+++ b/.readthedocs.yml
@@ -10,8 +10,8 @@ sphinx:
   configuration: docs/conf.py
 
 # Build documentation with MkDocs
-#mkdocs:
-#  configuration: mkdocs.yml
+# mkdocs:
+#   configuration: mkdocs.yml
 
 # Optionally build your docs in additional formats such as PDF
 formats:
@@ -22,10 +22,13 @@ formats:
 # conda:
 #   environment: environment.yml
 
-# Optionally set the version of Python and requirements required to build your docs
+build:  
+  os: ubuntu-22.04
+  tools:
+    python: "3.7"
+
 python:
-  version: 3.7
   install:
     - requirements: docs/requirements.txt
-    - method: setuptools
-      path: .
+    #- method: setuptools
+    #- path: .
diff --git a/Dockerfile b/Dockerfile
index 73856b21a..3f4a4d463 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -114,21 +114,26 @@ RUN apt-get update && \
     curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -o /tmp/miniconda.sh && \
     bash /tmp/miniconda.sh -bfp /opt/conda && \
     rm -rf /tmp/miniconda.sh && \
+    apt-get install -y gcc && \
     apt-get remove -y curl && \
     conda update conda && \
     conda clean --all --yes && \
-    rm -rf ~/.conda ~/.cache/pip/* && \
+    rm -rf ~/.conda ~/.cache//* && \
     apt-get clean && \
-    rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
+    rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/* 
+
 
 ## Create conda environment, including ANTs 2.2.0 and MRtrix 3.0.2
-ENV CONDA_ENV="py37cmp-core"
+ENV CONDA_ENV="py38cmp-core"
 COPY docker/spec-file.txt /app/spec-file.txt
 COPY docker/requirements.txt /app/requirements.txt
-RUN /bin/bash -c "conda config --set default_threads 4 &&\
-    conda create --name ${CONDA_ENV} --file /app/spec-file.txt &&\
+COPY docker/environment.yml /app/environment.yml
+RUN /bin/bash -c "conda config --set default_threads 3 &&\ 
+    conda env create --file /app/environment.yml &&\ 
     . activate ${CONDA_ENV} &&\
-    pip install -r /app/requirements.txt &&\
+    python --version &&\
+    pip install --upgrade pip setuptools wheel &&\
+    pip install -r /app/requirements.txt&&\
     conda clean -v --all --yes &&\
     rm -rf ~/.conda ~/.cache/pip/*"
 
@@ -140,7 +145,7 @@ RUN /bin/bash -c "conda config --set default_threads 4 &&\
 
 ##################################################################
 # Installation of Connectome Mapper 3 packages
-##################################################################
+################################################################t#
 FROM neurocondabuntu AS cmpbuntu
 
 # Docker build command arguments
@@ -183,7 +188,7 @@ RUN mkdir -p /cache/python-eggs && \
     chmod -R 777 /cache/python-eggs
 
 # Install cmp and cmtklib packages in the conda environment $CONDA_ENV
-ENV CONDA_ENV="py37cmp-core"
+ENV CONDA_ENV="py38cmp-core"
 RUN /bin/bash -c ". activate ${CONDA_ENV} &&\
     pip install ."
 
diff --git a/cmp/bidsappmanager/gui/handlers.py b/cmp/bidsappmanager/gui/handlers.py
index 10573d96c..8ccf515bc 100644
--- a/cmp/bidsappmanager/gui/handlers.py
+++ b/cmp/bidsappmanager/gui/handlers.py
@@ -858,15 +858,15 @@ def load_dataset(self, ui_info, debug=False):
                                     "DTI" in loaded_project.dmri_bids_acq
                                 ):
                                     loaded_project.diffusion_imaging_model = "DTI"
-                                elif ("hardi" in loaded_project.dmri_bids_acq) or (
-                                    "HARDI" in loaded_project.dmri_bids_acq
+                                elif ("HARDI" in loaded_project.dmri_bids_acq) or (
+                                    "HARDI_singleTissue" in loaded_project.dmri_bids_acq
                                 ):
-                                    loaded_project.diffusion_imaging_model = "HARDI"
-                                elif ("multishell" in loaded_project.dmri_bids_acq) or (
-                                    "MULTISHELL" in loaded_project.dmri_bids_acq
+                                    loaded_project.diffusion_imaging_model = "HARDI_singleTissue"
+                                elif ("HARDI_multiTissue" in loaded_project.dmri_bids_acq) or (
+                                    "multishell" in loaded_project.dmri_bids_acq
                                 ):
                                     loaded_project.diffusion_imaging_model = (
-                                        "multishell"
+                                        "HARDI_multiTissue"
                                     )
                                 else:
                                     loaded_project.diffusion_imaging_model = "DTI"
diff --git a/cmp/bidsappmanager/stages/diffusion/reconstruction.py b/cmp/bidsappmanager/stages/diffusion/reconstruction.py
index dce15c5fa..8b64deb67 100644
--- a/cmp/bidsappmanager/stages/diffusion/reconstruction.py
+++ b/cmp/bidsappmanager/stages/diffusion/reconstruction.py
@@ -105,12 +105,14 @@ class MRtrix_recon_configUI(MRtrixReconConfig):
     flip_table_axis = List(editor=CheckListEditor(values=["x", "y", "z"], cols=3))
 
     traits_view = View(
-        Item("flip_table_axis", style="custom", label="Flip gradient table:"),
-        Item("local_model", editor=EnumEditor(name="local_model_editor")),
-        Group(
-            Item("lmax_order"),
-            # Item('normalize_to_B0'),
-            Item("single_fib_thr", label="FA threshold"),
-            visible_when="local_model",
-        ),
-    )
+            Item("flip_table_axis", style="custom", label="Flip gradient table:"),
+            Item("local_model", editor=EnumEditor(name="local_model_editor")),
+            Group(
+                Group(
+                    Item("lmax_order"),
+                    # Item('normalize_to_B0'),
+                    Item("single_fib_thr", label="FA threshold"),
+                    visible_when="local_model",
+                ),
+            visible_when='imaging_model != "HARDI_multiTissue"',
+            ))
diff --git a/cmp/info.py b/cmp/info.py
index b5253aa15..b1104bca5 100644
--- a/cmp/info.py
+++ b/cmp/info.py
@@ -5,6 +5,7 @@
 #  This software is distributed under the open-source license Modified BSD.
 """This file contains cmp package information."""
 
+
 _version_major = 3
 _version_minor = 1
 _version_micro = 0
@@ -23,7 +24,7 @@
 # __current_year__ = datetime.datetime.now().strftime("%Y")
 __current_year__ = "2022"
 
-__author__ = "The CMP3 developers"
+__author__ = "The CMP3 developers"  
 
 __copyright__ = (
     "Copyright (C) 2009-{}, ".format(__current_year__)
diff --git a/cmp/pipelines/diffusion/diffusion.py b/cmp/pipelines/diffusion/diffusion.py
index d4618068e..1fb92726c 100644
--- a/cmp/pipelines/diffusion/diffusion.py
+++ b/cmp/pipelines/diffusion/diffusion.py
@@ -1102,12 +1102,14 @@ def remove_non_existing_scales(roi_volumes):
                 )
                 # fmt:on
 
+
         if self.stages["Diffusion"].enabled:
             diff_flow = self.create_stage_flow("Diffusion")
             # fmt:off
             diffusion_flow.connect(
                 [
-                    (preproc_flow, diff_flow, [("outputnode.diffusion_preproc", "inputnode.diffusion")]),
+                    (preproc_flow, diff_flow, [("outputnode.diffusion_preproc", "inputnode.diffusion"),
+                                               ("outputnode.in_5tt_orig", "inputnode.in_5tt_orig")],),
                     (reg_flow, diff_flow, [("outputnode.wm_mask_registered_crop", "inputnode.wm_mask_registered",),
                                            ("outputnode.brain_mask_registered_crop", "inputnode.brain_mask_registered",),
                                            ("outputnode.partial_volumes_registered_crop", "inputnode.partial_volumes",),
diff --git a/cmp/project.py b/cmp/project.py
index 14ab6f47c..85fb0100d 100644
--- a/cmp/project.py
+++ b/cmp/project.py
@@ -80,7 +80,7 @@ class ProjectInfo(HasTraits):
 
     diffusion_imaging_model : traits.Str
         Diffusion imaging model that can be
-        'DSI', 'DTI', 'HARDI' or 'multishell'
+        'DSI', 'DTI', 'HARDI_singleTissue' or 'HARDI_multiTissue'
 
     dmri_bids_acqs : traits.List
         List diffusion imaging models extracted from ``acq-<label>`` filename part.
@@ -228,7 +228,7 @@ class ProjectInfo(HasTraits):
     )
 
     # process_type = Enum('diffusion',['diffusion','fMRI'])
-    diffusion_imaging_model = Enum("DTI", ["DSI", "DTI", "HARDI", "multishell"])
+    diffusion_imaging_model = Enum("DTI", ["DSI", "DTI", "HARDI_singleTissue", "HARDI_multiTissue"])
     parcellation_scheme = Str("Lausanne2018")
     atlas_info = Dict()
     freesurfer_subjects_dir = Str("")
diff --git a/cmp/stages/diffusion/diffusion.py b/cmp/stages/diffusion/diffusion.py
index bea53a854..4c568e326 100644
--- a/cmp/stages/diffusion/diffusion.py
+++ b/cmp/stages/diffusion/diffusion.py
@@ -24,7 +24,7 @@ class DiffusionConfig(HasTraits):
 
     Attributes
     ----------
-    diffusion_imaging_model_editor : ['DSI', 'DTI', 'HARDI', 'multishell']
+    diffusion_imaging_model_editor : ['DSI', 'DTI', 'HARDI_singleTissue', 'HARDI_multiTissue']
         Available diffusion imaging models
 
     diffusion_imaging_model : traits.Str
@@ -89,7 +89,7 @@ class DiffusionConfig(HasTraits):
     cmp.stages.diffusion.diffusion.DiffusionStage
     """
 
-    diffusion_imaging_model_editor = List(["DSI", "DTI", "HARDI", "multishell"])
+    diffusion_imaging_model_editor = List(["DSI", "DTI", "HARDI_singleTissue", "HARDI_multiTissue"])
     diffusion_imaging_model = Str("DTI")
     dilate_rois = Bool(True)
     dilation_kernel = Enum(["Box", "Gauss", "Sphere"])
@@ -117,7 +117,8 @@ def __init__(self):
             tracking_processing_tool=self.tracking_processing_tool,
         )
         self.mrtrix_recon_config = MRtrixReconConfig(
-            imaging_model=self.diffusion_imaging_model, recon_mode=self.diffusion_model
+            imaging_model=self.diffusion_imaging_model, 
+            recon_mode=self.diffusion_model,
         )
         self.dipy_tracking_config = DipyTrackingConfig(
             imaging_model=self.diffusion_imaging_model,
@@ -336,6 +337,7 @@ def __init__(self, bids_dir, output_dir):
             "wm_mask_registered",
             "brain_mask_registered",
             "act_5tt_registered",
+            "in_5tt_orig",
             "gmwmi_registered",
             "roi_volumes",
             "grad",
@@ -460,6 +462,8 @@ def create_workflow(self, flow, inputnode, outputnode):
                     (inputnode, recon_flow, [("diffusion", "inputnode.diffusion")]),
                     (inputnode, recon_flow, [("grad", "inputnode.grad")]),
                     (inputnode, recon_flow, [("diffusion", "inputnode.diffusion_resampled")],),
+                    (inputnode, recon_flow, [("act_5tt_registered", "inputnode.act_5tt_registered")],),
+                    (inputnode, recon_flow, [("in_5tt_orig", "inputnode.in_5tt_orig")],),
                     (inputnode, recon_flow, [("brain_mask_registered", "inputnode.wm_mask_resampled")],),
                     (recon_flow, outputnode, [("outputnode.FA", "FA")]),
                     (recon_flow, outputnode, [("outputnode.ADC", "ADC")]),
diff --git a/cmp/stages/diffusion/reconstruction.py b/cmp/stages/diffusion/reconstruction.py
index e573ecc20..dcda37093 100644
--- a/cmp/stages/diffusion/reconstruction.py
+++ b/cmp/stages/diffusion/reconstruction.py
@@ -23,8 +23,10 @@
     MRtrix_mul,
     MRThreshold,
     MRConvert,
-    EstimateResponseForSH,
-    ConstrainedSphericalDeconvolution,
+    EstimateResponseForSHSingleTissue,
+    EstimateResponseForSHMultiTissue,
+    ConstrainedSphericalDeconvolutionSingleTissue,
+    ConstrainedSphericalDeconvolutionMultiTissue,
     DWI2Tensor,
     Tensor2Vector,
 )
@@ -38,7 +40,6 @@
 
 iflogger = logging.getLogger("nipype.interface")
 
-
 class DipyReconConfig(HasTraits):
     """Class used to store Dipy diffusion reconstruction sub-workflow configuration parameters.
 
@@ -133,7 +134,7 @@ class DipyReconConfig(HasTraits):
     flip_table_axis = List(["x", "y", "z"])
     # gradient_table = File
     local_model_editor = Dict(
-        {False: "1:Tensor", True: "2:Constrained Spherical Deconvolution"}
+        {False: "1:Tensor", True: "2:Constrained Spherical Deconvolution (single tissue)"}
     )
     local_model = Bool(True)
     lmax_order = Enum([2, 4, 6, 8, 10, 12, 14, 16])
@@ -208,12 +209,13 @@ def _imaging_model_changed(self, new):
         elif new == "DTI":
             self.local_model_editor = {
                 False: "1:Tensor",
-                True: "2:Constrained Spherical Deconvolution",
+                True: "2:Constrained Spherical Deconvolution (single tissue)",
             }
-        elif new == "multishell" or new == "HARDI":
-            self.local_model_editor = {True: "Constrained Spherical Deconvolution"}
+        elif new == "HARDI_singleTissue" or new == "HARDI_multiTissue":
+            self.local_model_editor = {True: "2:Constrained Spherical Deconvolution (single tissue)"}
             self.local_model = True
 
+
     def _recon_mode_changed(self, new):
         """Update ``local_model_editor`` and ``self.local_model`` when ``recon_mode`` is updated.
 
@@ -222,17 +224,15 @@ def _recon_mode_changed(self, new):
         new : string
             New value of ``recon_mode``
         """
-        if new == "Probabilistic" and self.imaging_model != "DSI":
-            self.local_model_editor = {True: "Constrained Spherical Deconvolution"}
-            self.local_model = True
-        elif new == "Probabilistic" and self.imaging_model == "DSI":
+
+        if new == "Probabilistic" and self.imaging_model == "DSI":
             pass
         else:
-            self.local_model_editor = {
-                False: "1:Tensor",
-                True: "2:Constrained Spherical Deconvolution",
-            }
-
+            pass
+            #self.local_model_editor = {
+            #    False: "1:Tensor",
+            #    True: "2:Constrained Spherical Deconvolution",
+            #}
 
 class MRtrixReconConfig(HasTraits):
     """Class used to store Dipy diffusion reconstruction sub-workflow configuration parameters.
@@ -259,10 +259,12 @@ class MRtrixReconConfig(HasTraits):
         Can be "Probabilistic" or "Deterministic"
     """
 
+    imaging_model = Str
     # gradient_table = File
     flip_table_axis = List(["x", "y", "z"])
     local_model_editor = Dict(
-        {False: "1:Tensor", True: "2:Constrained Spherical Deconvolution"}
+        {False: "1:Tensor", True: "2:Constrained Spherical Deconvolution (single tissue)",  False: "3:Constrained Spherical Deconvolution (multi tissues)"}
+
     )
     local_model = Bool(True)
     lmax_order = Enum([2, 4, 6, 8, 10, 12, 14, 16])
@@ -270,6 +272,7 @@ class MRtrixReconConfig(HasTraits):
     single_fib_thr = Float(0.7, min=0, max=1)
     recon_mode = Str
 
+
     def _imaging_model_changed(self, new):
         """Update ``local_model_editor`` and ``self.local_model`` when ``imaging_model`` is updated.
 
@@ -277,16 +280,24 @@ def _imaging_model_changed(self, new):
         ----------
         new : string
             New value of ``imaging_model``
-        """
+        """   
+
         if new == "DTI":
             self.local_model_editor = {
                 False: "1:Tensor",
-                True: "2:Constrained Spherical Deconvolution",
+                True: "2:Constrained Spherical Deconvolution (single tissue)",
             }
-        elif new == "multishell" or new == "HARDI":
-            self.local_model_editor = {True: "Constrained Spherical Deconvolution"}
+        elif new == "HARDI_singleTissue":
+            self.local_model_editor = {True: "2:Constrained Spherical Deconvolution (single tissue)"}
+            self.local_model = True
+            
+
+        elif new == "HARDI_multiTissue":
+            self.local_model_editor = {True: "3:Constrained Spherical Deconvolution (multi tissues)"}
             self.local_model = True
 
+    
+
     def _recon_mode_changed(self, new):
         """Update ``local_model_editor`` and ``self.local_model`` when ``recon_mode`` is updated.
 
@@ -295,14 +306,16 @@ def _recon_mode_changed(self, new):
         new : string
             New value of ``recon_mode``
         """
-        if new == "Probabilistic":
-            self.local_model_editor = {True: "Constrained Spherical Deconvolution"}
-            self.local_model = True
-        else:
-            self.local_model_editor = {
-                False: "1:Tensor",
-                True: "2:Constrained Spherical Deconvolution",
-            }
+        #if new == "Probabilistic":
+        #    self.local_model_editor = {True: "Constrained Spherical Deconvolution"}
+        #    self.local_model = True
+        #else:
+        #    self.local_model_editor = {
+        #        False: "1:Tensor",
+        #        True: "2:Constrained Spherical Deconvolution",
+        #    }
+        pass
+
 
 
 def create_dipy_recon_flow(config):
@@ -328,6 +341,8 @@ def create_dipy_recon_flow(config):
                 "wm_mask_resampled",
                 "bvals",
                 "bvecs",
+                "act_5tt_registered",
+                "in_5tt_orig",
             ]
         ),
         name="inputnode",
@@ -539,7 +554,7 @@ def create_mrtrix_recon_flow(config):
     flow = pe.Workflow(name="reconstruction")
     inputnode = pe.Node(
         interface=util.IdentityInterface(
-            fields=["diffusion", "diffusion_resampled", "wm_mask_resampled", "grad"]
+            fields=["diffusion", "diffusion_resampled", "wm_mask_resampled", "grad", "act_5tt_registered", "in_5tt_orig"]
         ),
         name="inputnode",
     )
@@ -619,74 +634,173 @@ def create_mrtrix_recon_flow(config):
 
     # Constrained Spherical Deconvolution
     if config.local_model:
-        print("CSD true")
-        # Compute single fiber voxel mask
-        mrtrix_erode = pe.Node(
-            interface=Erode(out_filename="wm_mask_res_eroded.nii.gz"),
-            name="mrtrix_erode",
-        )
-        mrtrix_erode.inputs.number_of_passes = 1
-        mrtrix_erode.inputs.filtertype = "erode"
-        mrtrix_mul_eroded_FA = pe.Node(
-            interface=MRtrix_mul(), name="mrtrix_mul_eroded_FA"
-        )
-        mrtrix_mul_eroded_FA.inputs.out_filename = "diffusion_resampled_tensor_FA_masked.mif"
-        mrtrix_thr_FA = pe.Node(
-            interface=MRThreshold(out_file="FA_th.mif"), name="mrtrix_thr"
-        )
-        mrtrix_thr_FA.inputs.abs_value = config.single_fib_thr
-        # fmt:off
-        flow.connect(
-            [
-                (inputnode, mrtrix_erode, [("wm_mask_resampled", "in_file")]),
-                (mrtrix_erode, mrtrix_mul_eroded_FA, [("out_file", "input2")]),
-                (mrtrix_tensor_metrics, mrtrix_mul_eroded_FA, [("out_fa", "input1")]),
-                (mrtrix_mul_eroded_FA, mrtrix_thr_FA, [("out_file", "in_file")]),
-            ]
-        )
-        # fmt:on
 
-        # Compute single fiber response function
-        mrtrix_rf = pe.Node(interface=EstimateResponseForSH(), name="mrtrix_rf")
-        mrtrix_rf.inputs.maximum_harmonic_order = int(config.lmax_order)
-        mrtrix_rf.inputs.algorithm = "tournier"
-        # mrtrix_rf.inputs.normalise = config.normalize_to_B0
-        # fmt:off
-        flow.connect(
+        ### Multi Tissue Spherical Deconvolution
+        if config.imaging_model=='HARDI_multiTissue':
+
+            print("CSD true")
+
+            # Compute single fiber voxel mask
+            mrtrix_erode = pe.Node(
+                interface=Erode(out_filename="wm_mask_res_eroded.nii.gz"),
+                name="mrtrix_erode",
+            )
+            mrtrix_erode.inputs.number_of_passes = 1
+            mrtrix_erode.inputs.filtertype = "erode"
+            mrtrix_mul_eroded_FA = pe.Node(
+                interface=MRtrix_mul(), name="mrtrix_mul_eroded_FA"
+            )
+            mrtrix_mul_eroded_FA.inputs.out_filename = "diffusion_resampled_tensor_FA_masked.mif"
+            mrtrix_thr_FA = pe.Node(
+                interface=MRThreshold(out_file="FA_th.mif"), name="mrtrix_thr"
+            )
+            mrtrix_thr_FA.inputs.abs_value = config.single_fib_thr
+            # fmt:off
+
+            flow.connect(
+                [
+                    (inputnode, mrtrix_erode, [("wm_mask_resampled", "in_file")]),
+                    (mrtrix_erode, mrtrix_mul_eroded_FA, [("out_file", "input2")]),
+                    (mrtrix_tensor_metrics, mrtrix_mul_eroded_FA, [("out_fa", "input1")]),
+                    (mrtrix_mul_eroded_FA, mrtrix_thr_FA, [("out_file", "in_file")]),
+                ]
+            )
+            # fmt:on
+
+            # Compute multi tissue  response function
+            mrtrix_rf = pe.Node(interface=EstimateResponseForSHMultiTissue(), name="mrtrix_rf")
+            #mrtrix_rf.inputs.maximum_harmonic_order = int(config.lmax_order)
+            mrtrix_rf.inputs.algorithm = "msmt_5tt"
+            # mrtrix_rf.inputs.normalise = config.normalize_to_B0
+            # fmt:off
+            flow.connect(
+                [
+                    (inputnode, mrtrix_rf, [("diffusion_resampled", "in_file")]),
+                    (inputnode, mrtrix_rf, [("in_5tt_orig", "in_5tt_file")]),
+                    #(mrtrix_thr_FA, mrtrix_rf, [("thresholded", "mask_image")]),
+
+                    (flip_table, mrtrix_rf, [("table", "encoding_file")]),
+                ]
+            )
+            # fmt:on
+
+        
+            # Perform multi tissue spherical deconvolution
+            mrtrix_CSD = pe.Node(
+                interface=ConstrainedSphericalDeconvolutionMultiTissue(), name="mrtrix_CSD"
+            )
+            mrtrix_CSD.inputs.algorithm = "msmt_csd"
+            #mrtrix_CSD.inputs.maximum_harmonic_order = int(config.lmax_order)
+            # mrtrix_CSD.inputs.normalise = config.normalize_to_B0
+
+            convert_CSD = pe.Node(
+                interface=MRConvert(out_filename="spherical_harmonics_image.nii.gz"),
+                name="convert_CSD",
+            )
+            convert_gm_CSD = pe.Node(
+                interface=MRConvert(out_filename="gm_spherical_harmonics_image.nii.gz"),
+                name="convert_gm_CSD",
+            )
+            convert_csf_CSD = pe.Node(
+                interface=MRConvert(out_filename="csf_spherical_harmonics_image.nii.gz"),
+                name="convert_csf_CSD",
+            )
+            # fmt:off
+            flow.connect(
             [
-                (inputnode, mrtrix_rf, [("diffusion_resampled", "in_file")]),
-                (mrtrix_thr_FA, mrtrix_rf, [("thresholded", "mask_image")]),
-                (flip_table, mrtrix_rf, [("table", "encoding_file")]),
+                    (inputnode, mrtrix_CSD, [("diffusion_resampled", "in_file")]),
+                    (mrtrix_rf, mrtrix_CSD, [("response", "response_file")]),
+                    (mrtrix_rf, mrtrix_CSD, [("response_gm", "response_gm_file")]),
+                    (mrtrix_rf, mrtrix_CSD, [("response_csf", "response_csf_file")]),
+                    (mrtrix_rf, outputnode, [("response", "RF")]),
+                    (flip_table, mrtrix_CSD, [("table", "encoding_file")]),
+
+                    #(mrtrix_CSD, convert_CSD, [("spherical_harmonics_image", "in_file")]),
+                    #(mrtrix_CSD, convert_CSD, [("CSD", "in_file")]),
+                    (mrtrix_CSD, convert_CSD, [("spherical_harmonics_image", "in_file")]),
+                    (mrtrix_CSD, convert_gm_CSD, [("gm_spherical_harmonics_image", "in_file")]),
+                    (mrtrix_CSD, convert_csf_CSD, [("csf_spherical_harmonics_image", "in_file")]),
+                    (convert_CSD, outputnode, [("converted", "DWI")])
+                    # (mrtrix_CSD,outputnode,[('spherical_harmonics_image','DWI')])
             ]
-        )
+            )
         # fmt:on
 
-        # Perform spherical deconvolution
-        mrtrix_CSD = pe.Node(
-            interface=ConstrainedSphericalDeconvolution(), name="mrtrix_CSD"
-        )
-        mrtrix_CSD.inputs.algorithm = "csd"
-        mrtrix_CSD.inputs.maximum_harmonic_order = int(config.lmax_order)
-        # mrtrix_CSD.inputs.normalise = config.normalize_to_B0
+            
 
-        convert_CSD = pe.Node(
-            interface=MRConvert(out_filename="spherical_harmonics_image.nii.gz"),
-            name="convert_CSD",
-        )
-        # fmt:off
-        flow.connect(
+        ### Single Tissue Spherical Deconvolution
+        else:
+            print("CSD true")
+            # Compute single fiber voxel mask
+            mrtrix_erode = pe.Node(
+                interface=Erode(out_filename="wm_mask_res_eroded.nii.gz"),
+                name="mrtrix_erode",
+            )
+            mrtrix_erode.inputs.number_of_passes = 1
+            mrtrix_erode.inputs.filtertype = "erode"
+            mrtrix_mul_eroded_FA = pe.Node(
+                interface=MRtrix_mul(), name="mrtrix_mul_eroded_FA"
+            )
+            mrtrix_mul_eroded_FA.inputs.out_filename = "diffusion_resampled_tensor_FA_masked.mif"
+            mrtrix_thr_FA = pe.Node(
+                interface=MRThreshold(out_file="FA_th.mif"), name="mrtrix_thr"
+            )
+            mrtrix_thr_FA.inputs.abs_value = config.single_fib_thr
+            # fmt:off
+            flow.connect(
+                [
+                    (inputnode, mrtrix_erode, [("wm_mask_resampled", "in_file")]),
+                    (mrtrix_erode, mrtrix_mul_eroded_FA, [("out_file", "input2")]),
+                    (mrtrix_tensor_metrics, mrtrix_mul_eroded_FA, [("out_fa", "input1")]),
+                    (mrtrix_mul_eroded_FA, mrtrix_thr_FA, [("out_file", "in_file")]),
+                ]
+            )
+            # fmt:on
+
+            # Compute single fiber response function
+            mrtrix_rf = pe.Node(interface=EstimateResponseForSHSingleTissue(), name="mrtrix_rf")
+            mrtrix_rf.inputs.maximum_harmonic_order = int(config.lmax_order)
+            mrtrix_rf.inputs.algorithm = "tournier"
+            # mrtrix_rf.inputs.normalise = config.normalize_to_B0
+            # fmt:off
+            flow.connect(
+                [
+                    (inputnode, mrtrix_rf, [("diffusion_resampled", "in_file")]),
+                    (mrtrix_thr_FA, mrtrix_rf, [("thresholded", "mask_image")]),
+                    (flip_table, mrtrix_rf, [("table", "encoding_file")]),
+                ]
+            )
+            # fmt:on
+
+        
+            # Perform spherical deconvolution
+            mrtrix_CSD = pe.Node(
+                interface=ConstrainedSphericalDeconvolutionSingleTissue(), name="mrtrix_CSD"
+            )
+            mrtrix_CSD.inputs.algorithm = "csd"
+            mrtrix_CSD.inputs.maximum_harmonic_order = int(config.lmax_order)
+            # mrtrix_CSD.inputs.normalise = config.normalize_to_B0
+
+            convert_CSD = pe.Node(
+                interface=MRConvert(out_filename="spherical_harmonics_image.nii.gz"),
+                name="convert_CSD",
+            )
+            # fmt:off
+            flow.connect(
             [
-                (inputnode, mrtrix_CSD, [("diffusion_resampled", "in_file")]),
-                (mrtrix_rf, mrtrix_CSD, [("response", "response_file")]),
-                (mrtrix_rf, outputnode, [("response", "RF")]),
-                (inputnode, mrtrix_CSD, [("wm_mask_resampled", "mask_image")]),
-                (flip_table, mrtrix_CSD, [("table", "encoding_file")]),
-                (mrtrix_CSD, convert_CSD, [("spherical_harmonics_image", "in_file")]),
-                (convert_CSD, outputnode, [("converted", "DWI")])
-                # (mrtrix_CSD,outputnode,[('spherical_harmonics_image','DWI')])
-            ]
-        )
-        # fmt:on
+                    (inputnode, mrtrix_CSD, [("diffusion_resampled", "in_file")]),
+                    (mrtrix_rf, mrtrix_CSD, [("response", "response_file")]),
+                    (mrtrix_rf, outputnode, [("response", "RF")]),
+                    (inputnode, mrtrix_CSD, [("wm_mask_resampled", "mask_image")]),
+                    (flip_table, mrtrix_CSD, [("table", "encoding_file")]),
+                    (mrtrix_CSD, convert_CSD, [("spherical_harmonics_image", "in_file")]),
+                    (convert_CSD, outputnode, [("converted", "DWI")])
+                    # (mrtrix_CSD,outputnode,[('spherical_harmonics_image','DWI')])
+                ]
+            )
+            # fmt:on
+
+
     else:
         # fmt:off
         flow.connect(
diff --git a/cmp/stages/preprocessing/preprocessing.py b/cmp/stages/preprocessing/preprocessing.py
index 2465ee8a2..fe2007d7e 100644
--- a/cmp/stages/preprocessing/preprocessing.py
+++ b/cmp/stages/preprocessing/preprocessing.py
@@ -189,6 +189,7 @@ def __init__(self, bids_dir, output_dir):
             "dwi_brain_mask",
             "T1",
             "act_5TT",
+            "in_5tt_orig",
             "gmwmi",
             "brain",
             "brain_mask",
@@ -987,6 +988,7 @@ def create_workflow(self, flow, inputnode, outputnode):
                 [
                     (processing_input, mrtrix_5tt, [("aparc_aseg", "in_file")]),
                     (mrtrix_5tt, fs_mriconvert_5tt, [("out_file", "in_file")]),
+                    (mrtrix_5tt, outputnode, [("out_file", "in_5tt_orig")]),
                     (fs_mriconvert_5tt, outputnode, [("out_file", "act_5TT")]),
                 ]
             )
diff --git a/cmtklib/data/parcellation/lausanne2018/README.txt b/cmtklib/data/parcellation/lausanne2018/README.txt
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors1.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors1.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors2.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors2.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors3.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors3.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors4.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors4.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors5.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors5.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors6.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/Priors2/priors6.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellum.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellum.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumExtractionMask.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumExtractionMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumMask.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumProbabilityMask.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumProbabilityMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumRegistrationMask.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_BrainCerebellumRegistrationMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_glm_4labelsJointFusion.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_glm_4labelsJointFusion.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_glm_6labelsJointFusion.nii.gz b/cmtklib/data/segmentation/ants_MICCAI2012_multi-atlas_challenge_data/T_template0_glm_6labelsJointFusion.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors1.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors1.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors2.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors2.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors3.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors3.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors4.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors4.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors5.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors5.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/Priors/priors6.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/Priors/priors6.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template0.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template0.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template1.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template1.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template2-malf.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template2-malf.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template2-malf_6Labels.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template2-malf_6Labels.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template2.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template2.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template2_BrainCerebellum.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template2_BrainCerebellum.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template3.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template3.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_templateExtractionMask.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_templateExtractionMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_templateProbabilityMask.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_templateProbabilityMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumExtractionMask.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumExtractionMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumMask.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumProbabilityMask.nii.gz b/cmtklib/data/segmentation/ants_template_IXI/T_template_BrainCerebellumProbabilityMask.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/thalamus2018/Thalamus_Nuclei-HCP-4DSPAMs.nii.gz b/cmtklib/data/segmentation/thalamus2018/Thalamus_Nuclei-HCP-4DSPAMs.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/data/segmentation/thalamus2018/mni_icbm152_t1_tal_nlin_sym_09b_hires_1.nii.gz b/cmtklib/data/segmentation/thalamus2018/mni_icbm152_t1_tal_nlin_sym_09b_hires_1.nii.gz
old mode 100755
new mode 100644
diff --git a/cmtklib/diffusion.py b/cmtklib/diffusion.py
index 6e8c1e2f0..053a45c77 100644
--- a/cmtklib/diffusion.py
+++ b/cmtklib/diffusion.py
@@ -8,7 +8,7 @@
 
 import os
 import subprocess
-
+import pandas
 import nibabel as nib
 import numpy as np
 import nibabel.trackvis as tv
diff --git a/cmtklib/interfaces/mrtrix3.py b/cmtklib/interfaces/mrtrix3.py
index 9bba12a02..980d87b56 100644
--- a/cmtklib/interfaces/mrtrix3.py
+++ b/cmtklib/interfaces/mrtrix3.py
@@ -1024,7 +1024,7 @@ def _gen_outfilename(self):
         return name + '_vector.mif'
 
 
-class EstimateResponseForSHInputSpec(CommandLineInputSpec):
+class EstimateResponseForSHSingleTissueInputSpec(CommandLineInputSpec):
     in_file = File(exists=True, argstr='%s', mandatory=True,
                    position=2, desc='Diffusion-weighted images')
 
@@ -1054,11 +1054,11 @@ class EstimateResponseForSHInputSpec(CommandLineInputSpec):
     debug = traits.Bool(argstr='-debug', desc='Display debugging messages.')
 
 
-class EstimateResponseForSHOutputSpec(TraitedSpec):
+class EstimateResponseForSHSingleTissueOutputSpec(TraitedSpec):
     response = File(exists=True, desc='Spherical harmonics image')
 
 
-class EstimateResponseForSH(CommandLine):
+class EstimateResponseForSHSingleTissue(CommandLine):
     """Estimates the fibre response function for use in spherical deconvolution using `dwi2response`.
 
     Example
@@ -1073,8 +1073,8 @@ class EstimateResponseForSH(CommandLine):
     """
 
     _cmd = 'dwi2response'
-    input_spec = EstimateResponseForSHInputSpec
-    output_spec = EstimateResponseForSHOutputSpec
+    input_spec = EstimateResponseForSHSingleTissueInputSpec
+    output_spec = EstimateResponseForSHSingleTissueOutputSpec
 
     def _list_outputs(self):
         outputs = self.output_spec().get()
@@ -1092,16 +1092,102 @@ def _gen_outfilename(self):
         return name + '_ER.mif'
 
 
-class ConstrainedSphericalDeconvolutionInputSpec(CommandLineInputSpec):
-    algorithm = traits.Enum('csd', argstr='%s', mandatory=True, position=-4,
-                            desc='use CSD algorithm for FOD estimation')
+class EstimateResponseForSHMultiTissueInputSpec(CommandLineInputSpec):
+    in_file = File(exists=True, argstr='%s', mandatory=True,
+                   position=2, desc='Diffusion-weighted images')
+    
+    in_5tt_file = File(exists=True, argstr='%s', mandatory=True,
+                   position=3, desc='Diffusion-weighted images')
+
+    algorithm = traits.Enum('dhollander', 'fa', 'manual', 'msmt_5tt', 'tax', 'tournier', argstr='%s', position=1,
+                            desc='Select the algorithm to be used to derive the response function; '
+                            'additional details and options become available once an algorithm is nominated. '
+                            'Options are: dhollander, fa, manual, msmt_5tt, tax, tournier')
+
+    out_filename = File(genfile=True, argstr='%s',
+                        position=4, desc='Output filename')
+    out_gm_filename = File(genfile=True, argstr='%s',
+                        position=5, desc='Output filename')
+    out_csf_filename = File(genfile=True, argstr='%s',
+                        position=6, desc='Output filename')
+
+    encoding_file = File(exists=True, argstr='-grad %s', mandatory=True, position=-2,
+                         desc='Gradient encoding, supplied as a 4xN text file with each line is in the format [ X Y Z b ], '
+                         'where [ X Y Z ] describe the direction of the applied gradient, and b gives the b-value in units (1000 s/mm^2). '
+                         'See FSL2MRTrix')
+    maximum_harmonic_order = traits.Int(argstr='-lmax %s', position=-3,
+                                        desc='set the maximum harmonic order for the output series. '
+                                        'By default, the program will use the highest possible lmax given the number of diffusion-weighted images.')
+    # normalise = traits.Bool(argstr='-normalise', desc='normalise the DW signal to the b=0 image')
+
+    quiet = traits.Bool(
+        argstr='-quiet', desc='Do not display information messages or progress status.')
+
+    debug = traits.Bool(argstr='-debug', desc='Display debugging messages.')
+
+
+class EstimateResponseForSHMultiTissueOutputSpec(TraitedSpec):
+    response = File(exists=True, desc='WM Spherical harmonics image')
+    response_gm = File(exists=True, desc='GM Spherical harmonics image')
+    response_csf = File(exists=True, desc='CSF Spherical harmonics image')
+
+
+class EstimateResponseForSHMultiTissue(CommandLine):
+    """Estimates the fibre response function for use in spherical deconvolution using `dwi2response`.
+
+    Example
+    -------
+    >>> import cmtklib.interfaces.mrtrix3 as mrt
+    >>> estresp = mrt.EstimateResponseForSH()
+    >>> estresp.inputs.in_file = 'dwi.mif'
+    >>> estresp.inputs.mask_image = 'dwi_WMProb.mif'
+    >>> estresp.inputs.encoding_file = 'encoding.txt'
+    >>> estresp.run()  # doctest: +SKIP
+
+    """
+
+    _cmd = 'dwi2response'
+    input_spec = EstimateResponseForSHMultiTissueInputSpec
+    output_spec = EstimateResponseForSHMultiTissueOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['response'] = op.abspath(self._gen_outfilename()[0])
+        outputs['response_gm'] = op.abspath(self._gen_outfilename()[1])
+        outputs['response_csf'] = op.abspath(self._gen_outfilename()[2])
+        return outputs
+
+    def _gen_filename(self, name):
+        if name is 'out_filename':
+            return self._gen_outfilename()[0]
+        elif name is 'out_gm_filename':
+            return self._gen_outfilename()[1] 
+        else:
+            return self._gen_outfilename()[2]
+
+    def _gen_outfilename(self):
+        _, name, _ = split_filename(self.inputs.in_file)
+        out_wm = name + '_ER.txt'
+        out_gm = name + '_ER_gm.txt'
+        out_csf = name + '_ER_csf.txt'
+        return  out_wm, out_gm, out_csf
+
+
+class ConstrainedSphericalDeconvolutionSingleTissueInputSpec(CommandLineInputSpec):
+    algorithm = traits.Enum('csd', 'msmt_csd', argstr='%s', mandatory=True, position=-4,
+                        desc='Select the CSD algorithm to be use for FOD estimation'
+                        'Options are: csd (single shell single tissue) or msmt_csd (multi-shell multi-tissues)')
+
+   # algorithm = traits.Enum('csd', argstr='%s', mandatory=True, position=-4,
+   #                         desc='Select the CSD algorithm to be use for FOD estimation'
+   #                         'Options are: csd (single shell single tissue) or msmt_csd (multi-shell multi-tissues)')
 
     in_file = File(exists=True, argstr='%s', mandatory=True,
                    position=-3, desc='diffusion-weighted image')
 
     response_file = File(exists=True, argstr='%s', mandatory=True, position=-2,
                          desc='the diffusion-weighted signal response function for a single fibre population (see EstimateResponse)')
-
+ 
     out_filename = File(genfile=True, argstr='%s',
                         position=-1, desc='Output filename')
 
@@ -1139,12 +1225,12 @@ class ConstrainedSphericalDeconvolutionInputSpec(CommandLineInputSpec):
     # normalise = traits.Bool(argstr='-normalise', position=3, desc="normalise the DW signal to the b=0 image")
 
 
-class ConstrainedSphericalDeconvolutionOutputSpec(TraitedSpec):
+class ConstrainedSphericalDeconvolutionSingleTissueOutputSpec(TraitedSpec):
     spherical_harmonics_image = File(
         exists=True, desc='Spherical harmonics image')
 
 
-class ConstrainedSphericalDeconvolution(CommandLine):
+class ConstrainedSphericalDeconvolutionSingleTissue(CommandLine):
     """Perform non-negativity constrained spherical deconvolution using `dwi2fod`.
 
     Note that this program makes use of implied symmetries in the diffusion profile.
@@ -1176,8 +1262,8 @@ class ConstrainedSphericalDeconvolution(CommandLine):
     """
 
     _cmd = 'dwi2fod'
-    input_spec = ConstrainedSphericalDeconvolutionInputSpec
-    output_spec = ConstrainedSphericalDeconvolutionOutputSpec
+    input_spec = ConstrainedSphericalDeconvolutionSingleTissueInputSpec
+    output_spec = ConstrainedSphericalDeconvolutionSingleTissueOutputSpec
 
     def _list_outputs(self):
         outputs = self.output_spec().get()
@@ -1196,6 +1282,143 @@ def _gen_outfilename(self):
         return name + '_CSD.mif'
 
 
+class ConstrainedSphericalDeconvolutionMultiTissueInputSpec(CommandLineInputSpec):
+    algorithm = traits.Enum('csd', 'msmt_csd', argstr='%s', mandatory=True, position=0,
+                        desc='Select the CSD algorithm to be use for FOD estimation'
+                        'Options are: csd (single shell single tissue) or msmt_csd (multi-shell multi-tissues)')
+
+   # algorithm = traits.Enum('csd', argstr='%s', mandatory=True, position=-4,
+   #                         desc='Select the CSD algorithm to be use for FOD estimation'
+   #                         'Options are: csd (single shell single tissue) or msmt_csd (multi-shell multi-tissues)')
+
+    in_file = File(exists=True, argstr='%s', mandatory=True,
+                   position=1, desc='diffusion-weighted image')
+
+    
+    response_file = File(exists=True, argstr='%s', mandatory=True, position=2,
+                         desc='the diffusion-weighted signal response function for a single fibre population WM (see EstimateResponse)')
+    
+    response_gm_file = File(exists=True, argstr='%s', mandatory=True, position=4,
+                         desc='the diffusion-weighted signal response function for a single fibre population GM (see EstimateResponse)')
+
+    response_csf_file = File(exists=True, argstr='%s', mandatory=True, position=6,
+                         desc='the diffusion-weighted signal response function for a single fibre population CSF (see EstimateResponse)')
+
+    out_filename = File(genfile=True, argstr='%s',
+                        position=3, desc='Output filename')
+
+    out_gm_filename = File(genfile=True, argstr='%s',
+                        position=5, desc='Output filename')
+    
+    out_csf_filename = File(genfile=True, argstr='%s',
+                        position=7, desc='Output filename')
+    
+    mask_image = File(exists=True, argstr='-mask %s', position=2,
+                      desc='only perform computation within the specified binary brain mask image')
+
+    encoding_file = File(exists=True, argstr='-grad %s', position=8,
+                         desc='Gradient encoding, supplied as a 4xN text file with each line is in the format [ X Y Z b ], '
+                         'where [ X Y Z ] describe the direction of the applied gradient, and b gives the b-value in units (1000 s/mm^2). '
+                         'See FSL2MRTrix')
+
+    filter_file = File(exists=True, argstr='-filter %s', position=-2,
+                       desc='a text file containing the filtering coefficients for each even harmonic order.'
+                            'the linear frequency filtering parameters used for the initial linear spherical deconvolution step (default = [ 1 1 1 0 0 ]).')
+
+    lambda_value = traits.Float(argstr='-norm_lambda %s',
+                                desc='the regularisation parameter lambda that controls the strength of the constraint (default = 1.0).')
+
+    maximum_harmonic_order = traits.Int(argstr='-lmax %s',
+                                        desc='set the maximum harmonic order for the output series. '
+                                        'By default, the program will use the highest possible lmax given the number of diffusion-weighted images.')
+
+    threshold_value = traits.Float(argstr='-threshold %s',
+                                   desc='the threshold below which the amplitude of the FOD is assumed to be zero, '
+                                   'expressed as a fraction of the mean value of the initial FOD (default = 0.1)')
+
+    iterations = traits.Int(argstr='-niter %s',
+                            desc='the maximum number of iterations to perform for each voxel (default = 50)')
+
+    directions_file = File(exists=True, argstr='-directions %s', position=-2,
+                           desc='a text file containing the [ el az ] pairs for the directions: '
+                           'Specify the directions over which to apply the non-negativity constraint '
+                           '(by default, the built-in 300 direction set is used)')
+
+    # normalise = traits.Bool(argstr='-normalise', position=3, desc="normalise the DW signal to the b=0 image")
+
+
+class ConstrainedSphericalDeconvolutionMultiTissueOutputSpec(TraitedSpec):
+    
+    #CSD = File(exists=True, desc='WM Spherical harmonics image')
+    #CSD_gm = File(exists=True, desc='GM Spherical harmonics image')
+    #CSD_csf = File(exists=True, desc='CSF Spherical harmonics image')
+    spherical_harmonics_image = File(exists=True, desc='WM Spherical harmonics image')
+    gm_spherical_harmonics_image = File(exists=True, desc='GM Spherical harmonics image')
+    csf_spherical_harmonics_image = File(exists=True, desc='CSF Spherical harmonics image')
+
+
+
+
+class ConstrainedSphericalDeconvolutionMultiTissue(CommandLine):
+    """Perform non-negativity constrained spherical deconvolution using `dwi2fod`.
+
+    Note that this program makes use of implied symmetries in the diffusion profile.
+    First, the fact the signal attenuation profile is real implies that it has conjugate symmetry,
+    i.e. Y(l,-m) = Y(l,m)* (where * denotes the complex conjugate). Second, the diffusion profile should be
+    antipodally symmetric (i.e. S(x) = S(-x)), implying that all odd l components should be zero.
+    Therefore, this program only computes the even elements. 	Note that the spherical harmonics equations used here
+    differ slightly from those conventionally used, in that the (-1)^m factor has been omitted. This should be taken
+    into account in all subsequent calculations. Each volume in the output image corresponds to a different spherical
+    harmonic component, according to the following convention:
+
+    * [0] Y(0,0)
+    * [1] Im {Y(2,2)}
+    * [2] Im {Y(2,1)}
+    * [3] Y(2,0)
+    * [4] Re {Y(2,1)}
+    * [5] Re {Y(2,2)}
+    * [6] Im {Y(4,4)}
+    * [7] Im {Y(4,3)}
+
+
+    Example
+    -------
+    >>> import cmtklib.interfaces.mrtrix3 as mrt
+    >>> csdeconv = mrt.ConstrainedSphericalDeconvolution()
+    >>> csdeconv.inputs.in_file = 'dwi.mif'
+    >>> csdeconv.inputs.encoding_file = 'encoding.txt'
+    >>> csdeconv.run()                                          # doctest: +SKIP
+    """
+
+    _cmd = 'dwi2fod'
+    input_spec = ConstrainedSphericalDeconvolutionMultiTissueInputSpec
+    output_spec = ConstrainedSphericalDeconvolutionMultiTissueOutputSpec
+
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs['spherical_harmonics_image'] = op.abspath(self._gen_outfilename()[0])
+        outputs['gm_spherical_harmonics_image'] = op.abspath(self._gen_outfilename()[1])
+        outputs['csf_spherical_harmonics_image'] = op.abspath(self._gen_outfilename()[2])
+        return outputs
+
+    def _gen_filename(self, name):
+        if name is 'out_filename':
+            return self._gen_outfilename()[0]
+        elif name is 'out_gm_filename':
+            return self._gen_outfilename()[1] 
+        else:
+            return self._gen_outfilename()[2]
+
+    def _gen_outfilename(self):
+        _, name, _ = split_filename(self.inputs.in_file)
+        out_wm = name + '_CSD.mif'
+        out_gm = name + '_CSD_gm.mif'
+        out_csf = name + '_CSD_csf.mif'
+        return  out_wm, out_gm, out_csf
+
+
+
 class Generate5ttInputSpec(CommandLineInputSpec):
     algorithm = traits.Enum(
         'fsl',
diff --git a/conda/environment.yml b/conda/environment.yml
index 831da4c3b..960c021f6 100644
--- a/conda/environment.yml
+++ b/conda/environment.yml
@@ -37,3 +37,4 @@ dependencies:
           - networkx==2.6.3
           - pydicom==2.2.2
           - pycartool==0.1.1
+          - pyface==7.4.4
diff --git a/docker/environment.yml b/docker/environment.yml
index 710fd1f37..bf9ecd33b 100644
--- a/docker/environment.yml
+++ b/docker/environment.yml
@@ -1,4 +1,4 @@
-name: py37cmp-core
+name: py38cmp-core
 
 channels:
 - conda-forge
@@ -6,10 +6,10 @@ channels:
 - mrtrix3
 
 dependencies:
-- python=3.7
+- python=3.8
 - pip=20.1.1
 - ants=2.2.0
-- numpy=1.19.2
+- numpy=1.21.0
 - matplotlib=3.2.2
 - graphviz=2.40.1
 - traits=6.2.0
@@ -24,7 +24,7 @@ dependencies:
 - mrtrix3=3.0.2
 
 - pip:
-    - coverage==5.2
+    - coverage==5.2.1
     - duecredit==0.8.0
     - mne==0.24.1
     - mne_connectivity==0.3
@@ -36,3 +36,4 @@ dependencies:
     - networkx==2.6.3
     - pydicom==2.2.2
     - pycartool==0.1.1
+    - pandas==1.4.0
diff --git a/docker/requirements.txt b/docker/requirements.txt
index c6597a894..2e159a48c 100644
--- a/docker/requirements.txt
+++ b/docker/requirements.txt
@@ -1,4 +1,4 @@
-coverage==5.2
+coverage==5.2.1
 duecredit==0.8.0
 nibabel==3.2.2
 mne==0.24.1
diff --git a/docs/LICENSE.rst b/docs/LICENSE.rst
index 98f827e87..e04326db3 100644
--- a/docs/LICENSE.rst
+++ b/docs/LICENSE.rst
@@ -5,7 +5,7 @@ BSD 3-Clause License
 
 Copyright (C) 2009-2022, Ecole Polytechnique Fédérale de Lausanne (EPFL) and
 Hospital Center and University of Lausanne (UNIL-CHUV), Switzerland, & Contributors,
-All rights reserved.
+All rights reserved.  
 
 Redistribution and use in source and binary forms, with or without
 modification, are permitted provided that the following conditions are met:
diff --git a/docs/images/anatomical_pipeline_graph.png b/docs/images/anatomical_pipeline_graph.png
old mode 100755
new mode 100644
diff --git a/docs/images/diffusion_pipeline_graph.png b/docs/images/diffusion_pipeline_graph.png
old mode 100755
new mode 100644
diff --git a/docs/images/registration_fs.png b/docs/images/registration_fs.png
deleted file mode 100644
index 5294d6bd8..000000000
Binary files a/docs/images/registration_fs.png and /dev/null differ
diff --git a/docs/make.bat b/docs/make.bat
old mode 100644
new mode 100755
diff --git a/docs/notebooks/EEG_tutorial_utils.py b/docs/notebooks/EEG_tutorial_utils.py
old mode 100755
new mode 100644
diff --git a/docs/support.rst b/docs/support.rst
index 65462cb40..be650d528 100644
--- a/docs/support.rst
+++ b/docs/support.rst
@@ -3,6 +3,4 @@
 Support, Bugs and New Feature Requests
 =======================================
 
-If you need any support or have any questions, you can post to the `CMTK-users group <http://groups.google.com/group/cmtk-users>`_.
-
-All bugs, concerns and enhancement requests for this software are managed on GitHub and can be submitted at `https://github.com/connectomicslab/connectomemapper3/issues <https://github.com/connectomicslab/connectomemapper3/issues>`_. (See :ref:`Contribute to Connectome Mapper <contributing>` for more details)
+If you need any support or have any questions, or for all bugs, concerns and enhancement requests for this software are managed on GitHub and can be submitted at `https://github.com/connectomicslab/connectomemapper3/issues <https://github.com/connectomicslab/connectomemapper3/issues>`_. (See :ref:`Contribute to Connectome Mapper <contributing>` for more details)
diff --git a/scripts/others/create_multishell_dwi.py b/scripts/others/create_multishell_dwi.py
old mode 100755
new mode 100644
diff --git a/scripts/utils/build_bidsapp.sh b/scripts/utils/build_bidsapp.sh
old mode 100755
new mode 100644
index cb78bf932..8c2ca836c
--- a/scripts/utils/build_bidsapp.sh
+++ b/scripts/utils/build_bidsapp.sh
@@ -16,6 +16,7 @@
 UTILSDIR=$(cd "$(dirname "$0")"; pwd)
 SCRIPTSDIR="$(dirname "$UTILSDIR")"
 BASEDIR="$(dirname "$SCRIPTSDIR")"
+cd $BASEDIR
 
 # Get the current date and time
 CMP_BUILD_DATE=$(date -u +"%Y-%m-%dT%H:%M:%SZ")


Test patch:
diff --git a/conda/environment_test.yml b/conda/environment_test.yml
new file mode 100644
index 000000000..812f8405c
--- /dev/null
+++ b/conda/environment_test.yml
@@ -0,0 +1,40 @@
+name: test_cmp
+
+channels:
+    - conda-forge
+    - mrtrix3
+
+dependencies:
+    - python=3.9
+    - pip=22.2
+    - graphviz=2.50.0
+    - dipy=1.3.0
+    - fsleyes=1.3.3
+    - indexed_gzip
+    - mrtrix3=3.0.3
+    - git-annex=10.20220724
+    - qt=5.15.4
+    - pyqt=5.15.4
+    - pyqt5-sip=12.9.0
+
+    - pip:
+          - numpy==1.21.6
+          - matplotlib==3.5.2
+          - configparser==5.2.0
+          - nibabel==3.2.2
+          - nipype==1.8.0
+          - traits==6.3.2
+          - traitsui==7.2.0
+          - duecredit==0.9.1
+          - obspy==1.2.2
+          - mne==0.24.1
+          - mne_connectivity==0.3
+          - datalad[full]==0.17.2
+          - datalad-container==1.1.6
+          - datalad-neuroimaging==0.3.1
+          - pybids==0.14.0
+          - statsmodels==0.13.1
+          - networkx==2.6.3
+          - pydicom==2.2.2
+          - pycartool==0.1.1
+          - pyface==7.4.4