Merge branch 'main' into arm

doc/changes/devel.rst

@@ -56,10 +56,11 @@ Bugs
 - Fix bug with axis clip box boundaries in :func:`mne.viz.plot_evoked_topo` and related functions (:gh:`11999` by `Eric Larson`_)
 - Fix bug with ``subject_info`` when loading data from and exporting to EDF file (:gh:`11952` by `Paul Roujansky`_)
 - Fix bug with delayed checking of :class:`info["bads"] <mne.Info>` (:gh:`12038` by `Eric Larson`_)
-- Fix handling of channel information in annotations when loading data from and exporting to EDF file (:gh:`11960` :gh:`12017` by `Paul Roujansky`_)
+- Fix handling of channel information in annotations when loading data from and exporting to EDF file (:gh:`11960` :gh:`12017` :gh:`12044` by `Paul Roujansky`_)
 - Add missing ``overwrite`` and ``verbose`` parameters to :meth:`Transform.save() <mne.transforms.Transform.save>` (:gh:`12004` by `Marijn van Vliet`_)
 - Fix parsing of eye-link :class:`~mne.Annotations` when ``apply_offsets=False`` is provided to :func:`~mne.io.read_raw_eyelink` (:gh:`12003` by `Mathieu Scheltienne`_)
 - Correctly prune channel-specific :class:`~mne.Annotations` when creating :class:`~mne.Epochs` without the channel(s) included in the channel specific annotations (:gh:`12010` by `Mathieu Scheltienne`_)
+- Fix :func:`~mne.viz.plot_volume_source_estimates` with :class:`~mne.VolSourceEstimate` which include a list of vertices (:gh:`12025` by `Mathieu Scheltienne`_)
 - Correctly handle passing ``"eyegaze"`` or ``"pupil"`` to :meth:`mne.io.Raw.pick` (:gh:`12019` by `Scott Huberty`_)
 
 API changes

mne/io/edf/edf.py

@@ -206,6 +206,7 @@ def __init__(
             )
             annotations = _read_annotations_edf(
                 tal_data[0],
+                ch_names=info["ch_names"],
                 encoding=encoding,
             )
             self.set_annotations(annotations, on_missing="warn")
@@ -1892,25 +1893,21 @@ def read_raw_gdf(
 
 
 @fill_doc
-def _read_annotations_edf(annotations, encoding="utf8"):
+def _read_annotations_edf(annotations, ch_names=None, encoding="utf8"):
     """Annotation File Reader.
 
     Parameters
     ----------
     annotations : ndarray (n_chans, n_samples) | str
         Channel data in EDF+ TAL format or path to annotation file.
+    ch_names : list of string
+        List of channels' names.
     %(encoding_edf)s
 
     Returns
     -------
-    onset : array of float, shape (n_annotations,)
-        The starting time of annotations in seconds after ``orig_time``.
-    duration : array of float, shape (n_annotations,)
-        Durations of the annotations in seconds.
-    description : array of str, shape (n_annotations,)
-        Array of strings containing description for each annotation. If a
-        string, all the annotations are given the same description. To reject
-        epochs, use description starting with keyword 'bad'. See example above.
+    annot : instance of Annotations
+        The annotations.
     """
     pat = "([+-]\\d+\\.?\\d*)(\x15(\\d+\\.?\\d*))?(\x14.*?)\x14\x00"
     if isinstance(annotations, str):
@@ -1949,7 +1946,11 @@ def _read_annotations_edf(annotations, encoding="utf8"):
         duration = float(ev[2]) if ev[2] else 0
         for description in ev[3].split("\x14")[1:]:
             if description:
-                if "@@" in description:
+                if (
+                    "@@" in description
+                    and ch_names is not None
+                    and description.split("@@")[1] in ch_names
+                ):
                     description, ch_name = description.split("@@")
                     key = f"{onset}_{duration}_{description}"
                 else:
@@ -1979,22 +1980,20 @@ def _read_annotations_edf(annotations, encoding="utf8"):
                 offset = -onset
 
     if events:
-        onset, duration, description, ch_names = zip(*events.values())
+        onset, duration, description, annot_ch_names = zip(*events.values())
     else:
-        onset, duration, description, ch_names = list(), list(), list(), list()
+        onset, duration, description, annot_ch_names = list(), list(), list(), list()
 
-    assert len(onset) == len(duration) == len(description) == len(ch_names)
+    assert len(onset) == len(duration) == len(description) == len(annot_ch_names)
 
-    annotations = Annotations(
+    return Annotations(
         onset=onset,
         duration=duration,
         description=description,
         orig_time=None,
-        ch_names=ch_names,
+        ch_names=annot_ch_names,
     )
 
-    return annotations
-
 
 def _get_annotations_gdf(edf_info, sfreq):
     onset, duration, desc = list(), list(), list()

mne/io/edf/tests/test_edf.py

@@ -24,7 +24,7 @@
 import pytest
 
 from mne import pick_types, Annotations
-from mne.annotations import events_from_annotations, read_annotations
+from mne.annotations import _ndarray_ch_names, events_from_annotations, read_annotations
 from mne.datasets import testing
 from mne.io import read_raw_edf, read_raw_bdf, read_raw_fif, edf, read_raw_gdf
 from mne.io.tests.test_raw import _test_raw_reader
@@ -504,6 +504,81 @@ def test_read_utf8_annotations():
     assert raw.annotations[1]["description"] == "仰卧"
 
 
+def test_read_annotations_edf(tmp_path):
+    """Test reading annotations from EDF file."""
+    annot = (
+        b"+1.1\x14Event A@@CH1\x14\x00\x00"
+        b"+1.2\x14Event A\x14\x00\x00"
+        b"+1.3\x14Event B@@CH1\x14\x00\x00"
+        b"+1.3\x14Event B@@CH2\x14\x00\x00"
+        b"+1.4\x14Event A@@CH3\x14\x00\x00"
+        b"+1.5\x14Event B\x14\x00\x00"
+    )
+    annot_file = tmp_path / "annotations.edf"
+    with open(annot_file, "wb") as f:
+        f.write(annot)
+
+    # Test reading annotations from channel data
+    with open(annot_file, "rb") as f:
+        tal_channel = _read_ch(
+            f,
+            subtype="EDF",
+            dtype="<i2",
+            samp=-1,
+            dtype_byte=None,
+        )
+
+    # Read annotations without input channel names: annotations are left untouched and
+    # assigned as global
+    annotations = _read_annotations_edf(tal_channel, ch_names=None, encoding="latin1")
+    assert_allclose(annotations.onset, [1.1, 1.2, 1.3, 1.3, 1.4, 1.5])
+    assert not any(annotations.duration)  # all durations are 0
+    assert_array_equal(
+        annotations.description,
+        [
+            "Event A@@CH1",
+            "Event A",
+            "Event B@@CH1",
+            "Event B@@CH2",
+            "Event A@@CH3",
+            "Event B",
+        ],
+    )
+    assert_array_equal(
+        annotations.ch_names, _ndarray_ch_names([(), (), (), (), (), ()])
+    )
+
+    # Read annotations with complete input channel names: each annotation is parsed and
+    # associated to a channel
+    annotations = _read_annotations_edf(
+        tal_channel, ch_names=["CH1", "CH2", "CH3"], encoding="latin1"
+    )
+    assert_allclose(annotations.onset, [1.1, 1.2, 1.3, 1.4, 1.5])
+    assert not any(annotations.duration)  # all durations are 0
+    assert_array_equal(
+        annotations.description, ["Event A", "Event A", "Event B", "Event A", "Event B"]
+    )
+    assert_array_equal(
+        annotations.ch_names,
+        _ndarray_ch_names([("CH1",), (), ("CH1", "CH2"), ("CH3",), ()]),
+    )
+
+    # Read annotations with incomplete input channel names: "CH3" is missing from input
+    # channels, turning the related annotation into a global one
+    annotations = _read_annotations_edf(
+        tal_channel, ch_names=["CH1", "CH2"], encoding="latin1"
+    )
+    assert_allclose(annotations.onset, [1.1, 1.2, 1.3, 1.4, 1.5])
+    assert not any(annotations.duration)  # all durations are 0
+    assert_array_equal(
+        annotations.description,
+        ["Event A", "Event A", "Event B", "Event A@@CH3", "Event B"],
+    )
+    assert_array_equal(
+        annotations.ch_names, _ndarray_ch_names([("CH1",), (), ("CH1", "CH2"), (), ()])
+    )
+
+
 def test_read_latin1_annotations(tmp_path):
     """Test if annotations encoded as Latin-1 can be read.
 

mne/minimum_norm/inverse.py

@@ -145,6 +145,16 @@ def _repr_html_(self):
         )
         return html
 
+    @property
+    def ch_names(self):
+        """Name of channels attached to the inverse operator."""
+        return self["info"].ch_names
+
+    @property
+    def info(self):
+        """:class:`~mne.Info` attached to the inverse operator."""
+        return self["info"]
+
 
 def _pick_channels_inverse_operator(ch_names, inv):
     """Return data channel indices to be used knowing an inverse operator.

mne/viz/_3d.py

@@ -2650,10 +2650,9 @@ def plot_volume_source_estimates(
     %(subject_none)s
         If ``None``, ``stc.subject`` will be used.
     %(subjects_dir)s
-    mode : str
-        The plotting mode to use. Either 'stat_map' (default) or 'glass_brain'.
-        For "glass_brain", activation absolute values are displayed
-        after being transformed to a standard MNI brain.
+    mode : ``'stat_map'`` | ``'glass_brain'``
+        The plotting mode to use. For ``'glass_brain'``, activation absolute values are
+        displayed after being transformed to a standard MNI brain.
     bg_img : instance of SpatialImage | str
         The background image used in the nilearn plotting function.
         Can also be a string to use the ``bg_img`` file in the subject's
@@ -2714,10 +2713,11 @@ def plot_volume_source_estimates(
     >>> morph = mne.compute_source_morph(src_sample, subject_to='fsaverage')  # doctest: +SKIP
     >>> fig = stc_vol_sample.plot(morph)  # doctest: +SKIP
     """  # noqa: E501
-    from matplotlib import pyplot as plt, colors
     import nibabel as nib
-    from ..source_estimate import VolSourceEstimate
+    from matplotlib import pyplot as plt, colors
+
     from ..morph import SourceMorph
+    from ..source_estimate import VolSourceEstimate
     from ..source_space._source_space import _ensure_src
 
     if not check_version("nilearn", "0.4"):
@@ -2745,8 +2745,9 @@ def plot_volume_source_estimates(
         level="debug",
     )
     subject = _check_subject(src_subject, subject, first_kind=kind)
-    stc_ijk = np.array(np.unravel_index(stc.vertices[0], img.shape[:3], order="F")).T
-    assert stc_ijk.shape == (len(stc.vertices[0]), 3)
+    vertices = np.hstack(stc.vertices)
+    stc_ijk = np.array(np.unravel_index(vertices, img.shape[:3], order="F")).T
+    assert stc_ijk.shape == (vertices.size, 3)
     del kind
 
     # XXX this assumes zooms are uniform, should probably mult by zooms...
@@ -2756,12 +2757,11 @@ def _cut_coords_to_idx(cut_coords, img):
         """Convert voxel coordinates to index in stc.data."""
         ijk = _cut_coords_to_ijk(cut_coords, img)
         del cut_coords
-        logger.debug("    Affine remapped cut coords to [%d, %d, %d] idx" % tuple(ijk))
+        logger.debug("    Affine remapped cut coords to [%d, %d, %d] idx", tuple(ijk))
         dist, loc_idx = dist_to_verts.query(ijk[np.newaxis])
         dist, loc_idx = dist[0], loc_idx[0]
         logger.debug(
-            "    Using vertex %d at a distance of %d voxels"
-            % (stc.vertices[0][loc_idx], dist)
+            "    Using vertex %d at a distance of %d voxels", (vertices[loc_idx], dist)
         )
         return loc_idx
 
@@ -2848,7 +2848,7 @@ def _update_timeslice(idx, params):
         plot_map_callback(params["img_idx"], title="", cut_coords=cut_coords)
 
     def _update_vertlabel(loc_idx):
-        vert_legend.get_texts()[0].set_text(f"{stc.vertices[0][loc_idx]}")
+        vert_legend.get_texts()[0].set_text(f"{vertices[loc_idx]}")
 
     @verbose_dec
     def _onclick(event, params, verbose=None):
@@ -2932,7 +2932,7 @@ def _onclick(event, params, verbose=None):
             (stc.times[time_idx],)
             + tuple(cut_coords)
             + tuple(ijk)
-            + (stc.vertices[0][loc_idx],)
+            + (vertices[loc_idx],)
         )
     )
     del ijk
@@ -3046,8 +3046,7 @@ def plot_and_correct(*args, **kwargs):
 
     plot_and_correct(stat_map_img=params["img_idx"], title="", cut_coords=cut_coords)
 
-    if show:
-        plt.show()
+    plt_show(show)
     fig.canvas.mpl_connect(
         "button_press_event", partial(_onclick, params=params, verbose=verbose)
     )

mne/viz/tests/test_3d_mpl.py

@@ -13,13 +13,21 @@
 import pytest
 
 from mne import (
+    compute_covariance,
+    compute_source_morph,
+    make_fixed_length_epochs,
+    make_forward_solution,
+    read_bem_solution,
     read_forward_solution,
-    VolSourceEstimate,
+    read_trans,
+    setup_volume_source_space,
     SourceEstimate,
+    VolSourceEstimate,
     VolVectorSourceEstimate,
-    compute_source_morph,
 )
 from mne.datasets import testing
+from mne.io import read_raw_fif
+from mne.minimum_norm import apply_inverse, make_inverse_operator
 from mne.utils import catch_logging, _record_warnings
 from mne.viz import plot_volume_source_estimates
 from mne.viz.utils import _fake_click, _fake_keypress
@@ -148,3 +156,58 @@ def test_plot_volume_source_estimates_morph():
         stc.plot(
             sample_src, "sample", subjects_dir, clim=dict(lims=[-1, 2, 3], kind="value")
         )
+
+
+@testing.requires_testing_data
+def test_plot_volume_source_estimates_on_vol_labels():
+    """Test plot of source estimate on srcs setup on 2 labels."""
+    pytest.importorskip("nibabel")
+    pytest.importorskip("dipy")
+    pytest.importorskip("nilearn")
+    raw = read_raw_fif(
+        data_dir / "MEG" / "sample" / "sample_audvis_trunc_raw.fif", preload=False
+    )
+    raw.pick("meg").crop(0, 10)
+    raw.pick(raw.ch_names[::2]).del_proj().load_data()
+    epochs = make_fixed_length_epochs(raw, preload=True).apply_baseline((None, None))
+    evoked = epochs.average()
+    subject = "sample"
+    bem = read_bem_solution(
+        subjects_dir / f"{subject}" / "bem" / "sample-320-bem-sol.fif"
+    )
+    pos = 25.0  # spacing in mm
+    volume_label = [
+        "Right-Cerebral-Cortex",
+        "Left-Cerebral-Cortex",
+    ]
+    src = setup_volume_source_space(
+        subject,
+        subjects_dir=subjects_dir,
+        pos=pos,
+        mri=subjects_dir / subject / "mri" / "aseg.mgz",
+        bem=bem,
+        volume_label=volume_label,
+        add_interpolator=False,
+    )
+    trans = read_trans(data_dir / "MEG" / "sample" / "sample_audvis_trunc-trans.fif")
+    fwd = make_forward_solution(
+        evoked.info,
+        trans,
+        src,
+        bem,
+        meg=True,
+        eeg=False,
+        mindist=0,
+        n_jobs=1,
+    )
+    cov = compute_covariance(
+        epochs,
+        tmin=None,
+        tmax=None,
+        method="empirical",
+    )
+    inverse_operator = make_inverse_operator(evoked.info, fwd, cov, loose=1, depth=0.8)
+    stc = apply_inverse(
+        evoked, inverse_operator, 1.0 / 3**2, method="sLORETA", pick_ori=None
+    )
+    stc.plot(src, subject, subjects_dir, initial_time=0.03)