Repository: nipreps/fmriprep

Issue:
Limit fieldmap estimation to ones that will be used to correct BOLD files
### What happened?

fMRIprep uses both fmap "intended for" fmri and dwi sequences for SDC (susceptibility distortion correction, due to "Inhomogeneities of the B0 field"), instead of using only the fmap relative to fmri. The nifti file used passed without errors the BIDS-validator  check. Note: the fmri dataset analyzed is multi-echo, made up of three echoes (with different TE). The html report end with the sentence: "Errors: No errors to report!". All work fine but fmap. 

Thank you in advance for your help
Giovanni


### What command did you use?

```shell
/opt/conda/bin/fmriprep /data /out participant --participant-label 00244356 --nprocs 16 --md-only-boilerplate --write-graph --me-output-echos -w /scratch --output-spaces MNI152NLin2009cAsym:res-native func
```


### What version of fMRIPrep are you running?

22.1.1

### How are you running fMRIPrep?

Docker

### Is your data BIDS valid?

Yes

### Are you reusing any previously computed results?

No

### Please copy and paste any relevant log output.

```shell
You are using fMRIPrep-22.1.1, and a newer version of fMRIPrep is available: 23.0.0.
Please check out our documentation about how and when to upgrade:
https://fmriprep.readthedocs.io/en/latest/faq.html#upgrading
--topup-max-vols input is no longer used and the flag will be removed in a future release.
bids-validator@1.8.0

	 BOLD series will be slice-timing corrected to an offset of 0.442s.
230314-14:22:54,796 nipype.workflow IMPORTANT:
	 Generating T2* map and optimally combined ME-EPI time series.
230314-14:22:55,289 nipype.workflow INFO:
	 Setting-up fieldmap "fmapFMRI" (EstimatorType.PEPOLAR) with <sub-00244356_ses-20230314_acq-SpinEchoFieldMapFMRI_dir-AP_epi.nii.gz, sub-00244356_ses-20230314_acq-SpinEchoFieldMapFMRI_dir-PA_epi.nii.gz>
230314-14:22:55,327 nipype.workflow INFO:
	 Setting-up fieldmap "fmapDWI" (EstimatorType.PEPOLAR) with <sub-00244356_ses-20230314_acq-SpinEchoFieldMapDWI_dir-AP_epi.nii.gz, sub-00244356_ses-20230314_acq-SpinEchoFieldMapDWI_dir-PA_epi.nii.gz>
230314-14:23:23,101 nipype.workflow INFO:
	 fMRIPrep workflow graph with 487 nodes built successfully.
```


### Additional information / screenshots

_No response_

Discussion:
effigies: @bpinsard Would you have a chance to try out this branch on the dataset that was causing you issues?
effigies: @mgxd I'm comfortable using this with the added tests. It doesn't filter out an `IntendedFor` in ds054: https://app.circleci.com/pipelines/github/nipreps/fmriprep/1396/workflows/c26e914d-93fd-4e1e-ba4c-87569723c69e/jobs/7660/artifacts

Patch:
diff --git a/fmriprep/workflows/base.py b/fmriprep/workflows/base.py
index 4d167248..a27a0415 100644
--- a/fmriprep/workflows/base.py
+++ b/fmriprep/workflows/base.py
@@ -36,12 +36,13 @@ from copy import deepcopy
 
 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 @@ It is released under the [CC0]\
             )
             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 "
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 247a9b98..cdde346a 100644
--- a/fmriprep/workflows/bold/base.py
+++ b/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


Test patch:
diff --git a/fmriprep/workflows/bold/tests/test_base.py b/fmriprep/workflows/bold/tests/test_base.py
new file mode 100644
index 00000000..96ecffff
--- /dev/null
+++ b/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()