ENH: Add a simplistic EPI masking algorithm

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"

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'

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

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,19 +117,19 @@ 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()
 
     with pytest.raises(ValueError):
         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,15 +138,15 @@ 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",
                     "IntendedFor": "file1.nii.gz",
                 },
             ),
             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,15 +174,15 @@ 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",
                     "IntendedFor": ["file1.nii.gz", "file2.nii.gz"],
                 },
             ),
             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:

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

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

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