Brainvision

bidsreader/_neurorad_algo.py

@@ -0,0 +1,81 @@
+"""Pair-location algorithms copied verbatim from the CML neurorad pipeline.
+
+Source of truth:
+    github.com/pennmem/neurorad_pipeline
+    File: localization.py :: Localization.get_pair_coordinate
+    Pinned revision (on rhino2): cf563def3087126b1aedb1dea274172200f3fa75
+    Path checked: /home2/iped/neurorad_pipeline/localization.py
+
+We copy rather than import to keep neurorad_pipeline off the runtime
+dependency list. A test in tests/test_pair_midpoint_parity.py imports
+neurorad_pipeline as a dev-only dependency and pins this copy against
+the upstream implementation.
+
+When the upstream algorithm changes, update this file AND the pinned
+revision in the parity test's docstring in the same commit.
+"""
+
+from __future__ import annotations
+
+from typing import Dict, Tuple
+
+import numpy as np
+
+
+# Mirrors Localization.VALID_COORDINATE_SPACES (localization.py:13-19).
+VALID_COORDINATE_SPACES: Tuple[str, ...] = (
+    "ct_voxel",
+    "fs",
+    "t1_mri",
+    "t2_mri",
+    "mni",
+)
+
+# Mirrors Localization.VALID_COORDINATE_TYPES (localization.py:21-24).
+VALID_COORDINATE_TYPES: Tuple[str, ...] = (
+    "raw",
+    "corrected",
+)
+
+# Mirrors Localization.VALID_ATLASES (localization.py:35-39).
+VALID_ATLASES: Tuple[str, ...] = (
+    "dk",
+    "whole_brain",
+    "mtl",
+)
+
+# BIDS space label -> (neurorad coordinate_space, coordinate_type).
+# Synthesized from event_creation/submission/neurorad_tasks.py
+# FIELD_NAMES_TABLE (lines 190-198) and the bids-convert CML_TO_BIDS_SPACE
+# map. Spaces not present here (Talairach, t1MRI, fsnativeDural) are
+# valid midpoint targets but have no neurorad atlas lookup.
+BIDS_SPACE_TO_NEURORAD: Dict[str, Tuple[str, str]] = {
+    "fsnative":            ("fs",        "raw"),
+    "fsnativeBrainshift":  ("fs",        "corrected"),
+    "fsaverage":           ("fsaverage", "raw"),
+    "fsaverageBrainshift": ("fsaverage", "corrected"),
+    "MNI152NLin6ASym":     ("mni",       "raw"),
+    "Pixels":              ("ct_voxel",  "raw"),
+}
+
+
+def pair_coordinate(coord1, coord2):
+    """Midpoint of two contact coordinates.
+
+    Verbatim behavior of
+    neurorad_pipeline.localization.Localization.get_pair_coordinate
+    (localization.py:274-289), stripped of its Localization-object
+    plumbing so callers can pass raw arrays.
+    """
+    if coord1 is None or coord2 is None:
+        return None
+    return (np.array(coord1) + np.array(coord2)) / 2
+
+
+def pair_coordinate_axis(a, b):
+    """Vectorized midpoint for one coordinate axis across N pairs.
+
+    NaN on either side propagates, matching the None-returning behavior
+    of scalar pair_coordinate. Intended for pandas/numpy column math.
+    """
+    return (np.asarray(a, dtype=float) + np.asarray(b, dtype=float)) / 2.0

bidsreader/cmlbidsreader.py

@@ -55,7 +55,7 @@ def _determine_device(self) -> Optional[str]:
         return None
 
     DEFAULT_SPACE = "MNI152NLin6ASym"
-    SPACE_PREFERENCE = ("MNI152NLin6ASym", "Talairach" ,"fsaverage", "Pixels", "fsaverageBrainshift" "fsnative", "fsnativeBrainshift", "fsnativeDural", "t1MRI")
+    SPACE_PREFERENCE = ("MNI152NLin6ASym", "Talairach", "fsaverage", "Pixels", "fsaverageBrainshift", "fsnative", "fsnativeBrainshift", "fsnativeDural", "t1MRI")
 
     def _coordsystem_dir(self) -> Path:
         subject_root = self._subject_root()
@@ -104,16 +104,24 @@ def _determine_space(self) -> str:
         if len(spaces) == 1:
             return spaces[0]
 
+        # Multiple spaces present: pick the highest-priority preferred space
+        # that is available.
         for preferred in self.SPACE_PREFERENCE:
             if preferred in spaces:
                 return preferred
 
-        raise AmbiguousMatchError(
+        # None of the preferred defaults is present. Rather than make the
+        # caller pass space= explicitly, fall back to a default: DEFAULT_SPACE
+        # if it happens to be available, otherwise the first space (sorted).
+        chosen = self.DEFAULT_SPACE if self.DEFAULT_SPACE in spaces else spaces[0]
+        warnings.warn(
             f"determine_space: multiple spaces found and none of the "
             f"preferred defaults {self.SPACE_PREFERENCE} present. "
-            f"Available spaces: {spaces}. "
-            f"Pass space=<one of these> when constructing CMLBIDSReader."
+            f"Available spaces: {spaces}. Defaulting to {chosen!r}. "
+            f"Pass space=<one of these> to choose explicitly.",
+            RuntimeWarning,
         )
+        return chosen
 
     def _validate_acq(self, acquisition: Optional[str]) -> Optional[str]:
         if not self.is_intracranial():
@@ -128,7 +136,7 @@ def _get_needed_fields(self):
     def _attach_bipolar_midpoint_montage(self, raw: mne.io.BaseRaw, space: Optional[str] = None) -> None:
         pairs_df = self.load_channels("bipolar")
         elec_df = self.load_electrodes(space=space)
-        combo = combine_bipolar_electrodes(pairs_df, elec_df)
+        combo = combine_bipolar_electrodes(pairs_df, elec_df, space=self.space)
 
         if not {"name", "x_mid", "y_mid", "z_mid"}.issubset(combo.columns):
             return
@@ -208,7 +216,7 @@ def load_combined_channels(self, acquisition: Optional[str] = None, space: Optio
         if acquisition == "monopolar" or acquisition is None:
             return channel_df.merge(elec_df, on="name", how="left", suffixes=("", "_elec"))
         if acquisition == "bipolar":
-            return combine_bipolar_electrodes(channel_df, elec_df)
+            return combine_bipolar_electrodes(channel_df, elec_df, space=self.space)
 
     @public_api
     def load_coordsystem_desc(self, space: Optional[str] = None) -> Dict:

bidsreader/helpers.py

@@ -1,8 +1,27 @@
+"""Utility helpers for bidsreader.
+
+Pair-coordinate math (``combine_bipolar_electrodes``) uses
+``pair_coordinate_axis`` from :mod:`._neurorad_algo`, which is a
+verbatim copy of the CML neurorad pipeline's
+``Localization.get_pair_coordinate``. The rule is the same in every
+BIDS coordinate space we support — ``MNI152NLin6ASym``, ``Talairach``,
+``fsaverage``, ``fsaverageBrainshift``, ``fsnative``,
+``fsnativeBrainshift``, ``fsnativeDural``, ``t1MRI``, ``Pixels`` —
+because those per-space contact coordinates have already been
+produced by the pipeline's space-specific transforms before BIDS
+export.
+
+For the BIDS space ↔ CML ``coordinate_space`` / ``coordinate_type``
+mapping see ``eeg_validation.preparers.montage.CML_TO_BIDS_SPACE`` or
+:attr:`bidsreader._neurorad_algo.BIDS_SPACE_TO_NEURORAD`.
+"""
+
 import numpy as np
 import pandas as pd
 from typing import Iterable, Any, Tuple, Sequence, Optional, Dict
 import re
 from .exc import InvalidOptionError
+from ._neurorad_algo import pair_coordinate_axis
 
 def validate_option(name: str, value: Any, allowed: Iterable[Any]) -> Any:
     if value is None:
@@ -95,7 +114,23 @@ def combine_bipolar_electrodes(
         pair_col: str = "name",
         elec_name_col: str = "name",
         region_cols: Sequence[str] = ("wb.region", "ind.region", "stein.region"),
+        space: Optional[str] = None,
     ) -> pd.DataFrame:
+    """Join bipolar pairs with electrode metadata and compute per-pair
+    coordinate midpoints, matching the neurorad pipeline's pair-location
+    algorithm (``Localization.get_pair_coordinate`` — see
+    :mod:`bidsreader._neurorad_algo`). Midpoint is taken in whatever BIDS
+    space ``elec_df`` carries; upstream brainshift / nonlinear-warp
+    corrections are expected to have already been applied per-contact.
+
+    Region columns from ``region_cols`` are carried through per-contact
+    (as ``{col}_ch1`` / ``{col}_ch2``) but no ``{col}_pair`` column is
+    synthesized. Neurorad assigns pair-region labels by independently
+    looking up the atlas at each pair's midpoint voxel, which cannot be
+    reproduced from contact-level agreement. Callers that need
+    pair-region labels should pull them from the upstream
+    ``pairs.json`` via ``enrich_pairs_with_cml_regions``.
+    """
     sep = "-"
     out = pairs_df.copy()
 
@@ -122,20 +157,18 @@ def combine_bipolar_electrodes(
     look2 = look.add_suffix("_ch2").rename(columns={f"{elec_name_col}_ch2": "ch2"})
     out = out.merge(look2, on="ch2", how="left")
 
-    # Region agreement
-    for rc in region_cols:
-        if f"{rc}_ch1" in out.columns and f"{rc}_ch2" in out.columns:
-            a = out[f"{rc}_ch1"]
-            b = out[f"{rc}_ch2"]
-            out[f"{rc}_pair"] = np.where(a.notna() & (a == b), a, np.nan)
-
-    # Midpoints for every detected coordinate triplet
+    # Midpoints for every detected coordinate triplet. pair_coordinate_axis
+    # is a verbatim copy of neurorad_pipeline.localization.Localization
+    # .get_pair_coordinate reduced to a single axis.
     for prefix, (xcol, ycol, zcol) in coord_triplets.items():
         for col in (xcol, ycol, zcol):
             a = out[f"{col}_ch1"]
             b = out[f"{col}_ch2"]
             mid_name = f"{col}_mid"  # e.g., "x_mid" or "tal.x_mid"
-            out[mid_name] = np.where(a.notna() & b.notna(), (a + b) / 2.0, np.nan)
+            mid = pair_coordinate_axis(a, b)
+            if space == "Pixels":
+                mid = np.floor(mid)
+            out[mid_name] = np.where(a.notna() & b.notna(), mid, np.nan)
 
     return out
 

bidsreader/region_enrichment.py

@@ -0,0 +1,94 @@
+"""Opt-in enrichment that attaches pair-level region labels from
+the upstream CML ``pairs.json`` to a bidsreader pairs DataFrame.
+
+Pure-BIDS users don't need this module. It exists only for the hybrid
+case where a subject's BIDS export sits alongside the original
+neurorad-pipeline output on rhino and the caller wants the pair-region
+labels that neurorad computed (via independent atlas lookup at each
+pair's midpoint voxel). The BIDS side can't reproduce those labels
+from contact-level data, so we pull them from the upstream artifact.
+
+See :mod:`._neurorad_algo` for why pair region labels differ from the
+contact-level agreement heuristic.
+"""
+
+from __future__ import annotations
+
+from typing import Iterable, Optional
+
+import pandas as pd
+
+
+_DEFAULT_REGION_COLS: tuple[str, ...] = (
+    "ind.region",
+    "ind.corrected.region",
+    "avg.region",
+    "avg.corrected.region",
+    "mni.region",
+    "hcp.region",
+    "stein.region",
+    "das.region",
+)
+
+
+def enrich_pairs_with_cml_regions(
+    pairs_df: pd.DataFrame,
+    *,
+    subject: Optional[str] = None,
+    experiment: Optional[str] = None,
+    session: Optional[int] = None,
+    localization: Optional[int] = 0,
+    montage: Optional[int] = 0,
+    reader=None,
+    label_col: str = "name",
+    region_cols: Optional[Iterable[str]] = None,
+) -> pd.DataFrame:
+    """Return a copy of ``pairs_df`` with CML ``*.region`` columns joined on.
+
+    ``pairs_df`` is the output of
+    :func:`bidsreader.helpers.combine_bipolar_electrodes`, which uses
+    ``name`` as the pair-label column by default.
+    """
+
+    try:
+        from cmlreaders import CMLReader
+    except ImportError as exc:
+        raise ImportError(
+            "enrich_pairs_with_cml_regions requires cmlreaders. Install "
+            "it (pip install cmlreaders) or load pairs.json manually."
+        ) from exc
+
+    if reader is None:
+        if subject is None or experiment is None or session is None:
+            raise ValueError(
+                "Provide either a ready CMLReader or "
+                "(subject, experiment, session)."
+            )
+        reader = CMLReader(
+            subject=subject,
+            experiment=experiment,
+            session=session,
+            localization=localization,
+            montage=montage,
+        )
+
+    cml_pairs = reader.load("pairs")
+
+    wanted = tuple(region_cols) if region_cols is not None else _DEFAULT_REGION_COLS
+    available = [c for c in wanted if c in cml_pairs.columns]
+
+    if "label" not in cml_pairs.columns:
+        raise KeyError(
+            "cml_pairs is missing the 'label' column — can't join on pair "
+            "name. Got columns: %r" % (list(cml_pairs.columns),)
+        )
+
+    right = cml_pairs[["label", *available]].copy()
+    right["label"] = right["label"].astype("string").str.strip()
+    if label_col != "label":
+        right = right.rename(columns={"label": label_col})
+
+    out = pairs_df.copy()
+    out[label_col] = out[label_col].astype("string").str.strip()
+
+    return out.merge(right, how="left", on=label_col)

tests/test_cmlbidsreader.py

@@ -142,13 +142,27 @@ def test_no_coordsystem_file_raises(self, tmp_root):
         with pytest.raises(FileNotFoundBIDSError, match="no.*coordsystem.json"):
             r._determine_space()
 
-    def test_multiple_coordsystem_files_raises(self, tmp_root):
+    def test_multiple_spaces_prefers_preferred_default(self, tmp_root):
+        """When several spaces are present, the highest-priority space in
+        SPACE_PREFERENCE wins (here MNI152NLin6ASym over fsnative)."""
+        data_dir = self._make_data_dir(tmp_root)
+        (data_dir / "sub-R1001P_ses-0_space-fsnative_coordsystem.json").touch()
+        (data_dir / "sub-R1001P_ses-0_space-MNI152NLin6ASym_coordsystem.json").touch()
+        r = CMLBIDSReader(root=tmp_root, subject="R1001P", task="FR1", session="0", device="ieeg")
+        assert r._determine_space() == "MNI152NLin6ASym"
+
+    def test_multiple_spaces_none_preferred_falls_back_with_warning(self, tmp_root):
+        """When multiple spaces are present and none is in SPACE_PREFERENCE,
+        _determine_space falls back to a default (warning) instead of raising.
+        Neither MNI nor TAL is a canonical preferred label, so it picks the
+        first available (sorted)."""
         data_dir = self._make_data_dir(tmp_root)
         (data_dir / "sub-R1001P_ses-0_space-MNI_coordsystem.json").touch()
         (data_dir / "sub-R1001P_ses-0_space-TAL_coordsystem.json").touch()
         r = CMLBIDSReader(root=tmp_root, subject="R1001P", task="FR1", session="0", device="ieeg")
-        with pytest.raises(AmbiguousMatchError, match="multiple coordsystem"):
-            r._determine_space()
+        with pytest.warns(RuntimeWarning, match="Defaulting to"):
+            chosen = r._determine_space()
+        assert chosen == "MNI"  # sorted(['MNI', 'TAL'])[0]
 
     def test_single_valid_match(self, tmp_root):
         data_dir = self._make_data_dir(tmp_root)