From eb78548a202d1019ddc71c4b88f80ff769379f69 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Thu, 23 Apr 2026 23:17:21 +0530
Subject: [PATCH 01/26] ENH: group triaxial OPM topomaps by orientation

---
 mne/viz/evoked.py             |  78 ++++++++---
 mne/viz/tests/test_ica.py     |   5 +-
 mne/viz/tests/test_topo.py    |   5 +-
 mne/viz/tests/test_topomap.py |  16 +++
 mne/viz/topomap.py            | 245 +++++++++++++++++++++++-----------
 5 files changed, 250 insertions(+), 99 deletions(-)

diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index a62d2379f03..57a657dff5a 100644
--- a/mne/viz/evoked.py
+++ b/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,32 @@ 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
+        from .topomap import _opm_coils
+
+        (
+            _,
+            _,
+            merge_channels,
+            _,
+            _,
+            _,
+            _,
+        ) = _prepare_topomap_plot(
+            evoked,
+            "mag",
+            sphere=topomap_args.get("sphere", None),
+        )
+        is_opm = any(ch["coil_type"] in _opm_coils for ch in evoked.info["chs"])
+        if is_opm and bool(merge_channels):
+            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 +2068,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:
diff --git a/mne/viz/tests/test_ica.py b/mne/viz/tests/test_ica.py
index 46d1145e851..34575e2eddc 100644
--- a/mne/viz/tests/test_ica.py
+++ b/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)
diff --git a/mne/viz/tests/test_topo.py b/mne/viz/tests/test_topo.py
index fe6c938d244..1ffe3188b32 100644
--- a/mne/viz/tests/test_topo.py
+++ b/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():
diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py
index d364d196f18..e1ddef78780 100644
--- a/mne/viz/tests/test_topomap.py
+++ b/mne/viz/tests/test_topomap.py
@@ -851,6 +851,22 @@ def test_split_opm_overlaps(triaxial_evoked):
     assert tangential == ["OPM002", "OPM003", "OPM005", "OPM006"]
 
 
+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()
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index fcbd213ce78..a78ef9ab43b 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -329,6 +329,36 @@ def _split_opm_overlaps(overlapping_channels):
     return radial, tangential
 
 
+def _compute_opm_orientation_topomap_data(data, ch_names, pos, overlapping_channels):
+    """Compute radial and tangential OPM topomap data from overlap sets."""
+    from ..channels.layout import _merge_ch_data
+
+    # Radial data matches the existing OPM merge behavior and position layout.
+    radial_data, radial_names = _merge_ch_data(
+        data.copy(), "mag", copy.copy(ch_names), modality="opm"
+    )
+    radial_pos = pos
+
+    name_lookup = [name.removesuffix("_MERGE-REMOVE") for name in ch_names]
+    tangential_data = []
+    tangential_names = []
+    tangential_pos = []
+    for overlap_set in overlapping_channels:
+        idx = [name_lookup.index(ch_name) for ch_name in overlap_set[1:]]
+        # Collapse multiple tangential channels at one location using RMS.
+        tangential_data.append(np.sqrt(np.mean(data[idx] ** 2, axis=0)))
+        tangential_names.append(f"{overlap_set[0]}t")
+        tangential_pos.append(radial_pos[radial_names.index(overlap_set[0])])
+
+    tangential_data = np.array(tangential_data)
+    tangential_pos = np.array(tangential_pos)
+
+    return [
+        ("radial", radial_data, radial_pos, radial_names),
+        ("tangential", tangential_data, tangential_pos, tangential_names),
+    ]
+
+
 def _plot_update_evoked_topomap(params, bools):
     """Update topomaps."""
     from ..channels.layout import _merge_ch_data
@@ -1729,7 +1759,6 @@ def plot_ica_components(
             clip_origin,
         ) = _prepare_topomap_plot(ica, ch_type, sphere=sphere)
         cmap = _setup_cmap(cmap, n_axes=len(picks))
-        disp_names = _prepare_sensor_names(names, show_names)
         outlines = _make_head_outlines(sphere, pos, outlines, clip_origin)
 
         data = np.dot(
@@ -1738,64 +1767,93 @@ def plot_ica_components(
         data = np.atleast_2d(data)
         data = data[:, data_picks]
 
+        use_opm_orientation_groups = (
+            ch_type == "mag"
+            and any(ch["coil_type"] in _opm_coils for ch in ica.info["chs"])
+            and bool(merge_channels)
+        )
+        n_group_axes = 2 if use_opm_orientation_groups else 1
+
         if title is None:
             title = "ICA components"
         user_passed_axes = _axes is not None
         if not user_passed_axes:
-            fig, _axes, _, _ = _prepare_trellis(len(data), ncols=ncols, nrows=nrows)
+            fig, _axes, _, _ = _prepare_trellis(
+                len(data) * n_group_axes, ncols=ncols, nrows=nrows
+            )
             fig.suptitle(title)
         else:
             _axes = [_axes] if isinstance(_axes, Axes) else _axes
+            if len(_axes) != len(data) * n_group_axes:
+                raise RuntimeError(
+                    "You must provide one axis per component and orientation "
+                    "group for colocated OPM data."
+                )
             fig = _axes[0].get_figure()
 
         subplot_titles = list()
-        for ii, data_, ax in zip(picks, data, _axes):
+        for comp_offset, (ii, data_) in enumerate(zip(picks, data)):
             kwargs = dict(color="gray") if ii in ica.exclude else dict()
             comp_title = ica._ica_names[ii]
             if len(set(ica.get_channel_types())) > 1:
                 comp_title += f" ({ch_type})"
-            subplot_titles.append(ax.set_title(comp_title, fontsize=12, **kwargs))
-            if merge_channels:
-                data_, names_ = _merge_ch_data(data_, ch_type, copy.copy(names))
-            # ↓↓↓ NOTE: we intentionally use the default norm=False here, so that
-            # ↓↓↓ we get vlims that are symmetric-about-zero, even if the data for
-            # ↓↓↓ a given component happens to be one-sided.
-            _vlim = _setup_vmin_vmax(data_, *vlim)
-            im = plot_topomap(
-                data_.flatten(),
-                pos,
-                ch_type=ch_type,
-                sensors=sensors,
-                names=disp_names,
-                contours=contours,
-                outlines=outlines,
-                sphere=sphere,
-                image_interp=image_interp,
-                extrapolate=extrapolate,
-                border=border,
-                res=res,
-                size=size,
-                cmap=cmap[0],
-                vlim=_vlim,
-                cnorm=cnorm,
-                axes=ax,
-                show=False,
-            )[0]
-
-            im.axes.set_label(ica._ica_names[ii])
-            if colorbar:
-                cbar, cax = _add_colorbar(
-                    ax,
-                    im,
-                    cmap,
-                    title="AU",
-                    format_=cbar_fmt,
-                    kind="ica_comp_topomap",
-                    ch_type=ch_type,
+
+            if use_opm_orientation_groups:
+                grouped_data = _compute_opm_orientation_topomap_data(
+                    data_[:, np.newaxis], names, pos, merge_channels
                 )
-                cbar.ax.tick_params(labelsize=12)
-                cbar.set_ticks(_vlim)
-            _hide_frame(ax)
+            else:
+                if merge_channels:
+                    data_, names_ = _merge_ch_data(data_, ch_type, copy.copy(names))
+                    grouped_data = [(None, data_[:, np.newaxis], pos, names_)]
+                else:
+                    grouped_data = [(None, data_[:, np.newaxis], pos, names)]
+
+            for group_idx, (group_label, group_data, group_pos, group_names) in enumerate(
+                grouped_data
+            ):
+                ax_idx = comp_offset * n_group_axes + group_idx
+                ax = _axes[ax_idx]
+                plot_title = comp_title
+                if group_label is not None:
+                    plot_title += f" [{group_label}]"
+                subplot_titles.append(ax.set_title(plot_title, fontsize=12, **kwargs))
+                _vlim = _setup_vmin_vmax(group_data[:, 0], *vlim)
+                im = plot_topomap(
+                    group_data[:, 0].flatten(),
+                    group_pos,
+                    ch_type=ch_type,
+                    sensors=sensors,
+                    names=_prepare_sensor_names(group_names, show_names),
+                    contours=contours,
+                    outlines=outlines,
+                    sphere=sphere,
+                    image_interp=image_interp,
+                    extrapolate=extrapolate,
+                    border=border,
+                    res=res,
+                    size=size,
+                    cmap=cmap[0],
+                    vlim=_vlim,
+                    cnorm=cnorm,
+                    axes=ax,
+                    show=False,
+                )[0]
+
+                im.axes.set_label(ica._ica_names[ii])
+                if colorbar:
+                    cbar, cax = _add_colorbar(
+                        ax,
+                        im,
+                        cmap,
+                        title="AU",
+                        format_=cbar_fmt,
+                        kind="ica_comp_topomap",
+                        ch_type=ch_type,
+                    )
+                    cbar.ax.tick_params(labelsize=12)
+                    cbar.set_ticks(_vlim)
+                _hide_frame(ax)
         del pos
         fig.canvas.draw()
 
@@ -2259,11 +2317,17 @@ def plot_evoked_topomap(
         clip_origin,
     ) = _prepare_topomap_plot(evoked, ch_type, sphere=sphere)
     outlines = _make_head_outlines(sphere, pos, outlines, clip_origin)
+    is_opm = any(ch["coil_type"] in _opm_coils for ch in evoked.info["chs"])
+    use_opm_orientation_groups = ch_type == "mag" and is_opm and bool(merge_channels)
     # check interactive
     axes_given = axes is not None
     interactive = isinstance(times, str) and times == "interactive"
     if interactive and axes_given:
         raise ValueError("User-provided axes not allowed when times='interactive'.")
+    if interactive and use_opm_orientation_groups:
+        raise NotImplementedError(
+            "times='interactive' is not supported for grouped OPM topomaps."
+        )
     # units, scalings
     key = "grad" if ch_type.startswith("planar") else ch_type
     default_scaling = _handle_default("scalings", None)[key]
@@ -2273,7 +2337,6 @@ def plot_evoked_topomap(
     unit = _handle_default("units", units)[key]
     # ch_names (required for NIRS)
     ch_names = names
-    names = _prepare_sensor_names(names, show_names)
     # apply projections before picking. NOTE: the `if proj is True`
     # anti-pattern is needed here to exclude proj='interactive'
     _check_option("proj", proj, (True, False, "interactive", "reconstruct"))
@@ -2299,7 +2362,8 @@ def plot_evoked_topomap(
             f"Times should be between {evoked.times[0]:0.3} and {evoked.times[-1]:0.3}."
         )
     # create axes
-    want_axes = n_times + int(colorbar)
+    n_groups = 2 if use_opm_orientation_groups else 1
+    want_axes = n_times * n_groups + int(colorbar)
     if interactive:
         height_ratios = [5, 1]
         nrows = 2
@@ -2313,7 +2377,7 @@ def plot_evoked_topomap(
             axes.append(plt.subplot(gs[0, ax_idx]))
     elif axes is None:
         fig, axes, ncols, nrows = _prepare_trellis(
-            n_times, ncols=ncols, nrows=nrows, size=size
+            n_times * n_groups, ncols=ncols, nrows=nrows, size=size
         )
     else:
         nrows, ncols = None, None  # Deactivate ncols when axes were passed
@@ -2386,19 +2450,26 @@ def plot_evoked_topomap(
 
     # apply scalings and merge channels
     data *= scaling
+    grouped_data = None
     if merge_channels:
         # check modality
-        if any(ch["coil_type"] in _opm_coils for ch in evoked.info["chs"]):
+        if is_opm:
             modality = "opm"
         elif ch_type in _fnirs_types:
             modality = "fnirs"
         else:
             modality = "other"
-        # merge data
-        data, ch_names = _merge_ch_data(data, ch_type, ch_names, modality=modality)
-        # if ch_type in _fnirs_types:
-        if modality != "other":
+        if modality == "opm" and use_opm_orientation_groups:
+            grouped_data = _compute_opm_orientation_topomap_data(
+                data, ch_names, pos, merge_channels
+            )
             merge_channels = False
+        else:
+            # merge data
+            data, ch_names = _merge_ch_data(data, ch_type, ch_names, modality=modality)
+            # if ch_type in _fnirs_types:
+            if modality != "other":
+                merge_channels = False
     # apply mask if requested
     if mask is not None:
         mask = mask.astype(bool, copy=False)
@@ -2409,8 +2480,14 @@ def plot_evoked_topomap(
         else:  # mag, eeg, planar1, planar2
             mask_ = mask[np.ix_(picks, time_idx)]
     # set up colormap
+    if grouped_data is None:
+        all_data = [data]
+    else:
+        all_data = [group_[1] for group_ in grouped_data]
     _vlim = [
-        _setup_vmin_vmax(data[:, i], *vlim, norm=merge_channels) for i in range(n_times)
+        _setup_vmin_vmax(group_data[:, i], *vlim, norm=merge_channels)
+        for group_data in all_data
+        for i in range(n_times)
     ]
     _vlim = [np.min(_vlim), np.max(_vlim)]
     cmap = _setup_cmap(cmap, n_axes=n_times, norm=_vlim[0] >= 0)
@@ -2427,7 +2504,6 @@ def plot_evoked_topomap(
     kwargs = dict(
         sensors=sensors,
         res=res,
-        names=names,
         cmap=cmap[0],
         cnorm=cnorm,
         mask_params=mask_params,
@@ -2441,33 +2517,42 @@ def plot_evoked_topomap(
         ch_type=ch_type,
     )
     images, contours_ = [], []
-    # loop over times
-    for average_idx, (time, this_average) in enumerate(zip(times, average)):
-        tp, cn, interp = _plot_topomap(
-            data[:, average_idx],
-            pos,
-            axes=axes[average_idx],
-            mask=mask_[:, average_idx] if mask is not None else None,
-            vmin=_vlim[0],
-            vmax=_vlim[1],
-            **kwargs,
-        )
+    if grouped_data is None:
+        grouped_data = [(None, data, pos, ch_names)]
 
-        images.append(tp)
-        if cn is not None:
-            contours_.append(cn)
-        if time_format != "":
-            if this_average is None:
-                axes_title = time_format % (time * scaling_time)
-            else:
-                tmin_ = averaged_times[average_idx][0]
-                tmax_ = averaged_times[average_idx][-1]
-                from_time = time_format % (tmin_ * scaling_time)
-                from_time = from_time.split(" ")[0]  # Remove unit
-                to_time = time_format % (tmax_ * scaling_time)
-                axes_title = f"{from_time} – {to_time}"
-                del from_time, to_time, tmin_, tmax_
-            axes[average_idx].set_title(axes_title)
+    for group_idx, (group_label, group_data, group_pos, group_names) in enumerate(
+        grouped_data
+    ):
+        kwargs["names"] = _prepare_sensor_names(group_names, show_names)
+        for average_idx, (time, this_average) in enumerate(zip(times, average)):
+            ax_idx = group_idx * n_times + average_idx
+            tp, cn, interp = _plot_topomap(
+                group_data[:, average_idx],
+                group_pos,
+                axes=axes[ax_idx],
+                mask=None,
+                vmin=_vlim[0],
+                vmax=_vlim[1],
+                **kwargs,
+            )
+
+            images.append(tp)
+            if cn is not None:
+                contours_.append(cn)
+            if time_format != "":
+                if this_average is None:
+                    axes_title = time_format % (time * scaling_time)
+                else:
+                    tmin_ = averaged_times[average_idx][0]
+                    tmax_ = averaged_times[average_idx][-1]
+                    from_time = time_format % (tmin_ * scaling_time)
+                    from_time = from_time.split(" ")[0]  # Remove unit
+                    to_time = time_format % (tmax_ * scaling_time)
+                    axes_title = f"{from_time} – {to_time}"
+                    del from_time, to_time, tmin_, tmax_
+                if group_label is not None:
+                    axes_title = f"{group_label}\n{axes_title}"
+                axes[ax_idx].set_title(axes_title)
 
     if interactive:
         # Add a slider to the figure and start publishing and subscribing to time_change

From be12b9e981e4960ced2731dc279abcf50188903c Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 24 Apr 2026 14:23:59 +0000
Subject: [PATCH 02/26] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 mne/viz/evoked.py  | 3 +--
 mne/viz/topomap.py | 9 ++++++---
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 57a657dff5a..ca079a6bcfe 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -1956,8 +1956,7 @@ def plot_evoked_joint(
     _, 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
-        from .topomap import _opm_coils
+        from .topomap import _opm_coils, _prepare_topomap_plot
 
         (
             _,
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index a78ef9ab43b..8a626379837 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -1809,9 +1809,12 @@ def plot_ica_components(
                 else:
                     grouped_data = [(None, data_[:, np.newaxis], pos, names)]
 
-            for group_idx, (group_label, group_data, group_pos, group_names) in enumerate(
-                grouped_data
-            ):
+            for group_idx, (
+                group_label,
+                group_data,
+                group_pos,
+                group_names,
+            ) in enumerate(grouped_data):
                 ax_idx = comp_offset * n_group_axes + group_idx
                 ax = _axes[ax_idx]
                 plot_title = comp_title

From d88772f4aa543e066717d5022c5ae816e7910d9b Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 24 Apr 2026 19:57:16 +0530
Subject: [PATCH 03/26] DOC: add changelog for #13866

---
 doc/changes/dev/13866.bugfix.rst | 1 +
 1 file changed, 1 insertion(+)
 create mode 100644 doc/changes/dev/13866.bugfix.rst

diff --git a/doc/changes/dev/13866.bugfix.rst b/doc/changes/dev/13866.bugfix.rst
new file mode 100644
index 00000000000..7c2e617e644
--- /dev/null
+++ b/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`_.

From 676c4c5e693af2d24a94dc4d9ab755237eff51e3 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 24 Apr 2026 20:10:40 +0530
Subject: [PATCH 04/26] FIX: remove dead mask assignment in topomap

---
 mne/viz/topomap.py | 9 ---------
 1 file changed, 9 deletions(-)

diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index 8a626379837..0998dd7f914 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -2473,15 +2473,6 @@ def plot_evoked_topomap(
             # if ch_type in _fnirs_types:
             if modality != "other":
                 merge_channels = False
-    # apply mask if requested
-    if mask is not None:
-        mask = mask.astype(bool, copy=False)
-        if ch_type == "grad":
-            mask_ = (
-                mask[np.ix_(picks[::2], time_idx)] | mask[np.ix_(picks[1::2], time_idx)]
-            )
-        else:  # mag, eeg, planar1, planar2
-            mask_ = mask[np.ix_(picks, time_idx)]
     # set up colormap
     if grouped_data is None:
         all_data = [data]

From 8b9c58d0a82bd71803af76ba647d7edd5794c914 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 24 Apr 2026 22:16:54 +0530
Subject: [PATCH 05/26] FIX: restrict OPM orientation grouping to triaxial
 overlaps

---
 mne/viz/tests/test_topomap.py | 26 ++++++++++++++++++++++++
 mne/viz/topomap.py            | 38 ++++++++++++++++++++++++++---------
 2 files changed, 54 insertions(+), 10 deletions(-)

diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py
index e1ddef78780..4b92d326dbd 100644
--- a/mne/viz/tests/test_topomap.py
+++ b/mne/viz/tests/test_topomap.py
@@ -851,6 +851,32 @@ 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(
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index 0998dd7f914..c29f6affbc6 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -359,6 +359,22 @@ 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).
+    """
+    if ch_type != "mag" or not merge_channels:
+        return False
+
+    pick_chs = [info["chs"][pick] for pick in picks]
+    if not pick_chs or not all(ch["coil_type"] in _opm_coils for ch in pick_chs):
+        return False
+
+    return any(len(overlap_set) >= 3 for overlap_set in merge_channels)
+
+
 def _plot_update_evoked_topomap(params, bools):
     """Update topomaps."""
     from ..channels.layout import _merge_ch_data
@@ -1744,9 +1760,9 @@ def plot_ica_components(
 
     axes = axes.flatten() if isinstance(axes, np.ndarray) else axes
     for k, picks in enumerate(pick_groups):
-        try:  # either an iterable, 1D numpy array or others
-            _axes = axes[k * max_subplots : (k + 1) * max_subplots]
-        except TypeError:  # None or Axes
+        if axes is None:
+            _axes = None
+        else:
             _axes = axes
 
         (
@@ -1767,10 +1783,8 @@ def plot_ica_components(
         data = np.atleast_2d(data)
         data = data[:, data_picks]
 
-        use_opm_orientation_groups = (
-            ch_type == "mag"
-            and any(ch["coil_type"] in _opm_coils for ch in ica.info["chs"])
-            and bool(merge_channels)
+        use_opm_orientation_groups = _should_use_opm_orientation_groups(
+            ica.info, data_picks, merge_channels, ch_type
         )
         n_group_axes = 2 if use_opm_orientation_groups else 1
 
@@ -2320,8 +2334,9 @@ def plot_evoked_topomap(
         clip_origin,
     ) = _prepare_topomap_plot(evoked, ch_type, sphere=sphere)
     outlines = _make_head_outlines(sphere, pos, outlines, clip_origin)
-    is_opm = any(ch["coil_type"] in _opm_coils for ch in evoked.info["chs"])
-    use_opm_orientation_groups = ch_type == "mag" and is_opm and bool(merge_channels)
+    use_opm_orientation_groups = _should_use_opm_orientation_groups(
+        evoked.info, picks, merge_channels, ch_type
+    )
     # check interactive
     axes_given = axes is not None
     interactive = isinstance(times, str) and times == "interactive"
@@ -2456,7 +2471,10 @@ def plot_evoked_topomap(
     grouped_data = None
     if merge_channels:
         # check modality
-        if is_opm:
+        is_opm_picks = len(picks) > 0 and all(
+            evoked.info["chs"][pick]["coil_type"] in _opm_coils for pick in picks
+        )
+        if is_opm_picks:
             modality = "opm"
         elif ch_type in _fnirs_types:
             modality = "fnirs"

From d695e8758a18858c8bf6ad6d1f82251756a37c03 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 24 Apr 2026 23:12:09 +0530
Subject: [PATCH 06/26] FIX: avoid stale pick indexing in OPM modality check

---
 mne/viz/topomap.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index c29f6affbc6..5b31a1e08d4 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -2471,8 +2471,8 @@ def plot_evoked_topomap(
     grouped_data = None
     if merge_channels:
         # check modality
-        is_opm_picks = len(picks) > 0 and all(
-            evoked.info["chs"][pick]["coil_type"] in _opm_coils for pick in picks
+        is_opm_picks = len(evoked.info["chs"]) > 0 and all(
+            ch["coil_type"] in _opm_coils for ch in evoked.info["chs"]
         )
         if is_opm_picks:
             modality = "opm"

From c24a2ea68ccef7480cb68e9639cd45256bb43f8f Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 01:00:47 +0530
Subject: [PATCH 07/26] DOC: add example showing grouped triaxial OPM topomaps

---
 .../plot_evoked_topomap_opm_grouped.py        | 41 +++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100644 examples/visualization/plot_evoked_topomap_opm_grouped.py

diff --git a/examples/visualization/plot_evoked_topomap_opm_grouped.py b/examples/visualization/plot_evoked_topomap_opm_grouped.py
new file mode 100644
index 00000000000..c577bfb80c8
--- /dev/null
+++ b/examples/visualization/plot_evoked_topomap_opm_grouped.py
@@ -0,0 +1,41 @@
+"""
+=====================================
+Plot grouped triaxial OPM topomaps
+=====================================
+
+This example demonstrates grouped radial/tangential topomap rendering for
+colocated triaxial OPM sensors using a small segment of the UCL OPM auditory
+dataset. The grouped rendering places radial maps alongside tangential maps
+so orientation information is explicit.
+
+"""
+# Authors: MNE contributors
+# License: BSD-3-Clause
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import mne
+
+# Load a small segment of the UCL OPM dataset
+subject = "sub-002"
+data_path = mne.datasets.ucl_opm_auditory.data_path()
+opm_file = (
+    data_path / subject / "ses-001" / "meg" / "sub-002_ses-001_task-aef_run-001_meg.bin"
+)
+
+# Read and crop for speed
+raw = mne.io.read_raw_fil(opm_file, verbose="error")
+raw.crop(120, 210).load_data()
+
+# Create epochs and average to get evoked
+events = mne.find_events(raw, min_duration=0.1)
+epochs = mne.Epochs(
+    raw, events, tmin=-0.1, tmax=0.4, baseline=(None, 0), verbose="error"
+)
+evoked = epochs.average()
+
+# Find a peak time and plot grouped topomap
+t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
+fig = evoked.plot_topomap(times=[float(t_peak)], ch_type="mag")
+plt.show()

From 0ae8b4b6cfc98ef42fcf8b3f4322ea51d4249182 Mon Sep 17 00:00:00 2001
From: "autofix-ci[bot]" <114827586+autofix-ci[bot]@users.noreply.github.com>
Date: Tue, 28 Apr 2026 19:36:01 +0000
Subject: [PATCH 08/26] [autofix.ci] apply automated fixes

---
 examples/visualization/plot_evoked_topomap_opm_grouped.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/examples/visualization/plot_evoked_topomap_opm_grouped.py b/examples/visualization/plot_evoked_topomap_opm_grouped.py
index c577bfb80c8..bf02b629002 100644
--- a/examples/visualization/plot_evoked_topomap_opm_grouped.py
+++ b/examples/visualization/plot_evoked_topomap_opm_grouped.py
@@ -11,6 +11,7 @@
 """
 # Authors: MNE contributors
 # License: BSD-3-Clause
+# Copyright the MNE-Python contributors.
 
 import matplotlib.pyplot as plt
 import numpy as np

From c6d31a97504de083808aa1ead27674243d76ee3a Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 01:25:37 +0530
Subject: [PATCH 09/26] FIX: use synthetic data in OPM topomap example to avoid
 dataset download timeout in CI

---
 .../plot_evoked_topomap_opm_grouped.py        | 51 ++++++++++---------
 1 file changed, 28 insertions(+), 23 deletions(-)

diff --git a/examples/visualization/plot_evoked_topomap_opm_grouped.py b/examples/visualization/plot_evoked_topomap_opm_grouped.py
index bf02b629002..c79703706fe 100644
--- a/examples/visualization/plot_evoked_topomap_opm_grouped.py
+++ b/examples/visualization/plot_evoked_topomap_opm_grouped.py
@@ -4,39 +4,44 @@
 =====================================
 
 This example demonstrates grouped radial/tangential topomap rendering for
-colocated triaxial OPM sensors using a small segment of the UCL OPM auditory
-dataset. The grouped rendering places radial maps alongside tangential maps
-so orientation information is explicit.
-
+colocated triaxial OPM sensors. The grouped rendering places radial maps
+alongside tangential maps so orientation information is explicit.
 """
 # Authors: MNE contributors
 # License: BSD-3-Clause
 # Copyright the MNE-Python contributors.
 
-import matplotlib.pyplot as plt
 import numpy as np
 
 import mne
 
-# Load a small segment of the UCL OPM dataset
-subject = "sub-002"
-data_path = mne.datasets.ucl_opm_auditory.data_path()
-opm_file = (
-    data_path / subject / "ses-001" / "meg" / "sub-002_ses-001_task-aef_run-001_meg.bin"
+# Create synthetic triaxial OPM data
+# Simulate three colocated OPM sensors with radial and tangential orientations
+n_channels = 6  # 3 locations × 2 orientations (radial + tangential)
+n_samples = 500
+sfreq = 100.0
+
+# Create synthetic channel info with triaxial OPM layout
+ch_names = [f"G{i // 2}-{['RAD', 'TAN'][i % 2]}" for i in range(n_channels)]
+ch_types = ["mag"] * n_channels
+info = mne.create_info(ch_names, sfreq, ch_types)
+
+# Set channel locations in a line to simulate colocated triplets
+locs = np.array([[0, i // 2 * 0.01, 0] for i in range(n_channels)])
+info.set_montage(
+    mne.channels.make_dig_montage(
+        ch_pos={name: loc for name, loc in zip(ch_names, locs)}
+    )
 )
 
-# Read and crop for speed
-raw = mne.io.read_raw_fil(opm_file, verbose="error")
-raw.crop(120, 210).load_data()
+# Create synthetic data
+data = np.random.randn(n_channels, n_samples) * 1e-12  # Tesla
+data[::2] += (
+    np.sin(2 * np.pi * 10 * np.arange(n_samples) / sfreq) * 1e-12
+)  # 10 Hz signal
 
-# Create epochs and average to get evoked
-events = mne.find_events(raw, min_duration=0.1)
-epochs = mne.Epochs(
-    raw, events, tmin=-0.1, tmax=0.4, baseline=(None, 0), verbose="error"
-)
-evoked = epochs.average()
+# Create evoked object
+evoked = mne.EvokedArray(data, info, tmin=-0.5)
 
-# Find a peak time and plot grouped topomap
-t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
-fig = evoked.plot_topomap(times=[float(t_peak)], ch_type="mag")
-plt.show()
+# Plot grouped topomap showing radial and tangential maps side-by-side
+fig = evoked.plot_topomap(times=[0.0], ch_type="mag")

From d491efc8acb7135774a6fc848a43d05fdb620f17 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 17:27:04 +0530
Subject: [PATCH 10/26] DOC: show grouped OPM topomaps in tutorial

---
 .../plot_evoked_topomap_opm_grouped.py        | 47 -------------------
 tutorials/preprocessing/80_opm_processing.py  |  7 +++
 2 files changed, 7 insertions(+), 47 deletions(-)
 delete mode 100644 examples/visualization/plot_evoked_topomap_opm_grouped.py

diff --git a/examples/visualization/plot_evoked_topomap_opm_grouped.py b/examples/visualization/plot_evoked_topomap_opm_grouped.py
deleted file mode 100644
index c79703706fe..00000000000
--- a/examples/visualization/plot_evoked_topomap_opm_grouped.py
+++ /dev/null
@@ -1,47 +0,0 @@
-"""
-=====================================
-Plot grouped triaxial OPM topomaps
-=====================================
-
-This example demonstrates grouped radial/tangential topomap rendering for
-colocated triaxial OPM sensors. The grouped rendering places radial maps
-alongside tangential maps so orientation information is explicit.
-"""
-# Authors: MNE contributors
-# License: BSD-3-Clause
-# Copyright the MNE-Python contributors.
-
-import numpy as np
-
-import mne
-
-# Create synthetic triaxial OPM data
-# Simulate three colocated OPM sensors with radial and tangential orientations
-n_channels = 6  # 3 locations × 2 orientations (radial + tangential)
-n_samples = 500
-sfreq = 100.0
-
-# Create synthetic channel info with triaxial OPM layout
-ch_names = [f"G{i // 2}-{['RAD', 'TAN'][i % 2]}" for i in range(n_channels)]
-ch_types = ["mag"] * n_channels
-info = mne.create_info(ch_names, sfreq, ch_types)
-
-# Set channel locations in a line to simulate colocated triplets
-locs = np.array([[0, i // 2 * 0.01, 0] for i in range(n_channels)])
-info.set_montage(
-    mne.channels.make_dig_montage(
-        ch_pos={name: loc for name, loc in zip(ch_names, locs)}
-    )
-)
-
-# Create synthetic data
-data = np.random.randn(n_channels, n_samples) * 1e-12  # Tesla
-data[::2] += (
-    np.sin(2 * np.pi * 10 * np.arange(n_samples) / sfreq) * 1e-12
-)  # 10 Hz signal
-
-# Create evoked object
-evoked = mne.EvokedArray(data, info, tmin=-0.5)
-
-# Plot grouped topomap showing radial and tangential maps side-by-side
-fig = evoked.plot_topomap(times=[0.0], ch_type="mag")
diff --git a/tutorials/preprocessing/80_opm_processing.py b/tutorials/preprocessing/80_opm_processing.py
index ba44a797a51..cc7180a1fab 100644
--- a/tutorials/preprocessing/80_opm_processing.py
+++ b/tutorials/preprocessing/80_opm_processing.py
@@ -245,6 +245,13 @@
 t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
 fig = evoked.plot_joint(picks="mag")
 
+# %%
+# Demonstrating grouped OPM topomaps
+# ----------------------------------
+# The UCL OPM auditory dataset has already been loaded above, so this extra
+# figure reuses the same real evoked object without any additional downloads.
+fig = evoked.plot_topomap(times=[t_peak], ch_type="mag")
+
 # %%
 # Visualizing coregistration
 # --------------------------

From da17857d8f93d5a3652e71e2b8b4e594c40098e1 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 18:01:50 +0530
Subject: [PATCH 11/26] FIX: avoid over-allocating axes for OPM joint plots

---
 mne/viz/evoked.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index ca079a6bcfe..5edfed8e756 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -1956,10 +1956,11 @@ def plot_evoked_joint(
     _, times_ts = _check_time_unit(ts_args["time_unit"], times_sec)
 
     if len(ch_types) == 1 and set(ch_types) == {"mag"}:
-        from .topomap import _opm_coils, _prepare_topomap_plot
+        from .topomap import _prepare_topomap_plot, _should_use_opm_orientation_groups
 
+        picks_topomap = None
         (
-            _,
+            picks_topomap,
             _,
             merge_channels,
             _,
@@ -1971,8 +1972,9 @@ def plot_evoked_joint(
             "mag",
             sphere=topomap_args.get("sphere", None),
         )
-        is_opm = any(ch["coil_type"] in _opm_coils for ch in evoked.info["chs"])
-        if is_opm and bool(merge_channels):
+        if _should_use_opm_orientation_groups(
+            evoked.info, picks_topomap, merge_channels, "mag"
+        ):
             opm_group_factor = 2
 
     # prepare axes for topomap

From c4e8b39a3ef217d0d0344af9b55201ac60d04fc0 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 19:41:29 +0530
Subject: [PATCH 12/26] DOC: move grouped OPM demo to kernel_phantom example

---
 examples/datasets/kernel_phantom.py          | 8 ++++++++
 tutorials/preprocessing/80_opm_processing.py | 8 --------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index da17f708454..9f6c7f23ee1 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -51,6 +51,14 @@
 t_peak = 0.016  # based on visual inspection of evoked
 fig.axes[0].axvline(t_peak, color="k", ls=":", lw=3, zorder=2)
 
+# %%
+# Demonstrating grouped OPM topomaps
+# The Kernel OPM phantom data is triaxial, so we can demonstrate the grouped
+# radial and tangential topomap rendering:
+
+fig = evoked.plot_topomap(times=[t_peak], ch_type="mag", merge_channels=True)
+fig = evoked.plot_joint(times=[t_peak], picks="mag")
+
 # %%
 # The data covariance has an interesting structure because of densely packed sensors:
 
diff --git a/tutorials/preprocessing/80_opm_processing.py b/tutorials/preprocessing/80_opm_processing.py
index cc7180a1fab..b650e8db7de 100644
--- a/tutorials/preprocessing/80_opm_processing.py
+++ b/tutorials/preprocessing/80_opm_processing.py
@@ -243,14 +243,6 @@
 )
 evoked = epochs.average()
 t_peak = evoked.times[np.argmax(np.std(evoked.copy().pick("meg").data, axis=0))]
-fig = evoked.plot_joint(picks="mag")
-
-# %%
-# Demonstrating grouped OPM topomaps
-# ----------------------------------
-# The UCL OPM auditory dataset has already been loaded above, so this extra
-# figure reuses the same real evoked object without any additional downloads.
-fig = evoked.plot_topomap(times=[t_peak], ch_type="mag")
 
 # %%
 # Visualizing coregistration

From 953cac8c1af10602c76662728891c4e162f9c909 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 19:58:07 +0530
Subject: [PATCH 13/26] DOC: fix kernel phantom grouped topomap example

---
 examples/datasets/kernel_phantom.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index 9f6c7f23ee1..23d8aaf6ee7 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -56,7 +56,7 @@
 # The Kernel OPM phantom data is triaxial, so we can demonstrate the grouped
 # radial and tangential topomap rendering:
 
-fig = evoked.plot_topomap(times=[t_peak], ch_type="mag", merge_channels=True)
+fig = evoked.plot_topomap(times=[t_peak], ch_type="mag")
 fig = evoked.plot_joint(times=[t_peak], picks="mag")
 
 # %%

From 40446920aeb13896cf4c3bc03b90c27088348b40 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Wed, 29 Apr 2026 20:59:43 +0530
Subject: [PATCH 14/26] DOC: remove non-functional grouped OPM demo from
 kernel_phantom example

---
 examples/datasets/kernel_phantom.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index 23d8aaf6ee7..da17f708454 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -51,14 +51,6 @@
 t_peak = 0.016  # based on visual inspection of evoked
 fig.axes[0].axvline(t_peak, color="k", ls=":", lw=3, zorder=2)
 
-# %%
-# Demonstrating grouped OPM topomaps
-# The Kernel OPM phantom data is triaxial, so we can demonstrate the grouped
-# radial and tangential topomap rendering:
-
-fig = evoked.plot_topomap(times=[t_peak], ch_type="mag")
-fig = evoked.plot_joint(times=[t_peak], picks="mag")
-
 # %%
 # The data covariance has an interesting structure because of densely packed sensors:
 

From 7468d765474df686b5a1e1d3472326614ef6381c Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Thu, 30 Apr 2026 23:25:15 +0530
Subject: [PATCH 15/26] FIX: Add type check for merge_channels in OPM grouping
 gate

---
 examples/datasets/kernel_phantom.py | 11 +++++++++++
 mne/viz/topomap.py                  |  4 ++++
 2 files changed, 15 insertions(+)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index da17f708454..db40b09632f 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -106,3 +106,14 @@
 )
 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. This is only done when merge_channels=True
+# and there are multiple colocated channels at each location:
+
+fig = evoked.plot_topomap(times=[t_peak], merge_channels=True, sphere=sphere)
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index 5b31a1e08d4..bb0f023c6ba 100644
--- a/mne/viz/topomap.py
+++ b/mne/viz/topomap.py
@@ -372,6 +372,10 @@ def _should_use_opm_orientation_groups(info, picks, merge_channels, ch_type):
     if not pick_chs or not all(ch["coil_type"] in _opm_coils for ch in pick_chs):
         return False
 
+    # merge_channels should be a list of overlap sets, not a boolean
+    if not isinstance(merge_channels, (list, tuple)):
+        return False
+
     return any(len(overlap_set) >= 3 for overlap_set in merge_channels)
 
 

From a650b2503dc4fdcf9e8ca9adcb1c068ef4548b55 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 1 May 2026 17:53:10 +0530
Subject: [PATCH 16/26] DOC: fix kernel phantom grouped topomap API usage

---
 examples/datasets/kernel_phantom.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index db40b09632f..1c749a2574d 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -116,4 +116,6 @@
 # components side-by-side. This is only done when merge_channels=True
 # and there are multiple colocated channels at each location:
 
-fig = evoked.plot_topomap(times=[t_peak], merge_channels=True, sphere=sphere)
+fig = evoked.plot_joint(
+    times=[t_peak], topomap_args=dict(merge_channels=True, sphere=sphere)
+)

From 8e81f8ddf4ac4074c0be17a2d58c72bede476661 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 1 May 2026 18:17:41 +0530
Subject: [PATCH 17/26] DOC: remove unsupported merge_channels from kernel
 phantom joint topomap

---
 examples/datasets/kernel_phantom.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index 1c749a2574d..e0e786db560 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -113,9 +113,8 @@
 #
 # 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. This is only done when merge_channels=True
-# and there are multiple colocated channels at each location:
+# components side-by-side when multiple colocated channels are detected:
 
 fig = evoked.plot_joint(
-    times=[t_peak], topomap_args=dict(merge_channels=True, sphere=sphere)
+    times=[t_peak], topomap_args=dict(sphere=sphere)
 )

From 8fd2600277b89e0e29b5dfbd414b18c67015e9c9 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 1 May 2026 12:48:03 +0000
Subject: [PATCH 18/26] [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci
---
 examples/datasets/kernel_phantom.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index e0e786db560..80d43f7fc34 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -115,6 +115,4 @@
 # 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)
-)
+fig = evoked.plot_joint(times=[t_peak], topomap_args=dict(sphere=sphere))

From f23049f69112d26f5d8a2bdbe60d0d6d932970e9 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Fri, 1 May 2026 22:52:08 +0530
Subject: [PATCH 19/26] 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.
---
 examples/datasets/kernel_phantom.py          | 6 +++---
 mne/viz/tests/test_topomap.py                | 7 +++----
 mne/viz/topomap.py                           | 9 +++++----
 tutorials/preprocessing/80_opm_processing.py | 1 +
 4 files changed, 12 insertions(+), 11 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index 80d43f7fc34..f08a2b36755 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/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))
diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py
index 4b92d326dbd..fb402b686a0 100644
--- a/mne/viz/tests/test_topomap.py
+++ b/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"
     )
diff --git a/mne/viz/topomap.py b/mne/viz/topomap.py
index bb0f023c6ba..8d491589588 100644
--- a/mne/viz/topomap.py
+++ b/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):
diff --git a/tutorials/preprocessing/80_opm_processing.py b/tutorials/preprocessing/80_opm_processing.py
index b650e8db7de..ba44a797a51 100644
--- a/tutorials/preprocessing/80_opm_processing.py
+++ b/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

From 37536ee44244602ad2bd4a37ade91501343c5967 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Sat, 2 May 2026 00:09:23 +0530
Subject: [PATCH 20/26] TEST: Update OPM test assertions for grouped biaxial
 rendering

Tests now expect doubled axis counts due to radial+tangential grouping
for biaxial OPM pairs. Kernel phantom data now correctly triggers grouped
visualization with 20 axes for 10 ICA components and 9 axes for 4-point
topomap (4 radial + 4 tangential + 1 colorbar).
---
 mne/viz/tests/test_ica.py     | 3 ++-
 mne/viz/tests/test_topomap.py | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/mne/viz/tests/test_ica.py b/mne/viz/tests/test_ica.py
index 34575e2eddc..c1b50340d17 100644
--- a/mne/viz/tests/test_ica.py
+++ b/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
diff --git a/mne/viz/tests/test_topomap.py b/mne/viz/tests/test_topomap.py
index fb402b686a0..3fbb31ce668 100644
--- a/mne/viz/tests/test_topomap.py
+++ b/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():

From c5b056a486ccf30653190c4c575a6669d980f4ad Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Mon, 4 May 2026 20:23:42 +0530
Subject: [PATCH 21/26] SCOPE: Remove plot_joint OPM grouping to narrow PR
 focus

---
 examples/datasets/kernel_phantom.py | 10 ++--
 mne/viz/evoked.py                   | 77 +++++++----------------------
 2 files changed, 20 insertions(+), 67 deletions(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index f08a2b36755..eaa954fa9ff 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -108,11 +108,7 @@
 mne.viz.set_3d_view(figure=fig, azimuth=30, elevation=70, distance=0.4)
 
 # %%
-# Grouped OPM topomap visualization
-# ==================================
+# For more information on OPM data visualization, see the OPM preprocessing
+# tutorial:
 #
-# 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))
+# - :ref:`ex-opm-processing`
diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 5edfed8e756..1dffe2ddabd 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -1868,7 +1868,6 @@ 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"}
@@ -1955,32 +1954,10 @@ 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) * opm_group_factor, figsize=(8.0, 4.2)
+            len(times_sec), figsize=(8.0, 4.2)
         )
         cbar_ax = None
     else:
@@ -2069,42 +2046,22 @@ def plot_evoked_joint(
 
     # connection lines
     # draw the connection lines between time series and topoplots
-    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)
+    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)
 
     # mark times in time series plot
     for timepoint in times_ts:

From 5edd9471ce00f635ad32d491df96fee95e132e53 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Mon, 4 May 2026 20:35:31 +0530
Subject: [PATCH 22/26] Minor: Format evoked.py line wrapping

---
 mne/viz/evoked.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/mne/viz/evoked.py b/mne/viz/evoked.py
index 1dffe2ddabd..923f80f8eda 100644
--- a/mne/viz/evoked.py
+++ b/mne/viz/evoked.py
@@ -1956,9 +1956,7 @@ def plot_evoked_joint(
 
     # 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), figsize=(8.0, 4.2))
         cbar_ax = None
     else:
         ts_ax = ts_args["axes"]

From 3a354cf20bc55e967c03028947802e66d6bd045a Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Mon, 4 May 2026 21:00:44 +0530
Subject: [PATCH 23/26] DOC: Update changelog entry as per the scope

---
 doc/changes/dev/13866.bugfix.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/changes/dev/13866.bugfix.rst b/doc/changes/dev/13866.bugfix.rst
index 7c2e617e644..9e0eff077d1 100644
--- a/doc/changes/dev/13866.bugfix.rst
+++ b/doc/changes/dev/13866.bugfix.rst
@@ -1 +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`_.
+Completed triaxial OPM topomap grouping by rendering separate radial and tangential maps in evoked topomap and ICA component plotting paths, by `Pragnya Khandelwal`_.

From 3f8bd1d814e1129bb2cd41e825d2da751b31facd Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Mon, 4 May 2026 21:25:06 +0530
Subject: [PATCH 24/26] TEST: Remove plot_joint OPM grouping test (feature
 removed from scope)

---
 mne/viz/tests/test_topo.py | 15 ---------------
 1 file changed, 15 deletions(-)

diff --git a/mne/viz/tests/test_topo.py b/mne/viz/tests/test_topo.py
index 1ffe3188b32..0cb01c03a4f 100644
--- a/mne/viz/tests/test_topo.py
+++ b/mne/viz/tests/test_topo.py
@@ -136,21 +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) >= 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():
     """Test plotting of ERP topography."""
     # Show topography

From e2b784e6c1b10946d5962f654c14958b2a96801c Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Mon, 4 May 2026 22:15:44 +0530
Subject: [PATCH 25/26] DOC: Fix OPM tutorial link in kernel phantom example

---
 examples/datasets/kernel_phantom.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index eaa954fa9ff..c516087e6d5 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/examples/datasets/kernel_phantom.py
@@ -111,4 +111,4 @@
 # For more information on OPM data visualization, see the OPM preprocessing
 # tutorial:
 #
-# - :ref:`ex-opm-processing`
+# - :ref:`tut-opm-processing`

From f97c89fe01ebb23a4eb20a2834fb49e8b3abf337 Mon Sep 17 00:00:00 2001
From: PragnyaKhandelwal <prag1704@gmail.com>
Date: Tue, 5 May 2026 10:03:12 +0530
Subject: [PATCH 26/26] DOC: Add grouped OPM topomap in kernel phantom example

---
 examples/datasets/kernel_phantom.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/examples/datasets/kernel_phantom.py b/examples/datasets/kernel_phantom.py
index c516087e6d5..da2297ba88f 100644
--- a/examples/datasets/kernel_phantom.py
+++ b/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: