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 and ICA component plotting paths, 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

@@ -2059,7 +2059,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

@@ -573,7 +573,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 pairs trigger grouped rendering (radial + tangential axes)
+    assert len(fig.axes) == 20
 
 
 @pytest.mark.slowtest
@@ -585,4 +586,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

@@ -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

@@ -794,7 +794,8 @@ 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 + 1 colorbar)
+    assert len(fig_evoked.axes) == 9
 
 
 def test_prepare_topomap_plot_opm_non_quspin_coils():
@@ -851,6 +852,47 @@ 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 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
+
+    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"]),
+    ]
+
+    # Both biaxial and triaxial overlaps should trigger grouping
+    assert 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()