Merge pull request #135 from GalKepler/diffusion

Diffusion workflow

pyproject.toml

@@ -37,7 +37,12 @@ sdcflows = "^2.5.1"
 tox = "^4.16.0"
 doc8 = "^1.1.1"
 pytest-cov = "^5.0.0"
+numpy = "^1.26.0"
 
+neuromaps = "^0.0.5"
+dipy = "^1.9.0"
+cvxpy = "^1.5.2"
+fury = "^0.10.0"
 [tool.poetry.dev-dependencies]
 coverage = "^7.5.4"  # testing
 mypy = "^1.10.0"  # linting

src/kepost/config.py

@@ -421,8 +421,12 @@ class workflow(_Config):
     """Threshold for the probabilistic segmentation of the gray matter."""
     atlases: list = ["all"]
     """Parcellation atlas(es) to use for the parcellation step. Available atlases are: `all`, `fan2016`, `huang2022`, `schaefer2018_{n_regions}_{n_networks}`."""
+    tensor_max_bval = 1000
+    """Maximum b-value to consider for tensor estimation."""
     dipy_reconstruction_method = "NLLS"
     """Reconstruction method to use for the estimation of tensor-derived parameters using dipy."""
+    dipy_reconstruction_sigma = None
+    """Sigma parameter for the RESTORE algorithm. If none provided, sigma will be estimated."""
 
 
 class loggers:

src/kepost/interfaces/bids/utils.py

@@ -50,3 +50,20 @@ def collect_data(
         )
 
     return subj_data, session_data
+
+
+def gen_acq_label(max_bval: int) -> str:
+    """
+    Generate the acquisition label
+
+    Parameters
+    ----------
+    max_bval : int
+        The maximum bval
+
+    Returns
+    -------
+    str
+        The acquisition label
+    """
+    return f"shell{int(max_bval)}"

src/kepost/interfaces/dipy/__init__.py

@@ -0,0 +1 @@
+from kepost.interfaces.dipy.reconst import ReconstDTI  # noqa: F401

src/kepost/interfaces/dipy/reconst.py

@@ -0,0 +1,84 @@
+from nipype.interfaces.base import File, TraitedSpec, isdefined, traits
+from nipype.interfaces.dipy.base import (
+    DipyBaseInterfaceInputSpec,
+    DipyDiffusionInterface,
+)
+
+
+class ReconstDTIInputSpec(DipyBaseInterfaceInputSpec):
+    mask_file = File(exists=True, desc="An optional white matter mask")
+    fit_method = traits.Str(desc="The method to fit the tensor", default="WLS")
+    sigma = traits.Float(desc="The standard deviation of the noise")
+
+
+class ReconstDTIOutputSpec(TraitedSpec):
+    tensor_file = File(exists=True, desc="The output tensor file")
+    fa_file = File(exists=True, desc="The output fractional anisotropy file")
+    ga_file = File(exists=True, desc="The output geodesic anisotropy file")
+    rgb_file = File(exists=True, desc="The output RGB file")
+    md_file = File(exists=True, desc="The output mean diffusivity file")
+    ad_file = File(exists=True, desc="The output axial diffusivity file")
+    rd_file = File(exists=True, desc="The output radial diffusivity file")
+    mode_file = File(exists=True, desc="The output mode file")
+    evec_file = File(exists=True, desc="The output eigenvectors file")
+    eval_file = File(exists=True, desc="The output eigenvalues file")
+
+
+class ReconstDTI(DipyDiffusionInterface):
+    """
+    Calculates the diffusion tensor model parameters
+    Example
+    -------
+    >>> import nipype.interfaces.dipy as dipy
+    >>> dti = dipy.DTI()
+    >>> dti.inputs.in_file = 'diffusion.nii'
+    >>> dti.inputs.in_bvec = 'bvecs'
+    >>> dti.inputs.in_bval = 'bvals'
+    >>> dti.run()                                   # doctest: +SKIP
+    """
+
+    input_spec = ReconstDTIInputSpec
+    output_spec = ReconstDTIOutputSpec
+
+    def _run_interface(self, runtime):
+        from dipy.workflows.reconst import ReconstDtiFlow
+
+        flow = ReconstDtiFlow()
+        flow.run(
+            input_files=self.inputs.in_file,
+            bvalues_files=self.inputs.in_bval,
+            bvectors_files=self.inputs.in_bvec,
+            mask_files=self.inputs.mask_file,
+            fit_method=self.inputs.fit_method,
+            sigma=self.inputs.sigma if isdefined(self.inputs.sigma) else None,
+            out_tensor=self._gen_filename("tensor"),
+            out_fa=self._gen_filename("fa"),
+            out_ga=self._gen_filename("ga"),
+            out_rgb=self._gen_filename("rgb"),
+            out_md=self._gen_filename("md"),
+            out_ad=self._gen_filename("ad"),
+            out_rd=self._gen_filename("rd"),
+            out_mode=self._gen_filename("mode"),
+            out_evec=self._gen_filename("evec"),
+            out_eval=self._gen_filename("eval"),
+        )
+
+        return runtime
+
+    def _list_outputs(self):
+        outputs = self._outputs().get()
+        for metric in [
+            "tensor",
+            "fa",
+            "ga",
+            "rgb",
+            "md",
+            "ad",
+            "rd",
+            "mode",
+            "evec",
+            "eval",
+        ]:
+            outputs[f"{metric}_file"] = self._gen_filename(metric)
+
+        return outputs

src/kepost/interfaces/mrtrix3/__init__.py

@@ -0,0 +1 @@
+from kepost.interfaces.mrtrix3.utils import MRConvert  # noqa: F401

src/kepost/interfaces/mrtrix3/utils.py

@@ -0,0 +1,111 @@
+import os.path as op
+
+from nipype.interfaces.base import File, TraitedSpec, traits
+from nipype.interfaces.mrtrix3.base import MRTrix3Base, MRTrix3BaseInputSpec
+
+
+class MRConvertInputSpec(MRTrix3BaseInputSpec):
+    in_file = File(
+        exists=True,
+        argstr="%s",
+        mandatory=True,
+        position=-2,
+        desc="input image",
+    )
+    out_file = File(
+        "dwi.mif",
+        argstr="%s",
+        mandatory=True,
+        position=-1,
+        usedefault=True,
+        desc="output image",
+    )
+    coord = traits.List(
+        traits.Int,
+        sep=" ",
+        argstr="-coord %s",
+        desc="extract data at the specified coordinates",
+    )
+    vox = traits.List(
+        traits.Float,
+        sep=",",
+        argstr="-vox %s",
+        desc="change the voxel dimensions",
+    )
+    axes = traits.List(
+        traits.Int,
+        sep=",",
+        argstr="-axes %s",
+        desc="specify the axes that will be used",
+    )
+    scaling = traits.List(
+        traits.Float,
+        sep=",",
+        argstr="-scaling %s",
+        desc="specify the data scaling parameter",
+    )
+    json_import = File(
+        exists=True,
+        argstr="-json_import %s",
+        mandatory=False,
+        desc="import data from a JSON file into header key-value pairs",
+    )
+    json_export = File(
+        exists=False,
+        argstr="-json_export %s",
+        mandatory=False,
+        desc="export data from an image header key-value pairs into a JSON file",
+    )
+    out_mrtrix_grad = File(
+        exists=False,
+        argstr="-export_grad_mrtrix %s",
+        mandatory=False,
+        desc="export the gradient table in MRtrix format",
+    )
+
+
+class MRConvertOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="output image")
+    json_export = File(
+        exists=True,
+        desc="exported data from an image header key-value pairs in a JSON file",
+    )
+    out_bvec = File(exists=True, desc="export bvec file in FSL format")
+    out_bval = File(exists=True, desc="export bvec file in FSL format")
+    out_mrtrix_grad = File(
+        exists=True, desc="export the gradient table in MRtrix format"
+    )
+
+
+class MRConvert(MRTrix3Base):  # pylint: disable=abstract-method
+    """
+    Perform conversion between different file types and optionally extract a subset of the input image.
+
+    Example
+    -------
+
+    >>> import nipype.interfaces.mrtrix3 as mrt
+    >>> mrconvert = mrt.MRConvert()
+    >>> mrconvert.inputs.in_file = 'dwi.nii.gz'
+    >>> mrconvert.inputs.grad_fsl = ('bvecs', 'bvals')
+    >>> mrconvert.cmdline                             # doctest: +ELLIPSIS
+    'mrconvert -fslgrad bvecs bvals dwi.nii.gz dwi.mif'
+    >>> mrconvert.run()                               # doctest: +SKIP
+    """
+
+    _cmd = "mrconvert"
+    input_spec = MRConvertInputSpec
+    output_spec = MRConvertOutputSpec
+
+    def _list_outputs(self):
+        outputs = self.output_spec().get()
+        outputs["out_file"] = op.abspath(self.inputs.out_file)
+        if self.inputs.json_export:
+            outputs["json_export"] = op.abspath(self.inputs.json_export)
+        if self.inputs.out_bvec:
+            outputs["out_bvec"] = op.abspath(self.inputs.out_bvec)
+        if self.inputs.out_bval:
+            outputs["out_bval"] = op.abspath(self.inputs.out_bval)
+        if self.inputs.out_mrtrix_grad:
+            outputs["out_mrtrix_grad"] = op.abspath(self.inputs.out_mrtrix_grad)
+        return outputs

src/kepost/workflows/anatomical/__init__.py

@@ -0,0 +1 @@
+from kepost.workflows.anatomical.anatomical import init_anatomical_wf  # noqa: F401

src/kepost/workflows/anatomical/anatomical.py

@@ -7,12 +7,12 @@
 
 from kepost import config
 from kepost.atlases.utils import get_atlas_properties
-from kepost.workflows.anatomical.procedures.crop_to_gm import init_gm_cropping_wf
-from kepost.workflows.anatomical.procedures.derivatives import init_derivatives_wf
-from kepost.workflows.anatomical.procedures.five_tissue_type import (
+from kepost.workflows.anatomical.procedures import (
+    init_derivatives_wf,
     init_five_tissue_type_wf,
+    init_gm_cropping_wf,
+    init_registration_wf,
 )
-from kepost.workflows.anatomical.procedures.register_atlas import init_registration_wf
 
 
 def init_anatomical_wf(
@@ -48,6 +48,7 @@ def init_anatomical_wf(
     outputnode = pe.Node(
         interface=niu.IdentityInterface(
             fields=[
+                "atlas_name",
                 "whole_brain_parcellation",
                 "gm_cropped_parcellation",
             ]
@@ -64,6 +65,13 @@ def init_anatomical_wf(
     )
     workflow.connect(
         [
+            (
+                inputnode,
+                outputnode,
+                [
+                    ("atlas_name", "atlas_name"),
+                ],
+            ),
             (
                 atlases_node,
                 get_atlas_info_node,

src/kepost/workflows/anatomical/procedures/__init__.py

@@ -0,0 +1,12 @@
+from kepost.workflows.anatomical.procedures.crop_to_gm import (  # noqa: F401
+    init_gm_cropping_wf,
+)
+from kepost.workflows.anatomical.procedures.derivatives import (  # noqa: F401
+    init_derivatives_wf,
+)
+from kepost.workflows.anatomical.procedures.five_tissue_type import (  # noqa: F401
+    init_five_tissue_type_wf,
+)
+from kepost.workflows.anatomical.procedures.register_atlas import (  # noqa: F401
+    init_registration_wf,
+)