rf: Allow callers to pass a pre-computed epi2anat transform for SyN

Author: effigies

sdcflows/workflows/fit/syn.py

@@ -325,6 +325,7 @@ def init_syn_preprocessing_wf(
     debug=False,
     name="syn_preprocessing_wf",
     omp_nthreads=1,
+    coregister=True,
     auto_bold_nss=False,
     t1w_inversion=None,
     sd_prior=True,
@@ -348,6 +349,9 @@ def init_syn_preprocessing_wf(
         Name for this workflow
     omp_nthreads : :obj:`int`
         Parallelize internal tasks across the number of CPUs given by this option.
+    coregister: :class:`bool`
+        Run BOLD-to-Anat coregistration. If set to ``False``, ``epi2anat_xfm`` must be
+        provided.
     auto_bold_nss : :obj:`bool`
         Set up the reference workflow to automatically execute nonsteady states detection
         of BOLD images.
@@ -417,6 +421,7 @@ def init_syn_preprocessing_wf(
                 "in_anat",
                 "mask_anat",
                 "std2anat_xfm",
+                "epi2anat_xfm",
             ]
         ),
         name="inputnode",
@@ -464,28 +469,44 @@ def init_syn_preprocessing_wf(
         DenoiseImage(copy_header=True), name="ref_anat", n_procs=omp_nthreads
     )
 
-    epi2anat = pe.Node(
-        Registration(from_file=data.load("affine.json")),
-        name="epi2anat",
-        n_procs=omp_nthreads,
-    )
-    epi2anat.inputs.output_warped_image = debug
-    epi2anat.inputs.output_inverse_warped_image = debug
-    if debug:
-        epi2anat.inputs.args = "--write-interval-volumes 5"
-
-    def _remove_first_mask(in_file):
-        if not isinstance(in_file, list):
-            in_file = [in_file]
-
-        in_file.insert(0, "NULL")
-        return in_file
-
     anat_dilmsk = pe.Node(BinaryDilation(), name="anat_dilmsk")
     epi_dilmsk = pe.Node(BinaryDilation(), name="epi_dilmsk")
 
     sampling_ref = pe.Node(GenerateSamplingReference(), name="sampling_ref")
 
+    if coregister:
+        epi2anat = pe.Node(
+            Registration(from_file=data.load("affine.json")),
+            name="epi2anat",
+            n_procs=omp_nthreads,
+        )
+        epi2anat.inputs.output_warped_image = debug
+        epi2anat.inputs.output_inverse_warped_image = debug
+        if debug:
+            epi2anat.inputs.args = "--write-interval-volumes 5"
+
+        def _remove_first_mask(in_file):
+            if not isinstance(in_file, list):
+                in_file = [in_file]
+
+            in_file.insert(0, "NULL")
+            return in_file
+
+        workflow.connect([
+            (ref_anat, epi2anat, [("output_image", "fixed_image")]),
+            (anat_dilmsk, epi2anat, [("out_file", "fixed_image_masks")]),
+            (deob_epi, epi2anat, [("out_file", "moving_image")]),
+            (epi_dilmsk, epi2anat, [
+                (("out_file", _remove_first_mask), "moving_image_masks")]),
+            (epi2anat, anat2epi, [("forward_transforms", "transforms")]),
+            (epi2anat, mask2epi, [("forward_transforms", "transforms")]),
+        ])  # fmt:skip
+    else:
+        workflow.connect([
+            (inputnode, anat2epi, [("epi2anat_xfm", "transforms")]),
+            (inputnode, mask2epi, [("epi2anat_xfm", "transforms")]),
+        ])
+
     if sd_prior:
         # Mapping & preparing prior knowledge
         # Concatenate transform files:
@@ -511,12 +532,20 @@ def _remove_first_mask(in_file):
 
         workflow.connect([
             (inputnode, transform_list, [("std2anat_xfm", "in2")]),
-            (epi2anat, transform_list, [("forward_transforms", "in1")]),
             (transform_list, prior2epi, [("out", "transforms")]),
             (sampling_ref, prior2epi, [("out_file", "reference_image")]),
             (prior2epi, outputnode, [("output_image", "sd_prior")]),
         ])  # fmt:skip
 
+        if coregister:
+            workflow.connect([
+                (epi2anat, transform_list, [("forward_transforms", "in1")]),
+            ])  # fmt:skip
+        else:
+            workflow.connect([
+                (inputnode, transform_list, [("epi2anat_xfm", "in1")]),
+            ])
+
     else:
         # no prior to be used
         # MG: Future goal is to allow using alternative mappings
@@ -536,16 +565,9 @@ def _remove_first_mask(in_file):
         (clip_anat, ref_anat, [("out_file", "input_image")]),
         (deob_epi, epi_brain, [("out_file", "in_file")]),
         (epi_brain, epi_dilmsk, [("out_mask", "in_file")]),
-        (ref_anat, epi2anat, [("output_image", "fixed_image")]),
-        (anat_dilmsk, epi2anat, [("out_file", "fixed_image_masks")]),
-        (deob_epi, epi2anat, [("out_file", "moving_image")]),
-        (epi_dilmsk, epi2anat, [
-            (("out_file", _remove_first_mask), "moving_image_masks")]),
         (deob_epi, sampling_ref, [("out_file", "fixed_image")]),
         (ref_anat, anat2epi, [("output_image", "input_image")]),
-        (epi2anat, anat2epi, [("forward_transforms", "transforms")]),
         (sampling_ref, anat2epi, [("out_file", "reference_image")]),
-        (epi2anat, mask2epi, [("forward_transforms", "transforms")]),
         (sampling_ref, mask2epi, [("out_file", "reference_image")]),
         (mask2epi, mask_dtype, [("output_image", "in_file")]),
         (anat2epi, outputnode, [("output_image", "anat_ref")]),

sdcflows/workflows/fit/tests/test_syn.py

@@ -22,6 +22,8 @@
 #
 """Test fieldmap-less SDC-SyN."""
 import json
+
+import acres
 import pytest
 from nipype.pipeline import engine as pe
 
@@ -30,8 +32,15 @@
 
 @pytest.mark.veryslow
 @pytest.mark.slow
-@pytest.mark.parametrize("sd_prior", [True, False])
-def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
+@pytest.mark.parametrize(
+    ("n_bold", "coregister", "sd_prior"),
+    [
+        (1, True, True),
+        # Switch to False once we have a transform in tests/data
+        (2, True, False),
+    ]
+)
+def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, n_bold, coregister, sd_prior):
     """Build and run an SDC-SyN workflow."""
     derivs_path = datadir / "ds000054" / "derivatives"
     smriprep = derivs_path / "smriprep-0.6" / "sub-100185" / "anat"
@@ -43,6 +52,7 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
         debug=sloppy_mode,
         auto_bold_nss=True,
         sd_prior=sd_prior,
+        coregister=coregister,
     )
     prep_wf.inputs.inputnode.in_epis = [
         str(
@@ -59,10 +69,10 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
             / "func"
             / "sub-100185_task-machinegame_run-02_bold.nii.gz"
         ),
-    ]
+    ][:n_bold]
     prep_wf.inputs.inputnode.in_meta = [
         json.loads((datadir / "ds000054" / "task-machinegame_bold.json").read_text()),
-    ] * 2
+    ] * n_bold
     prep_wf.inputs.inputnode.std2anat_xfm = str(
         smriprep / "sub-100185_from-MNI152NLin2009cAsym_to-T1w_mode-image_xfm.h5"
     )
@@ -72,6 +82,11 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
     prep_wf.inputs.inputnode.mask_anat = str(
         smriprep / "sub-100185_desc-brain_mask.nii.gz"
     )
+    if not coregister:
+        test_data = acres.Loader('sdcflows.tests')
+        prep_wf.inputs.inputnode.epi_ref = str(
+            test_data('data/anat2epi_xfm.txt')
+        )
 
     syn_wf = init_syn_sdc_wf(
         debug=sloppy_mode,