Colocated topomaps for OPM (#13144)

Co-authored-by: Eric Larson <[email protected]>

Author: harrisonritz

doc/changes/devel/13144.newfeature.rst

@@ -0,0 +1 @@
+Allow for ``topomap`` plotting of optically pumped MEG (OPM) sensors with overlapping channel locations. When channel locations overlap, plot the most radially oriented channel. By :newcontrib:`Harrison Ritz`.

doc/changes/names.inc

@@ -107,6 +107,7 @@
 .. _Hamid Maymandi: https://github.com/HamidMandi
 .. _Hamza Abdelhedi: https://github.com/BabaSanfour
 .. _Hari Bharadwaj: https://github.com/haribharadwaj
+.. _Harrison Ritz: https://github.com/harrisonritz
 .. _Hasrat Ali Arzoo: https://github.com/hasrat17
 .. _Henrich Kolkhorst: https://github.com/hekolk
 .. _Hongjiang Ye: https://github.com/hongjiang-ye

mne/channels/layout.py

@@ -902,7 +902,7 @@ def _auto_topomap_coords(info, picks, ignore_overlap, to_sphere, sphere):
     # Use channel locations if available
     locs3d = np.array([ch["loc"][:3] for ch in chs])
 
-    # If electrode locations are not available, use digization points
+    # If electrode locations are not available, use digitization points
     if not _check_ch_locs(info=info, picks=picks):
         logging.warning(
             "Did not find any electrode locations (in the info "
@@ -1089,7 +1089,7 @@ def _pair_grad_sensors(
         return picks
 
 
-def _merge_ch_data(data, ch_type, names, method="rms"):
+def _merge_ch_data(data, ch_type, names, method="rms", *, modality="opm"):
     """Merge data from channel pairs.
 
     Parameters
@@ -1102,6 +1102,8 @@ def _merge_ch_data(data, ch_type, names, method="rms"):
         List of channel names.
     method : str
         Can be 'rms' or 'mean'.
+    modality : str
+        The modality of the data, either 'grad', 'fnirs', or 'opm'
 
     Returns
     -------
@@ -1112,9 +1114,13 @@ def _merge_ch_data(data, ch_type, names, method="rms"):
     """
     if ch_type == "grad":
         data = _merge_grad_data(data, method)
-    else:
-        assert ch_type in _FNIRS_CH_TYPES_SPLIT
+    elif modality == "fnirs" or ch_type in _FNIRS_CH_TYPES_SPLIT:
         data, names = _merge_nirs_data(data, names)
+    elif modality == "opm" and ch_type == "mag":
+        data, names = _merge_opm_data(data, names)
+    else:
+        raise ValueError(f"Unknown modality {modality} for channel type {ch_type}")
+
     return data, names
 
 
@@ -1180,6 +1186,37 @@ def _merge_nirs_data(data, merged_names):
     return data, merged_names
 
 
+def _merge_opm_data(data, merged_names):
+    """Merge data from multiple opm channel by just using the radial component.
+
+    Channel names that end in "MERGE_REMOVE" (ie non-radial channels) will be
+    removed. Only the the radial channel is kept.
+
+    Parameters
+    ----------
+    data : array, shape = (n_channels, ..., n_times)
+        Data for channels.
+    merged_names : list
+        List of strings containing the channel names. Channels that are to be
+        removed end in "MERGE_REMOVE".
+
+    Returns
+    -------
+    data : array
+        Data for channels with requested channels merged. Channels used in the
+        merge are removed from the array.
+    """
+    to_remove = np.empty(0, dtype=np.int32)
+    for idx, ch in enumerate(merged_names):
+        if ch.endswith("MERGE-REMOVE"):
+            to_remove = np.append(to_remove, idx)
+    to_remove = np.unique(to_remove)
+    for rem in sorted(to_remove, reverse=True):
+        del merged_names[rem]
+    data = np.delete(data, to_remove, axis=0)
+    return data, merged_names
+
+
 def generate_2d_layout(
     xy,
     w=0.07,

mne/datasets/config.py

@@ -87,7 +87,7 @@
 # update the checksum in the MNE_DATASETS dict below, and change version
 # here: ↓↓↓↓↓↓↓↓
 RELEASES = dict(
-    testing="0.156",
+    testing="0.161",
     misc="0.27",
     phantom_kit="0.2",
     ucl_opm_auditory="0.2",
@@ -115,7 +115,7 @@
 # Testing and misc are at the top as they're updated most often
 MNE_DATASETS["testing"] = dict(
     archive_name=f"{TESTING_VERSIONED}.tar.gz",
-    hash="md5:d94fe9f3abe949a507eaeb865fb84a3f",
+    hash="md5:a32cfb9e098dec39a5f3ed6c0833580d",
     url=(
         "https://codeload.github.com/mne-tools/mne-testing-data/"
         f"tar.gz/{RELEASES['testing']}"

mne/preprocessing/tests/test_ica.py

@@ -1027,7 +1027,6 @@ def f(x, y):
 
 def test_get_explained_variance_ratio(tmp_path, short_raw_epochs):
     """Test ICA.get_explained_variance_ratio()."""
-    pytest.importorskip("sklearn")
     raw, epochs, _ = short_raw_epochs
     ica = ICA(max_iter=1)
 

mne/viz/tests/test_ica.py

@@ -17,8 +17,10 @@
     pick_types,
     read_cov,
     read_events,
+    read_evokeds,
 )
-from mne.io import read_raw_fif
+from mne.datasets import testing
+from mne.io import RawArray, read_raw_fif
 from mne.preprocessing import ICA, create_ecg_epochs, create_eog_epochs
 from mne.utils import _record_warnings, catch_logging
 from mne.viz.ica import _create_properties_layout, plot_ica_properties
@@ -32,6 +34,9 @@
 event_id, tmin, tmax = 1, -0.1, 0.2
 raw_ctf_fname = base_dir / "test_ctf_raw.fif"
 
+testing_path = testing.data_path(download=False)
+opm_fname = testing_path / "OPM" / "opm-evoked-ave.fif"
+
 pytest.importorskip("sklearn")
 
 
@@ -526,3 +531,15 @@ def test_plot_instance_components(browser_backend):
     fig._fake_click((x, y), xform="data")
     fig._click_ch_name(ch_index=0, button=1)
     fig._fake_keypress("escape")
+
+
+@pytest.mark.slowtest
+@pytest.mark.filterwarnings("ignore:.*did not converge.*:")
+@testing.requires_testing_data
+def test_plot_components_opm():
+    """Test for gh-12934."""
+    evoked = read_evokeds(opm_fname, kind="average")[0]
+    ica = ICA(max_iter=1, random_state=0, n_components=10)
+    ica.fit(RawArray(evoked.data, evoked.info), picks="mag", verbose="error")
+    fig = ica.plot_components()
+    assert len(fig.axes) == 10

mne/viz/tests/test_topomap.py

@@ -62,6 +62,8 @@
 subjects_dir = data_dir / "subjects"
 ecg_fname = data_dir / "MEG" / "sample" / "sample_audvis_ecg-proj.fif"
 triux_fname = data_dir / "SSS" / "TRIUX" / "triux_bmlhus_erm_raw.fif"
+opm_fname = data_dir / "OPM" / "opm-evoked-ave.fif"
+
 
 base_dir = Path(__file__).parents[2] / "io" / "tests" / "data"
 evoked_fname = base_dir / "test-ave.fif"
@@ -776,6 +778,19 @@ def test_plot_topomap_bads_grad():
     plot_topomap(data, info, res=8)
 
 
+@testing.requires_testing_data
+def test_plot_topomap_opm():
+    """Test plotting topomap with OPM data."""
+    # load data
+    evoked = read_evokeds(opm_fname, kind="average")[0]
+
+    # plot evoked topomap
+    fig_evoked = evoked.plot_topomap(
+        times=[-0.1, 0, 0.1, 0.2], ch_type="mag", show=False
+    )
+    assert len(fig_evoked.axes) == 5
+
+
 def test_plot_topomap_nirs_overlap(fnirs_epochs):
     """Test plotting nirs topomap with overlapping channels (gh-7414)."""
     fig = fnirs_epochs["A"].average(picks="hbo").plot_topomap()