Add more options to sphere argument in topomap plots. (#13400)

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

Author: wmvanvliet

doc/changes/dev/13400.newfeature.rst

@@ -0,0 +1 @@
+Add more options for the ``sphere`` parameter of :func:`mne.viz.plot_sensors`, by `Marijn van Vliet`_

mne/bem.py

@@ -1090,10 +1090,17 @@ def _fit_sphere_to_headshape(info, dig_kinds, *, verbose=None):
     o_mm = origin_head * 1e3
     o_d = origin_device * 1e3
     if np.linalg.norm(origin_head[:2]) > 0.02:
-        warn(
+        msg = (
             f"(X, Y) fit ({o_mm[0]:0.1f}, {o_mm[1]:0.1f}) "
             "more than 20 mm from head frame origin"
         )
+        if dig_kinds == "auto":
+            logger.info(msg)
+            logger.info("Trying again with all digitization points.")
+            return _fit_sphere_to_headshape(
+                info, dig_kinds=("extra", "eeg", "hpi", "cardinal"), verbose=verbose
+            )
+        warn(msg)
     logger.info(
         "Origin head coordinates:".ljust(30)
         + f"{o_mm[0]:0.1f} {o_mm[1]:0.1f} {o_mm[2]:0.1f} mm"

mne/source_space/tests/test_source_space.py

@@ -371,7 +371,7 @@ def test_volume_source_space(tmp_path):
     src = setup_volume_source_space(pos=10, sphere=(0.0, 0.0, 0.0, 0.09))
     src_new = setup_volume_source_space(pos=10, sphere=sphere)
     _compare_source_spaces(src, src_new, mode="exact")
-    with pytest.raises(ValueError, match="sphere, if str"):
+    with pytest.raises(ValueError, match="Invalid value for the 'sphere' parameter"):
         setup_volume_source_space(sphere="foo")
     # Need a radius
     sphere = make_sphere_model(head_radius=None)

mne/utils/check.py

@@ -1040,90 +1040,90 @@ def _check_sphere(sphere, info=None, sphere_units="m"):
                 sphere = "auto"
 
     if isinstance(sphere, str):
-        if sphere not in ("auto", "eeglab"):
-            raise ValueError(
-                f'sphere, if str, must be "auto" or "eeglab", got {sphere}'
-            )
+        _check_option(
+            "sphere", sphere, ("auto", "eeglab", "extra", "eeg", "cardinal", "hpi")
+        )
+
+    if isinstance(sphere, str) and sphere == "eeglab":
         assert info is not None
 
-        if sphere == "auto":
-            R, r0, _ = fit_sphere_to_headshape(
-                info, verbose=_verbose_safe_false(), units="m"
+        # We need coordinates for the 2D plane formed by
+        # Fpz<->Oz and T7<->T8, as this plane will be the horizon (i.e. it
+        # will determine the location of the head circle).
+        #
+        # We implement some special-handling in case Fpz is missing, as this seems to be
+        # a quite common situation in numerous EEG labs.
+        montage = info.get_montage()
+        if montage is None:
+            raise ValueError(
+                'No montage was set on your data, but sphere="eeglab" can only work if '
+                "digitization points for the EEG channels are available. Consider "
+                "calling set_montage() to apply a montage."
             )
-            sphere = tuple(r0) + (R,)
-            sphere_units = "m"
-        elif sphere == "eeglab":
-            # We need coordinates for the 2D plane formed by
-            # Fpz<->Oz and T7<->T8, as this plane will be the horizon (i.e. it
-            # will determine the location of the head circle).
-            #
-            # We implement some special-handling in case Fpz is missing, as
-            # this seems to be a quite common situation in numerous EEG labs.
-            montage = info.get_montage()
-            if montage is None:
-                raise ValueError(
-                    'No montage was set on your data, but sphere="eeglab" '
-                    "can only work if digitization points for the EEG "
-                    "channels are available. Consider calling set_montage() "
-                    "to apply a montage."
-                )
-            ch_pos = montage.get_positions()["ch_pos"]
-            horizon_ch_names = ("Fpz", "Oz", "T7", "T8")
-
-            if "FPz" in ch_pos:  # "fix" naming
-                ch_pos["Fpz"] = ch_pos["FPz"]
-                del ch_pos["FPz"]
-            elif "Fpz" not in ch_pos and "Oz" in ch_pos:
-                logger.info(
-                    "Approximating Fpz location by mirroring Oz along the X and Y axes."
+        ch_pos = montage.get_positions()["ch_pos"]
+        horizon_ch_names = ("Fpz", "Oz", "T7", "T8")
+
+        if "FPz" in ch_pos:  # "fix" naming
+            ch_pos["Fpz"] = ch_pos["FPz"]
+            del ch_pos["FPz"]
+        elif "Fpz" not in ch_pos and "Oz" in ch_pos:
+            logger.info(
+                "Approximating Fpz location by mirroring Oz along the X and Y axes."
+            )
+            # This assumes Fpz and Oz have the same Z coordinate
+            ch_pos["Fpz"] = ch_pos["Oz"] * [-1, -1, 1]
+
+        for ch_name in horizon_ch_names:
+            if ch_name not in ch_pos:
+                msg = (
+                    f'sphere="eeglab" requires digitization points of the following '
+                    f"electrode locations in the data: {', '.join(horizon_ch_names)}, "
+                    f"but could not find: {ch_name}"
                 )
-                # This assumes Fpz and Oz have the same Z coordinate
-                ch_pos["Fpz"] = ch_pos["Oz"] * [-1, -1, 1]
-
-            for ch_name in horizon_ch_names:
-                if ch_name not in ch_pos:
-                    msg = (
-                        f'sphere="eeglab" requires digitization points of '
-                        f"the following electrode locations in the data: "
-                        f"{', '.join(horizon_ch_names)}, but could not find: "
-                        f"{ch_name}"
-                    )
-                    if ch_name == "Fpz":
-                        msg += ", and was unable to approximate its location from Oz"
-                    raise ValueError(msg)
-
-            # Calculate the radius from: T7<->T8, Fpz<->Oz
-            radius = np.abs(
+                if ch_name == "Fpz":
+                    msg += ", and was unable to approximate its location from Oz"
+                raise ValueError(msg)
+
+        # Calculate the radius from: T7<->T8, Fpz<->Oz
+        radius = np.abs(
+            [
+                ch_pos["T7"][0],  # X axis
+                ch_pos["T8"][0],  # X axis
+                ch_pos["Fpz"][1],  # Y axis
+                ch_pos["Oz"][1],  # Y axis
+            ]
+        ).mean()
+
+        # Calculate the center of the head sphere.
+        # Use 4 digpoints for each of the 3 axes to hopefully get a better approximation
+        # than when using just 2 digpoints.
+        sphere_locs = dict()
+        for idx, axis in enumerate(("X", "Y", "Z")):
+            sphere_locs[axis] = np.mean(
                 [
-                    ch_pos["T7"][0],  # X axis
-                    ch_pos["T8"][0],  # X axis
-                    ch_pos["Fpz"][1],  # Y axis
-                    ch_pos["Oz"][1],  # Y axis
+                    ch_pos["T7"][idx],
+                    ch_pos["T8"][idx],
+                    ch_pos["Fpz"][idx],
+                    ch_pos["Oz"][idx],
                 ]
-            ).mean()
-
-            # Calculate the center of the head sphere
-            # Use 4 digpoints for each of the 3 axes to hopefully get a better
-            # approximation than when using just 2 digpoints.
-            sphere_locs = dict()
-            for idx, axis in enumerate(("X", "Y", "Z")):
-                sphere_locs[axis] = np.mean(
-                    [
-                        ch_pos["T7"][idx],
-                        ch_pos["T8"][idx],
-                        ch_pos["Fpz"][idx],
-                        ch_pos["Oz"][idx],
-                    ]
-                )
-            sphere = (sphere_locs["X"], sphere_locs["Y"], sphere_locs["Z"], radius)
-            sphere_units = "m"
-            del sphere_locs, radius, montage, ch_pos
+            )
+        sphere = (sphere_locs["X"], sphere_locs["Y"], sphere_locs["Z"], radius)
+        sphere_units = "m"
+        del sphere_locs, radius, montage, ch_pos
+    elif isinstance(sphere, str) or (
+        isinstance(sphere, list) and all(isinstance(s, str) for s in sphere)
+    ):
+        # Fit a sphere to the head points.
+        R, r0, _ = fit_sphere_to_headshape(
+            info, dig_kinds=sphere, verbose=_verbose_safe_false(), units="m"
+        )
+        sphere = tuple(r0) + (R,)
+        sphere_units = "m"
     elif isinstance(sphere, ConductorModel):
         if not sphere["is_sphere"] or len(sphere["layers"]) == 0:
             raise ValueError(
-                "sphere, if a ConductorModel, must be spherical "
-                "with multiple layers, not a BEM or single-layer "
-                f"sphere (got {sphere})"
+                "sphere, if a ConductorModel, must be spherical with multiple layers, "
+                f"not a BEM or single-layer sphere (got {sphere})"
             )
         sphere = tuple(sphere["r0"]) + (sphere["layers"][0]["rad"],)
         sphere_units = "m"
@@ -1132,8 +1132,8 @@ def _check_sphere(sphere, info=None, sphere_units="m"):
         sphere = np.concatenate([[0.0] * 3, [sphere]])
     if sphere.shape != (4,):
         raise ValueError(
-            "sphere must be float or 1D array of shape (4,), got "
-            f"array-like of shape {sphere.shape}"
+            "sphere must be float or 1D array of shape (4,), "
+            f"got array-like of shape {sphere.shape}"
         )
     _check_option("sphere_units", sphere_units, ("m", "mm"))
     if sphere_units == "mm":

mne/utils/docs.py

@@ -4217,21 +4217,40 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
 """
 
 docdict["sphere_topomap_auto"] = f"""\
-sphere : float | array-like | instance of ConductorModel | None  | 'auto' | 'eeglab'
-    The sphere parameters to use for the head outline. Can be array-like of
-    shape (4,) to give the X/Y/Z origin and radius in meters, or a single float
-    to give just the radius (origin assumed 0, 0, 0). Can also be an instance
-    of a spherical :class:`~mne.bem.ConductorModel` to use the origin and
-    radius from that object. If ``'auto'`` the sphere is fit to digitization
-    points. If ``'eeglab'`` the head circle is defined by EEG electrodes
-    ``'Fpz'``, ``'Oz'``, ``'T7'``, and ``'T8'`` (if ``'Fpz'`` is not present,
-    it will be approximated from the coordinates of ``'Oz'``). ``None`` (the
-    default) is equivalent to ``'auto'`` when enough extra digitization points
-    are available, and (0, 0, 0, {HEAD_SIZE_DEFAULT}) otherwise.
+sphere : float | array-like of float | instance of ConductorModel | str | list of str | None
+    The sphere parameters to use for the head outline.
+    Can be array-like of shape (4,) to give the X/Y/Z origin and radius in meters, or a
+    single float to give just the radius (origin assumed 0, 0, 0).
+    Can also be an instance of a spherical :class:`~mne.bem.ConductorModel` to use the
+    origin and radius from that object.
+    Can also be a ``str``, in which case:
+
+    - ``'auto'``: the sphere is fit to external digitization points first, and to
+      external + EEG digitization points if the former fails.
+
+    - ``'eeglab'``: the head circle is defined by EEG electrodes ``'Fpz'``, ``'Oz'``,
+      ``'T7'``, and ``'T8'`` (if ``'Fpz'`` is not present, it will be approximated from
+      the coordinates of ``'Oz'``).
+
+      - ``'extra'``: the sphere is fit to external digitization points.
+
+      - ``'eeg'``: the sphere is fit to EEG digitization points.
+
+      - ``'cardinal'``: the sphere is fit to cardinal digitization points.
+
+      - ``'hpi'``: the sphere is fit to HPI coil digitization points.
+
+    Can also be a list of ``str``, in which case the sphere is fit to the specified
+    digitization points, which can be any combination of ``'extra'``, ``'eeg'``,
+    ``'cardinal'``, and ``'hpi'``, as specified above.
+    ``None`` (the default) is equivalent to ``'auto'`` when enough extra digitization
+    points are available, and (0, 0, 0, {HEAD_SIZE_DEFAULT}) otherwise.
 
     .. versionadded:: 0.20
     .. versionchanged:: 1.1 Added ``'eeglab'`` option.
-"""
+    .. versionchanged:: 1.11 Added ``'extra'``, ``'eeg'``, ``'cardinal'``, ``'hpi'`` and
+       list of ``str`` options.
+"""  # noqa E501
 
 docdict["splash"] = """
 splash : bool

mne/utils/tests/test_check.py

@@ -10,9 +10,10 @@
 
 import numpy as np
 import pytest
+from numpy.testing import assert_allclose, assert_equal
 
 import mne
-from mne import pick_channels_cov, read_vectorview_selection
+from mne import create_info, pick_channels_cov, read_vectorview_selection
 from mne._fiff.pick import _picks_to_idx
 from mne.datasets import testing
 from mne.utils import (
@@ -364,15 +365,50 @@ def looks_stable(version):
 
 
 @testing.requires_testing_data
-def test_check_sphere_verbose():
-    """Test that verbose is handled properly in _check_sphere."""
+def test_check_sphere():
+    """Test the _check_sphere function."""
     info = mne.io.read_info(fname_raw)
-    with info._unlock():
-        info["dig"] = info["dig"][:20]
+    info_eeglab = create_info(
+        ch_names=["Fpz", "Oz", "T7", "T8"], sfreq=100, ch_types="eeg"
+    )
+    info_eeglab.set_montage("biosemi64")
+
+    # Test passing None.
+    assert_equal(_check_sphere(None), [0, 0, 0, 0.095])  # default head pos
+    assert not np.any(_check_sphere(None, info) == 0)  # fit to dig points
+
+    # Test passing a 4-element array-like as sphere parameter.
+    assert_equal(_check_sphere([1, 2, 3, 4], info), [1, 2, 3, 4])
+    assert_equal(_check_sphere([1, 2, 3, 4], info=None), [1, 2, 3, 4])
+    with pytest.raises(ValueError, match=r"1D array of shape \(4,\)"):
+        _check_sphere([1, 2, 3], info)
+
+    # Test passing various string values for `sphere`.
+    sphere_auto = _check_sphere("auto", info)
+    sphere_eeglab = _check_sphere("eeglab", info_eeglab)
+    sphere_extra = _check_sphere("extra", info)
+    sphere_eeg = _check_sphere("eeg", info)
+    with _record_warnings(), pytest.warns(RuntimeWarning, match="may be inaccurate"):
+        sphere_hpi = _check_sphere("hpi", info)
+    sphere_all = _check_sphere(["extra", "eeg", "cardinal", "hpi"], info)
+
+    assert_allclose(sphere_auto, sphere_extra)
+    assert not np.allclose(sphere_auto, sphere_eeglab, rtol=1e-4, atol=1e-4)
+    assert not np.allclose(sphere_auto, sphere_eeg, rtol=1e-4, atol=1e-4)
+    assert not np.allclose(sphere_auto, sphere_hpi, rtol=1e-4, atol=1e-4)
+    assert not np.allclose(sphere_auto, sphere_all, rtol=1e-4, atol=1e-4)
+
+    with pytest.raises(TypeError, match="Item must be an instance of Info"):
+        _check_sphere("auto", info=None)
+
+    # Test that verbose is handled properly in _check_sphere.
+    info_trunc = info.copy()
+    with info_trunc._unlock():
+        info_trunc["dig"] = info_trunc["dig"][:20]
     with _record_warnings(), pytest.warns(RuntimeWarning, match="may be inaccurate"):
-        _check_sphere("auto", info)
+        _check_sphere("auto", info_trunc)
     with mne.use_log_level("error"):
-        _check_sphere("auto", info)
+        _check_sphere("auto", info_trunc)
 
 
 def test_soft_import():