RF: Filter fieldmaps based on whether they will be used to correct a BOLD series (#3025)

This PR takes a chunk out of the `init_func_preproc_wf` function and
pre-runs it in `init_single_subject_wf` to filter out fieldmaps that are
not needed to correct BOLD files.

Closes #2968.
Closes nipreps/sdcflows#359.

Author: effigies

fmriprep/workflows/base.py

@@ -36,12 +36,13 @@
 
 from nipype.interfaces import utility as niu
 from nipype.pipeline import engine as pe
+from niworkflows.utils.connections import listify
 from packaging.version import Version
 
 from .. import config
 from ..interfaces import DerivativesDataSink
 from ..interfaces.reports import AboutSummary, SubjectSummary
-from .bold import init_func_preproc_wf
+from .bold.base import get_estimator, init_func_preproc_wf
 
 
 def init_fmriprep_wf():
@@ -418,6 +419,16 @@ def init_single_subject_wf(subject_id: str):
             )
             fmap_estimators = [f for f in fmap_estimators if f.method == fm.EstimatorType.ANAT]
 
+        # Do not calculate fieldmaps that we will not use
+        if fmap_estimators:
+            used_estimators = {
+                key
+                for bold_file in subject_data['bold']
+                for key in get_estimator(config.execution.layout, listify(bold_file)[0])
+            }
+
+            fmap_estimators = [fmap for fmap in fmap_estimators if fmap.bids_id in used_estimators]
+
         if fmap_estimators:
             config.loggers.workflow.info(
                 "B0 field inhomogeneity map will be estimated with "

fmriprep/workflows/bold/base.py

@@ -283,22 +283,7 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     config.loggers.workflow.info(sbref_msg)
 
     if has_fieldmap:
-        # First check if specified via B0FieldSource
-        estimator_key = listify(metadata.get("B0FieldSource"))
-
-        if not estimator_key:
-            import re
-            from pathlib import Path
-
-            from sdcflows.fieldmaps import get_identifier
-
-            # Fallback to IntendedFor
-            intended_rel = re.sub(
-                r"^sub-[a-zA-Z0-9]*/",
-                "",
-                str(Path(bold_file if not multiecho else bold_file[0]).relative_to(layout.root)),
-            )
-            estimator_key = get_identifier(intended_rel)
+        estimator_key = get_estimator(layout, bold_file if not multiecho else bold_file[0])
 
         if not estimator_key:
             has_fieldmap = False
@@ -1356,3 +1341,21 @@ def get_img_orientation(imgf):
     """Return the image orientation as a string"""
     img = nb.load(imgf)
     return "".join(nb.aff2axcodes(img.affine))
+
+
+def get_estimator(layout, fname):
+    field_source = layout.get_metadata(fname).get("B0FieldSource")
+    if isinstance(field_source, str):
+        field_source = (field_source,)
+
+    if field_source is None:
+        import re
+        from pathlib import Path
+
+        from sdcflows.fieldmaps import get_identifier
+
+        # Fallback to IntendedFor
+        intended_rel = re.sub(r"^sub-[a-zA-Z0-9]*/", "", str(Path(fname).relative_to(layout.root)))
+        field_source = get_identifier(intended_rel)
+
+    return field_source

fmriprep/workflows/bold/tests/test_base.py

@@ -0,0 +1,125 @@
+from copy import deepcopy
+
+import bids
+import pytest
+from niworkflows.utils.testing import generate_bids_skeleton
+from sdcflows.fieldmaps import clear_registry
+from sdcflows.utils.wrangler import find_estimators
+
+from ..base import get_estimator
+
+BASE_LAYOUT = {
+    "01": {
+        "anat": [{"suffix": "T1w"}],
+        "func": [
+            {
+                "task": "rest",
+                "run": i,
+                "suffix": "bold",
+                "metadata": {"PhaseEncodingDirection": "j", "TotalReadoutTime": 0.6},
+            }
+            for i in range(1, 3)
+        ],
+        "fmap": [
+            {"suffix": "phasediff", "metadata": {"EchoTime1": 0.005, "EchoTime2": 0.007}},
+            {"suffix": "magnitude1", "metadata": {"EchoTime": 0.005}},
+            {
+                "suffix": "epi",
+                "direction": "PA",
+                "metadata": {"PhaseEncodingDirection": "j", "TotalReadoutTime": 0.6},
+            },
+            {
+                "suffix": "epi",
+                "direction": "AP",
+                "metadata": {"PhaseEncodingDirection": "j-", "TotalReadoutTime": 0.6},
+            },
+        ],
+    },
+}
+
+
+def test_get_estimator_none(tmp_path):
+    bids_dir = tmp_path / "bids"
+
+    # No IntendedFors/B0Fields
+    generate_bids_skeleton(bids_dir, BASE_LAYOUT)
+    layout = bids.BIDSLayout(bids_dir)
+    bold_files = sorted(layout.get(suffix='bold', extension='.nii.gz', return_type='file'))
+
+    assert get_estimator(layout, bold_files[0]) == ()
+    assert get_estimator(layout, bold_files[1]) == ()
+
+
+def test_get_estimator_b0field_and_intendedfor(tmp_path):
+    bids_dir = tmp_path / "bids"
+
+    # Set B0FieldSource for run 1
+    spec = deepcopy(BASE_LAYOUT)
+    spec['01']['func'][0]['metadata']['B0FieldSource'] = 'epi'
+    spec['01']['fmap'][2]['metadata']['B0FieldIdentifier'] = 'epi'
+    spec['01']['fmap'][3]['metadata']['B0FieldIdentifier'] = 'epi'
+
+    # Set IntendedFor for run 2
+    spec['01']['fmap'][0]['metadata']['IntendedFor'] = 'func/sub-01_task-rest_run-2_bold.nii.gz'
+
+    generate_bids_skeleton(bids_dir, spec)
+    layout = bids.BIDSLayout(bids_dir)
+    estimators = find_estimators(layout=layout, subject='01')
+
+    bold_files = sorted(layout.get(suffix='bold', extension='.nii.gz', return_type='file'))
+
+    assert get_estimator(layout, bold_files[0]) == ('epi',)
+    assert get_estimator(layout, bold_files[1]) == ('auto_00000',)
+    clear_registry()
+
+
+def test_get_estimator_overlapping_specs(tmp_path):
+    bids_dir = tmp_path / "bids"
+
+    # Set B0FieldSource for both runs
+    spec = deepcopy(BASE_LAYOUT)
+    spec['01']['func'][0]['metadata']['B0FieldSource'] = 'epi'
+    spec['01']['func'][1]['metadata']['B0FieldSource'] = 'epi'
+    spec['01']['fmap'][2]['metadata']['B0FieldIdentifier'] = 'epi'
+    spec['01']['fmap'][3]['metadata']['B0FieldIdentifier'] = 'epi'
+
+    # Set IntendedFor for both runs
+    spec['01']['fmap'][0]['metadata']['IntendedFor'] = [
+        'func/sub-01_task-rest_run-1_bold.nii.gz',
+        'func/sub-01_task-rest_run-2_bold.nii.gz',
+    ]
+
+    generate_bids_skeleton(bids_dir, spec)
+    layout = bids.BIDSLayout(bids_dir)
+    estimators = find_estimators(layout=layout, subject='01')
+
+    bold_files = sorted(layout.get(suffix='bold', extension='.nii.gz', return_type='file'))
+
+    # B0Fields take precedence
+    assert get_estimator(layout, bold_files[0]) == ('epi',)
+    assert get_estimator(layout, bold_files[1]) == ('epi',)
+    clear_registry()
+
+
+def test_get_estimator_multiple_b0fields(tmp_path):
+    bids_dir = tmp_path / "bids"
+
+    # Set B0FieldSource for both runs
+    spec = deepcopy(BASE_LAYOUT)
+    spec['01']['func'][0]['metadata']['B0FieldSource'] = ('epi', 'phasediff')
+    spec['01']['func'][1]['metadata']['B0FieldSource'] = 'epi'
+    spec['01']['fmap'][0]['metadata']['B0FieldIdentifier'] = 'phasediff'
+    spec['01']['fmap'][1]['metadata']['B0FieldIdentifier'] = 'phasediff'
+    spec['01']['fmap'][2]['metadata']['B0FieldIdentifier'] = 'epi'
+    spec['01']['fmap'][3]['metadata']['B0FieldIdentifier'] = 'epi'
+
+    generate_bids_skeleton(bids_dir, spec)
+    layout = bids.BIDSLayout(bids_dir)
+    estimators = find_estimators(layout=layout, subject='01')
+
+    bold_files = sorted(layout.get(suffix='bold', extension='.nii.gz', return_type='file'))
+
+    # Always get an iterable; don't care if it's a list or tuple
+    assert get_estimator(layout, bold_files[0]) == ['epi', 'phasediff']
+    assert get_estimator(layout, bold_files[1]) == ('epi',)
+    clear_registry()