diff --git a/sdcflows/conftest.py b/sdcflows/conftest.py index e837168818..180e4addbb 100644 --- a/sdcflows/conftest.py +++ b/sdcflows/conftest.py @@ -16,9 +16,9 @@ if p.is_dir() } -data_dir = Path(__file__).parent / "tests" / "data" / "dsA" +data_dir = Path(__file__).parent / "tests" / "data" -layouts["dsA"] = BIDSLayout(data_dir, validate=False, derivatives=False) +layouts["dsA"] = BIDSLayout(data_dir / "dsA", validate=False, derivatives=False) def pytest_report_header(config): @@ -40,7 +40,7 @@ def add_np(doctest_namespace): for key, val in list(layouts.items()): doctest_namespace[key] = Path(val.root) - doctest_namespace["testdata_dir"] = data_dir + doctest_namespace["dsA_dir"] = data_dir / "dsA" @pytest.fixture @@ -68,3 +68,8 @@ def datadir(): @pytest.fixture def testdata_dir(): return data_dir + + +@pytest.fixture +def dsA_dir(): + return data_dir / "dsA" diff --git a/sdcflows/fieldmaps.py b/sdcflows/fieldmaps.py index 67592a7ce9..e62995eed6 100644 --- a/sdcflows/fieldmaps.py +++ b/sdcflows/fieldmaps.py @@ -82,65 +82,65 @@ class FieldmapFile: Examples -------- - >>> f = FieldmapFile(testdata_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz") + >>> f = FieldmapFile(dsA_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz") >>> f.suffix 'T1w' >>> FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", ... find_meta=False ... ) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): MetadataError: >>> f = FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", ... ) >>> f.metadata['TotalReadoutTime'] 0.005 >>> f = FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_dir-LR_epi.nii.gz", ... metadata={'TotalReadoutTime': 0.006} ... ) >>> f.metadata['TotalReadoutTime'] 0.006 >>> FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz", ... find_meta=False ... ) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): MetadataError: >>> f = FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz" + ... dsA_dir / "sub-01" / "fmap" / "sub-01_phasediff.nii.gz" ... ) >>> f.metadata['EchoTime2'] 0.00746 >>> FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz", ... find_meta=False ... ) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): MetadataError: >>> f = FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz" + ... dsA_dir / "sub-01" / "fmap" / "sub-01_phase2.nii.gz" ... ) >>> f.metadata['EchoTime'] 0.00746 >>> FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz", + ... dsA_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz", ... find_meta=False ... ) # doctest: +IGNORE_EXCEPTION_DETAIL Traceback (most recent call last): MetadataError: >>> f = FieldmapFile( - ... testdata_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz" + ... dsA_dir / "sub-01" / "fmap" / "sub-01_fieldmap.nii.gz" ... ) >>> f.metadata['Units'] 'rad/s' diff --git a/sdcflows/interfaces/epi.py b/sdcflows/interfaces/epi.py index eaad8ea1e1..65214048ec 100644 --- a/sdcflows/interfaces/epi.py +++ b/sdcflows/interfaces/epi.py @@ -8,6 +8,7 @@ traits, SimpleInterface, ) +from nipype.utils.filemanip import fname_presuffix class _GetReadoutTimeInputSpec(BaseInterfaceInputSpec): @@ -35,3 +36,27 @@ def _run_interface(self, runtime): ) self._results["pe_direction"] = self.inputs.metadata["PhaseEncodingDirection"] return runtime + + +class _EPIMaskInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, desc="EPI image") + + +class _EPIMaskOutputSpec(TraitedSpec): + out_file = File(exists=True) + + +class EPIMask(SimpleInterface): + """Calculate the readout time from available metadata.""" + + input_spec = _EPIMaskInputSpec + output_spec = _EPIMaskOutputSpec + + def _run_interface(self, runtime): + from ..utils.epimanip import epi_mask + + self._results["out_file"] = epi_mask( + self.inputs.in_file, + fname_presuffix(self.inputs.in_file, suffix="_mask", newpath=runtime.cwd), + ) + return runtime diff --git a/sdcflows/tests/data/epi.nii.gz b/sdcflows/tests/data/epi.nii.gz new file mode 100644 index 0000000000..26956cb439 Binary files /dev/null and b/sdcflows/tests/data/epi.nii.gz differ diff --git a/sdcflows/tests/test_fieldmaps.py b/sdcflows/tests/test_fieldmaps.py index be26fd4e16..820e0cb97c 100644 --- a/sdcflows/tests/test_fieldmaps.py +++ b/sdcflows/tests/test_fieldmaps.py @@ -3,9 +3,9 @@ from .. import fieldmaps as fm -def test_FieldmapFile(testdata_dir): +def test_FieldmapFile(dsA_dir): """Test one existing file.""" - fm.FieldmapFile(testdata_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz") + fm.FieldmapFile(dsA_dir / "sub-01" / "anat" / "sub-01_T1w.nii.gz") @pytest.mark.parametrize( @@ -56,9 +56,9 @@ def test_FieldmapFile(testdata_dir): ), ], ) -def test_FieldmapEstimation(testdata_dir, inputfiles, method, nsources, raises): +def test_FieldmapEstimation(dsA_dir, inputfiles, method, nsources, raises): """Test errors.""" - sub_dir = testdata_dir / "sub-01" + sub_dir = dsA_dir / "sub-01" sources = [sub_dir / f for f in inputfiles] @@ -105,9 +105,9 @@ def test_FieldmapEstimation(testdata_dir, inputfiles, method, nsources, raises): (("anat/sub-01_T1w.nii.gz", "fmap/sub-01_phase2.nii.gz"), TypeError), ], ) -def test_FieldmapEstimationError(testdata_dir, inputfiles, errortype): +def test_FieldmapEstimationError(dsA_dir, inputfiles, errortype): """Test errors.""" - sub_dir = testdata_dir / "sub-01" + sub_dir = dsA_dir / "sub-01" fm.clear_registry() @@ -117,7 +117,7 @@ def test_FieldmapEstimationError(testdata_dir, inputfiles, errortype): fm.clear_registry() -def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): +def test_FieldmapEstimationIdentifier(monkeypatch, dsA_dir): """Check some use cases of B0FieldIdentifier.""" fm.clear_registry() @@ -125,11 +125,11 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): fm.FieldmapEstimation( [ fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_0"}, ), fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_1"}, ), ] @@ -138,7 +138,7 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): fe = fm.FieldmapEstimation( [ fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", metadata={ "Units": "Hz", "B0FieldIdentifier": "fmap_0", @@ -146,7 +146,7 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): }, ), fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_0"}, ), ] @@ -159,11 +159,11 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): fm.FieldmapEstimation( [ fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_0"}, ), fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_0"}, ), ] @@ -174,7 +174,7 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): fe = fm.FieldmapEstimation( [ fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_fieldmap.nii.gz", metadata={ "Units": "Hz", "B0FieldIdentifier": "fmap_1", @@ -182,7 +182,7 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): }, ), fm.FieldmapFile( - testdata_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", + dsA_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz", metadata={"Units": "Hz", "B0FieldIdentifier": "fmap_1"}, ), ] @@ -192,7 +192,7 @@ def test_FieldmapEstimationIdentifier(monkeypatch, testdata_dir): assert fm.get_identifier("file2.nii.gz") == ("fmap_1",) assert not fm.get_identifier("file3.nii.gz") assert fm.get_identifier( - str(testdata_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz"), by="sources" + str(dsA_dir / "sub-01" / "fmap/sub-01_magnitude.nii.gz"), by="sources" ) == ("fmap_1",) with monkeypatch.context() as m: diff --git a/sdcflows/utils/epimanip.py b/sdcflows/utils/epimanip.py index 25d9d983b0..dab5213939 100644 --- a/sdcflows/utils/epimanip.py +++ b/sdcflows/utils/epimanip.py @@ -200,3 +200,48 @@ def get_trt(in_meta, in_file=None): return ees * (npe - 1) raise ValueError("Unknown total-readout time specification") + + +def epi_mask(in_file, out_file=None): + """Use grayscale morphological operations to obtain a quick mask of EPI data.""" + from pathlib import Path + import nibabel as nb + import numpy as np + from scipy import ndimage + from skimage.morphology import ball + + if out_file is None: + out_file = Path("mask.nii.gz").absolute() + + img = nb.load(in_file) + data = img.get_fdata(dtype="float32") + # First open to blur out the skull around the brain + opened = ndimage.grey_opening(data, structure=ball(3)) + # Second, close large vessels and the ventricles + closed = ndimage.grey_closing(opened, structure=ball(2)) + + # Window filter on percentile 30 + closed -= np.percentile(closed, 30) + # Window filter on percentile 90 of data + 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)) + # Erode the picture of the brain by a lot + eroded = ndimage.grey_erosion(closed, structure=ball(5)) + # Calculate the residual + wshed = opened - eroded + wshed -= wshed.min() + wshed = np.round(1e3 * wshed / wshed.max()).astype(np.uint16) + markers = np.zeros_like(wshed, dtype=int) + markers[cm] = 2 + markers[0, 0, -1] = -1 + # Run watershed + labels = ndimage.watershed_ift(wshed, markers) + + hdr = img.header.copy() + hdr.set_data_dtype("uint8") + nb.Nifti1Image( + 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..b11727620c --- /dev/null +++ b/sdcflows/utils/tests/test_epimanip.py @@ -0,0 +1,11 @@ +"""Test EPI manipulation routines.""" +import numpy as np +import nibabel as nb +from ..epimanip import epi_mask + + +def test_epi_mask(tmpdir, testdata_dir): + """Check mask algorithm.""" + tmpdir.chdir() + mask = epi_mask(testdata_dir / "epi.nii.gz") + assert abs(np.asanyarray(nb.load(mask).dataobj).sum() - 189052) < 10 diff --git a/sdcflows/workflows/base.py b/sdcflows/workflows/base.py index 04e26d2057..374698c12b 100644 --- a/sdcflows/workflows/base.py +++ b/sdcflows/workflows/base.py @@ -123,16 +123,6 @@ def init_fmap_preproc_wf( (inputnode, est_wf, [(f, f"inputnode.{f}") for f in fields]) ]) # fmt:on - else: - # PEPOLAR and ANAT do not produce masks - # fmt:off - workflow.connect([ - (est_wf, fmap_reports_wf, [ - ("outputnode.fmap_mask", "inputnode.fmap_mask"), - ]), - (est_wf, out_map, [("outputnode.fmap_mask", "fmap_mask")]), - ]) - # fmt:on # fmt:off workflow.connect([ @@ -144,11 +134,13 @@ def init_fmap_preproc_wf( (est_wf, fmap_reports_wf, [ ("outputnode.fmap", "inputnode.fieldmap"), ("outputnode.fmap_ref", "inputnode.fmap_ref"), + ("outputnode.fmap_mask", "inputnode.fmap_mask"), ]), (est_wf, out_map, [ ("outputnode.fmap", "fmap"), ("outputnode.fmap_ref", "fmap_ref"), ("outputnode.fmap_coeff", "fmap_coeff"), + ("outputnode.fmap_mask", "fmap_mask"), ]), ]) # fmt:on diff --git a/sdcflows/workflows/fit/pepolar.py b/sdcflows/workflows/fit/pepolar.py index 4b728f712b..6ec148c349 100644 --- a/sdcflows/workflows/fit/pepolar.py +++ b/sdcflows/workflows/fit/pepolar.py @@ -57,6 +57,8 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"): The path of the estimated fieldmap. fmap_ref : :obj:`str` 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. @@ -65,7 +67,7 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"): from niworkflows.interfaces.nibabel import MergeSeries from niworkflows.interfaces.images import IntraModalMerge - from ...interfaces.epi import GetReadoutTime + from ...interfaces.epi import GetReadoutTime, EPIMask from ...interfaces.utils import Flatten workflow = Workflow(name=name) @@ -78,7 +80,15 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"): ) outputnode = pe.Node( niu.IdentityInterface( - fields=["fmap", "fmap_ref", "fmap_coeff", "jacobians", "xfms", "out_warps"] + fields=[ + "fmap", + "fmap_ref", + "fmap_coeff", + "fmap_mask", + "jacobians", + "xfms", + "out_warps", + ] ), name="outputnode", ) @@ -102,6 +112,7 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"): IntraModalMerge(hmc=False, to_ras=False), name="merge_corrected" ) + epi_mask = pe.Node(EPIMask(), name="epi_mask") # fmt: off workflow.connect([ (inputnode, flatten, [("in_data", "in_data"), @@ -118,7 +129,9 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"): ("out_jacs", "jacobians"), ("out_mats", "xfms"), ("out_warps", "out_warps")]), + (merge_corrected, epi_mask, [("out_avg", "in_file")]), (merge_corrected, outputnode, [("out_avg", "fmap_ref")]), + (epi_mask, outputnode, [("out_file", "fmap_mask")]), ]) # fmt: on diff --git a/sdcflows/workflows/fit/tests/test_pepolar.py b/sdcflows/workflows/fit/tests/test_pepolar.py index 8a27da2131..3499013880 100644 --- a/sdcflows/workflows/fit/tests/test_pepolar.py +++ b/sdcflows/workflows/fit/tests/test_pepolar.py @@ -44,7 +44,6 @@ def test_topup_wf(tmpdir, datadir, workdir, outdir, epi_file): topup_wf.inputs.inputnode.metadata = metadata if outdir: - from nipype.interfaces.afni import Automask from ...outputs import init_fmap_derivatives_wf, init_fmap_reports_wf outdir = outdir / "unittests" / epi_file[0].split("/")[0] @@ -62,14 +61,11 @@ def test_topup_wf(tmpdir, datadir, workdir, outdir, epi_file): ) fmap_reports_wf.inputs.inputnode.source_files = in_data - pre_mask = pe.Node(Automask(dilate=1, outputtype="NIFTI_GZ"), name="pre_mask") - # fmt: off wf.connect([ - (topup_wf, pre_mask, [("outputnode.fmap_ref", "in_file")]), - (pre_mask, fmap_reports_wf, [("out_file", "inputnode.fmap_mask")]), (topup_wf, fmap_reports_wf, [("outputnode.fmap", "inputnode.fieldmap"), - ("outputnode.fmap_ref", "inputnode.fmap_ref")]), + ("outputnode.fmap_ref", "inputnode.fmap_ref"), + ("outputnode.fmap_mask", "inputnode.fmap_mask")]), (topup_wf, fmap_derivatives_wf, [ ("outputnode.fmap", "inputnode.fieldmap"), ("outputnode.fmap_ref", "inputnode.fmap_ref"), diff --git a/sdcflows/workflows/tests/test_base.py b/sdcflows/workflows/tests/test_base.py index 8c3126775f..64bc82e6da 100644 --- a/sdcflows/workflows/tests/test_base.py +++ b/sdcflows/workflows/tests/test_base.py @@ -1,4 +1,5 @@ """Test the base workflow.""" +from pathlib import Path import os import pytest from ... import fieldmaps as fm @@ -13,6 +14,9 @@ ) def test_fmap_wf(tmpdir, workdir, outdir, bids_layouts, dataset, subject): """Test the encompassing of the wrangler and the workflow creator.""" + if outdir is None: + outdir = Path(str(tmpdir)) + outdir = outdir / "test_base" / dataset fm._estimators.clear() estimators = find_estimators(bids_layouts[dataset], subject=subject)