From 19677a1c44284ec49242285ef5a70471dbac651c Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sat, 6 Jul 2024 21:18:32 -0400 Subject: [PATCH 01/15] Repair 3dQwarp workflow --- sdcflows/workflows/fit/pepolar.py | 133 ++++++++++++++----- sdcflows/workflows/fit/tests/test_pepolar.py | 28 ++-- 2 files changed, 118 insertions(+), 43 deletions(-) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 9bda683a41..79b81d334f 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -21,9 +21,8 @@ # https://www.nipreps.org/community/licensing/ # """Datasets with multiple phase encoded directions.""" -from nipype.pipeline import engine as pe from nipype.interfaces import utility as niu - +from nipype.pipeline import engine as pe from niworkflows.engine.workflows import LiterateWorkflow as Workflow from ... import data @@ -32,7 +31,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( grid_reference=0, @@ -88,13 +87,13 @@ def init_topup_wf( """ from nipype.interfaces.fsl.epi import TOPUP - from niworkflows.interfaces.nibabel import MergeSeries, ReorientImage from niworkflows.interfaces.images import RobustAverage + from niworkflows.interfaces.nibabel import MergeSeries, ReorientImage - from ...utils.misc import front as _front - from ...interfaces.epi import GetReadoutTime, SortPEBlips - from ...interfaces.utils import UniformGrid, PadSlices, ReorientImageAndMetadata from ...interfaces.bspline import TOPUPCoeffReorient + from ...interfaces.epi import GetReadoutTime, SortPEBlips + from ...interfaces.utils import PadSlices, ReorientImageAndMetadata, UniformGrid + from ...utils.misc import front as _front from ..ancillary import init_brainextraction_wf workflow = Workflow(name=name) @@ -118,7 +117,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 +237,23 @@ 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,46 @@ 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.func.util import init_enhance_and_skullstrip_bold_wf + from niworkflows.interfaces.fixes import FixHeaderRegistration as Registration + from niworkflows.interfaces.freesurfer import StructuralReference from niworkflows.interfaces.header import CopyHeader - from niworkflows.interfaces.fixes import ( - FixHeaderRegistration as Registration, - FixHeaderApplyTransforms as ApplyTransforms, + + from ...interfaces.bspline import ( + DEFAULT_HF_ZOOMS_MM, + DEFAULT_ZOOMS_MM, + ApplyCoeffsField, + BSplineApprox, ) - 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 ...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,7 +327,22 @@ 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") @@ -355,14 +400,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 +426,32 @@ 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")]), + (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 diff --git a/sdcflows/workflows/fit/tests/test_pepolar.py b/sdcflows/workflows/fit/tests/test_pepolar.py index 995c867f96..9cfbab1f1a 100644 --- a/sdcflows/workflows/fit/tests/test_pepolar.py +++ b/sdcflows/workflows/fit/tests/test_pepolar.py @@ -25,13 +25,14 @@ 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.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 +41,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 +75,10 @@ def test_topup_wf(tmpdir, bids_layouts, workdir, outdir, ds): # fmt: off wf.connect([ - (topup_wf, fmap_reports_wf, [("outputnode.fmap", "inputnode.fieldmap"), + (pepolar_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_derivatives_wf, [ ("outputnode.fmap", "inputnode.fieldmap"), ("outputnode.fmap_ref", "inputnode.fmap_ref"), ("outputnode.fmap_coeff", "inputnode.fmap_coeff"), @@ -78,7 +86,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) From dc296de65d11d05bc6a0a0bc0e907e13ac83db9e Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sat, 6 Jul 2024 21:47:12 -0400 Subject: [PATCH 02/15] Lint --- sdcflows/workflows/fit/pepolar.py | 81 ++++++++++++-------- sdcflows/workflows/fit/tests/test_pepolar.py | 8 +- 2 files changed, 55 insertions(+), 34 deletions(-) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 79b81d334f..46d379f796 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -33,6 +33,7 @@ echo-planar imaging (EPI) references """ _PEPOLAR_METHOD = "PEB/PEPOLAR (phase-encoding based / PE-POLARity)" + def init_topup_wf( grid_reference=0, omp_nthreads=1, @@ -92,7 +93,11 @@ def init_topup_wf( from ...interfaces.bspline import TOPUPCoeffReorient from ...interfaces.epi import GetReadoutTime, SortPEBlips - from ...interfaces.utils import PadSlices, ReorientImageAndMetadata, UniformGrid + from ...interfaces.utils import ( + PadSlices, + ReorientImageAndMetadata, + UniformGrid, + ) from ...utils.misc import front as _front from ..ancillary import init_brainextraction_wf @@ -101,7 +106,9 @@ def init_topup_wf( {_PEPOLAR_DESC} with `topup` (@topup; FSL {TOPUP().version}). """ - inputnode = pe.Node(niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode") + inputnode = pe.Node( + niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode" + ) outputnode = pe.Node( niu.IdentityInterface( fields=[ @@ -139,19 +146,23 @@ def init_topup_wf( SortPEBlips(), name="sort_pe_blips", run_without_submitting=True ) # Merge into one 4D file - concat_blips = pe.Node(MergeSeries(affine_tolerance=1e-4), name="concat_blips") + concat_blips = pe.Node( + MergeSeries(affine_tolerance=1e-4), name="concat_blips" + ) # Pad dimensions so that they meet TOPUP's expectations pad_blip_slices = pe.Node(PadSlices(), name="pad_blip_slices") # Run 3dVolReg between runs: uses RobustAverage for consistency and to generate # debugging artifacts (typically, one wants to look at the average across uncorrected runs) - setwise_avg = pe.Node(RobustAverage(num_threads=omp_nthreads), name="setwise_avg") + setwise_avg = pe.Node( + RobustAverage(num_threads=omp_nthreads), name="setwise_avg" + ) # The core of the implementation # Feed the input images in LAS orientation, so FSL does not run funky reorientations - to_las = pe.Node(ReorientImageAndMetadata(target_orientation="LAS"), name="to_las") + to_las = pe.Node( + ReorientImageAndMetadata(target_orientation="LAS"), name="to_las" + ) topup = pe.Node( - TOPUP( - config=str(data.load(f"flirtsch/b02b0{'_quick' * sloppy}.cnf")) - ), + TOPUP(config=str(data.load(f"flirtsch/b02b0{'_quick' * sloppy}.cnf"))), name="topup", ) # "Generalize" topup coefficients and store them in a spatially-correct NIfTI file @@ -160,7 +171,9 @@ def init_topup_wf( ) # Average the output - ref_average = pe.Node(RobustAverage(num_threads=omp_nthreads), name="ref_average") + ref_average = pe.Node( + RobustAverage(num_threads=omp_nthreads), name="ref_average" + ) # Sophisticated brain extraction of fMRIPrep brainextraction_wf = init_brainextraction_wf() @@ -238,10 +251,8 @@ def init_topup_wf( def init_3dQwarp_wf( - omp_nthreads=1, - debug=False, - sloppy=False, - name="pepolar_estimate_wf"): + omp_nthreads=1, debug=False, sloppy=False, name="pepolar_estimate_wf" +): """ Create the PEPOLAR field estimation workflow based on AFNI's ``3dQwarp``. @@ -251,8 +262,8 @@ def init_3dQwarp_wf( ``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 + 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 @@ -276,11 +287,11 @@ def init_3dQwarp_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 + 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 + 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 @@ -301,7 +312,9 @@ def init_3dQwarp_wf( """ from nipype.interfaces import afni from niworkflows.func.util import init_enhance_and_skullstrip_bold_wf - from niworkflows.interfaces.fixes import FixHeaderRegistration as Registration + from niworkflows.interfaces.fixes import ( + FixHeaderRegistration as Registration, + ) from niworkflows.interfaces.freesurfer import StructuralReference from niworkflows.interfaces.header import CopyHeader @@ -316,7 +329,6 @@ def init_3dQwarp_wf( from ...interfaces.utils import ConvertWarp, Flatten from ...utils.misc import front, last - workflow = Workflow(name=name) workflow.__desc__ = f"""{_PEPOLAR_DESC} \ with `3dQwarp` (@afni; AFNI {''.join(['%02d' % v for v in afni.Info().version() or []])}). @@ -329,13 +341,15 @@ def init_3dQwarp_wf( outputnode = pe.Node( niu.IdentityInterface( fields=[ - "fmap", - "fmap_ref", - "fmap_mask", - "fmap_coeff", + "fmap", + "fmap_ref", + "fmap_mask", + "fmap_coeff", "method", - "out_warps"]), - name="outputnode" + "out_warps", + ] + ), + name="outputnode", ) outputnode.inputs.method = _PEPOLAR_METHOD @@ -347,7 +361,9 @@ def init_3dQwarp_wf( 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, ) @@ -499,11 +515,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_pepolar.py b/sdcflows/workflows/fit/tests/test_pepolar.py index 9cfbab1f1a..00792a62c7 100644 --- a/sdcflows/workflows/fit/tests/test_pepolar.py +++ b/sdcflows/workflows/fit/tests/test_pepolar.py @@ -29,9 +29,11 @@ @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")) -@pytest.mark.parametrize("workflow", ("topup", "3dQwarp")) +@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] From d7f58a5687272ef993384e31c3d5c8479b1f52f6 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sat, 6 Jul 2024 21:50:26 -0400 Subject: [PATCH 03/15] More linting --- sdcflows/workflows/fit/tests/test_pepolar.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sdcflows/workflows/fit/tests/test_pepolar.py b/sdcflows/workflows/fit/tests/test_pepolar.py index 00792a62c7..3091fa6aa8 100644 --- a/sdcflows/workflows/fit/tests/test_pepolar.py +++ b/sdcflows/workflows/fit/tests/test_pepolar.py @@ -78,8 +78,8 @@ def test_pepolar_wf(tmpdir, bids_layouts, workdir, outdir, ds, workflow): # fmt: off wf.connect([ (pepolar_wf, fmap_reports_wf, [("outputnode.fmap", "inputnode.fieldmap"), - ("outputnode.fmap_ref", "inputnode.fmap_ref"), - ("outputnode.fmap_mask", "inputnode.fmap_mask")]), + ("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"), From f8f601c39b2ef1b282b593308fa51adb66e8e87a Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sun, 7 Jul 2024 07:45:44 -0400 Subject: [PATCH 04/15] Add missing fsl dependency for testing --- .github/workflows/unittests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From afc1a7de22576c41db276e685fca7fdacec6e44d Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sun, 7 Jul 2024 08:09:12 -0400 Subject: [PATCH 05/15] Add fsl-beet2 to env.yml for docker --- env.yml | 1 + 1 file changed, 1 insertion(+) 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 From 48ad74cb5ba88a74a7de3f81872e1347943a3a88 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Sun, 7 Jul 2024 16:24:15 -0400 Subject: [PATCH 06/15] Keep 3dQwarp in docker image --- Dockerfile | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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 From 1a118591bc81cdf520322b277061913a0bb6594d Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Mon, 8 Jul 2024 08:15:27 -0400 Subject: [PATCH 07/15] Allow getting trt when TotalReadoutTime not in metadata --- sdcflows/workflows/fit/pepolar.py | 1 + 1 file changed, 1 insertion(+) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 46d379f796..c1c5b136b3 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -454,6 +454,7 @@ def init_3dQwarp_wf( (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")]), From 7b72302e8ed8a545ac3a58f64a5d456e7113cf73 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Mon, 8 Jul 2024 08:16:02 -0400 Subject: [PATCH 08/15] Add informative error when TotalReadoutTime not in metadata and in_file not provided --- sdcflows/utils/epimanip.py | 14 +++++++++--- sdcflows/utils/tests/test_epimanip.py | 33 +++++++++++++++++++++++++++ 2 files changed, 44 insertions(+), 3 deletions(-) create mode 100644 sdcflows/utils/tests/test_epimanip.py diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index 6af5be7f21..836f327b9e 100644 --- a/sdcflows/utils/epimanip.py +++ b/sdcflows/utils/epimanip.py @@ -188,10 +188,14 @@ def get_trt(in_meta, in_file=None): raise ValueError(f"'{trt}'") return trt + elif in_file is None: + msg = "Unable to find TotalReadoutTime in metadata and in_file \ + not defined." + raise AssertionError(msg) # npe = N voxels PE direction pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0]) - npe = nb.load(in_file).shape[pe_index] + npe = nb.loadsave.load(in_file).shape[pe_index] # Use case 2: EES is defined ees = in_meta.get("EffectiveEchoSpacing") @@ -252,7 +256,9 @@ def epi_mask(in_file, out_file=None): maxnorm = np.percentile(closed[closed > 0], 90) closed = np.clip(closed, a_min=0.0, a_max=maxnorm) # Calculate index of center of masses - cm = tuple(np.round(ndimage.measurements.center_of_mass(closed)).astype(int)) + cm = tuple( + np.round(ndimage.measurements.center_of_mass(closed)).astype(int) + ) # Erode the picture of the brain by a lot eroded = ndimage.grey_erosion(closed, structure=ball(5)) # Calculate the residual @@ -268,6 +274,8 @@ def epi_mask(in_file, out_file=None): hdr = img.header.copy() hdr.set_data_dtype("uint8") nb.Nifti1Image( - ndimage.binary_dilation(labels == 2, ball(2)).astype("uint8"), img.affine, hdr + ndimage.binary_dilation(labels == 2, ball(2)).astype("uint8"), + img.affine, + hdr, ).to_filename(out_file) return out_file diff --git a/sdcflows/utils/tests/test_epimanip.py b/sdcflows/utils/tests/test_epimanip.py new file mode 100644 index 0000000000..d6d08a7eea --- /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(AssertionError): + get_trt(in_meta={}) From dd5f251ccacb9c3e8cc290ef95ac98ddc5902f35 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Mon, 8 Jul 2024 10:34:55 -0400 Subject: [PATCH 09/15] Improve test coverage --- sdcflows/workflows/fit/tests/test_fit.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) 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}), ), From 9b40a9d8cb66c4648cfa3f42bf7c3937cf3ba199 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 18:34:03 -0400 Subject: [PATCH 10/15] Improve exception Co-authored-by: Chris Markiewicz --- sdcflows/utils/epimanip.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index 836f327b9e..4461f10a5d 100644 --- a/sdcflows/utils/epimanip.py +++ b/sdcflows/utils/epimanip.py @@ -189,9 +189,8 @@ def get_trt(in_meta, in_file=None): return trt elif in_file is None: - msg = "Unable to find TotalReadoutTime in metadata and in_file \ - not defined." - raise AssertionError(msg) + 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]) From 923ab4a091e1f91c177c5f589f69bf080ed0ef40 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 18:39:14 -0400 Subject: [PATCH 11/15] manually undo accidental linting --- sdcflows/utils/epimanip.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index 4461f10a5d..d5f9844795 100644 --- a/sdcflows/utils/epimanip.py +++ b/sdcflows/utils/epimanip.py @@ -194,7 +194,7 @@ def get_trt(in_meta, in_file=None): # npe = N voxels PE direction pe_index = "ijk".index(in_meta["PhaseEncodingDirection"][0]) - npe = nb.loadsave.load(in_file).shape[pe_index] + npe = nb.load(in_file).shape[pe_index] # Use case 2: EES is defined ees = in_meta.get("EffectiveEchoSpacing") @@ -255,9 +255,7 @@ def epi_mask(in_file, out_file=None): maxnorm = np.percentile(closed[closed > 0], 90) closed = np.clip(closed, a_min=0.0, a_max=maxnorm) # Calculate index of center of masses - cm = tuple( - np.round(ndimage.measurements.center_of_mass(closed)).astype(int) - ) + cm = tuple(np.round(ndimage.measurements.center_of_mass(closed)).astype(int)) # Erode the picture of the brain by a lot eroded = ndimage.grey_erosion(closed, structure=ball(5)) # Calculate the residual @@ -273,8 +271,6 @@ def epi_mask(in_file, out_file=None): hdr = img.header.copy() hdr.set_data_dtype("uint8") nb.Nifti1Image( - ndimage.binary_dilation(labels == 2, ball(2)).astype("uint8"), - img.affine, - hdr, + ndimage.binary_dilation(labels == 2, ball(2)).astype("uint8"), img.affine, hdr ).to_filename(out_file) return out_file From d66c86efd9bd462e873f6368832ce3e89d67f402 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 18:39:39 -0400 Subject: [PATCH 12/15] Update sdcflows/utils/tests/test_epimanip.py Co-authored-by: Chris Markiewicz --- sdcflows/utils/tests/test_epimanip.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdcflows/utils/tests/test_epimanip.py b/sdcflows/utils/tests/test_epimanip.py index d6d08a7eea..1ec1b949bb 100644 --- a/sdcflows/utils/tests/test_epimanip.py +++ b/sdcflows/utils/tests/test_epimanip.py @@ -29,5 +29,5 @@ 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(AssertionError): + with pytest.raises(ValueError): get_trt(in_meta={}) From 0275bb27af86509b850a08b43fb9e503c54587ff Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 18:42:46 -0400 Subject: [PATCH 13/15] manually undo accidental linting --- sdcflows/workflows/fit/pepolar.py | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index c1c5b136b3..0d4114d386 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -88,17 +88,13 @@ def init_topup_wf( """ from nipype.interfaces.fsl.epi import TOPUP - from niworkflows.interfaces.images import RobustAverage from niworkflows.interfaces.nibabel import MergeSeries, ReorientImage + from niworkflows.interfaces.images import RobustAverage - from ...interfaces.bspline import TOPUPCoeffReorient - from ...interfaces.epi import GetReadoutTime, SortPEBlips - from ...interfaces.utils import ( - PadSlices, - ReorientImageAndMetadata, - UniformGrid, - ) from ...utils.misc import front as _front + from ...interfaces.epi import GetReadoutTime, SortPEBlips + from ...interfaces.utils import UniformGrid, PadSlices, ReorientImageAndMetadata + from ...interfaces.bspline import TOPUPCoeffReorient from ..ancillary import init_brainextraction_wf workflow = Workflow(name=name) @@ -106,9 +102,7 @@ def init_topup_wf( {_PEPOLAR_DESC} with `topup` (@topup; FSL {TOPUP().version}). """ - inputnode = pe.Node( - niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode" - ) + inputnode = pe.Node(niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode") outputnode = pe.Node( niu.IdentityInterface( fields=[ From 166be74dd229c3e2487b3f6a4dff9b5bddbcfc2b Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 18:43:41 -0400 Subject: [PATCH 14/15] manually undo accidental linting --- sdcflows/workflows/fit/pepolar.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 0d4114d386..ef892e33d8 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -21,8 +21,9 @@ # https://www.nipreps.org/community/licensing/ # """Datasets with multiple phase encoded directions.""" -from nipype.interfaces import utility as niu from nipype.pipeline import engine as pe +from nipype.interfaces import utility as niu + from niworkflows.engine.workflows import LiterateWorkflow as Workflow from ... import data From ac5f87aa704c29fd368f640851d39bc2f1330334 Mon Sep 17 00:00:00 2001 From: Patrick Sadil Date: Wed, 16 Oct 2024 20:11:26 -0400 Subject: [PATCH 15/15] more accidental linting --- sdcflows/workflows/fit/pepolar.py | 20 +++++++------------- 1 file changed, 7 insertions(+), 13 deletions(-) diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index ef892e33d8..928a2b5e5b 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -141,23 +141,19 @@ def init_topup_wf( SortPEBlips(), name="sort_pe_blips", run_without_submitting=True ) # Merge into one 4D file - concat_blips = pe.Node( - MergeSeries(affine_tolerance=1e-4), name="concat_blips" - ) + concat_blips = pe.Node(MergeSeries(affine_tolerance=1e-4), name="concat_blips") # Pad dimensions so that they meet TOPUP's expectations pad_blip_slices = pe.Node(PadSlices(), name="pad_blip_slices") # Run 3dVolReg between runs: uses RobustAverage for consistency and to generate # debugging artifacts (typically, one wants to look at the average across uncorrected runs) - setwise_avg = pe.Node( - RobustAverage(num_threads=omp_nthreads), name="setwise_avg" - ) + setwise_avg = pe.Node(RobustAverage(num_threads=omp_nthreads), name="setwise_avg") # The core of the implementation # Feed the input images in LAS orientation, so FSL does not run funky reorientations - to_las = pe.Node( - ReorientImageAndMetadata(target_orientation="LAS"), name="to_las" - ) + to_las = pe.Node(ReorientImageAndMetadata(target_orientation="LAS"), name="to_las") topup = pe.Node( - TOPUP(config=str(data.load(f"flirtsch/b02b0{'_quick' * sloppy}.cnf"))), + TOPUP( + config=str(data.load(f"flirtsch/b02b0{'_quick' * sloppy}.cnf")) + ), name="topup", ) # "Generalize" topup coefficients and store them in a spatially-correct NIfTI file @@ -166,9 +162,7 @@ def init_topup_wf( ) # Average the output - ref_average = pe.Node( - RobustAverage(num_threads=omp_nthreads), name="ref_average" - ) + ref_average = pe.Node(RobustAverage(num_threads=omp_nthreads), name="ref_average") # Sophisticated brain extraction of fMRIPrep brainextraction_wf = init_brainextraction_wf()