ENH: insert gradunwarp in bold workflow

fmriprep/cli/parser.py

@@ -279,6 +279,11 @@ def _slice_time_ref(value, parser):
         help="Ignore selected aspects of the input dataset to disable corresponding "
         "parts of the workflow (a space delimited list)",
     )
+    g_conf.add_argument(
+        "--gradunwarp-file",
+        metavar="PATH",
+        help="Path to vendor file for gradunwarp gradient distortion  correction.",
+    )
     g_conf.add_argument(
         "--output-spaces",
         nargs="*",

fmriprep/config.py

@@ -536,7 +536,9 @@ class workflow(_Config):
     ignore = None
     """Ignore particular steps for *fMRIPrep*."""
     longitudinal = False
-    """Run FreeSurfer ``recon-all`` with the ``-logitudinal`` flag."""
+    """Run FreeSurfer ``recon-all`` with the ``-longitudinal`` flag."""
+    gradunwarp_file = None
+    """Run GradUnwarp using the provided gradient or coefficient file"""
     medial_surface_nan = None
     """Fill medial surface with :abbr:`NaNs (not-a-number)` when sampling."""
     project_goodvoxels = False

fmriprep/workflows/base.py

@@ -319,6 +319,8 @@ def init_single_subject_wf(subject_id: str):
         spaces=spaces,
         t1w=subject_data['t1w'],
         t2w=subject_data['t2w'],
+        flair=subject_data['flair'],
+        gradunwarp_file=config.workflow.gradunwarp_file,
         cifti_output=config.workflow.cifti_output,
     )
     # fmt:off

fmriprep/workflows/bold/base.py

@@ -111,6 +111,8 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space
     fsnative2t1w_xfm
         LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w
+    gradunwarp_file
+        Vendor provided .coeff or .grad file for gradient distortion correction.
 
     Outputs
     -------
@@ -187,6 +189,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
 
     """
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
+    from niworkflows.workflows.gradunwarp import init_gradunwarp_wf
     from niworkflows.func.util import init_bold_reference_wf
     from niworkflows.interfaces.nibabel import ApplyMask
     from niworkflows.interfaces.reportlets.registration import (
@@ -478,6 +481,17 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     # Top-level BOLD splitter
     bold_split = pe.Node(FSLSplit(dimension="t"), name="bold_split", mem_gb=mem_gb["filesize"] * 3)
 
+    # Gradient unwarping
+    if config.workflow.gradunwarp_file:
+        gradunwarp_wf = init_gradunwarp_wf()
+        gradunwarp_wf.inputs.inputnode.grad_file = config.workflow.gradunwarp_file
+        # input corrected file, raw not necessary, only uses acquisition matrix size
+        # unless mask/ref are used for sdcflows
+        workflow.connect([
+            (initial_boldref_wf, gradunwarp_wf, [
+                ('outputnode.raw_ref_image', 'inputnode.input_file')]),
+        ])
+
     # HMC on the BOLD
     bold_hmc_wf = init_bold_hmc_wf(
         name="bold_hmc_wf", mem_gb=mem_gb["filesize"], omp_nthreads=omp_nthreads
@@ -505,6 +519,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
         use_compression=False,
     )
     bold_t1_trans_wf.inputs.inputnode.fieldwarp = "identity"
+    bold_t1_trans_wf.inputs.inputnode.gradient_warp = "identity"
 
     # get confounds
     bold_confounds_wf = init_bold_confs_wf(
@@ -637,6 +652,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
             ("mask", "inputnode.ref_bold_mask"),
             ("boldref", "inputnode.ref_bold_brain"),
         ]),
+        (gradunwarp_wf, bold_t1_trans_wf, [('outputnode.warp_file', 'inputnode.gradient_warp')]),
         (bold_t1_trans_wf, outputnode, [
             ("outputnode.bold_t1", "bold_t1"),
             ("outputnode.bold_t1_ref", "bold_t1_ref"),
@@ -779,6 +795,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
             use_compression=not config.execution.low_mem,
         )
         bold_std_trans_wf.inputs.inputnode.fieldwarp = "identity"
+        bold_std_trans_wf.inputs.inputnode.gradient_warp = "identity"
 
         # fmt:off
         workflow.connect([
@@ -801,6 +818,9 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
                 ("outputnode.bold_std_ref", "bold_std_ref"),
                 ("outputnode.bold_mask_std", "bold_mask_std"),
             ]),
+            (gradunwarp_wf, bold_std_trans_wf, [
+                ('outputnode.warp_file', 'inputnode.gradient_warp')
+            ]),
         ])
         # fmt:on
 
@@ -1010,6 +1030,7 @@ def _last(inlist):
             name="bold_bold_trans_wf",
         )
         bold_bold_trans_wf.inputs.inputnode.fieldwarp = "identity"
+        bold_bold_trans_wf.inputs.inputnode.gradient_warp = "identity"
 
         # fmt:off
         workflow.connect([
@@ -1019,6 +1040,7 @@ def _last(inlist):
             (bold_hmc_wf, bold_bold_trans_wf, [
                 ("outputnode.xforms", "inputnode.hmc_xforms"),
             ]),
+            (gradunwarp_wf, bold_bold_trans_wf, [('outputnode.warp_file', 'inputnode.gradient_warp')]),
         ])
 
         workflow.connect([

fmriprep/workflows/bold/registration.py

@@ -325,6 +325,7 @@ def init_bold_t1_trans_wf(
                 't1w_aparc',
                 'bold_split',
                 'fieldwarp',
+                'gradient_warp',
                 'hmc_xforms',
                 'itk_bold_to_t1',
             ]
@@ -400,7 +401,7 @@ def init_bold_t1_trans_wf(
 
     # Merge transforms placing the head motion correction last
     merge_xforms = pe.Node(
-        niu.Merge(3),
+        niu.Merge(4),
         name='merge_xforms',
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
@@ -409,6 +410,7 @@ def init_bold_t1_trans_wf(
     workflow.connect([
         (inputnode, merge, [('name_source', 'header_source')]),
         (inputnode, merge_xforms, [
+            ('fieldwarp', 'in4'),  # May be 'identity' if no gradunwarp applied
             ('hmc_xforms', 'in3'),  # May be 'identity' if HMC already applied
             ('fieldwarp', 'in2'),   # May be 'identity' if SDC already applied
             ('itk_bold_to_t1', 'in1')]),

fmriprep/workflows/bold/resampling.py

@@ -912,6 +912,7 @@ def init_bold_std_trans_wf(
                 "bold_split",
                 "t2star",
                 "fieldwarp",
+                "gradient_warp",
                 "hmc_xforms",
                 "itk_bold_to_t1",
                 "name_source",
@@ -964,7 +965,7 @@ def init_bold_std_trans_wf(
     )
 
     merge_xforms = pe.Node(
-        niu.Merge(4),
+        niu.Merge(5),
         name="merge_xforms",
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
@@ -994,7 +995,8 @@ def init_bold_std_trans_wf(
                                  ("templates", "keys")]),
         (inputnode, mask_std_tfm, [("bold_mask", "input_image")]),
         (inputnode, gen_ref, [(("bold_split", _first), "moving_image")]),
-        (inputnode, merge_xforms, [("hmc_xforms", "in4"),
+        (inputnode, merge_xforms, [("gradient_warp", "in5"),
+                                   ("hmc_xforms", "in4"),
                                    ("fieldwarp", "in3"),
                                    (("itk_bold_to_t1", _aslist), "in2")]),
         (inputnode, merge, [("name_source", "header_source")]),
@@ -1170,7 +1172,9 @@ def init_bold_preproc_trans_wf(
     )
 
     inputnode = pe.Node(
-        niu.IdentityInterface(fields=["name_source", "bold_file", "hmc_xforms", "fieldwarp"]),
+        niu.IdentityInterface(
+            fields=["name_source", "bold_file", "hmc_xforms", "fieldwarp", "gradient_warp"]
+        ),
         name="inputnode",
     )
 
@@ -1180,7 +1184,7 @@ def init_bold_preproc_trans_wf(
     )
 
     merge_xforms = pe.Node(
-        niu.Merge(2),
+        niu.Merge(3),
         name="merge_xforms",
         run_without_submitting=True,
         mem_gb=DEFAULT_MEMORY_MIN_GB,
@@ -1203,7 +1207,8 @@ def init_bold_preproc_trans_wf(
     # fmt:off
     workflow.connect([
         (inputnode, merge_xforms, [("fieldwarp", "in1"),
-                                   ("hmc_xforms", "in2")]),
+                                   ("hmc_xforms", "in2"),
+                                   ("gradient_warp", "in3")]),
         (inputnode, bold_transform, [("bold_file", "input_image"),
                                      (("bold_file", _first), "reference_image")]),
         (inputnode, merge, [("name_source", "header_source")]),