Merge pull request #2338 from matteomancini/enh-mrtrix3

[ENH] Extended MRtrix3 interface

Author: mgxd

.zenodo.json

@@ -563,6 +563,11 @@
       "affiliation": "MIT, HMS",
       "name": "Ghosh, Satrajit",
       "orcid": "0000-0002-5312-6729"
+    },
+    {
+      "affiliation": "University College London",
+      "name": "Mancini, Matteo",
+      "orcid": "0000-0001-7194-4568"
     }
   ],
   "keywords": [

nipype/interfaces/mrtrix3/__init__.py

@@ -4,7 +4,7 @@
 # -*- coding: utf-8 -*-
 
 from .utils import (Mesh2PVE, Generate5tt, BrainMask, TensorMetrics,
-                    ComputeTDI, TCK2VTK)
+                    ComputeTDI, TCK2VTK, MRMath, MRConvert, DWIExtract)
 from .preprocess import ResponseSD, ACTPrepareFSL, ReplaceFSwithFIRST
 from .tracking import Tractography
 from .reconst import FitTensor, EstimateFOD

nipype/interfaces/mrtrix3/preprocess.py

@@ -16,87 +16,47 @@
 import os.path as op
 
 from ..base import (CommandLineInputSpec, CommandLine, traits, TraitedSpec,
-                    File, isdefined)
+                    File, isdefined, Undefined)
 from .base import MRTrix3BaseInputSpec, MRTrix3Base
 
 
 class ResponseSDInputSpec(MRTrix3BaseInputSpec):
-    in_file = File(exists=True, argstr='%s', mandatory=True, position=-2,
-                   desc='input diffusion weighted images')
-
-    out_file = File(
-        'response.txt', argstr='%s', mandatory=True, position=-1,
-        usedefault=True, desc='output file containing SH coefficients')
-
-    # DW Shell selection options
-    shell = traits.List(traits.Float, sep=',', argstr='-shell %s',
-                        desc='specify one or more dw gradient shells')
+    algorithm = traits.Enum('msmt_5tt','dhollander','tournier','tax', argstr='%s', position=-6,
+        mandatory=True, desc='response estimation algorithm (multi-tissue)')
+    in_file = File(exists=True, argstr='%s', position=-5,
+        mandatory=True, desc='input DWI image')
+    mtt_file = File(argstr='%s', position=-4, desc='input 5tt image')
+    wm_file = File('wm.txt', argstr='%s', position=-3, usedefault=True,
+        desc='output WM response text file')
+    gm_file = File(argstr='%s', position=-2, desc='output GM response text file')
+    csf_file = File(argstr='%s', position=-1, desc='output CSF response text file')
     in_mask = File(exists=True, argstr='-mask %s',
-                   desc='provide initial mask image')
+        desc='provide initial mask image')
     max_sh = traits.Int(8, argstr='-lmax %d',
-                        desc='maximum harmonic degree of response function')
-    out_sf = File('sf_mask.nii.gz', argstr='-sf %s',
-                  desc='write a mask containing single-fibre voxels')
-    test_all = traits.Bool(False, argstr='-test_all',
-                           desc='re-test all voxels at every iteration')
-
-    # Optimization
-    iterations = traits.Int(0, argstr='-max_iters %d',
-                            desc='maximum number of iterations per pass')
-    max_change = traits.Float(
-        argstr='-max_change %f',
-        desc=('maximum percentile change in any response function coefficient;'
-              ' if no individual coefficient changes by more than this '
-              'fraction, the algorithm is terminated.'))
-
-    # Thresholds
-    vol_ratio = traits.Float(
-        .15, argstr='-volume_ratio %f',
-        desc=('maximal volume ratio between the sum of all other positive'
-              ' lobes in the voxel and the largest FOD lobe'))
-    disp_mult = traits.Float(
-        1., argstr='-dispersion_multiplier %f',
-        desc=('dispersion of FOD lobe must not exceed some threshold as '
-              'determined by this multiplier and the FOD dispersion in other '
-              'single-fibre voxels. The threshold is: (mean + (multiplier * '
-              '(mean - min))); default = 1.0. Criterion is only applied in '
-              'second pass of RF estimation.'))
-    int_mult = traits.Float(
-        2., argstr='-integral_multiplier %f',
-        desc=('integral of FOD lobe must not be outside some range as '
-              'determined by this multiplier and FOD lobe integral in other'
-              ' single-fibre voxels. The range is: (mean +- (multiplier * '
-              'stdev)); default = 2.0. Criterion is only applied in second '
-              'pass of RF estimation.'))
+        desc='maximum harmonic degree of response function')
 
 
 class ResponseSDOutputSpec(TraitedSpec):
-    out_file = File(exists=True, desc='the output response file')
-    out_sf = File(desc=('mask containing single-fibre voxels'))
+    wm_file = File(argstr='%s', desc='output WM response text file')
+    gm_file = File(argstr='%s', desc='output GM response text file')
+    csf_file = File(argstr='%s', desc='output CSF response text file')
 
 
 class ResponseSD(MRTrix3Base):
 
     """
-    Generate an appropriate response function from the image data for
-    spherical deconvolution.
-
-    .. [1] Tax, C. M.; Jeurissen, B.; Vos, S. B.; Viergever, M. A. and
-      Leemans, A., Recursive calibration of the fiber response function
-      for spherical deconvolution of diffusion MRI data. NeuroImage,
-      2014, 86, 67-80
-
+    Estimate response function(s) for spherical deconvolution using the specified algorithm.
 
     Example
     -------
 
     >>> import nipype.interfaces.mrtrix3 as mrt
     >>> resp = mrt.ResponseSD()
     >>> resp.inputs.in_file = 'dwi.mif'
-    >>> resp.inputs.in_mask = 'mask.nii.gz'
+    >>> resp.inputs.algorithm = 'tournier'
     >>> resp.inputs.grad_fsl = ('bvecs', 'bvals')
     >>> resp.cmdline                               # doctest: +ELLIPSIS
-    'dwi2response -fslgrad bvecs bvals -mask mask.nii.gz dwi.mif response.txt'
+    'dwi2response -fslgrad bvecs bvals tournier dwi.mif wm.txt'
     >>> resp.run()                                 # doctest: +SKIP
     """
 
@@ -106,10 +66,11 @@ class ResponseSD(MRTrix3Base):
 
     def _list_outputs(self):
         outputs = self.output_spec().get()
-        outputs['out_file'] = op.abspath(self.inputs.out_file)
-
-        if isdefined(self.inputs.out_sf):
-            outputs['out_sf'] = op.abspath(self.inputs.out_sf)
+        outputs['wm_file'] = op.abspath(self.inputs.wm_file)
+        if self.inputs.gm_file != Undefined:
+            outputs['gm_file'] = op.abspath(self.inputs.gm_file)
+        if self.inputs.csf_file != Undefined:
+            outputs['csf_file'] = op.abspath(self.inputs.csf_file)
         return outputs
 
 

nipype/interfaces/mrtrix3/reconst.py

@@ -15,7 +15,7 @@
 
 import os.path as op
 
-from ..base import traits, TraitedSpec, File
+from ..base import traits, TraitedSpec, File, Undefined
 from .base import MRTrix3BaseInputSpec, MRTrix3Base
 
 
@@ -74,108 +74,55 @@ def _list_outputs(self):
 
 
 class EstimateFODInputSpec(MRTrix3BaseInputSpec):
-    in_file = File(exists=True, argstr='%s', mandatory=True, position=-3,
-                   desc='input diffusion weighted images')
-    response = File(
-        exists=True, argstr='%s', mandatory=True, position=-2,
-        desc=('a text file containing the diffusion-weighted signal response '
-              'function coefficients for a single fibre population'))
-    out_file = File(
-        'fods.mif', argstr='%s', mandatory=True, position=-1,
-        usedefault=True, desc=('the output spherical harmonics coefficients'
-                               ' image'))
+    algorithm = traits.Enum('csd','msmt_csd', argstr='%s', position=-8,
+        mandatory=True, desc='FOD algorithm')
+    in_file = File(exists=True, argstr='%s', position=-7,
+        mandatory=True, desc='input DWI image')
+    wm_txt = File(argstr='%s', position=-6,
+        mandatory=True, desc='WM response text file')
+    wm_odf = File('wm.mif', argstr='%s', position=-5, usedefault=True,
+        mandatory=True, desc='output WM ODF')
+    gm_txt = File(argstr='%s', position=-4, desc='GM response text file')
+    gm_odf = File('gm.mif', argstr='%s', position=-3, desc='output GM ODF')
+    csf_txt = File(argstr='%s', position=-2, desc='CSF response text file')
+    csf_odf = File('csf.mif', argstr='%s', position=-1, desc='output CSF ODF')
+    mask_file = File(exists=True, argstr='-mask %s', desc='mask image')
 
     # DW Shell selection options
     shell = traits.List(traits.Float, sep=',', argstr='-shell %s',
                         desc='specify one or more dw gradient shells')
-
-    # Spherical deconvolution options
     max_sh = traits.Int(8, argstr='-lmax %d',
                         desc='maximum harmonic degree of response function')
-    in_mask = File(exists=True, argstr='-mask %s',
-                   desc='provide initial mask image')
     in_dirs = File(
         exists=True, argstr='-directions %s',
         desc=('specify the directions over which to apply the non-negativity '
               'constraint (by default, the built-in 300 direction set is '
               'used). These should be supplied as a text file containing the '
               '[ az el ] pairs for the directions.'))
-    sh_filter = File(
-        exists=True, argstr='-filter %s',
-        desc=('the linear frequency filtering parameters used for the initial '
-              'linear spherical deconvolution step (default = [ 1 1 1 0 0 ]). '
-              'These should be supplied as a text file containing the '
-              'filtering coefficients for each even harmonic order.'))
-
-    neg_lambda = traits.Float(
-        1.0, argstr='-neg_lambda %f',
-        desc=('the regularisation parameter lambda that controls the strength'
-              ' of the non-negativity constraint'))
-    thres = traits.Float(
-        0.0, argstr='-threshold %f',
-        desc=('the threshold below which the amplitude of the FOD is assumed '
-              'to be zero, expressed as an absolute amplitude'))
-
-    n_iter = traits.Int(
-        50, argstr='-niter %d', desc=('the maximum number of iterations '
-                                      'to perform for each voxel'))
 
 
 class EstimateFODOutputSpec(TraitedSpec):
-    out_file = File(exists=True, desc='the output response file')
+    wm_odf = File(argstr='%s', desc='output WM ODF')
+    gm_odf = File(argstr='%s', desc='output GM ODF')
+    csf_odf = File(argstr='%s', desc='output CSF ODF')
 
 
 class EstimateFOD(MRTrix3Base):
 
     """
-    Convert diffusion-weighted images to tensor images
-
-    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.
-    The spherical harmonic coefficients are stored as follows. First, since
-    the signal attenuation profile is real, 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, only the
-    even elements are computed.
-
-    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. Each volume will correspond to the following:
-
-    volume 0: l = 0, m = 0
-    volume 1: l = 2, m = -2 (imaginary part of m=2 SH)
-    volume 2: l = 2, m = -1 (imaginary part of m=1 SH)
-    volume 3: l = 2, m = 0
-    volume 4: l = 2, m = 1 (real part of m=1 SH)
-    volume 5: l = 2, m = 2 (real part of m=2 SH)
-    etc...
-
-
+    Estimate fibre orientation distributions from diffusion data using spherical deconvolution
 
     Example
     -------
 
     >>> import nipype.interfaces.mrtrix3 as mrt
     >>> fod = mrt.EstimateFOD()
+    >>> fod.inputs.algorithm = 'csd'
     >>> fod.inputs.in_file = 'dwi.mif'
-    >>> fod.inputs.response = 'response.txt'
-    >>> fod.inputs.in_mask = 'mask.nii.gz'
+    >>> fod.inputs.wm_txt = 'wm.txt'
     >>> fod.inputs.grad_fsl = ('bvecs', 'bvals')
     >>> fod.cmdline                               # doctest: +ELLIPSIS
-    'dwi2fod -fslgrad bvecs bvals -mask mask.nii.gz dwi.mif response.txt\
- fods.mif'
+    'dwi2fod -fslgrad bvecs bvals csd dwi.mif wm.txt wm.mif'
     >>> fod.run()                                 # doctest: +SKIP
     """
 
@@ -185,5 +132,12 @@ class EstimateFOD(MRTrix3Base):
 
     def _list_outputs(self):
         outputs = self.output_spec().get()
-        outputs['out_file'] = op.abspath(self.inputs.out_file)
+        outputs['wm_odf'] = op.abspath(self.inputs.wm_odf)
+        if self.inputs.gm_odf != Undefined:
+            outputs['gm_odf'] = op.abspath(self.inputs.gm_odf)
+        if self.inputs.csf_odf != Undefined:
+            outputs['csf_odf'] = op.abspath(self.inputs.csf_odf)
         return outputs
+
+
+

nipype/interfaces/mrtrix3/tests/test_auto_DWIExtract.py

@@ -0,0 +1,63 @@
+# AUTO-GENERATED by tools/checkspecs.py - DO NOT EDIT
+from __future__ import unicode_literals
+from ..utils import DWIExtract
+
+
+def test_DWIExtract_inputs():
+    input_map = dict(args=dict(argstr='%s',
+    ),
+    bval_scale=dict(argstr='-bvalue_scaling %s',
+    ),
+    bzero=dict(argstr='-bzero',
+    ),
+    environ=dict(nohash=True,
+    usedefault=True,
+    ),
+    grad_file=dict(argstr='-grad %s',
+    ),
+    grad_fsl=dict(argstr='-fslgrad %s %s',
+    ),
+    ignore_exception=dict(deprecated='1.0.0',
+    nohash=True,
+    usedefault=True,
+    ),
+    in_bval=dict(),
+    in_bvec=dict(argstr='-fslgrad %s %s',
+    ),
+    in_file=dict(argstr='%s',
+    mandatory=True,
+    position=-2,
+    ),
+    nobzero=dict(argstr='-nobzero',
+    ),
+    nthreads=dict(argstr='-nthreads %d',
+    nohash=True,
+    ),
+    out_file=dict(argstr='%s',
+    mandatory=True,
+    position=-1,
+    ),
+    shell=dict(argstr='-shell %s',
+    sep=',',
+    ),
+    singleshell=dict(argstr='-singleshell',
+    ),
+    terminal_output=dict(deprecated='1.0.0',
+    nohash=True,
+    ),
+    )
+    inputs = DWIExtract.input_spec()
+
+    for key, metadata in list(input_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(inputs.traits()[key], metakey) == value
+
+
+def test_DWIExtract_outputs():
+    output_map = dict(out_file=dict(),
+    )
+    outputs = DWIExtract.output_spec()
+
+    for key, metadata in list(output_map.items()):
+        for metakey, value in list(metadata.items()):
+            assert getattr(outputs.traits()[key], metakey) == value