ENH: group triaxial OPM topomaps by orientation (#13866)

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com>

Author: 3 people

doc/changes/dev/13866.bugfix.rst

@@ -0,0 +1 @@
+Completed triaxial OPM topomap grouping by rendering separate radial and tangential maps in :func:`~mne.viz.plot_evoked_topomap` and :meth:`~mne.preprocessing.ICA.plot_components`, and fixed shape-mismatch crashes in :func:`~mne.viz.plot_projs_topomap` and :meth:`~mne.Evoked.animate_topomap` with colocated OPM channels, by `Pragnya Khandelwal`_.

examples/datasets/kernel_phantom.py

@@ -51,6 +51,11 @@
 t_peak = 0.016  # based on visual inspection of evoked
 fig.axes[0].axvline(t_peak, color="k", ls=":", lw=3, zorder=2)
 
+# %%
+# Because these OPM sensors are colocated in biaxial pairs, topomaps are
+# grouped into radial and tangential components.
+evoked.plot_topomap(times=[t_peak], ch_type="mag", show=True)
+
 # %%
 # The data covariance has an interesting structure because of densely packed sensors:
 
@@ -106,3 +111,9 @@
 )
 mne.viz.plot_dipole_locations(dipoles=dip, mode="arrow", color=(0.2, 1.0, 0.5), fig=fig)
 mne.viz.set_3d_view(figure=fig, azimuth=30, elevation=70, distance=0.4)
+
+# %%
+# For more information on OPM data visualization, see the OPM preprocessing
+# tutorial:
+#
+# - :ref:`tut-opm-processing`

mne/viz/evoked.py

@@ -1954,9 +1954,24 @@ def plot_evoked_joint(
     del times
     _, times_ts = _check_time_unit(ts_args["time_unit"], times_sec)
 
-    # prepare axes for topomap
+    ch_type = ch_types.pop()  # set should only contain one element
+    use_opm_orientation_groups = False
+    if ch_type == "mag":
+        from .topomap import _should_use_opm_orientation_groups
+
+        _, _, merge_channels, *_ = _prepare_topomap_plot(
+            evoked, ch_type, sphere=topomap_args.get("sphere")
+        )
+        use_opm_orientation_groups = _should_use_opm_orientation_groups(
+            merge_channels, ch_type
+        )
+    n_group_axes = 2 if use_opm_orientation_groups else 1
+
+    # prepare axes for topomap and butterfly plots
     if not got_axes:
-        fig, ts_ax, map_ax = _prepare_joint_axes(len(times_sec), figsize=(8.0, 4.2))
+        fig, ts_ax, map_ax = _prepare_joint_axes(
+            len(times_sec) * n_group_axes, figsize=(8.0, 4.2)
+        )
         cbar_ax = None
     else:
         ts_ax = ts_args["axes"]
@@ -2002,7 +2017,7 @@ def plot_evoked_joint(
 
     # topomap
     contours = topomap_args.get("contours", 6)
-    ch_type = ch_types.pop()  # set should only contain one element
+
     # Since the data has all the ch_types, we get the limits from the plot.
     vmin, vmax = (None, None)
     norm = ch_type == "grad"
@@ -2059,7 +2074,7 @@ def plot_evoked_joint(
             zorder=1,
             clip_on=False,
         )
-        fig.add_artist(con)
+        ts_ax.add_artist(con)
 
     # mark times in time series plot
     for timepoint in times_ts:

mne/viz/tests/test_ica.py

@@ -616,7 +616,8 @@ def test_plot_components_opm():
     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
+    # Biaxial OPM overlaps render grouped radial+tangential maps.
+    assert len(fig.axes) == 20
 
 
 @pytest.mark.slowtest
@@ -628,4 +629,4 @@ def test_plot_components_opm_triaxial(triaxial_raw):
     ica = ICA(max_iter=1, random_state=0, n_components=3)
     ica.fit(triaxial_raw, picks="mag", verbose="error")
     fig = ica.plot_components()
-    assert len(fig.axes) == 3
+    assert len(fig.axes) == 6

mne/viz/tests/test_topo.py

@@ -136,18 +136,6 @@ def return_inds(d):  # to test function kwarg to zorder arg of evoked.plot
     plt.close("all")
 
 
-def test_plot_joint_opm_triaxial(triaxial_evoked):
-    """Test joint plot with triaxial colocated OPM channels."""
-    fig = triaxial_evoked.plot_joint(
-        times=[0.0],
-        picks="mag",
-        show=False,
-        ts_args=dict(time_unit="s"),
-        topomap_args=dict(time_unit="s", contours=0, res=8, sensors=False),
-    )
-    assert len(fig.axes) >= 2
-
-
 def test_plot_topo():
     """Test plotting of ERP topography."""
     # Show topography

mne/viz/tests/test_topomap.py

@@ -812,7 +812,9 @@ def test_plot_topomap_opm():
     fig_evoked = evoked.plot_topomap(
         times=[-0.1, 0, 0.1, 0.2], ch_type="mag", show=False
     )
-    assert len(fig_evoked.axes) == 5
+    # Biaxial OPM pairs trigger grouped rendering
+    # (4 radial + 4 tangential + 2 colorbars)
+    assert len(fig_evoked.axes) == 10
 
 
 def test_prepare_topomap_plot_opm_non_quspin_coils():
@@ -869,6 +871,102 @@ def test_split_opm_overlaps(triaxial_evoked):
     assert tangential == ["OPM002", "OPM003", "OPM005", "OPM006"]
 
 
+def test_opm_tangential_rms_unsigned(triaxial_evoked):
+    """Test that tangential OPM data is RMS magnitude and unsigned."""
+    picks, pos, merge_channels, names, *_ = topomap._prepare_topomap_plot(
+        triaxial_evoked, "mag"
+    )
+    data = triaxial_evoked.data[picks]
+    grouped = topomap._compute_orientation_group_data(
+        data,
+        names,
+        pos,
+        ch_type="mag",
+        modality="opm",
+        merge_channels=merge_channels,
+        use_opm_orientation_groups=True,
+    )
+    tangential = [group for group in grouped if group[0] == "tangential"][0]
+    assert np.all(tangential[1] >= 0)
+    assert tangential[4]
+
+
+def test_should_use_opm_orientation_groups_only_for_triaxial():
+    """Test that OPM orientation grouping works for biaxial and triaxial overlaps."""
+    ch_names = [f"OPM{k:03}" for k in range(1, 7)]
+    info = create_info(ch_names, 1000.0, ch_types="mag")
+    with info._unlock():
+        for ch in info["chs"]:
+            ch["coil_type"] = FIFF.FIFFV_COIL_FIELDLINE_OPM_MAG_GEN1
+
+    pair_overlaps = [
+        np.array(["OPM001", "OPM002"]),
+        np.array(["OPM003", "OPM004"]),
+    ]
+    triax_overlaps = [
+        np.array(["OPM001", "OPM002", "OPM003"]),
+        np.array(["OPM004", "OPM005", "OPM006"]),
+    ]
+
+    # Both biaxial and triaxial overlaps should trigger grouping
+    assert topomap._should_use_opm_orientation_groups(pair_overlaps, "mag")
+    assert topomap._should_use_opm_orientation_groups(triax_overlaps, "mag")
+
+
+def test_plot_evoked_topomap_opm_triaxial_groups(triaxial_evoked):
+    """Test grouped radial/tangential topomap rendering for triaxial OPM."""
+    fig = triaxial_evoked.plot_topomap(
+        times=[0.0],
+        ch_type="mag",
+        contours=0,
+        res=8,
+        sensors=False,
+        show=False,
+    )
+    assert len(fig.axes) == 4
+    titles = [ax.get_title() for ax in fig.axes]
+    assert any("radial" in title for title in titles)
+    assert any("tangential" in title for title in titles)
+
+
+def test_plot_projs_topomap_opm(triaxial_evoked):
+    """Test plot_projs_topomap does not crash on colocated OPM channels (gh-13866)."""
+    from mne import compute_proj_evoked
+
+    projs = compute_proj_evoked(triaxial_evoked, n_mag=2)
+    # Should not raise a shape mismatch between data and pos
+    fig = plot_projs_topomap(projs, triaxial_evoked.info, show=False)
+    assert len(fig.axes) >= 1
+
+
+@pytest.mark.filterwarnings("ignore:.*No contour levels.*:UserWarning")
+def test_animate_topomap_opm(triaxial_evoked):
+    """Test animate_topomap does not crash on colocated OPM channels (gh-13866)."""
+    fig, anim = triaxial_evoked.animate_topomap(ch_type="mag", times=[0.0], show=False)
+    anim._func(0)
+    assert len(fig.axes) >= 1
+
+
+def test_plot_arrowmap_opm():
+    """Test plot_arrowmap does not crash on colocated OPM channels (gh-13866)."""
+    from mne.viz import plot_arrowmap
+
+    # Need at least 3 unique sensor locations for Delaunay triangulation
+    ch_names = [f"OPM{k:03d}" for k in range(1, 10)]
+    info = create_info(ch_names, 1000.0, ch_types="mag")
+    positions = np.array(
+        [[0.03, 0.0, 0.05]] * 3 + [[-0.03, 0.0, 0.05]] * 3 + [[0.0, 0.03, 0.05]] * 3
+    )
+    with info._unlock():
+        for idx, ch in enumerate(info["chs"]):
+            ch["coil_type"] = FIFF.FIFFV_COIL_FIELDLINE_OPM_MAG_GEN1
+            ch["loc"][:3] = positions[idx]
+    rng = np.random.default_rng(0)
+    data_snap = rng.standard_normal(9)
+    fig = plot_arrowmap(data_snap, info, show=False)
+    assert len(fig.axes) == 1
+
+
 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()