ENH: Expand BOLD coregistration targets

.circleci/bcp_full_outputs.txt

@@ -81,8 +81,8 @@ sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-boldref_to-T1w_
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-boldref_to-T1w_mode-image_desc-coreg_xfm.txt
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-boldref_to-auto00000_mode-image_desc-fmap_xfm.json
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-boldref_to-auto00000_mode-image_desc-fmap_xfm.txt
-sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-orig_to-boldref_mode-image_desc-hmc_xfm.json
-sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-orig_to-boldref_mode-image_desc-hmc_xfm.txt
+sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-orig_to-run_mode-image_desc-hmc_xfm.json
+sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_from-orig_to-run_mode-image_desc-hmc_xfm.txt
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_space-MNIInfant+1_boldref.json
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_space-MNIInfant+1_boldref.nii.gz
 sub-01/ses-1mo/func/sub-01_ses-1mo_task-rest_acq-PA_run-001_space-MNIInfant+1_desc-brain_mask.json

docs/outputs.md

@@ -249,17 +249,40 @@ The mask file is part of the *minimal* processing level. The BOLD series
 is only generated at the *full* processing level.
 :::
 
+#### BOLD reference spaces
+
+NiBabies defines a consistent hierarchy of BOLD reference spaces, reflected in
+the `space` and `from`/`to` entities of functional outputs:
+
+- **`orig`** — the native BOLD acquisition space.
+- **`run`** — the per-run boldref space. This is the target of HMC and is derived
+  from the sbref or a robust average of the BOLD series.
+- **`boldref`** — the final coregistration reference space registered to anatomy.
+  When `--bold-coreg-level run` (default), `boldref` is identical to `run`.
+  When `--bold-coreg-level session`, `boldref` is a template image built from all
+  run-level boldrefs.
+- **`anat`** — the anatomical reference space (T1w or T2w).
+
+The full transform chain applied during resampling is therefore:
+
+```
+orig → run      (from-orig_to-run, HMC)
+run  → boldref  (from-run_to-boldref, identity when --bold-coreg-level run)
+boldref → anat  (from-boldref_to-anat)
+anat → std      (anat2std)
+```
+
 #### Motion correction outputs
 
 Head-motion correction (HMC) produces a reference image to which all volumes
 are aligned, and a corresponding transform that maps the original BOLD series
-to the reference image::
+to the run-level boldref:
 
 ```
 sub-<subject_label>/[ses-<session_label>/]
   func/
     sub-<subject_label>_[specifiers]_desc-hmc_boldref.nii.gz
-    sub-<subject_label>_[specifiers]_from-orig_to_boldref_mode-image_desc-hmc_xfm.nii.gz
+    sub-<subject_label>_[specifiers]_from-orig_to-run_mode-image_desc-hmc_xfm.txt
 ```
 
 :::{note}
@@ -268,8 +291,15 @@ Motion correction outputs are part of the *minimal* processing level.
 
 #### Coregistration outputs
 
-Registration of the BOLD series to the anatomical reference image generates a further reference
-image and affine transform
+Registration of the BOLD series to the anatomical image generates a further reference
+image and affine transform:
+
+:::{tip}
+The coregistration reference used for anatomy registration depends on
+`--bold-coreg-level`. See [BOLD coregistration level](usage.md#bold-coregistration-level---bold-coreg-level)
+for details.
+:::
+
 
 ```
 sub-<subject_label>/[ses-<session_label>/]
@@ -291,7 +321,7 @@ is identified with `"B0FieldIdentifier": "TOPUP"`, the generated transform will
 ```
 sub-<subject_label>/[ses-<session_label>/]
   func/
-    sub-<subject_label>_[specifiers]_from-boldref_to-TOPUP_mode-image_xfm.nii.gz
+    sub-<subject_label>_[specifiers]_from-boldref_to-TOPUP_mode-image_xfm.txt
 ```
 
 If the association is discovered through the `IntendedFor` field of the

docs/usage.md

@@ -86,6 +86,33 @@ The segmentation directory layout should consist of one or more template directo
 
 :::
 
+### BOLD coregistration level (`--bold-coreg-level`)
+
+By default (`run`), each BOLD run is independently registered to the anatomical image.
+The `session` option instead builds a session-wide mean BOLD template from all runs before
+registering that template to the anatomical image.
+This can improve BOLD-to-anatomical registration quality by averaging signal across runs
+before coregistration.
+
+The session-level pipeline follows four steps:
+
+1. Head-motion correction (HMC) is applied per-volume within each run.
+2. Susceptibility distortion correction (SDC) is applied to each run's BOLD reference image.
+3. Each run's SDC-corrected boldref is registered to a session-wide BOLD template using FreeSurfer's `mri_robust_template`.
+4. The session template is registered to the anatomical image.
+
+The full transform chain applied during resampling is therefore:
+`bold volume → run boldref (HMC) → boldref (run2boldref) → T1w (boldref2anat) → standard space (anat2std)`.
+
+If `--bold-coreg-level session` is selected, NiBabies will first validate whether all runs can contribute to a common session template. If the data are incompatible (for example, when only some runs have SDC applied, or when all runs are SDC-less but have mixed phase-encoding directions), NiBabies falls back to `run`-level coregistration instead.
+
+:::{admonition} Output spaces and template-space BOLD
+:class: tip
+
+Adding `boldref` to `--output-spaces` when using `--bold-coreg-level session` saves each run
+resampled into the boldref space.
+:::
+
 ## Using the nibabies wrapper
 
 The wrapper will generate a Docker or Singularity command line for you, print it out for reporting purposes, and then execute it without further action needed.

nibabies/cli/parser.py

@@ -786,6 +786,12 @@ def _str_none(val):
         help='For certain adult templates (MNI152NLin6Asym), perform two step '
         'registrations (native -> MNIInfant -> template) and concatenate into a single xfm',
     )
+    g_baby.add_argument(
+        '--bold-coreg-level',
+        choices=('run', 'session'),
+        default='run',
+        help='Determine BOLD reference image for coregistration at run or session level.',
+    )
     return parser
 
 

nibabies/config.py

@@ -560,6 +560,8 @@ class workflow(_Config):
     bold2anat_init = 'auto'
     """Whether to use standard coregistration ('register') or to initialize coregistration from the
     BOLD image-header ('header')."""
+    bold_coreg_level = 'run'
+    """Level to coregister BOLD runs to a common reference, either 'run' (default) or 'session'."""
     cifti_output = None
     """Generate HCP Grayordinates, accepts either ``'91k'`` (default) or ``'170k'``."""
     dummy_scans = None
@@ -824,6 +826,7 @@ def dismiss_entities(entities: list | None = None) -> list:
         'surface_recon_method',
         'bold2anat_dof',
         'bold2anat_init',
+        'bold_coreg_level',
         'dummy_scans',
         'fd_radius',
         'fmap_bspline',

nibabies/data/io_spec_func.json

@@ -26,22 +26,31 @@
         "hmc": {
           "datatype": "func",
           "from": "orig",
-          "to": "boldref",
+          "to": "run",
           "mode": "image",
           "suffix": "xfm",
           "extension": ".txt"
         },
         "boldref2anat": {
           "datatype": "func",
-          "from": "orig",
+          "from": "boldref",
           "to": "anat",
           "mode": "image",
           "suffix": "xfm",
           "extension": ".txt"
         },
         "boldref2fmap": {
           "datatype": "func",
-          "from": "orig",
+          "from": "boldref",
+          "mode": "image",
+          "suffix": "xfm",
+          "extension": ".txt"
+        },
+        "run2boldref": {
+          "datatype": "func",
+          "desc": "coreg",
+          "from": "run",
+          "to": "boldref",
           "mode": "image",
           "suffix": "xfm",
           "extension": ".txt"

nibabies/interfaces/reports.py

@@ -51,6 +51,7 @@
 \t\t\t<li>Slice timing correction: {stc}</li>
 \t\t\t<li>Susceptibility distortion correction: {sdc}</li>
 \t\t\t<li>Registration: {registration}</li>
+\t\t\t<li>Registration level: {bold_coreg_level}</li>
 \t\t\t<li>Non-steady-state volumes: {dummy_scan_desc}</li>
 \t\t</ul>
 \t\t</details>
@@ -239,6 +240,9 @@ class FunctionalSummaryInputSpec(TraitedSpec):
         'FSL', 'FreeSurfer', mandatory=True, desc='Functional/anatomical registration method'
     )
     fallback = traits.Bool(desc='Boundary-based registration rejected')
+    bold_coreg_level = traits.Enum(
+        'run', 'session', desc='BOLD registration level (run or session)'
+    )
     registration_dof = traits.Enum(
         6, 9, 12, desc='Registration degrees of freedom', mandatory=True
     )
@@ -316,6 +320,7 @@ def _generate_segment(self):
             stc=stc,
             sdc=self.inputs.distortion_correction,
             registration=reg,
+            bold_coreg_level=self.inputs.bold_coreg_level,
             tr=self.inputs.tr,
             dummy_scan_desc=dummy_scan_msg,
             multiecho=multiecho,

nibabies/tests/test_config.py

@@ -150,7 +150,7 @@ def _load_spaces(age):
 
 def test_hash_config():
     # This may change with changes to config defaults / new attributes!
-    expected = 'cfee5aaf'
+    expected = '84226605'
     assert config.hash_config(config.get()) == expected
     _reset_config()
 

nibabies/utils/bids.py

@@ -362,6 +362,56 @@ def age_to_months(age: int | float, units: ty.Literal['weeks', 'months', 'years'
     return round(age)
 
 
+def is_valid_bold_template(
+    bold_runs: list[list[str]],
+    estimator_map: dict,
+    layout,
+) -> bool:
+    """Return True if BOLD template creation is possible.
+
+    Three independent conditions are checked:
+
+    1. **Insufficient runs** – fewer than two runs cannot form a meaningful
+       session template.
+
+    2. **Mixed SDC coverage** – if only some runs are distortion corrected,
+       the resulting boldrefs live in different geometric spaces even when their
+       phase-encoding directions agree.
+
+    3. **PE-direction heterogeneity** – when all runs are SDC-less, if more
+       than one distinct phase-encoding direction is present (including ``None``
+       for runs that lack the metadata), the distortions oppose each other in a
+       way that rigid/affine registration cannot recover.
+
+    Parameters
+    ----------
+    bold_runs
+        List of bold series (each series is a list of file paths).
+    estimator_map
+        Mapping of bold file → fieldmap estimator ID.  An absent key or None
+        value means no fieldmap / no SDC for that run.
+    layout
+        A :class:`bids.layout.BIDSLayout` instance used to read file metadata.
+
+    Returns
+    -------
+    compatible : bool
+        ``True`` if there are at least two runs, all share the same SDC status,
+        and SDC-less runs span exactly one phase-encoding direction.
+    """
+    if len(bold_runs) < 2:
+        return False
+
+    sdc_corrected = [bool(estimator_map.get(series[0])) for series in bold_runs]
+    if any(sdc_corrected):
+        return all(sdc_corrected)
+
+    pe_dirs = {
+        layout.get_metadata(series[0]).get('PhaseEncodingDirection') for series in bold_runs
+    }
+    return len(pe_dirs) == 1
+
+
 def combine_space(space, cohort) -> str:
     """Combine space and cohort into a single string.
 

nibabies/utils/derivatives.py

@@ -119,14 +119,31 @@ def collect_functional_derivatives(
             continue
         derivs_cache[f'{key}_boldref'] = item[0] if len(item) == 1 else item
 
+    # coreg_boldref is the legacy derivative name for run_boldref
+    if 'coreg_boldref' in derivs_cache and 'run_boldref' not in derivs_cache:
+        derivs_cache['run_boldref'] = derivs_cache['coreg_boldref']
+
     for xfm, qry in spec['transforms'].items():
         query = {**qry, **entities}
         if xfm == 'boldref2fmap':
             query['to'] = fieldmap_id
         item = layout.get(return_type='filename', **query)
+        if not item and xfm == 'hmc':
+            # legacy: from-orig_to-boldref (now from-orig_to-run)
+            item = layout.get(return_type='filename', **{**query, 'to': 'boldref'})
         if not item:
             continue
         derivs_cache[xfm] = item[0] if len(item) == 1 else item
+
+    # Session queries use {**entities, **qry} so spec null-values (run, task) override
+    # per-run entity values, ensuring only session-level files (lacking those entities) match.
+    for key, qry in spec.get('session', {}).items():
+        query = {**entities, **qry}
+        item = layout.get(return_type='filename', **query)
+        if not item:
+            continue
+        derivs_cache[f'session_{key}'] = item[0] if len(item) == 1 else item
+
     return derivs_cache
 
 

nibabies/utils/tests/test_bids.py

@@ -2,7 +2,46 @@
 
 import pytest
 
-from nibabies.utils.bids import _get_age_from_tsv
+from nibabies.utils.bids import _get_age_from_tsv, is_valid_bold_template
+
+
+class _MockLayout:
+    """Minimal BIDSLayout stub for bold_coreg_compat tests."""
+
+    def __init__(self, pe_map: dict):
+        self._pe_map = pe_map
+
+    def get_metadata(self, f):
+        pe = self._pe_map.get(f)
+        return {'PhaseEncodingDirection': pe} if pe is not None else {}
+
+
+@pytest.mark.parametrize(
+    ('bold_runs', 'estimator_map', 'pe_map', 'expected'),
+    [
+        ([], {}, {}, False),
+        ([['a.nii']], {}, {'a.nii': 'j'}, False),
+        ([['a.nii'], ['b.nii']], {'a.nii': 'fmap1', 'b.nii': 'fmap2'}, {}, True),
+        ([['a.nii'], ['b.nii']], {'a.nii': 'fmap1'}, {'a.nii': 'j', 'b.nii': 'j'}, False),
+        ([['a.nii'], ['b.nii']], {}, {'a.nii': 'j', 'b.nii': 'j'}, True),
+        ([['a.nii'], ['b.nii']], {}, {'a.nii': 'j', 'b.nii': 'j-'}, False),
+        ([['a.nii'], ['b.nii']], {}, {}, True),
+        ([['a.nii'], ['b.nii']], {}, {'a.nii': 'j'}, False),
+    ],
+    ids=[
+        'false-no_runs',
+        'false-one_run',
+        'true-all_sdc',
+        'false-mixed_sdc',
+        'true-no_sdc_single_pe',
+        'false-no_sdc_opposing_pe',
+        'true-no_sdc_no_pe',
+        'false-no_sdc_missing_pe',
+    ],
+)
+def test_is_valid_bold_template(bold_runs, estimator_map, pe_map, expected):
+    layout = _MockLayout(pe_map)
+    assert is_valid_bold_template(bold_runs, estimator_map, layout) is expected
 
 
 def create_tsv(data: dict, out_file: Path) -> None: