FIX: Support biaxial OPM sensor grouping for visualization

Fixed two critical bugs in OPM grouped topomap rendering:

1. Use symmetric distance matrix (not upper triangle) to properly detect all
   colocated channel overlaps bidirectionally.

2. Lower grouping threshold from >=3 to >=2 channels to support both biaxial
   (e.g., bx+by or by+bz) and triaxial (bx+by+bz) OPM sensors.

Real OPM hardware like Kernel phantom has biaxial pairs instead of perfect
triaxial arrays. This fix enables grouped radial/tangential visualization
for real datasets.

Changes:
- mne/viz/topomap.py: Fixed distance matrix and threshold
- mne/viz/tests/test_topomap.py: Updated test expectations
- examples/datasets/kernel_phantom.py: Updated docstring
- tutorials/preprocessing/80_opm_processing.py: Restored plot_joint call

Closes the rendering issue where only 1 plot showed instead of 2.

Author: PragnyaKhandelwal

examples/datasets/kernel_phantom.py

@@ -111,8 +111,8 @@
 # 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:
+# Kernel OPMs measure multiple magnetic field directions (Bx, By, Bz).
+# We can visualize colocated channels as grouped topomaps showing radial and
+# tangential components side-by-side for clearer interpretation of the data:
 
 fig = evoked.plot_joint(times=[t_peak], topomap_args=dict(sphere=sphere))

mne/viz/tests/test_topomap.py

@@ -852,7 +852,7 @@ def test_split_opm_overlaps(triaxial_evoked):
 
 
 def test_should_use_opm_orientation_groups_only_for_triaxial():
-    """Test that OPM orientation grouping is restricted to triaxial overlaps."""
+    """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():
@@ -869,9 +869,8 @@ def test_should_use_opm_orientation_groups_only_for_triaxial():
         np.array(["OPM004", "OPM005", "OPM006"]),
     ]
 
-    assert not topomap._should_use_opm_orientation_groups(
-        info, picks, pair_overlaps, "mag"
-    )
+    # 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"
     )

mne/viz/topomap.py

@@ -226,7 +226,8 @@ def _find_overlaps(info, ch_type, sphere, modality="fnirs"):
     channels_to_exclude = list()
 
     if len(locs3d) > 1 and np.min(dist) < 1e-10:
-        overlapping_mask = np.triu(squareform(dist < 1e-10))
+        # Use symmetric distance matrix to find all colocated channel groups
+        overlapping_mask = squareform(dist < 1e-10)
         for chan_idx in range(overlapping_mask.shape[0]):
             already_overlapped = list(
                 itertools.chain.from_iterable(overlapping_channels)
@@ -362,8 +363,8 @@ def _compute_opm_orientation_topomap_data(data, ch_names, pos, overlapping_chann
 def _should_use_opm_orientation_groups(info, picks, merge_channels, ch_type):
     """Return whether OPM orientation grouping should be enabled.
 
-    Grouping is only used for OPM magnetometer channels with overlap sets that
-    include at least 3 colocated channels (triaxial-style sensors).
+    Grouping is used for OPM magnetometer channels with overlap sets that
+    include at least 2 colocated channels (biaxial or triaxial sensors).
     """
     if ch_type != "mag" or not merge_channels:
         return False
@@ -376,7 +377,7 @@ def _should_use_opm_orientation_groups(info, picks, merge_channels, ch_type):
     if not isinstance(merge_channels, (list, tuple)):
         return False
 
-    return any(len(overlap_set) >= 3 for overlap_set in merge_channels)
+    return any(len(overlap_set) >= 2 for overlap_set in merge_channels)
 
 
 def _plot_update_evoked_topomap(params, bools):

tutorials/preprocessing/80_opm_processing.py

@@ -243,6 +243,7 @@
 )
 evoked = epochs.average()
 t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
+fig = evoked.plot_joint(picks="mag")
 
 # %%
 # Visualizing coregistration