Repair 3dQwarp workflow

.github/workflows/unittests.yml

@@ -107,7 +107,7 @@ jobs:
         pip install datalad datalad-osf
     - name: Install fsl
       run: |
-        conda install fsl-fugue fsl-topup
+        conda install fsl-fugue fsl-topup fsl-bet2
     - uses: actions/checkout@v4
     - name: Install dependencies
       timeout-minutes: 5

Dockerfile

@@ -72,7 +72,8 @@ RUN mkdir -p /opt/afni-latest \
         -name "3dTshift" -or \
         -name "3dUnifize" -or \
         -name "3dAutomask" -or \
-        -name "3dvolreg" \) -delete
+        -name "3dvolreg" -or \
+        -name "3dQwarp" \) -delete
 
 # Convert3d 1.4.0
 FROM downloader as c3d

env.yml

@@ -23,3 +23,4 @@ dependencies:
   # Workflow dependencies: FSL (versions pinned in 6.0.6.2)
   - fsl-fugue=2201.2
   - fsl-topup=2203.1
+  - fsl-bet2=2111.8

sdcflows/utils/epimanip.py

@@ -188,6 +188,9 @@ def get_trt(in_meta, in_file=None):
             raise ValueError(f"'{trt}'")
 
         return trt
+    elif in_file is None:
+        msg = "`in_file` must be defined if TotalReadoutTime does not appear in `in_meta`."
+        raise ValueError(msg)
 
     # npe = N voxels PE direction
     pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0])

sdcflows/utils/tests/test_epimanip.py

@@ -0,0 +1,33 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2021 The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Test EPI manipulation routines."""
+import pytest
+from sdcflows.utils.epimanip import get_trt
+
+
+def test_get_trt_err_wo_trt_and_in_file():
+    """Test that calling get_trt with dict that does not have TotalReadoutTime \
+    and no in_file raises AssertionError.
+    """
+    with pytest.raises(ValueError):
+        get_trt(in_meta={})

sdcflows/workflows/fit/pepolar.py

@@ -32,6 +32,7 @@
 _PEPOLAR_DESC = """\
 A *B<sub>0</sub>*-nonuniformity map (or *fieldmap*) was estimated based on two (or more)
 echo-planar imaging (EPI) references """
+_PEPOLAR_METHOD = "PEB/PEPOLAR (phase-encoding based / PE-POLARity)"
 
 
 def init_topup_wf(
@@ -118,7 +119,7 @@ def init_topup_wf(
         ),
         name="outputnode",
     )
-    outputnode.inputs.method = "PEB/PEPOLAR (phase-encoding based / PE-POLARity)"
+    outputnode.inputs.method = _PEPOLAR_METHOD
 
     # Calculate the total readout time of each run
     readout_time = pe.MapNode(
@@ -238,13 +239,21 @@ def init_topup_wf(
     return workflow
 
 
-def init_3dQwarp_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
+def init_3dQwarp_wf(
+    omp_nthreads=1, debug=False, sloppy=False, name="pepolar_estimate_wf"
+):
     """
     Create the PEPOLAR field estimation workflow based on AFNI's ``3dQwarp``.
 
     This workflow takes in two EPI files that MUST have opposed
     :abbr:`PE (phase-encoding)` direction.
     Therefore, EPIs with orthogonal PE directions are not supported.
+    ``3dQwarp`` is used to generate a displacement field and correct
+    the reference image. The workflow also returns an estimated fieldmap,
+    which is the result of converting the displacement field to a fieldmap
+    and then regularizing it with a bspline field. This means that the unwarped
+    image is in general not what one would get by reconstructing the fieldmap
+    from fmap_coeff and warping the in_data.
 
     Workflow Graph
         .. workflow ::
@@ -267,25 +276,47 @@ def init_3dQwarp_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
     ------
     in_data : :obj:`list` of :obj:`str`
         A list of two EPI files, the first of which will be taken as reference.
+        The reference PhaseEncodingDirection should match the value for
+        in_reference.
+    metadata : :obj:`list` of :obj:`dict`
+        A list with length matching the length of in_data. Each element should be a
+        dict with keys that are strings and values of any type. One key should be
+        PhaseEncodingDirection and the values should be BIDS-valid codings.
 
     Outputs
     -------
     fmap : :obj:`str`
         The path of the estimated fieldmap.
     fmap_ref : :obj:`str`
-        The path of an unwarped conversion of the first element of ``in_data``.
+        The path of an unwarped conversion of files in ``in_data``.
+    fmap_mask : :obj:`str`
+        The path of mask corresponding to the ``fmap_ref`` output.
+    fmap_coeff : :obj:`str` or :obj:`list` of :obj:`str`
+        The path(s) of the B-Spline coefficients supporting the fieldmap.
+    method: :obj:`str`
+        Short description of the estimation method that was run.
+    out_warps: :obj:`str`
+        The displacement field from 3dQwarp, in ANTS format.
 
     """
     from nipype.interfaces import afni
-    from niworkflows.interfaces.header import CopyHeader
+    from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf
     from niworkflows.interfaces.fixes import (
         FixHeaderRegistration as Registration,
-        FixHeaderApplyTransforms as ApplyTransforms,
     )
     from niworkflows.interfaces.freesurfer import StructuralReference
-    from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf
-    from ...utils.misc import front as _front, last as _last
-    from ...interfaces.utils import Flatten, ConvertWarp
+    from niworkflows.interfaces.header import CopyHeader
+
+    from ...interfaces.bspline import (
+        DEFAULT_HF_ZOOMS_MM,
+        DEFAULT_ZOOMS_MM,
+        ApplyCoeffsField,
+        BSplineApprox,
+    )
+    from ...interfaces.epi import GetReadoutTime
+    from ...interfaces.fmap import DisplacementsField2Fieldmap
+    from ...interfaces.utils import ConvertWarp, Flatten
+    from ...utils.misc import front, last
 
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""{_PEPOLAR_DESC} \
@@ -297,12 +328,31 @@ def init_3dQwarp_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
     )
 
     outputnode = pe.Node(
-        niu.IdentityInterface(fields=["fmap", "fmap_ref"]), name="outputnode"
+        niu.IdentityInterface(
+            fields=[
+                "fmap",
+                "fmap_ref",
+                "fmap_mask",
+                "fmap_coeff",
+                "method",
+                "out_warps",
+            ]
+        ),
+        name="outputnode",
+    )
+    outputnode.inputs.method = _PEPOLAR_METHOD
+
+    readout_time = pe.Node(
+        GetReadoutTime(),
+        name="readout_time",
+        run_without_submitting=True,
     )
 
     flatten = pe.Node(Flatten(), name="flatten")
     sort_pe = pe.Node(
-        niu.Function(function=_sorted_pe, output_names=["sorted", "qwarp_args"]),
+        niu.Function(
+            function=_sorted_pe, output_names=["sorted", "qwarp_args"]
+        ),
         name="sort_pe",
         run_without_submitting=True,
     )
@@ -355,14 +405,24 @@ def init_3dQwarp_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
 
     cphdr_warp = pe.Node(CopyHeader(), name="cphdr_warp", mem_gb=0.01)
 
-    unwarp_reference = pe.Node(
-        ApplyTransforms(
-            dimension=3,
-            float=True,
-            interpolation="LanczosWindowedSinc",
-        ),
-        name="unwarp_reference",
+    # Extract the corresponding fieldmap in Hz
+    extract_field = pe.Node(
+        DisplacementsField2Fieldmap(), name="extract_field"
+    )
+
+    # Regularize with B-Splines
+    bs_filter = pe.Node(
+        BSplineApprox(debug=debug, extrapolate=not debug),
+        name="bs_filter",
     )
+    bs_filter.interface._always_run = debug
+    bs_filter.inputs.bs_spacing = (
+        [DEFAULT_HF_ZOOMS_MM] if not sloppy else [DEFAULT_ZOOMS_MM]
+    )
+    if sloppy:
+        bs_filter.inputs.zooms_min = 4.0
+
+    unwarp = pe.Node(ApplyCoeffsField(), name="unwarp")
 
     # fmt: off
     workflow.connect([
@@ -371,20 +431,33 @@ def init_3dQwarp_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
         (flatten, sort_pe, [("out_list", "inlist")]),
         (sort_pe, qwarp, [("qwarp_args", "args")]),
         (sort_pe, merge_pes, [("sorted", "in_files")]),
-        (merge_pes, pe0_wf, [(("out_file", _front), "inputnode.in_file")]),
-        (merge_pes, pe1_wf, [(("out_file", _last), "inputnode.in_file")]),
+        (merge_pes, pe0_wf, [(("out_file", front), "inputnode.in_file")]),
+        (merge_pes, pe1_wf, [(("out_file", last), "inputnode.in_file")]),
         (pe0_wf, align_pes, [("outputnode.skull_stripped_file", "fixed_image")]),
         (pe1_wf, align_pes, [("outputnode.skull_stripped_file", "moving_image")]),
         (pe0_wf, qwarp, [("outputnode.skull_stripped_file", "in_file")]),
         (align_pes, qwarp, [("warped_image", "base_file")]),
-        (inputnode, cphdr_warp, [(("in_data", _front), "hdr_file")]),
+        (inputnode, cphdr_warp, [(("in_data", front), "hdr_file")]),
         (qwarp, cphdr_warp, [("source_warp", "in_file")]),
         (cphdr_warp, to_ants, [("out_file", "in_file")]),
-        (to_ants, unwarp_reference, [("out_file", "transforms")]),
-        (inputnode, unwarp_reference, [("in_reference", "reference_image"),
-                                       ("in_reference", "input_image")]),
-        (unwarp_reference, outputnode, [("output_image", "fmap_ref")]),
-        (to_ants, outputnode, [("out_file", "fmap")]),
+        (pe0_wf, extract_field, [("outputnode.skull_stripped_file", "epi")]),
+        (to_ants, extract_field, [("out_file", "transform")]),
+        (inputnode, readout_time, [(("metadata", front), "metadata")]),
+        (pe0_wf, readout_time, [("outputnode.skull_stripped_file", "in_file")]),
+        (readout_time, extract_field, [("readout_time", "ro_time"),
+                                       ("pe_direction", "pe_dir")]),
+        (pe1_wf, unwarp, [("outputnode.skull_stripped_file", "in_data")]),
+        (pe0_wf, bs_filter, [("outputnode.mask_file", "in_mask")]),
+        (extract_field, bs_filter, [("out_file", "in_data")]),
+        (bs_filter, unwarp, [("out_coeff", "in_coeff")]),
+        (readout_time, unwarp, [("readout_time", "ro_time"),
+                                ("pe_direction", "pe_dir")]),
+        (bs_filter, outputnode, [("out_coeff", "fmap_coeff")]),
+        (qwarp, outputnode, [("warped_source", "fmap_ref")]),
+        (unwarp, outputnode, [("out_field", "fmap")]),
+        (pe0_wf, outputnode, [("outputnode.mask_file", "fmap_mask")]),
+        (to_ants, outputnode, [("out_file", "out_warps")])
+
     ])
     # fmt: on
     return workflow
@@ -432,11 +505,14 @@ def _sorted_pe(inlist):
         elif pe[0] == ref_pe[0]:
             out_opp.append(d)
         else:
-            raise ValueError("Cannot handle orthogonal PE encodings.")
+            msg = "Cannot handle orthogonal PE encodings."
+            raise ValueError(msg)
 
     return (
         [out_ref, out_opp],
-        {"i": "-noYdis -noZdis", "j": "-noXdis -noZdis", "k": "-noXdis -noYdis"}[
-            ref_pe[0]
-        ],
+        {
+            "i": "-noYdis -noZdis",
+            "j": "-noXdis -noZdis",
+            "k": "-noXdis -noYdis",
+        }[ref_pe[0]],
     )

sdcflows/workflows/fit/tests/test_fit.py

@@ -34,8 +34,12 @@
     (
         ("sdcflows.workflows.fit.fieldmap.init_fmap_wf", {"mode": "mapped"}),
         ("sdcflows.workflows.fit.fieldmap.init_fmap_wf", {}),
-        ("sdcflows.workflows.fit.fieldmap.init_phdiff_wf", {"omp_nthreads": 1}),
+        (
+            "sdcflows.workflows.fit.fieldmap.init_phdiff_wf",
+            {"omp_nthreads": 1},
+        ),
         ("sdcflows.workflows.fit.pepolar.init_3dQwarp_wf", {}),
+        ("sdcflows.workflows.fit.pepolar.init_3dQwarp_wf", {"sloppy": True}),
         ("sdcflows.workflows.fit.pepolar.init_topup_wf", {}),
         ("sdcflows.workflows.fit.syn.init_syn_sdc_wf", {"omp_nthreads": 1}),
     ),

sdcflows/workflows/fit/tests/test_pepolar.py

@@ -25,13 +25,16 @@
 import pytest
 from nipype.pipeline import engine as pe
 
-from ..pepolar import init_topup_wf
+from ..pepolar import init_3dQwarp_wf, init_topup_wf
 
 
 @pytest.mark.skipif(os.getenv("TRAVIS") == "true", reason="this is TravisCI")
-@pytest.mark.skipif(os.getenv("GITHUB_ACTIONS") == "true", reason="this is GH Actions")
-@pytest.mark.parametrize("ds", ("ds001771", "HCP101006"))
-def test_topup_wf(tmpdir, bids_layouts, workdir, outdir, ds):
+@pytest.mark.skipif(
+    os.getenv("GITHUB_ACTIONS") == "true", reason="this is GH Actions"
+)
+@pytest.mark.parametrize("ds", ["ds001771", "HCP101006"])
+@pytest.mark.parametrize("workflow", ["topup", "3dQwarp"])
+def test_pepolar_wf(tmpdir, bids_layouts, workdir, outdir, ds, workflow):
     """Test preparation workflow."""
     layout = bids_layouts[ds]
     epi_path = sorted(
@@ -40,17 +43,24 @@ def test_topup_wf(tmpdir, bids_layouts, workdir, outdir, ds):
     )
     in_data = [f.path for f in epi_path]
 
-    wf = pe.Workflow(name=f"topup_{ds}")
-    topup_wf = init_topup_wf(omp_nthreads=2, debug=True, sloppy=True)
+    wf = pe.Workflow(name=f"{workflow}_{ds}")
+    if workflow == "topup":
+        init_pepolar = init_topup_wf
+    elif workflow == "3dQwarp":
+        init_pepolar = init_3dQwarp_wf
+    else:
+        msg = f"Unknown workflow: {workflow}"
+        raise ValueError(msg)
+    pepolar_wf = init_pepolar(omp_nthreads=2, debug=True, sloppy=True)
     metadata = [layout.get_metadata(f.path) for f in epi_path]
 
-    topup_wf.inputs.inputnode.in_data = in_data
-    topup_wf.inputs.inputnode.metadata = metadata
+    pepolar_wf.inputs.inputnode.in_data = in_data
+    pepolar_wf.inputs.inputnode.metadata = metadata
 
     if outdir:
         from ...outputs import init_fmap_derivatives_wf, init_fmap_reports_wf
 
-        outdir = outdir / "unittests" / f"topup_{ds}"
+        outdir = outdir / "unittests" / f"{workflow}_{ds}"
         fmap_derivatives_wf = init_fmap_derivatives_wf(
             output_dir=str(outdir),
             write_coeff=True,
@@ -67,18 +77,18 @@ def test_topup_wf(tmpdir, bids_layouts, workdir, outdir, ds):
 
         # fmt: off
         wf.connect([
-            (topup_wf, fmap_reports_wf, [("outputnode.fmap", "inputnode.fieldmap"),
-                                         ("outputnode.fmap_ref", "inputnode.fmap_ref"),
-                                         ("outputnode.fmap_mask", "inputnode.fmap_mask")]),
-            (topup_wf, fmap_derivatives_wf, [
+            (pepolar_wf, fmap_reports_wf, [("outputnode.fmap", "inputnode.fieldmap"),
+                                           ("outputnode.fmap_ref", "inputnode.fmap_ref"),
+                                           ("outputnode.fmap_mask", "inputnode.fmap_mask")]),
+            (pepolar_wf, fmap_derivatives_wf, [
                 ("outputnode.fmap", "inputnode.fieldmap"),
                 ("outputnode.fmap_ref", "inputnode.fmap_ref"),
                 ("outputnode.fmap_coeff", "inputnode.fmap_coeff"),
             ]),
         ])
         # fmt: on
     else:
-        wf.add_nodes([topup_wf])
+        wf.add_nodes([pepolar_wf])
 
     if workdir:
         wf.base_dir = str(workdir)