ENH: group triaxial OPM topomaps by orientation

doc/changes/dev/13866.bugfix.rst

@@ -0,0 +1 @@
+Completed triaxial OPM topomap grouping by rendering separate radial and tangential maps in evoked topomap, joint plot, and ICA component plotting paths, by `Pragnya Khandelwal`_.

examples/datasets/kernel_phantom.py

@@ -106,3 +106,13 @@
 )
 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)
+
+# %%
+# Grouped OPM topomap visualization
+# ==================================
+#
+# Since Kernel OPMs are triaxial sensors (measuring Bx, By, Bz directions),
+# we can visualize them as grouped topomaps showing radial and tangential
+# components side-by-side when multiple colocated channels are detected:
+
+fig = evoked.plot_joint(times=[t_peak], topomap_args=dict(sphere=sphere))

mne/viz/evoked.py

@@ -1868,6 +1868,7 @@ def plot_evoked_joint(
         ts_args.get("time_unit", "s"), evoked.times
     )
     topomap_args = dict() if topomap_args is None else topomap_args.copy()
+    opm_group_factor = 1
 
     got_axes = False
     illegal_args = {"show", "times", "exclude"}
@@ -1954,9 +1955,33 @@ def plot_evoked_joint(
     del times
     _, times_ts = _check_time_unit(ts_args["time_unit"], times_sec)
 
+    if len(ch_types) == 1 and set(ch_types) == {"mag"}:
+        from .topomap import _prepare_topomap_plot, _should_use_opm_orientation_groups
+
+        picks_topomap = None
+        (
+            picks_topomap,
+            _,
+            merge_channels,
+            _,
+            _,
+            _,
+            _,
+        ) = _prepare_topomap_plot(
+            evoked,
+            "mag",
+            sphere=topomap_args.get("sphere", None),
+        )
+        if _should_use_opm_orientation_groups(
+            evoked.info, picks_topomap, merge_channels, "mag"
+        ):
+            opm_group_factor = 2
+
     # prepare axes for topomap
     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) * opm_group_factor, figsize=(8.0, 4.2)
+        )
         cbar_ax = None
     else:
         ts_ax = ts_args["axes"]
@@ -2044,22 +2069,42 @@ def plot_evoked_joint(
 
     # connection lines
     # draw the connection lines between time series and topoplots
-    for timepoint, map_ax_ in zip(times_ts, map_ax):
-        con = ConnectionPatch(
-            xyA=[timepoint, ts_ax.get_ylim()[1]],
-            xyB=[0.5, 0],
-            coordsA="data",
-            coordsB="axes fraction",
-            axesA=ts_ax,
-            axesB=map_ax_,
-            color="grey",
-            linestyle="-",
-            linewidth=1.5,
-            alpha=0.66,
-            zorder=1,
-            clip_on=False,
-        )
-        fig.add_artist(con)
+    if opm_group_factor == 1:
+        for timepoint, map_ax_ in zip(times_ts, map_ax):
+            con = ConnectionPatch(
+                xyA=[timepoint, ts_ax.get_ylim()[1]],
+                xyB=[0.5, 0],
+                coordsA="data",
+                coordsB="axes fraction",
+                axesA=ts_ax,
+                axesB=map_ax_,
+                color="grey",
+                linestyle="-",
+                linewidth=1.5,
+                alpha=0.66,
+                zorder=1,
+                clip_on=False,
+            )
+            ts_ax.add_artist(con)
+    else:
+        for time_idx, timepoint in enumerate(times_ts):
+            for group_idx in range(opm_group_factor):
+                map_ax_ = map_ax[time_idx + group_idx * len(times_ts)]
+                con = ConnectionPatch(
+                    xyA=[timepoint, ts_ax.get_ylim()[1]],
+                    xyB=[0.5, 0],
+                    coordsA="data",
+                    coordsB="axes fraction",
+                    axesA=ts_ax,
+                    axesB=map_ax_,
+                    color="grey",
+                    linestyle="-",
+                    linewidth=1.0,
+                    alpha=0.5,
+                    zorder=1,
+                    clip_on=False,
+                )
+                ts_ax.add_artist(con)
 
     # mark times in time series plot
     for timepoint in times_ts:

mne/viz/tests/test_ica.py

@@ -585,4 +585,7 @@ 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
+    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)

mne/viz/tests/test_topo.py

@@ -145,7 +145,10 @@ def test_plot_joint_opm_triaxial(triaxial_evoked):
         ts_args=dict(time_unit="s"),
         topomap_args=dict(time_unit="s", contours=0, res=8, sensors=False),
     )
-    assert len(fig.axes) >= 2
+    assert len(fig.axes) >= 3
+    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_topo():

mne/viz/tests/test_topomap.py

@@ -851,6 +851,48 @@ def test_split_opm_overlaps(triaxial_evoked):
     assert tangential == ["OPM002", "OPM003", "OPM005", "OPM006"]
 
 
+def test_should_use_opm_orientation_groups_only_for_triaxial():
+    """Test that OPM orientation grouping is restricted to 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
+
+    picks = np.arange(len(ch_names))
+    pair_overlaps = [
+        np.array(["OPM001", "OPM002"]),
+        np.array(["OPM003", "OPM004"]),
+    ]
+    triax_overlaps = [
+        np.array(["OPM001", "OPM002", "OPM003"]),
+        np.array(["OPM004", "OPM005", "OPM006"]),
+    ]
+
+    assert not topomap._should_use_opm_orientation_groups(
+        info, picks, pair_overlaps, "mag"
+    )
+    assert topomap._should_use_opm_orientation_groups(
+        info, picks, 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) == 3
+    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_topomap_nirs_overlap(fnirs_epochs):
     """Test plotting nirs topomap with overlapping channels (gh-7414)."""
     fig = fnirs_epochs["A"].average(picks="hbo").plot_topomap()