Merge pull request #247 from nipreps/fix/b0field-identifier

FIX: Implement ``B0FieldIdentifier`` / ``B0FieldSource``

Author: oesteban

sdcflows/conftest.py

@@ -40,9 +40,10 @@
 }
 
 data_dir = Path(__file__).parent / "tests" / "data"
-
-layouts["dsA"] = BIDSLayout(data_dir / "dsA", validate=False, derivatives=False)
-layouts["dsB"] = BIDSLayout(data_dir / "dsB", validate=False, derivatives=False)
+layouts.update({
+    folder.name: BIDSLayout(folder, validate=False, derivatives=False)
+    for folder in data_dir.glob("ds*") if folder.is_dir()
+})
 
 
 def pytest_report_header(config):
@@ -66,6 +67,7 @@ def add_np(doctest_namespace):
 
     doctest_namespace["dsA_dir"] = data_dir / "dsA"
     doctest_namespace["dsB_dir"] = data_dir / "dsB"
+    doctest_namespace["dsC_dir"] = data_dir / "dsC"
 
 
 @pytest.fixture

sdcflows/tests/data/dsC/dataset_description.json

@@ -0,0 +1,11 @@
+{
+    "Name": "Test Dataset A, only empty files",
+    "BIDSVersion": "",
+    "License": "CC0",
+    "Authors": ["Esteban O."],
+    "Acknowledgements": "",
+    "HowToAcknowledge": "",
+    "Funding": "",
+    "ReferencesAndLinks": [""],
+    "DatasetDOI": ""
+}

sdcflows/tests/data/dsC/sub-01/anat/sub-01_FLAIR.nii.gz

@@ -0,0 +1,4 @@
+{
+    "PhaseEncodingDirection": "j-",
+    "TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/anat/sub-01_T1w.nii.gz

@@ -0,0 +1,4 @@
+{
+    "PhaseEncodingDirection": "i",
+    "TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/anat/sub-01_T2w.nii.gz

@@ -0,0 +1,4 @@
+{
+    "PhaseEncodingDirection": "j",
+    "TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-AP_dwi.nii.gz

@@ -0,0 +1,4 @@
+{
+    "PhaseEncodingDirection": "i-",
+    "TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-AP_sbref.json

@@ -0,0 +1,10 @@
+{
+    "PhaseEncodingDirection": "j",
+    "TotalReadoutTime": 0.005,
+    "IntendedFor": [
+        "dwi/sub-01_dir-AP_dwi.nii.gz",
+        "dwi/sub-01_dir-AP_sbref.nii.gz",
+        "func/sub-01_task-rest_bold.nii.gz",
+        "func/sub-01_task-rest_sbref.nii.gz"
+    ]
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-AP_sbref.nii.gz

@@ -0,0 +1,5 @@
+{
+	"B0FieldIdentifier": "pepolar4pe",
+	"PhaseEncodingDirection": "j-",
+	"TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-LR_dwi.nii.gz

@@ -0,0 +1,5 @@
+{
+	"B0FieldIdentifier": "pepolar4pe",
+	"PhaseEncodingDirection": "i",
+	"TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-LR_sbref.json

@@ -0,0 +1,5 @@
+{
+	"B0FieldIdentifier": "pepolar4pe",
+	"PhaseEncodingDirection": "j",
+	"TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-LR_sbref.nii.gz

@@ -0,0 +1,5 @@
+{
+	"B0FieldIdentifier": "pepolar4pe",
+	"PhaseEncodingDirection": "i-",
+	"TotalReadoutTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-PA_dwi.nii.gz

@@ -0,0 +1,3 @@
+{
+	"Units": "rad/s"
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-PA_sbref.json

@@ -0,0 +1,3 @@
+{
+	"EchoTime": 0.005
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-PA_sbref.nii.gz

@@ -0,0 +1,3 @@
+{
+	"EchoTime": 0.00746
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-RL_dwi.nii.gz

@@ -0,0 +1,4 @@
+{
+	"EchoTime1": 0.00500,
+	"EchoTime2": 0.00746
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-RL_sbref.json

@@ -0,0 +1,5 @@
+{
+	"TaskName": "Rest",
+	"PhaseEncodingDirection": "j-",
+	"TotalReadoutTime": 0.05
+}

sdcflows/tests/data/dsC/sub-01/dwi/sub-01_dir-RL_sbref.nii.gz

@@ -214,9 +214,35 @@ def find_estimators(*, layout, subject, fmapless=True, force_fmapless=False):
     FieldmapEstimation(sources=<2 files>, method=<EstimatorType.ANAT: 5>,
                        bids_id='auto_00012')]
 
+    When the ``B0FieldIdentifier`` metadata is set for one or more fieldmaps, then
+    the heuristics that use ``IntendedFor`` are dismissed:
+
+    >>> find_estimators(
+    ...     layout=layouts['dsC'],
+    ...     subject="01",
+    ... )  # doctest: +ELLIPSIS
+    [FieldmapEstimation(sources=<5 files>, method=<EstimatorType.PEPOLAR: 2>,
+                        bids_id='pepolar4pe')]
+
+    The only exception to the priority of ``B0FieldIdentifier`` is when fieldmaps
+    are searched with the ``force_fmapless`` argument on:
+
+    >>> fm.clear_registry()  # Necessary as `pepolar4pe` is not changing.
+    >>> find_estimators(
+    ...     layout=layouts['dsC'],
+    ...     subject="01",
+    ...     fmapless=True,
+    ...     force_fmapless=True,
+    ... )  # doctest: +ELLIPSIS
+    [FieldmapEstimation(sources=<5 files>, method=<EstimatorType.PEPOLAR: 2>,
+                        bids_id='pepolar4pe'),
+    FieldmapEstimation(sources=<2 files>, method=<EstimatorType.ANAT: 5>,
+                       bids_id='auto_00000')]
+
     """
     from .. import fieldmaps as fm
     from bids.layout import Query
+    from bids.exceptions import BIDSEntityError
 
     base_entities = {
         "subject": subject,
@@ -232,57 +258,74 @@ def find_estimators(*, layout, subject, fmapless=True, force_fmapless=False):
 
     estimators = []
 
-    # Set up B0 fieldmap strategies:
-    for fmap in layout.get(suffix=["fieldmap", "phasediff", "phase1"], **base_entities):
-        e = fm.FieldmapEstimation(
-            fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
-        )
-        estimators.append(e)
+    b0_ids = tuple()
+    with suppress(BIDSEntityError):
+        b0_ids = layout.get_B0FieldIdentifiers(**base_entities)
 
-    # A bunch of heuristics to select EPI fieldmaps
-    sessions = layout.get_sessions()
-    acqs = tuple(layout.get_acquisitions(suffix="epi") + [None])
-    contrasts = tuple(layout.get_ceagents(suffix="epi") + [None])
-
-    for ses, acq, ce in product(sessions or (None,), acqs, contrasts):
+    for b0_id in b0_ids:
+        # Found B0FieldIdentifier metadata entries
         entities = base_entities.copy()
-        entities.update(
-            {"suffix": "epi", "session": ses, "acquisition": acq, "ceagent": ce}
-        )
-        dirs = layout.get_directions(**entities)
-        if len(dirs) > 1:
+        entities["B0FieldIdentifier"] = b0_id
+
+        _e = fm.FieldmapEstimation([
+            fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
+            for fmap in layout.get(**entities)
+        ])
+        estimators.append(_e)
+
+    # As no B0FieldIdentifier(s) were found, try several heuristics
+    if not estimators:
+        # Set up B0 fieldmap strategies:
+        for fmap in layout.get(suffix=["fieldmap", "phasediff", "phase1"], **base_entities):
             e = fm.FieldmapEstimation(
-                [
-                    fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
-                    for fmap in layout.get(direction=dirs, **entities)
-                ]
+                fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
             )
             estimators.append(e)
 
-    # At this point, only single-PE _epi files WITH ``IntendedFor`` can be automatically processed
-    # (this will be easier with bids-standard/bids-specification#622 in).
-    has_intended = tuple()
-    with suppress(ValueError):
-        has_intended = layout.get(suffix="epi", IntendedFor=Query.ANY, **base_entities)
-
-    for epi_fmap in has_intended:
-        if epi_fmap.path in fm._estimators.sources:
-            continue  # skip EPI images already considered above
-
-        targets = [epi_fmap] + [
-            layout.get_file(str(subject_root / intent))
-            for intent in epi_fmap.get_metadata()["IntendedFor"]
-        ]
-
-        epi_sources = []
-        for fmap in targets:
-            with suppress(fm.MetadataError):
-                epi_sources.append(
-                    fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
+        # A bunch of heuristics to select EPI fieldmaps
+        sessions = layout.get_sessions(subject=subject)
+        acqs = tuple(layout.get_acquisitions(subject=subject, suffix="epi") + [None])
+        contrasts = tuple(layout.get_ceagents(subject=subject, suffix="epi") + [None])
+
+        for ses, acq, ce in product(sessions or (None,), acqs, contrasts):
+            entities = base_entities.copy()
+            entities.update(
+                {"suffix": "epi", "session": ses, "acquisition": acq, "ceagent": ce}
+            )
+            dirs = layout.get_directions(**entities)
+            if len(dirs) > 1:
+                e = fm.FieldmapEstimation(
+                    [
+                        fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
+                        for fmap in layout.get(direction=dirs, **entities)
+                    ]
                 )
+                estimators.append(e)
+
+        # At this point, only single-PE _epi files WITH ``IntendedFor`` can
+        # be automatically processed.
+        has_intended = tuple()
+        with suppress(ValueError):
+            has_intended = layout.get(suffix="epi", IntendedFor=Query.ANY, **base_entities)
+
+        for epi_fmap in has_intended:
+            if epi_fmap.path in fm._estimators.sources:
+                continue  # skip EPI images already considered above
+
+            targets = [epi_fmap] + [
+                layout.get_file(str(subject_root / intent))
+                for intent in epi_fmap.get_metadata()["IntendedFor"]
+            ]
+
+            epi_sources = []
+            for fmap in targets:
+                with suppress(fm.MetadataError):
+                    epi_sources.append(
+                        fm.FieldmapFile(fmap.path, metadata=fmap.get_metadata())
+                    )
 
-        with suppress(ValueError, TypeError):
-            estimators.append(fm.FieldmapEstimation(epi_sources))
+            with suppress(ValueError, TypeError):
+                estimators.append(fm.FieldmapEstimation(epi_sources))
 
     if estimators and not force_fmapless:
         fmapless = False
@@ -295,6 +338,8 @@ def find_estimators(*, layout, subject, fmapless=True, force_fmapless=False):
 
     from .epimanip import get_trt
 
+    # Sessions may not be defined at this point if some id was found.
+    sessions = layout.get_sessions(subject=subject)
     for ses, suffix in sorted(product(sessions or (None,), fmapless)):
         candidates = layout.get(suffix=suffix, session=ses, **base_entities)
 
@@ -306,6 +351,7 @@ def find_estimators(*, layout, subject, fmapless=True, force_fmapless=False):
         for candidate in candidates:
             meta = candidate.get_metadata()
             pe_dir = meta.get("PhaseEncodingDirection")
+
             if not pe_dir:
                 continue