ENH: Leverage T2w if available for BOLD -> anat coregistration (#3208)

This PR adds a new parameter `use_t2w` to `init_bold_reg_wf`, which
enables/disables using a T2w reference as an intermediate for
bbregister.

This PR only changes the behavior of the `bbreg_wf` fork, the FSL
workflow remains using the T1w reference.

cc @madisoth as this was inspired based on some comparison work from him

Author: mgxd

fmriprep/cli/parser.py

@@ -31,7 +31,7 @@ def _build_parser(**kwargs):
 
     ``kwargs`` are passed to ``argparse.ArgumentParser`` (mainly useful for debugging).
     """
-    from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
+    from argparse import Action, ArgumentDefaultsHelpFormatter, ArgumentParser
     from functools import partial
     from pathlib import Path
 
@@ -40,6 +40,27 @@ def _build_parser(**kwargs):
 
     from .version import check_latest, is_flagged
 
+    deprecations = {
+        # parser attribute name: (replacement flag, version slated to be removed in)
+        'use_aroma': (None, '24.0.0'),
+        'aroma_melodic_dim': (None, '24.0.0'),
+        'aroma_err_on_warn': (None, '24.0.0'),
+        'bold2t1w_init': ('--bold2anat-init', '24.2.0'),
+        'bold2t1w_dof': ('--bold2anat-dof', '24.2.0'),
+    }
+
+    class DeprecatedAction(Action):
+        def __call__(self, parser, namespace, values, option_string=None):
+            new_opt, rem_vers = deprecations.get(self.dest, (None, None))
+            msg = (
+                f"{self.option_strings} has been deprecated and will be removed in "
+                f"{rem_vers or 'a later version'}."
+            )
+            if new_opt:
+                msg += f" Please use `{new_opt}` instead."
+            print(msg, file=sys.stderr)
+            delattr(namespace, self.dest)
+
     def _path_exists(path, parser):
         """Ensure a given path exists."""
         if path is None or not Path(path).exists():
@@ -313,19 +334,32 @@ def _slice_time_ref(value, parser):
     )
     g_conf.add_argument(
         "--bold2t1w-init",
-        action="store",
-        default="register",
+        action=DeprecatedAction,
         choices=["register", "header"],
-        help='Either "register" (the default) to initialize volumes at center or "header"'
-        " to use the header information when coregistering BOLD to T1w images.",
+        help="Deprecated - use `--bold2anat-init` instead.",
     )
     g_conf.add_argument(
         "--bold2t1w-dof",
+        action=DeprecatedAction,
+        choices=[6, 9, 12],
+        type=int,
+        help="Deprecated - use `--bold2anat-dof` instead.",
+    )
+    g_conf.add_argument(
+        "--bold2anat-init",
+        choices=["auto", "t1w", "t2w", "header"],
+        default="auto",
+        help="Method of initial BOLD to anatomical coregistration. If `auto`, a T2w image is used "
+        "if available, otherwise the T1w image. `t1w` forces use of the T1w, `t2w` forces use of "
+        "the T2w, and `header` uses the BOLD header information without an initial registration.",
+    )
+    g_conf.add_argument(
+        "--bold2anat-dof",
         action="store",
         default=6,
         choices=[6, 9, 12],
         type=int,
-        help="Degrees of freedom when registering BOLD to T1w images. "
+        help="Degrees of freedom when registering BOLD to anatomical images. "
         "6 degrees (rotation and translation) are used by default.",
     )
     g_conf.add_argument(
@@ -445,23 +479,20 @@ def _slice_time_ref(value, parser):
     g_aroma = parser.add_argument_group("[DEPRECATED] Options for running ICA_AROMA")
     g_aroma.add_argument(
         "--use-aroma",
-        action="store_true",
-        default=False,
+        action=DeprecatedAction,
         help="Deprecated. Will raise an error in 24.0.",
     )
     g_aroma.add_argument(
         "--aroma-melodic-dimensionality",
         dest="aroma_melodic_dim",
-        action="store",
-        default=0,
+        action=DeprecatedAction,
         type=int,
         help="Deprecated. Will raise an error in 24.0.",
     )
     g_aroma.add_argument(
         "--error-on-aroma-warnings",
-        action="store_true",
+        action=DeprecatedAction,
         dest="aroma_err_on_warn",
-        default=False,
         help="Deprecated. Will raise an error in 24.0.",
     )
 

fmriprep/config.py

@@ -525,11 +525,12 @@ class workflow(_Config):
     aroma_melodic_dim = None
     """Number of ICA components to be estimated by MELODIC
     (positive = exact, negative = maximum)."""
-    bold2t1w_dof = None
-    """Degrees of freedom of the BOLD-to-T1w registration steps."""
-    bold2t1w_init = "register"
-    """Whether to use standard coregistration ('register') or to initialize coregistration from the
-    BOLD image-header ('header')."""
+    bold2anat_dof = None
+    """Degrees of freedom of the BOLD-to-anatomical registration steps."""
+    bold2anat_init = "auto"
+    """Method of initial BOLD to anatomical coregistration. If `auto`, a T2w image is used
+    if available, otherwise the T1w image. `t1w` forces use of the T1w, `t2w` forces use of
+    the T2w, and `header` uses the BOLD header information without an initial registration."""
     cifti_output = None
     """Generate HCP Grayordinates, accepts either ``'91k'`` (default) or ``'170k'``."""
     dummy_scans = None

fmriprep/data/tests/config.toml

@@ -32,7 +32,7 @@ write_graph = false
 anat_only = false
 aroma_err_on_warn = false
 aroma_melodic_dim = -200
-bold2t1w_dof = 6
+bold2anat_dof = 6
 fmap_bspline = false
 force_syn = false
 hires = true

fmriprep/interfaces/patches.py

@@ -29,7 +29,10 @@
 from random import randint
 from time import sleep
 
+import nipype.interfaces.freesurfer as fs
+import nipype.interfaces.io as nio
 from nipype.algorithms import confounds as nac
+from nipype.interfaces.base import File, traits
 from numpy.linalg.linalg import LinAlgError
 
 
@@ -77,3 +80,38 @@ def _run_interface(self, runtime):
                 sleep(randint(start + 4, start + 10))
 
         return runtime
+
+
+class _MRICoregInputSpec(fs.registration.MRICoregInputSpec):
+    reference_file = File(
+        argstr="--ref %s",
+        desc="reference (target) file",
+        copyfile=False,
+    )
+    subject_id = traits.Str(
+        argstr="--s %s",
+        position=1,
+        requires=["subjects_dir"],
+        desc="freesurfer subject ID (implies ``reference_mask == "
+        "aparc+aseg.mgz`` unless otherwise specified)",
+    )
+
+
+class MRICoreg(fs.MRICoreg):
+    """
+    Patched that allows setting both a reference file and the subjects dir.
+    """
+
+    input_spec = _MRICoregInputSpec
+
+
+class _FSSourceOutputSpec(nio.FSSourceOutputSpec):
+    T2 = File(desc="Intensity normalized whole-head volume", loc="mri")
+
+
+class FreeSurferSource(nio.FreeSurferSource):
+    """
+    Patch to allow grabbing the T2 volume, if available
+    """
+
+    output_spec = _FSSourceOutputSpec

fmriprep/interfaces/reports.py

@@ -213,11 +213,12 @@ class FunctionalSummaryInputSpec(TraitedSpec):
         6, 9, 12, desc='Registration degrees of freedom', mandatory=True
     )
     registration_init = traits.Enum(
-        'register',
+        't1w',
+        't2w',
         'header',
         mandatory=True,
         desc='Whether to initialize registration with the "header"'
-        ' or by centering the volumes ("register")',
+        ' or by centering the volumes ("t1w" or "t2w")',
     )
     tr = traits.Float(desc='Repetition time', mandatory=True)
     dummy_scans = traits.Either(traits.Int(), None, desc='number of dummy scans specified by user')

fmriprep/workflows/base.py

@@ -637,6 +637,15 @@ def init_single_subject_wf(subject_id: str):
         num_bold=len(bold_runs)
     )
 
+    # Before initializing BOLD workflow, select/verify anatomical target for coregistration
+    if config.workflow.bold2anat_init in ('auto', 't2w'):
+        has_t2w = subject_data['t2w'] or 't2w_preproc' in anatomical_cache
+        if config.workflow.bold2anat_init == 't2w' and not has_t2w:
+            raise OSError(
+                "A T2w image is expected for BOLD-to-anatomical coregistration and was not found"
+            )
+        config.workflow.bold2anat_init = 't2w' if has_t2w else 't1w'
+
     for bold_series in bold_runs:
         bold_file = bold_series[0]
         fieldmap_id = estimator_map.get(bold_file)

fmriprep/workflows/bold/fit.py

@@ -203,6 +203,7 @@ def init_bold_fit_wf(
     from fmriprep.utils.misc import estimate_bold_mem_usage
 
     layout = config.execution.layout
+    bids_filters = config.execution.get().get('bids_filters', {})
 
     # Fitting operates on the shortest echo
     # This could become more complicated in the future
@@ -211,7 +212,7 @@ def init_bold_fit_wf(
     # Collect sbref files, sorted by EchoTime
     sbref_files = get_sbrefs(
         bold_series,
-        entity_overrides=config.execution.get().get('bids_filters', {}).get('sbref', {}),
+        entity_overrides=bids_filters.get('sbref', {}),
         layout=layout,
     )
 
@@ -312,8 +313,8 @@ def init_bold_fit_wf(
         FunctionalSummary(
             distortion_correction="None",  # Can override with connection
             registration=("FSL", "FreeSurfer")[config.workflow.run_reconall],
-            registration_dof=config.workflow.bold2t1w_dof,
-            registration_init=config.workflow.bold2t1w_init,
+            registration_dof=config.workflow.bold2anat_dof,
+            registration_init=config.workflow.bold2anat_init,
             pe_direction=metadata.get("PhaseEncodingDirection"),
             echo_idx=entities.get("echo", []),
             tr=metadata["RepetitionTime"],
@@ -586,8 +587,8 @@ def init_bold_fit_wf(
     if not boldref2anat_xform:
         # calculate BOLD registration to T1w
         bold_reg_wf = init_bold_reg_wf(
-            bold2t1w_dof=config.workflow.bold2t1w_dof,
-            bold2t1w_init=config.workflow.bold2t1w_init,
+            bold2anat_dof=config.workflow.bold2anat_dof,
+            bold2anat_init=config.workflow.bold2anat_init,
             use_bbr=config.workflow.use_bbr,
             freesurfer=config.workflow.run_reconall,
             omp_nthreads=omp_nthreads,