add ants and brainsuite cli options and move eddy to separate workflow

Author: josephmje

dmriprep/cli.py

@@ -128,6 +128,20 @@
     type=click.Choice(['denoising', 'unringing', 'fieldmaps']),
     multiple=True
 )
+@click.option(
+    '--use_ants',
+    is_flag=True,
+    default=False,
+    help='Use ANTs SDC',
+    type=bool
+)
+@click.option(
+    '--use_brainsuite',
+    is_flag=True,
+    default=False,
+    help='Use BrainSuite SDC',
+    type=bool
+)
 @click.option(
     '--work_dir',
     '-w',
@@ -163,6 +177,8 @@ def main(
     nthreads,
     omp_nthreads,
     ignore,
+    use_ants,
+    use_syn,
     work_dir,
     synb0_dir,
     write_graph
@@ -247,6 +263,8 @@ def main(
         acqp_file=acqp_file,
         omp_nthreads=omp_nthreads,
         ignore=list(ignore),
+        use_ants=use_ants,
+        use_brainsuite=use_brainsuite,
         work_dir=work_dir,
         synb0_dir=synb0_dir
     )

dmriprep/workflows/base.py

@@ -31,6 +31,8 @@ def init_dmriprep_wf(
     acqp_file,
     omp_nthreads,
     ignore,
+    use_ants,
+    use_brainsuite,
     work_dir,
     synb0_dir
 ):
@@ -58,6 +60,8 @@ def init_dmriprep_wf(
         acqp_file='',
         omp_nthreads=1,
         ignore=[],
+        use_ants=False,
+        use_brainsuite=False,
         work_dir='.',
         synb0_dir=''
     )
@@ -114,6 +118,8 @@ def init_dmriprep_wf(
             acqp_file=acqp_file,
             omp_nthreads=omp_nthreads,
             ignore=ignore,
+            use_ants=use_ants,
+            use_brainsuite=use_brainsuite,
             work_dir=work_dir,
             synb0_dir=synb0_dir
         )
@@ -144,6 +150,8 @@ def init_single_subject_wf(
     omp_nthreads,
     acqp_file,
     ignore,
+    use_ants,
+    use_brainsuite,
     work_dir,
     synb0_dir
 ):
@@ -176,6 +184,8 @@ def init_single_subject_wf(
         acqp_file='',
         omp_nthreads=1,
         ignore=[],
+        use_ants=False,
+        use_brainsuite=False,
         work_dir='.',
         synb0_dir=''
     )
@@ -263,6 +273,8 @@ def init_single_subject_wf(
                 omp_nthreads=omp_nthreads,
                 acqp_file=acqp_file,
                 ignore=ignore,
+                use_ants=use_ants,
+                use_brainsuite=use_brainsuite,
                 synb0_dir=synb0_dir
             )
 
@@ -276,9 +288,10 @@ def init_single_subject_wf(
             inputspec.inputs.out_dir = os.path.abspath(output_dir)
 
             subject_wf.connect([
-                (dwi_concat_wf, dwi_preproc_wf, [('outputnode.dwi_file', 'inputnode.dwi_file'),
-                                                 ('outputnode.bvec_file', 'inputnode.bvec_file'),
-                                                 ('outputnode.bval_file', 'inputnode.bval_file')])
+                (dwi_concat_wf, dwi_preproc_wf,
+                    [('outputnode.dwi_file', 'inputnode.dwi_file'),
+                     ('outputnode.bvec_file', 'inputnode.bvec_file'),
+                     ('outputnode.bval_file', 'inputnode.bval_file')])
             ])
 
         else:
@@ -299,6 +312,8 @@ def init_single_subject_wf(
                 omp_nthreads=omp_nthreads,
                 acqp_file=acqp_file,
                 ignore=ignore,
+                use_ants=use_ants,
+                use_brainsuite=use_brainsuite,
                 synb0_dir=synb0_dir
             )
 
@@ -332,37 +347,22 @@ def init_single_subject_wf(
         ds_inputspec.inputs.output_folder = output_dir
         ds_inputspec.inputs.metadata = metadata
 
-        if session_id:
-            wf_name = 'sub_' + subject_id + '_ses_' + session_id + '_preproc_wf'
-        else:
-            wf_name = 'sub_' + subject_id + '_preproc_wf'
-        full_wf = pe.Workflow(name=wf_name)
-
-        full_wf.connect(
-            [
-                (
-                    dwi_preproc_wf,
-                    datasink_wf,
-                    [
-                        ('outputnode.out_dwi', 'inputnode.dwi'),
-                        ('outputnode.out_bval', 'inputnode.bval'),
-                        ('outputnode.out_bvec', 'inputnode.bvec'),
-                        ('outputnode.index', 'inputnode.index'),
-                        ('outputnode.acq_params', 'inputnode.acq_params'),
-                        ('outputnode.out_mask', 'inputnode.mask'),
-                        ('outputnode.out_b0_pre', 'inputnode.b0'),
-                        ('outputnode.out_b0_mask_pre', 'inputnode.b0_mask'),
-                        ('outputnode.out_eddy_quad_json', 'inputnode.eddy_quad_json'),
-                        ('outputnode.out_eddy_quad_pdf', 'inputnode.eddy_quad_pdf'),
-                        ('outputnode.out_dtifit_FA', 'inputnode.dtifit_FA'),
-                        ('outputnode.out_dtifit_V1', 'inputnode.dtifit_V1'),
-                        ('outputnode.out_dtifit_sse', 'inputnode.dtifit_sse'),
-                        ('outputnode.out_noise', 'inputnode.noise')
-                    ]
-                )
-            ]
-        )
-
-        subject_wf.add_nodes([full_wf])
+        subject_wf.connect([
+            (dwi_preproc_wf, datasink_wf,
+                [('outputnode.out_dwi', 'inputnode.dwi'),
+                 ('outputnode.out_bval', 'inputnode.bval'),
+                 ('outputnode.out_bvec', 'inputnode.bvec'),
+                 ('outputnode.index', 'inputnode.index'),
+                 ('outputnode.acq_params', 'inputnode.acq_params'),
+                 ('outputnode.out_mask', 'inputnode.mask'),
+                 ('outputnode.out_b0_pre', 'inputnode.b0'),
+                 ('outputnode.out_b0_mask_pre', 'inputnode.b0_mask'),
+                 ('outputnode.out_eddy_quad_json', 'inputnode.eddy_quad_json'),
+                 ('outputnode.out_eddy_quad_pdf', 'inputnode.eddy_quad_pdf'),
+                 ('outputnode.out_dtifit_FA', 'inputnode.dtifit_FA'),
+                 ('outputnode.out_dtifit_V1', 'inputnode.dtifit_V1'),
+                 ('outputnode.out_dtifit_sse', 'inputnode.dtifit_sse'),
+                 ('outputnode.out_noise', 'inputnode.noise')])
+        ])
 
     return subject_wf

dmriprep/workflows/dwi/base.py

@@ -10,11 +10,11 @@
 
 from nipype.pipeline import engine as pe
 from nipype.interfaces import fsl, mrtrix3, utility as niu
-from numba import cuda
 
 from .artifacts import init_dwi_artifacts_wf
+from .eddy import init_dwi_eddy_wf
 from .tensor import init_dwi_tensor_wf
-from ..fieldmap.base import init_sdc_prep_wf
+from ..fieldmap.base import init_sdc_wf
 
 FMAP_PRIORITY = {'epi': 0, 'fieldmap': 1, 'phasediff': 2, 'phase': 3, 'syn': 4}
 
@@ -32,6 +32,8 @@ def init_dwi_preproc_wf(
     acqp_file,
     omp_nthreads,
     ignore,
+    use_ants,
+    use_brainsuite,
     synb0_dir
 ):
     """
@@ -50,6 +52,37 @@ def init_dwi_preproc_wf(
 
     wf_name = _get_wf_name(dwi_file)
 
+    if isinstance(dwi_file, list):
+        ref_file = dwi_file[0]
+    else:
+        ref_file = dwi_file
+
+    # For doc building purposes
+    if not hasattr(layout, 'parse_file_entities'):
+        metadata = {
+            'PhaseEncodingDirection': 'j-',
+            'TotalReadoutTime': 0.05
+        }
+        fmaps = [{
+            'suffix': 'phasediff',
+            'phasediff': 'sub-01/fmap/sub-01_phasediff.nii.gz',
+            'magnitude1': 'sub-01/fmap/sub-01_magnitude1.nii.gz',
+            'magnitude2': 'sub-01/fmap/sub-01_magnitude2.nii.gz'
+        }]
+    else:
+        metadata = layout.get_metadata(ref_file)
+        # Find fieldmaps. Options: (epi|fieldmap|phasediff|phase1|phase2|syn)
+        fmaps = []
+        if 'fieldmaps' not in ignore:
+            for fmap in layout.get_fieldmap(ref_file, return_list=True):
+                fmap['metadata'] = layout.get_metadata(fmap[fmap['suffix']])
+                fmaps.append(fmap)
+
+        if (use_ants and not fmaps):
+            fmaps.append({'suffix': 'ants'})
+        if (use_brainsuite and not fmaps):
+            fmaps.append({'suffix': 'brainsuite'})
+
     dwi_wf = pe.Workflow(name=wf_name)
 
     # # If use_synb0 set, get synb0 from files
@@ -69,7 +102,7 @@ def init_dwi_preproc_wf(
             )
             fmaps.append(fmap)
 
-    sdc_wf = init_sdc_prep_wf(
+    dwi_sdc_wf = init_sdc_wf(
         subject_id,
         fmaps,
         metadata,
@@ -246,14 +279,16 @@ def b0_average(in_dwi, in_bval, b0_thresh, out_file=None):
         fsl.BET(frac=bet_dwi, mask=True, robust=True), name='bet_dwi_pre'
     )
 
-    ecc = pe.Node(
-        fsl.Eddy(num_threads=omp_nthreads, repol=True, cnr_maps=True, residuals=True),
-        name='fsl_eddy')
-    try:
-        if cuda.gpus:
-            ecc.inputs.use_cuda = True
-    except:
-        ecc.inputs.use_cuda = False
+    dwi_eddy_wf = init_dwi_eddy_wf()
+
+    # ecc = pe.Node(
+    #     fsl.Eddy(num_threads=omp_nthreads, repol=True, cnr_maps=True, residuals=True),
+    #     name='fsl_eddy')
+    # try:
+    #     if cuda.gpus:
+    #         ecc.inputs.use_cuda = True
+    # except:
+    #     ecc.inputs.use_cuda = False
 
     denoise_eddy = pe.Node(mrtrix3.DWIDenoise(), name='denoise_eddy')
 
@@ -286,44 +321,45 @@ def get_b0_mask_fn(b0_file):
             input_names=['b0_file'], output_names=['mask_file'], function=get_b0_mask_fn),
         name='getB0Mask')
 
-    tensor_wf = init_dwi_tensor_wf()
+    dwi_tensor_wf = init_dwi_tensor_wf()
 
     dwi_wf.connect([
         (inputnode, gen_idx, [('dwi_file', 'in_file')]),
         (inputnode, acqp, [('dwi_file', 'in_file'), ('dwi_meta', 'metadata')]),
         (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_artifacts_wf, dwi_eddy_wf, [('outputnode.out_file', 'inputnode.dwi_file')]),
         (inputnode, avg_b0_0, [('bval_file', 'in_bval')]),
         (avg_b0_0, bet_dwi0, [('out_file', 'in_file')]),
-        (inputnode, ecc, [('bval_file', 'in_bval'),
-                          ('bvec_file', 'in_bvec')]),
-        (bet_dwi0, ecc, [('mask_file', 'in_mask')]),
-        (gen_idx, ecc, [('out_file', 'in_index')]),
-        (ecc, denoise_eddy, [('out_corrected', 'in_file')]),
-        (ecc, bias_correct, [('out_corrected', 'in_file'),
-                             ('out_rotated_bvecs', 'in_bvec')]),
+        (inputnode, dwi_eddy_wf, [('bval_file', 'inputnode.bval_file'),
+                                  ('bvec_file', 'inputnode.bvec_file')]),
+        (bet_dwi0, dwi_eddy_wf, [('mask_file', 'inputnode.mask_file')]),
+        (gen_idx, dwi_eddy_wf, [('out_file', 'inputnode.index')]),
+        (dwi_eddy_wf, denoise_eddy, [('outputnode.out_file', 'in_file')]),
+        (dwi_eddy_wf, bias_correct, [('outputnode.out_file', 'in_file'),
+                                 ('outputnode.out_bvec', 'in_bvec')]),
         (inputnode, bias_correct, [('bval_file', 'in_bval')]),
         (bias_correct, fslroi, [('out_file', 'in_file')]),
         (fslroi, b0mask_node, [('roi_file', 'b0_file')]),
-        (ecc, eddy_quad, [(('out_corrected', get_path), 'base_name'),
-                          (('out_corrected', get_qc_path), 'output_dir'),
-                          ('out_rotated_bvecs', 'bvec_file')]),
+        (dwi_eddy_wf, eddy_quad, [(('outputnode.out_file', get_path), 'base_name'),
+                                  (('outputnode.out_file', get_qc_path), 'output_dir'),
+                                  ('outputnode.out_bvec', 'bvec_file')]),
         (inputnode, eddy_quad, [('bval_file', 'bval_file')]),
         (b0mask_node, eddy_quad, [('mask_file', 'mask_file')]),
         (gen_idx, eddy_quad, [('out_file', 'idx_file')]),
-        (inputnode, tensor_wf, [('bval_file', 'inputnode.bval_file')]),
-        (b0mask_node, tensor_wf, [('mask_file', 'inputnode.mask_file')]),
-        (ecc, tensor_wf, [('out_corrected', 'inputnode.dwi_file'),
-                          ('out_rotated_bvecs', 'inputnode.bvec_file')])
+        (inputnode, dwi_tensor_wf, [('bval_file', 'inputnode.bval_file')]),
+        (b0mask_node, dwi_tensor_wf, [('mask_file', 'inputnode.mask_file')]),
+        (dwi_eddy_wf, dwi_tensor_wf, [('outputnode.out_file', 'inputnode.dwi_file'),
+                                      ('outputnode.out_bvec', 'inputnode.bvec_file')])
     ])
 
     if acqp_file:
-        ecc.inputs.in_acqp = acqp_file
+        eddy_inputspec = dwi_eddy_wf.get_node('inputnode')
+        eddy_inputspec.inputs.acqp = acqp_file
         eddy_quad.inputs.param_file = acqp_file
     else:
         dwi_wf.connect([
-            (acqp, ecc, [('out_file', 'in_acqp')]),
+            (acqp, dwi_eddy_wf, [('out_file', 'inputnode.acqp')]),
             (acqp, eddy_quad, [('out_file', 'param_file')])
         ])
 
@@ -349,38 +385,39 @@ def get_b0_mask_fn(b0_file):
     #     ecc.inputs.in_acqp = acqp_file
     # else:
     # Decide what ecc will take: topup or fmap
-    fmaps.sort(key=lambda fmap: FMAP_PRIORITY[fmap['suffix']])
-    fmap = fmaps[0]
-    # Else If epi files detected
-    if fmap['suffix'] == 'epi':
-        dwi_wf.connect([
-            (sdc_wf, ecc, [('outputnode.out_topup', 'in_topup_fieldcoef'),
-                           ('outputnode.out_movpar', 'in_topup_movpar')])
-        ])
-    # Otherwise (fieldmaps)
-    else:
-        dwi_wf.connect([
-            (bet_dwi0, sdc_wf, [('out_file', 'inputnode.b0_stripped')]),
-            (sdc_wf, ecc, [('outputnode.out_fmap', 'field')])
-        ])
+    if fmaps:
+        fmaps.sort(key=lambda fmap: FMAP_PRIORITY[fmap['suffix']])
+        fmap = fmaps[0]
+        # Else If epi files detected
+        if fmap['suffix'] == 'epi':
+            dwi_wf.connect([
+                (dwi_sdc_wf, dwi_eddy_wf, [('outputnode.out_topup', 'inputnode.topup_fieldcoef'),
+                                           ('outputnode.out_movpar', 'inputnode.topup_movpar')])
+            ])
+        # Otherwise (fieldmaps)
+        else:
+            dwi_wf.connect([
+                (bet_dwi0, dwi_sdc_wf, [('out_file', 'inputnode.b0_stripped')]),
+                (dwi_sdc_wf, dwi_eddy_wf, [('outputnode.out_fmap', 'inputnode.fieldmap_file')])
+            ])
 
     dwi_wf.connect([
-        (ecc, outputnode, [('out_corrected', 'out_dwi')]),
+        (dwi_eddy_wf, outputnode, [('outputnode.out_file', 'out_dwi')]),
         (inputnode, outputnode, [('bval_file', 'out_bval')]),
-        (ecc, outputnode, [('out_rotated_bvecs', 'out_bvec')]),
+        (dwi_eddy_wf, outputnode, [('outputnode.out_bvec', 'out_bvec')]),
         (gen_idx, outputnode, [('out_file', 'index')]),
         (acqp, outputnode, [('out_file', 'acq_params')]),
         (b0mask_node, outputnode, [('mask_file', 'out_mask')]),
         (avg_b0_0, outputnode, [('out_file', 'out_b0_pre')]),
         (bet_dwi0, outputnode, [('mask_file', 'out_b0_mask_pre')]),
         (eddy_quad, outputnode, [('qc_json', 'out_eddy_quad_json'),
                                  ('qc_pdf', 'out_eddy_quad_pdf')]),
-        (tensor_wf, outputnode, [('outputnode.FA_file', 'out_dtifit_FA'),
-                                 ('outputnode.MD_file', 'out_dtifit_MD'),
-                                 ('outputnode.AD_file', 'out_dtifit_AD'),
-                                 ('outputnode.RD_file', 'out_dtifit_RD'),
-                                 ('outputnode.V1_file', 'out_dtifit_V1'),
-                                 ('outputnode.sse_file', 'out_dtifit_sse')])
+        (dwi_tensor_wf, outputnode, [('outputnode.FA_file', 'out_dtifit_FA'),
+                                     ('outputnode.MD_file', 'out_dtifit_MD'),
+                                     ('outputnode.AD_file', 'out_dtifit_AD'),
+                                     ('outputnode.RD_file', 'out_dtifit_RD'),
+                                     ('outputnode.V1_file', 'out_dtifit_V1'),
+                                     ('outputnode.sse_file', 'out_dtifit_sse')])
     ])
 
     return dwi_wf

dmriprep/workflows/fieldmap/__init__.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-from .base import init_sdc_prep_wf
+from .base import init_sdc_wf
 from .pepolar import init_pepolar_wf
 from .fmap import init_fmap_wf
 from .phasediff import init_phase_wf, init_phdiff_wf