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

Dockerfile

@@ -74,12 +74,6 @@ RUN mkdir -p /opt/afni-latest \
         -name "3dAutomask" -or \
         -name "3dvolreg" \) -delete
 
-# Convert3d 1.4.0
-FROM downloader as c3d
-RUN mkdir /opt/convert3d && \
-    curl -fsSL --retry 5 https://sourceforge.net/projects/c3d/files/c3d/Experimental/c3d-1.4.0-Linux-gcc64.tar.gz/download \
-    | tar -xz -C /opt/convert3d --strip-components 1
-
 # Micromamba
 FROM downloader as micromamba
 WORKDIR /
@@ -100,7 +94,7 @@ RUN /opt/conda/envs/sdcflows/bin/pip install --no-cache-dir -r /tmp/requirements
 #
 # Main stage
 #
-FROM --platform=linux/amd64 ${BASE_IMAGE} as sdcflows
+FROM --platform=linux/amd64 ${BASE_IMAGE} as dev
 
 # Configure apt
 ENV DEBIAN_FRONTEND="noninteractive" \
@@ -155,7 +149,6 @@ RUN apt-get update -qq \
 
 # Install files from stages
 COPY --from=afni /opt/afni-latest /opt/afni-latest
-COPY --from=c3d /opt/convert3d/bin/c3d_affine_tool /usr/bin/c3d_affine_tool
 
 # AFNI config
 ENV PATH="/opt/afni-latest:$PATH" \
@@ -207,6 +200,8 @@ ENV SUBJECTS_DIR="$FREESURFER_HOME/subjects" \
 ENV MKL_NUM_THREADS=1 \
     OMP_NUM_THREADS=1
 
+FROM dev as sdcflows
+
 # Installing SDCFlows
 COPY --from=src /src/dist/*.whl .
 RUN pip install --no-cache-dir $( ls *.whl )[all]

sdcflows/workflows/fit/syn.py

@@ -342,8 +342,9 @@
     debug=False,
     name="syn_preprocessing_wf",
     omp_nthreads=1,
+    coregister=True,
     auto_bold_nss=False,
-    t1w_inversion=False,
+    t1w_inversion=None,
     sd_prior=True,
 ):
     """
@@ -365,11 +366,15 @@
         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.
     t1w_inversion : :obj:`bool`
         Run T1w intensity inversion so that it looks more like a T2 contrast.
+        (DEPRECATED. Does nothing.)
     sd_prior : :obj:`bool`
         Enable using a prior map to regularize the SyN cost function.
 
@@ -418,6 +423,13 @@
     from ...interfaces.utils import Deoblique, DenoiseImage
     from ...interfaces.brainmask import BrainExtraction, BinaryDilation
 
+    if t1w_inversion is not None:
+        import warnings
+        warnings.warn(
+            "The `t1w_inversion` argument is deprecated and does nothing.",
+            DeprecationWarning,
+        )
+
     workflow = Workflow(name=name)
 
     inputnode = pe.Node(
@@ -429,6 +441,7 @@
                 "in_anat",
                 "mask_anat",
                 "std2anat_xfm",
+                "epi2anat_xfm",
             ]
         ),
         name="inputnode",
@@ -476,28 +489,44 @@
         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:
@@ -523,12 +552,20 @@
 
         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
@@ -548,16 +585,9 @@
         (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"
@@ -42,8 +51,8 @@ def test_syn_wf(tmpdir, datadir, workdir, outdir, sloppy_mode, sd_prior):
         omp_nthreads=4,
         debug=sloppy_mode,
         auto_bold_nss=True,
-        t1w_inversion=True,
         sd_prior=sd_prior,
+        coregister=coregister,
     )
     prep_wf.inputs.inputnode.in_epis = [
         str(
@@ -60,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"
     )
@@ -73,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,