- add nthreads option

- update workflow docstrings
- move resizing back to artifact workflow
- update eddy workflow

Author: josephmje

dmriprep/cli.py

@@ -62,11 +62,18 @@
     multiple=True
 )
 # specific options for eddy
+@click.option(
+    '--nthreads',
+    default=1,
+    show_default=True,
+    help='Maximum number of threads across all processes',
+    type=int,
+)
 @click.option(
     '--omp_nthreads',
     default=1,
     show_default=True,
-    help='Maximum number of threads for eddy',
+    help='Maximum number of threads per process',
     type=int,
 )
 @click.option(
@@ -136,6 +143,7 @@ def main(
     output_resolution,
     bet_dwi,
     bet_mag,
+    nthreads,
     omp_nthreads,
     eddy_niter,
     synb0_dir,

dmriprep/workflows/base.py

@@ -31,7 +31,7 @@ def init_dmriprep_wf(
     synb0_dir
 ):
     """
-    This workflow organizes the execusion of dMRIPrep, with a sub-workflow for
+    This workflow organizes the execusion of dMRIprep, with a sub-workflow for
     each subject.
 
     .. workflow::
@@ -41,7 +41,19 @@ def init_dmriprep_wf(
     from collections import namedtuple
     from dmriprep.workflows.base import init_dmriprep_wf
     BIDSLayout = namedtuple('BIDSLayout', ['root'])
-    wf = init_dmriprep_wf(subject_list=['dmripreptest'], layout=BIDSLayout('.'), output_dir='.', work_dir='.', ignore=[], b0_thresh=5, output_resolution=(1, 1, 1), bet_dwi=0.3, bet_mag=0.3, omp_nthreads=1, synb0_dir='')
+    wf = init_dmriprep_wf(
+        subject_list=['dmripreptest'],
+        layout=BIDSLayout('.'),
+        output_dir='.',
+        work_dir='.',
+        ignore=[],
+        b0_thresh=5,
+        output_resolution=(1, 1, 1),
+        bet_dwi=0.3,
+        bet_mag=0.3,
+        omp_nthreads=1,
+        synb0_dir=''
+    )
 
     Parameters
 
@@ -137,7 +149,20 @@ def init_single_subject_wf(
     from collections import namedtuple
     from dmriprep.workflows.base import init_single_subject_wf
     BIDSLayout = namedtuple('BIDSLayout', ['root'])
-    wf = init_single_subject_wf(subject_id='test', name='single_subject_wf', layout=BIDSLayout, output_dir='.', work_dir='.', ignore=[], b0_thresh=5, output_resolution=(1, 1, 1), bet_dwi=0.3, bet_mag=0.3, omp_nthreads=1, synb0_dir='')
+    wf = init_single_subject_wf(
+        subject_id='test',
+        name='single_subject_wf',
+        layout=BIDSLayout,
+        output_dir='.',
+        work_dir='.',
+        ignore=[],
+        b0_thresh=5,
+        output_resolution=(1, 1, 1),
+        bet_dwi=0.3,
+        bet_mag=0.3,
+        omp_nthreads=1,
+        synb0_dir=''
+    )
 
     Parameters
 
@@ -167,14 +192,19 @@ def init_single_subject_wf(
             Direction in which synb0 derivatives are saved
 
     """
-
-    dwi_files = layout.get(
-        subject=subject_id,
-        datatype='dwi',
-        suffix='dwi',
-        extensions=['.nii', '.nii.gz'],
-        return_type='filename'
-    )
+    # for documentation purposes
+    if name in ('single_subject_wf', 'single_subject_dmripreptest_wf'):
+        dwi_files = ['/madeup/path/sub-01_dwi.nii.gz']
+    else:
+        sessions = layout.get_sessions(subject=subject_id)
+
+        dwi_files = layout.get(
+            subject=subject_id,
+            datatype='dwi',
+            suffix='dwi',
+            extensions=['.nii', '.nii.gz'],
+            return_type='filename'
+        )
 
     if not dwi_files:
         raise Exception(
@@ -184,7 +214,7 @@ def init_single_subject_wf(
 
     subject_wf = pe.Workflow(name=name)
 
-    dwi_files = group_dwis(dwi_files, concat_dwis)
+    dwi_files = group_dwis(dwi_files, sessions, concat_dwis)
 
     for dwi_file in dwi_files:
         entities = layout.parse_file_entities(dwi_file)
@@ -268,7 +298,7 @@ def init_single_subject_wf(
     return subject_wf
 
 
-def group_dwis(dwi_files, concat_dwis):
+def group_dwis(dwi_files, sessions, concat_dwis):
 
     all_dwis = []
 

dmriprep/workflows/dwi/__init__.py

@@ -7,7 +7,6 @@
 
 .. automodule:: dmriprep.workflows.dwi.base
 .. automodule:: dmriprep.workflows.dwi.artifacts
-.. automodule:: dmriprep.workflows.dwi.util
 .. automodule:: dmriprep.workflows.dwi.eddy
 .. automodule:: dmriprep.workflows.dwi.tensor
 .. automodule:: dmriprep.workflows.dwi.outputs
@@ -16,15 +15,13 @@
 
 from .base import init_dwi_preproc_wf
 from .artifacts import init_dwi_artifacts_wf
-from .util import init_dwi_resize_wf
 from .eddy import init_dwi_eddy_wf
 from .tensor import init_dwi_tensor_wf
 from .outputs import init_dwi_derivatives_wf
 
 __all__ = [
     'init_dwi_preproc_wf',
     'init_dwi_artifacts_wf',
-    'init_dwi_resize_wf',
     'init_dwi_eddy_wf',
     'init_dwi_tensor_wf',
     'init_dwi_derivatives_wf'

dmriprep/workflows/dwi/artifacts.py

@@ -1,8 +1,8 @@
 # -*- coding: utf-8 -*-
 
 """
-Artefact removal
-^^^^^^^^^^^^^^^^
+Artefact removal and resizing
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 .. autofunction:: init_dwi_artifacts_wf
 
@@ -12,7 +12,7 @@
 from nipype.interfaces import mrtrix3, utility as niu
 
 
-def init_dwi_artifacts_wf(ignore):
+def init_dwi_artifacts_wf(ignore, output_resolution):
     """
     This workflow performs denoising and unringing on the input dwi image.
 
@@ -25,7 +25,7 @@ def init_dwi_artifacts_wf(ignore):
         :simple_form: yes
 
         from dmriprep.workflows.dwi import init_dwi_artifacts_wf
-        wf = init_dwi_artifacts_wf(ignore=[])
+        wf = init_dwi_artifacts_wf(ignore=[], output_resolution=(1, 1, 1))
 
     **Parameters**
 
@@ -50,10 +50,14 @@ def init_dwi_artifacts_wf(ignore):
 
     outputnode = pe.Node(niu.IdentityInterface(fields=['out_file']), name='outputnode')
 
+    dwibuffer = pe.Node(niu.IdentityInterface(fields=['dwi_file']), name='dwibuffer')
+
     denoise = pe.Node(mrtrix3.DWIDenoise(), name='denoise')
 
     unring = pe.Node(mrtrix3.MRDeGibbs(), name='unring')
 
+    resize = pe.Node(mrtrix3.MRResize(voxel_size=output_resolution), name='resize')
+
     if ignore == ['denoise']:
         wf.connect([
             (inputnode, unring, [('dwi_file', 'in_file')]),

dmriprep/workflows/dwi/base.py

@@ -14,7 +14,6 @@
 from numba import cuda
 
 from .artifacts import init_dwi_artifacts_wf
-from .util import init_dwi_resize_wf
 from .tensor import init_dwi_tensor_wf
 
 from ..fieldmap.base import init_sdc_prep_wf
@@ -46,8 +45,8 @@ def init_dwi_preproc_wf(
         from dmriprep.workflows.dwi import init_dwi_preproc_wf
         BIDSLayout = namedtuple('BIDSLayout', ['root'])
         wf = init_dwi_preproc_wf(
-            subjectid=,
-            dwi_file=,
+            subjectid='dmripreptest',
+            dwi_file='/madeup/path/sub-01_dwi.nii.gz',
             metadata=,
             layout=BIDSLayout('.'),
             ignore=[],
@@ -56,7 +55,7 @@ def init_dwi_preproc_wf(
             bet_dwi=0.3,
             bet_mag=0.3,
             omp_nthreads=1,
-            synb0_dir='.'
+            synb0_dir=''
         )
 
     """
@@ -140,26 +139,13 @@ def init_dwi_preproc_wf(
 
     dwi_artifacts_wf = init_dwi_artifacts_wf(ignore)
 
-    if output_resolution:
-        dwi_resize_wf = init_dwi_resize_wf(output_resolution)
-
-        dwi_wf.connect(
-            [
-                (inputnode, dwi_artifacts_wf, [("dwi_file", "inputnode.dwi_file")]),
-                (dwi_artifacts_wf, dwi_resize_wf, [("outputnode.out_file", "inputnode.in_file")]),
-                (dwi_resize_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
-                (dwi_resize_wf, ecc, [("outputnode.out_file", "in_file")])
-            ]
-        )
-
-    else:
-        dwi_wf.connect(
-            [
-                (inputnode, dwi_artifacts_wf, [("dwi_file", "inputnode.dwi_file")]),
-                (dwi_artifacts_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
-                (dwi_artifacts_wf, ecc, [("outputnode.out_file", "in_file")])
-            ]
-        )
+    dwi_wf.connect(
+        [
+            (inputnode, dwi_artifacts_wf, [("dwi_file", "inputnode.dwi_file")]),
+            (dwi_artifacts_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
+            (dwi_artifacts_wf, ecc, [("outputnode.out_file", "in_file")])
+        ]
+    )
 
     def gen_index(in_file):
         """

dmriprep/workflows/dwi/eddy.py

@@ -13,46 +13,99 @@
 from numba import cuda
 
 
-def init_dwi_eddy_wf(omp_nthreads):
+def init_dwi_eddy_wf(omp_nthreads, fmap_type):
+    """
+    This workflow runs eddy on the input dwi image.
 
-    wf = pe.Workflow(name="dwi_eddy_wf")
+    .. workflow::
+        :graph2use: orig
+        :simple_form: yes
 
-    inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_file",
-                                                      "bvec_file",
-                                                      "bval_file",
-                                                      "mask_file",
-                                                      "acqp",
-                                                      "index"]),
-                       name="inputnode")
+        from dmriprep.workflows.dwi import init_dwi_eddy_wf
+        wf = init_dwi_eddy_wf(omp_nthreads=1)
 
-    outputnode = pe.Node(niu.IdentityInterface(fields=["out_file", "out_bvec"]),
-                         name="outputnode")
+    **Parameters**
+
+        omp_nthreads: int
+            Number of threads to run eddy
+
+    **Inputs**
+
+        dwi_file
+            dwi NIfTI file
+        bvec_file
+            bvec file
+        bval_file
+            bval file
+        mask_file
+            brain mask file
+        fieldmap_file
+            fieldmap file
+        topup_fieldcoef
+            topup file containing field coefficients
+        topup_movpar
+            topup movpar.txt file
+        acqp
+            acquisition parameters file
+        index
+            index file
+
+    **Outputs**
+
+        out_file
+            output eddy-corrected dwi image
+        out_bvec
+            output rotated bvecs after eddy correction
+
+    """
+
+    wf = pe.Workflow(name='dwi_eddy_wf')
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=['dwi_file',
+                                                      'bvec_file',
+                                                      'bval_file',
+                                                      'mask_file',
+                                                      'fieldmap_file',
+                                                      'topup_fieldcoef',
+                                                      'topup_movpar',
+                                                      'acqp',
+                                                      'index']),
+                       name='inputnode')
+
+    outputnode = pe.Node(niu.IdentityInterface(fields=['out_file', 'out_bvec']),
+                         name='outputnode')
 
     ecc = pe.Node(
-        fsl.Eddy(num_threads=1, repol=True, cnr_maps=True, residuals=True),
-        name="fsl_eddy",
+        fsl.Eddy(num_threads=omp_nthreads, repol=True, cnr_maps=True, residuals=True),
+        name='fsl_eddy',
     )
 
-    if omp_nthreads:
-        ecc.inputs.num_threads = omp_nthreads
-
-    # this doesn't work well
-    # recognizes that a computer is cuda capable but failed if cuda isn't loaded
     try:
         if cuda.gpus:
             ecc.inputs.use_cuda = True
     except:
         ecc.inputs.use_cuda = False
 
     wf.connect([
-        (inputnode, ecc, [("dwi_file", "in_file"),
-                          ("bval_file", "in_bval"),
-                          ("bvec_file", "in_bvec"),
-                          ("acqp", "in_acqp"),
-                          ("index", "in_index"),
-                          ("mask_file", "in_mask")]),
-        (ecc, outputnode, [("out_corrected", "out_file"),
-                           ("out_rotated_bvecs", "out_bvec")])
+        (inputnode, ecc, [('dwi_file', 'in_file'),
+                          ('bval_file', 'in_bval'),
+                          ('bvec_file', 'in_bvec'),
+                          ('acqp', 'in_acqp'),
+                          ('index', 'in_index'),
+                          ('mask_file', 'in_mask')]),
+        (ecc, outputnode, [('out_corrected', 'out_file'),
+                           ('out_rotated_bvecs', 'out_bvec')])
+        ])
+
+    if fmap_type == 'fieldmap':
+        wf.connect([
+            (inputnode, ecc, [('fieldmap_file', 'field')])
+        ])
+
+    if fmap_type == 'topup':
+        wf.connect([
+            (inputnode, ecc, [('topup_fieldcoef', 'in_topup_fieldcoef'),
+                              ('topup_movpar', 'in_topup_movpar')])
         ])
 
     return wf