diff --git a/.github/workflows/unittests.yml b/.github/workflows/unittests.yml index 06ec4b298f..5a5e46ccf0 100644 --- a/.github/workflows/unittests.yml +++ b/.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 diff --git a/Dockerfile b/Dockerfile index fd226a8f75..d9b7ed3062 100644 --- a/Dockerfile +++ b/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 diff --git a/env.yml b/env.yml index 28b9e9be11..32571c396b 100644 --- a/env.yml +++ b/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 diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index 6af5be7f21..d5f9844795 100644 --- a/sdcflows/utils/epimanip.py +++ b/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]) diff --git a/sdcflows/utils/tests/test_epimanip.py b/sdcflows/utils/tests/test_epimanip.py new file mode 100644 index 0000000000..1ec1b949bb --- /dev/null +++ b/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 +# +# 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={}) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 9bda683a41..928a2b5e5b 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -32,6 +32,7 @@ _PEPOLAR_DESC = """\ A *B0*-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]], ) diff --git a/sdcflows/workflows/fit/tests/test_fit.py b/sdcflows/workflows/fit/tests/test_fit.py index 6b7b862594..32af447fef 100644 --- a/sdcflows/workflows/fit/tests/test_fit.py +++ b/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}), ), diff --git a/sdcflows/workflows/fit/tests/test_pepolar.py b/sdcflows/workflows/fit/tests/test_pepolar.py index 995c867f96..3091fa6aa8 100644 --- a/sdcflows/workflows/fit/tests/test_pepolar.py +++ b/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,10 +77,10 @@ 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"), @@ -78,7 +88,7 @@ def test_topup_wf(tmpdir, bids_layouts, workdir, outdir, ds): ]) # fmt: on else: - wf.add_nodes([topup_wf]) + wf.add_nodes([pepolar_wf]) if workdir: wf.base_dir = str(workdir)