move mrtrix3 commands into separate workflow

Author: josephmje

dmriprep/cli.py

@@ -34,7 +34,9 @@ def __init__(
         omp_nthreads,
         eddy_niter,
         synb0_dir,
-        acqp_file
+        acqp_file,
+        tbss,
+        diffusivity_meas,
     ):
 
         self.layout = layout
@@ -45,13 +47,15 @@ def __init__(
         self.work_dir = work_dir
         self.ignore = ignore
         self.b0_thresh = b0_thresh
-        self.output_resolution = resize_scale
+        self.output_resolution = output_resolution
         self.bet_dwi = bet_dwi
         self.bet_mag = bet_mag
         self.omp_nthreads = omp_nthreads
         self.eddy_niter = eddy_niter
         self.synb0_dir = synb0_dir
         self.acqp_file = acqp_file
+        self.tbss = (tbss,)
+        self.diffusivity_meas = diffusivity_meas
 
 
 @click.command()
@@ -61,15 +65,11 @@ def __init__(
     "output_dir", type=click.Path(exists=True, file_okay=False, writable=True)
 )
 @click.argument(
-    "analysis_level",
-    default="participant",
-    type=click.Choice(["participant", "group"]),
+    "analysis_level", default="participant", type=click.Choice(["participant", "group"])
 )
 # optional arguments
 # options for filtering BIDS queries
-@click.option(
-    "--skip_bids_validation", help="Skip BIDS validation", is_flag=True
-)
+@click.option("--skip_bids_validation", help="Skip BIDS validation", is_flag=True)
 @click.option(
     "--participant_label",
     default=None,
@@ -96,7 +96,7 @@ def __init__(
 @click.option(
     "--output_resolution",
     help="The isotropic voxel size in mm the data will be resampled to before eddy.",
-    type=float
+    type=float,
 )
 # specific options for eddy
 @click.option(
@@ -154,17 +154,13 @@ def __init__(
     help="working directory",
     type=click.Path(exists=True, file_okay=False, writable=True),
 )
-@click.option(
-    "--tbss",
-    help="Run TBSS"
-    is_flag=True
-)
+@click.option("--tbss", help="Run TBSS", is_flag=True)
 @click.option(
     "--diffusivity_meas",
     help="Specify which measures to calculate.",
     default=("FA",),
     type=click.Choice(["FA", "MD", "AD", "RD"]),
-    multiple=True
+    multiple=True,
 )
 @click.option(
     "--synb0_dir",
@@ -188,6 +184,8 @@ def main(
     eddy_niter,
     synb0_dir,
     acqp_file,
+    tbss,
+    diffusivity_meas,
 ):
     """
     BIDS_DIR: The directory with the input dataset formatted according to the
@@ -237,10 +235,12 @@ def main(
         eddy_niter=eddy_niter,
         synb0_dir=synb0_dir,
         acqp_file=acqp_file,
+        tbss=tbss,
+        diffusivity_meas=list(diffusivity_meas),
     )
 
     wf = init_dmriprep_wf(parameters)
-    wf.write_graph(graph2use="colored")
+    wf.write_graph()
     wf.config["execution"]["remove_unnecessary_outputs"] = False
     wf.config["execution"]["keep_inputs"] = True
     wf.config["execution"]["crashfile_format"] = "txt"

dmriprep/utils/hemisphere.py

@@ -5,26 +5,26 @@
 
 
 def is_hemispherical(vecs):
-    """Test whether all points on a unit sphere lie in the same hemisphere.
+    """
+    Test whether all points on a unit sphere lie in the same hemisphere.
+
+    **Inputs**
 
-    Parameters
-    ----------
     vecs : numpy.ndarray
         2D numpy array with shape (N, 3) where N is the number of points.
         All points must lie on the unit sphere.
 
-    Returns
-    -------
+    **Outputs**
+
     is_hemi : bool
         If True, one can find a hemisphere that contains all the points.
         If False, then the points do not lie in any hemisphere
-
     pole : numpy.ndarray
         If `is_hemi == True`, then pole is the "central" pole of the
         input vectors. Otherwise, pole is the zero vector.
 
-    References
-    ----------
+    **References**
+
     https://rstudio-pubs-static.s3.amazonaws.com/27121_a22e51b47c544980bad594d5e0bb2d04.html  # noqa
     """
     if vecs.shape[1] != 3:

dmriprep/workflows/dwi/__init__.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
 from .base import init_dwi_preproc_wf
-from .remove_artefacts import init_remove_artefacts_wf
+from .prep_dwi import init_prep_dwi_wf
 from .tensor import init_tensor_wf
 from .outputs import init_output_wf

dmriprep/workflows/dwi/base.py

@@ -11,7 +11,7 @@
 from nipype.interfaces import ants, fsl, mrtrix3, utility as niu
 from numba import cuda
 
-from .remove_artefacts import init_remove_artefacts_wf
+from .prep_dwi import init_prep_dwi_wf
 from .tensor import init_tensor_wf
 
 # from ..fieldmap.base import init_sdc_prep_wf
@@ -91,13 +91,7 @@ def init_dwi_preproc_wf(subject_id, dwi_file, metadata, parameters):
         name="outputnode",
     )
 
-    denoise = pe.Node(mrtrix3.DWIDenoise(), name="denoise")
-
-    unring = pe.Node(mrtrix3.MRDeGibbs(), name="unring")
-
-    resample = pe.Node(
-        mrtrix3.MRResize(voxel_size=parameters.output_resolution), name="resample"
-    )
+    dwi_prep_wf = init_prep_dwi_wf(parameters.ignore, parameters.output_resolution)
 
     def gen_index(in_file):
         import os
@@ -242,51 +236,9 @@ def b0_average(in_dwi, in_bval, b0_thresh, out_file=None):
             ]
         )
 
-    if "denoise" in parameters.ignore:
-        unring_wf = init_unring_wf()
-
-        dwi_wf.connect(
-            [
-                (inputnode, unring_wf, [("dwi_file", "inputnode.dwi_file")]),
-                (unring_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
-                (unring_wf, ecc, [("outputnode.out_file", "in_file")]),
-            ]
-        )
-
-    elif "unring" in parameters.ignore:
-        denoise_wf = init_denoise_wf()
-
-        dwi_wf.connect(
-            [
-                (inputnode, denoise_wf, [("dwi_file", "inputnode.dwi_file")]),
-                (denoise_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
-                (denoise_wf, ecc, [("outputnode.out_file", "in_file")]),
-            ]
-        )
-
-    else:
-        denoise_wf = init_denoise_wf()
-        unring_wf = init_unring_wf()
-
-        dwi_wf.connect(
-            [
-                (inputnode, denoise_wf, [("dwi_file", "inputnode.dwi_file")]),
-                (
-                    denoise_wf,
-                    unring_wf,
-                    [("outputnode.out_file", "inputnode.dwi_file")],
-                ),
-                (unring_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
-                (unring_wf, ecc, [("outputnode.out_file", "in_file")]),
-            ]
-        )
-
     fslroi = pe.Node(fsl.ExtractROI(t_min=0, t_size=1), name="fslroi")
 
-    bias_correct = pe.Node(
-        ants.N4BiasFieldCorrection(save_bias=True, copy_header=True, dimension=3),
-        name="bias_correct",
-    )
+    bias_correct = pe.Node(mrtrix3.DWIBiasCorrect(use_ants=True), name="bias_correct")
 
     def get_b0_mask_fn(b0_file):
         import os
@@ -312,16 +264,24 @@ def get_b0_mask_fn(b0_file):
 
     dwi_wf.connect(
         [
+            (inputnode, dwi_prep_wf, [("dwi_file", "inputnode.dwi_file")]),
+            (dwi_prep_wf, avg_b0_0, [("outputnode.out_file", "in_dwi")]),
             (inputnode, avg_b0_0, [("bval_file", "in_bval")]),
             (avg_b0_0, bet_dwi0, [("out_file", "in_file")]),
+            (dwi_prep_wf, ecc, [("outputnode.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")]),
             (acqp, ecc, [("out_file", "in_acqp")]),
             (ecc, denoise_eddy, [("out_corrected", "in_file")]),
-            (ecc, fslroi, [("out_corrected", "in_file")]),
-            (fslroi, bias_correct, [("roi_file", "input_image")]),
-            (bias_correct, b0mask_node, [("output_image", "b0_file")]),
+            (
+                ecc,
+                bias_correct,
+                [("out_corrected", "in_file"), ("out_rotated_bvecs", "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,

dmriprep/workflows/dwi/denoise.py

@@ -10,49 +10,50 @@
 from nipype.interfaces import mrtrix3, utility as niu
 
 
-def init_remove_artefacts_wf(parameters):
+def init_prep_dwi_wf(ignore, output_resolution):
 
-    dwi_prep_wf = pe.Workflow(name="remove_artefacts_wf")
+    prep_dwi_wf = pe.Workflow(name="prep_dwi_wf")
+
+    inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_file"]), name="inputnode")
+
+    outputnode = pe.Node(niu.IdentityInterface(fields=["out_file"]), name="outputnode")
 
     denoise = pe.Node(mrtrix3.DWIDenoise(), name="denoise")
 
     unring = pe.Node(mrtrix3.MRDeGibbs(), name="unring")
 
-    resample = pe.Node(
-        mrtrix3.MRResize(voxel_size=parameters.output_resolution), name="resample"
-    )
+    resample = pe.Node(mrtrix3.MRResize(voxel_size=[1, 1, 1]), name="resample")
 
-    inputnode = pe.Node(niu.IdentityInterface(fields=["dwi_file"]), name="inputnode")
-
-    outputnode = pe.Node(niu.IdentityInterface(fields=["out_file"]), name="outputnode")
+    prep_full = ["denoise", "unring", "resample"]
 
-    prep_full = ["denoise", "unring"]
+    if not output_resolution:
+        ignore.append("resample")
 
-    prep_wanted_str = [node for node in prep_full if not (node in ignore_nodes)]
+    prep_wanted_str = [node for node in prep_full if not (node in ignore)]
 
     # Translate string input to node names
-    str2node = {"denoise": denoise, "unring": unring}
+    str2node = {"denoise": denoise, "unring": unring, "resample": resample}
 
     prep_wanted = [str2node[str_node] for str_node in prep_wanted_str]
 
     # If no steps selected, just connect input to output
     if not (prep_wanted):
-        dwi_prep_wf.connect([(inputnode, outputnode, [("dwi_file", "out_file")])])
+        prep_dwi_wf.connect([(inputnode, outputnode, [("dwi_file", "out_file")])])
 
     # If there are steps
     else:
         # Must at least connect input node to first node
         first_node = prep_wanted[0]
-        dwi_prep_wf.connect([(inputnode, first_node, [("dwi_file", "in_file")])])
+        prep_dwi_wf.connect([(inputnode, first_node, [("dwi_file", "in_file")])])
         # Loop through the prep order
         # Note: only works if each node has in_file and out_file
         # Can work around this by wrapping in node/workflow with in_file+out_file
         prev = first_node
         for curr_node in prep_wanted[1:]:
-            dwi_prep_wf.connect([(prev, curr_node, [("out_file", "in_file")])])
+            prep_dwi_wf.connect([(prev, curr_node, [("out_file", "in_file")])])
             prev = curr_node
         # Connect last node to output node
         last_node = prep_wanted[-1]
-        dwi_prep_wf.connect([(last_node, outputnode, [("out_file", "out_file")])])
+        prep_dwi_wf.connect([(last_node, outputnode, [("out_file", "out_file")])])
 
-    return dwi_prep_wf
+    return prep_dwi_wf