_3d.py

"""Functions to make 3D plots with M/EEG data."""

# Authors: Alexandre Gramfort <alexandre.gramfort@inria.fr>
#          Denis Engemann <denis.engemann@gmail.com>
#          Martin Luessi <mluessi@nmr.mgh.harvard.edu>
#          Eric Larson <larson.eric.d@gmail.com>
#          Mainak Jas <mainak@neuro.hut.fi>
#          Mark Wronkiewicz <wronk.mark@gmail.com>
#
# License: Simplified BSD

import os
import os.path as op
import warnings
from collections.abc import Iterable
from dataclasses import dataclass
from functools import partial
from itertools import cycle
from typing import Optional
from pathlib import Path

import numpy as np

from ._dipole import _check_concat_dipoles, _plot_dipole_mri_outlines, _plot_dipole_3d
from ..defaults import DEFAULTS
from ..fixes import _crop_colorbar, _get_img_fdata
from .._freesurfer import (
    _read_mri_info,
    _check_mri,
    _get_head_surface,
    _get_skull_surface,
)
from ..io import _loc_to_coil_trans
from ..io.pick import (
    pick_types,
    channel_type,
    pick_info,
    _FNIRS_CH_TYPES_SPLIT,
    _MEG_CH_TYPES_SPLIT,
)
from ..io.constants import FIFF
from ..io.meas_info import read_fiducials, create_info
from ..source_space import (
    _ensure_src,
    _create_surf_spacing,
    _check_spacing,
    SourceSpaces,
    read_freesurfer_lut,
)

from ..surface import (
    get_meg_helmet_surf,
    _read_mri_surface,
    _DistanceQuery,
    _project_onto_surface,
    _reorder_ccw,
    _CheckInside,
)
from ..transforms import (
    apply_trans,
    rot_to_quat,
    combine_transforms,
    _get_trans,
    _ensure_trans,
    Transform,
    rotation,
    read_ras_mni_t,
    _print_coord_trans,
    _find_trans,
    transform_surface_to,
    _frame_to_str,
    _get_transforms_to_coord_frame,
)
from ..utils import (
    get_subjects_dir,
    logger,
    _check_subject,
    verbose,
    warn,
    check_version,
    fill_doc,
    _pl,
    get_config,
    _import_nibabel,
    _ensure_int,
    _validate_type,
    _check_option,
    _to_rgb,
)
from ._3d_overlay import _LayeredMesh
from .utils import (
    mne_analyze_colormap,
    _get_color_list,
    _get_cmap,
    plt_show,
    tight_layout,
    figure_nobar,
    _check_time_unit,
)
from ..bem import ConductorModel, _bem_find_surface, _ensure_bem_surfaces


verbose_dec = verbose
FIDUCIAL_ORDER = (FIFF.FIFFV_POINT_LPA, FIFF.FIFFV_POINT_NASION, FIFF.FIFFV_POINT_RPA)


# XXX: to unify with digitization
def _fiducial_coords(points, coord_frame=None):
    """Generate 3x3 array of fiducial coordinates."""
    points = points or []  # None -> list
    if coord_frame is not None:
        points = [p for p in points if p["coord_frame"] == coord_frame]
    points_ = {p["ident"]: p for p in points if p["kind"] == FIFF.FIFFV_POINT_CARDINAL}
    if points_:
        return np.array([points_[i]["r"] for i in FIDUCIAL_ORDER])
    else:
        # XXX eventually this should probably live in montage.py
        if coord_frame is None or coord_frame == FIFF.FIFFV_COORD_HEAD:
            # Try converting CTF HPI coils to fiducials
            out = np.empty((3, 3))
            out.fill(np.nan)
            for p in points:
                if p["kind"] == FIFF.FIFFV_POINT_HPI:
                    if np.isclose(p["r"][1:], 0, atol=1e-6).all():
                        out[0 if p["r"][0] < 0 else 2] = p["r"]
                    elif np.isclose(p["r"][::2], 0, atol=1e-6).all():
                        out[1] = p["r"]
            if np.isfinite(out).all():
                return out
        return np.array([])


@fill_doc
def plot_head_positions(
    pos,
    mode="traces",
    cmap="viridis",
    direction="z",
    show=True,
    destination=None,
    info=None,
    color="k",
    axes=None,
):
    """Plot head positions.

    Parameters
    ----------
    pos : ndarray, shape (n_pos, 10) | list of ndarray
        The head position data. Can also be a list to treat as a
        concatenation of runs.
    mode : str
        Can be 'traces' (default) to show position and quaternion traces,
        or 'field' to show the position as a vector field over time.
    cmap : colormap
        Colormap to use for the trace plot, default is "viridis".
    direction : str
        Can be any combination of "x", "y", or "z" (default: "z") to show
        directional axes in "field" mode.
    show : bool
        Show figure if True. Defaults to True.
    destination : str | array-like, shape (3,) | None
        The destination location for the head, assumed to be in head
        coordinates. See :func:`mne.preprocessing.maxwell_filter` for
        details.

        .. versionadded:: 0.16
    %(info)s If provided, will be used to show the destination position when
        ``destination is None``, and for showing the MEG sensors.

        .. versionadded:: 0.16
    color : color object
        The color to use for lines in ``mode == 'traces'`` and quiver
        arrows in ``mode == 'field'``.

        .. versionadded:: 0.16
    axes : array-like, shape (3, 2)
        The matplotlib axes to use. Only used for ``mode == 'traces'``.

        .. versionadded:: 0.16

    Returns
    -------
    fig : instance of matplotlib.figure.Figure
        The figure.
    """
    from ..chpi import head_pos_to_trans_rot_t
    from ..preprocessing.maxwell import _check_destination
    import matplotlib.pyplot as plt

    _check_option("mode", mode, ["traces", "field"])
    dest_info = dict(dev_head_t=None) if info is None else info
    destination = _check_destination(destination, dest_info, head_frame=True)
    if destination is not None:
        destination = _ensure_trans(destination, "head", "meg")  # probably inv
        destination = destination["trans"][:3].copy()
        destination[:, 3] *= 1000

    if not isinstance(pos, (list, tuple)):
        pos = [pos]
    for ii, p in enumerate(pos):
        p = np.array(p, float)
        if p.ndim != 2 or p.shape[1] != 10:
            raise ValueError(
                "pos (or each entry in pos if a list) must be "
                "dimension (N, 10), got %s" % (p.shape,)
            )
        if ii > 0:  # concatenation
            p[:, 0] += pos[ii - 1][-1, 0] - p[0, 0]
        pos[ii] = p
    borders = np.cumsum([len(pp) for pp in pos])
    pos = np.concatenate(pos, axis=0)
    trans, rot, t = head_pos_to_trans_rot_t(pos)  # also ensures pos is okay
    # trans, rot, and t are for dev_head_t, but what we really want
    # is head_dev_t (i.e., where the head origin is in device coords)
    use_trans = (
        np.einsum("ijk,ik->ij", rot[:, :3, :3].transpose([0, 2, 1]), -trans) * 1000
    )
    use_rot = rot.transpose([0, 2, 1])
    use_quats = -pos[:, 1:4]  # inverse (like doing rot.T)
    surf = rrs = lims = None
    if info is not None:
        meg_picks = pick_types(info, meg=True, ref_meg=False, exclude=())
        if len(meg_picks) > 0:
            rrs = 1000 * np.array(
                [info["chs"][pick]["loc"][:3] for pick in meg_picks], float
            )
            if mode == "traces":
                lims = np.array((rrs.min(0), rrs.max(0))).T
            else:  # mode == 'field'
                surf = get_meg_helmet_surf(info)
                transform_surface_to(surf, "meg", info["dev_head_t"], copy=False)
                surf["rr"] *= 1000.0
    helmet_color = DEFAULTS["coreg"]["helmet_color"]
    if mode == "traces":
        if axes is None:
            axes = plt.subplots(3, 2, sharex=True)[1]
        else:
            axes = np.array(axes)
        if axes.shape != (3, 2):
            raise ValueError("axes must have shape (3, 2), got %s" % (axes.shape,))
        fig = axes[0, 0].figure

        labels = ["xyz", ("$q_1$", "$q_2$", "$q_3$")]
        for ii, (quat, coord) in enumerate(zip(use_quats.T, use_trans.T)):
            axes[ii, 0].plot(t, coord, color, lw=1.0, zorder=3)
            axes[ii, 0].set(ylabel=labels[0][ii], xlim=t[[0, -1]])
            axes[ii, 1].plot(t, quat, color, lw=1.0, zorder=3)
            axes[ii, 1].set(ylabel=labels[1][ii], xlim=t[[0, -1]])
            for b in borders[:-1]:
                for jj in range(2):
                    axes[ii, jj].axvline(t[b], color="r")
        for ii, title in enumerate(("Position (mm)", "Rotation (quat)")):
            axes[0, ii].set(title=title)
            axes[-1, ii].set(xlabel="Time (s)")
        if rrs is not None:
            pos_bads = np.any(
                [
                    (use_trans[:, ii] <= lims[ii, 0])
                    | (use_trans[:, ii] >= lims[ii, 1])
                    for ii in range(3)
                ],
                axis=0,
            )
            for ii in range(3):
                oidx = list(range(ii)) + list(range(ii + 1, 3))
                # knowing it will generally be spherical, we can approximate
                # how far away we are along the axis line by taking the
                # point to the left and right with the smallest distance
                from scipy.spatial.distance import cdist

                dists = cdist(rrs[:, oidx], use_trans[:, oidx])
                left = rrs[:, [ii]] < use_trans[:, ii]
                left_dists_all = dists.copy()
                left_dists_all[~left] = np.inf
                # Don't show negative Z direction
                if ii != 2 and np.isfinite(left_dists_all).any():
                    idx = np.argmin(left_dists_all, axis=0)
                    left_dists = rrs[idx, ii]
                    bads = (
                        ~np.isfinite(left_dists_all[idx, np.arange(len(idx))])
                        | pos_bads
                    )
                    left_dists[bads] = np.nan
                    axes[ii, 0].plot(
                        t, left_dists, color=helmet_color, ls="-", lw=0.5, zorder=2
                    )
                else:
                    axes[ii, 0].axhline(
                        lims[ii][0], color=helmet_color, ls="-", lw=0.5, zorder=2
                    )
                right_dists_all = dists
                right_dists_all[left] = np.inf
                if np.isfinite(right_dists_all).any():
                    idx = np.argmin(right_dists_all, axis=0)
                    right_dists = rrs[idx, ii]
                    bads = (
                        ~np.isfinite(right_dists_all[idx, np.arange(len(idx))])
                        | pos_bads
                    )
                    right_dists[bads] = np.nan
                    axes[ii, 0].plot(
                        t, right_dists, color=helmet_color, ls="-", lw=0.5, zorder=2
                    )
                else:
                    axes[ii, 0].axhline(
                        lims[ii][1], color=helmet_color, ls="-", lw=0.5, zorder=2
                    )

        for ii in range(3):
            axes[ii, 1].set(ylim=[-1, 1])

        if destination is not None:
            vals = np.array(
                [destination[:, 3], rot_to_quat(destination[:, :3])]
            ).T.ravel()
            for ax, val in zip(fig.axes, vals):
                ax.axhline(val, color="r", ls=":", zorder=2, lw=1.0)

    else:  # mode == 'field':
        from matplotlib.colors import Normalize
        from mpl_toolkits.mplot3d.art3d import Line3DCollection
        from mpl_toolkits.mplot3d import Axes3D  # noqa: F401, analysis:ignore

        fig, ax = plt.subplots(1, subplot_kw=dict(projection="3d"))

        # First plot the trajectory as a colormap:
        # http://matplotlib.org/examples/pylab_examples/multicolored_line.html
        pts = use_trans[:, np.newaxis]
        segments = np.concatenate([pts[:-1], pts[1:]], axis=1)
        norm = Normalize(t[0], t[-2])
        lc = Line3DCollection(segments, cmap=cmap, norm=norm)
        lc.set_array(t[:-1])
        ax.add_collection(lc)
        # now plot the head directions as a quiver
        dir_idx = dict(x=0, y=1, z=2)
        kwargs = dict(pivot="tail")
        for d, length in zip(direction, [5.0, 2.5, 1.0]):
            use_dir = use_rot[:, :, dir_idx[d]]
            # draws stems, then heads
            array = np.concatenate((t, np.repeat(t, 2)))
            ax.quiver(
                use_trans[:, 0],
                use_trans[:, 1],
                use_trans[:, 2],
                use_dir[:, 0],
                use_dir[:, 1],
                use_dir[:, 2],
                norm=norm,
                cmap=cmap,
                array=array,
                length=length,
                **kwargs,
            )
            if destination is not None:
                ax.quiver(
                    destination[0, 3],
                    destination[1, 3],
                    destination[2, 3],
                    destination[dir_idx[d], 0],
                    destination[dir_idx[d], 1],
                    destination[dir_idx[d], 2],
                    color=color,
                    length=length,
                    **kwargs,
                )
        mins = use_trans.min(0)
        maxs = use_trans.max(0)
        if surf is not None:
            ax.plot_trisurf(
                *surf["rr"].T,
                triangles=surf["tris"],
                color=helmet_color,
                alpha=0.1,
                shade=False,
            )
            ax.scatter(*rrs.T, s=1, color=helmet_color)
            mins = np.minimum(mins, rrs.min(0))
            maxs = np.maximum(maxs, rrs.max(0))
        scale = (maxs - mins).max() / 2.0
        xlim, ylim, zlim = (maxs + mins)[:, np.newaxis] / 2.0 + [-scale, scale]
        ax.set(xlabel="x", ylabel="y", zlabel="z", xlim=xlim, ylim=ylim, zlim=zlim)
        _set_aspect_equal(ax)
        ax.view_init(30, 45)
    tight_layout(fig=fig)
    plt_show(show)
    return fig


def _set_aspect_equal(ax):
    # XXX recent MPL throws an error for 3D axis aspect setting, not much
    # we can do about it at this point
    try:
        ax.set_aspect("equal")
    except NotImplementedError:
        pass


@verbose
def plot_evoked_field(
    evoked,
    surf_maps,
    time=None,
    time_label="t = %0.0f ms",
    n_jobs=None,
    fig=None,
    vmax=None,
    n_contours=21,
    *,
    interaction="terrain",
    verbose=None,
):
    """Plot MEG/EEG fields on head surface and helmet in 3D.

    Parameters
    ----------
    evoked : instance of mne.Evoked
        The evoked object.
    surf_maps : list
        The surface mapping information obtained with make_field_map.
    time : float | None
        The time point at which the field map shall be displayed. If None,
        the average peak latency (across sensor types) is used.
    time_label : str | None
        How to print info about the time instant visualized.
    %(n_jobs)s
    fig : instance of Figure3D | None
        If None (default), a new figure will be created, otherwise it will
        plot into the given figure.

        .. versionadded:: 0.20
    vmax : float | None
        Maximum intensity. Can be None to use the max(abs(data)).

        .. versionadded:: 0.21
    n_contours : int
        The number of contours.

        .. versionadded:: 0.21
    %(interaction_scene)s
        Defaults to ``'terrain'``.

        .. versionadded:: 1.1
    %(verbose)s

    Returns
    -------
    fig : instance of Figure3D
        The figure.
    """
    # Update the backend
    from .backends.renderer import _get_renderer

    types = [t for t in ["eeg", "grad", "mag"] if t in evoked]
    _validate_type(vmax, (None, "numeric"), "vmax")
    n_contours = _ensure_int(n_contours, "n_contours")
    _check_option("interaction", interaction, ["trackball", "terrain"])

    time_idx = None
    if time is None:
        time = np.mean([evoked.get_peak(ch_type=t)[1] for t in types])
    del types

    if not evoked.times[0] <= time <= evoked.times[-1]:
        raise ValueError("`time` (%0.3f) must be inside `evoked.times`" % time)
    time_idx = np.argmin(np.abs(evoked.times - time))

    # Plot them
    alphas = [1.0, 0.5]
    colors = [(0.6, 0.6, 0.6), (1.0, 1.0, 1.0)]
    colormap = mne_analyze_colormap(format="vtk")
    colormap_lines = np.concatenate(
        [
            np.tile([0.0, 0.0, 255.0, 255.0], (127, 1)),
            np.tile([0.0, 0.0, 0.0, 255.0], (2, 1)),
            np.tile([255.0, 0.0, 0.0, 255.0], (127, 1)),
        ]
    )

    renderer = _get_renderer(fig, bgcolor=(0.0, 0.0, 0.0), size=(600, 600))
    renderer.set_interaction(interaction)

    for ii, this_map in enumerate(surf_maps):
        surf = this_map["surf"]
        map_data = this_map["data"]
        map_type = this_map["kind"]
        map_ch_names = this_map["ch_names"]

        if map_type == "eeg":
            pick = pick_types(evoked.info, meg=False, eeg=True)
        else:
            pick = pick_types(evoked.info, meg=True, eeg=False, ref_meg=False)

        ch_names = [evoked.ch_names[k] for k in pick]

        set_ch_names = set(ch_names)
        set_map_ch_names = set(map_ch_names)
        if set_ch_names != set_map_ch_names:
            message = ["Channels in map and data do not match."]
            diff = set_map_ch_names - set_ch_names
            if len(diff):
                message += ["%s not in data file. " % list(diff)]
            diff = set_ch_names - set_map_ch_names
            if len(diff):
                message += ["%s not in map file." % list(diff)]
            raise RuntimeError(" ".join(message))

        data = np.dot(map_data, evoked.data[pick, time_idx])

        # Make a solid surface
        if vmax is None:
            vmax = np.max(np.abs(data))
        vmax = float(vmax)
        alpha = alphas[ii]
        mesh = _LayeredMesh(
            renderer=renderer,
            vertices=surf["rr"],
            triangles=surf["tris"],
            normals=surf["nn"],
        )
        mesh.map()
        color = _to_rgb(colors[ii], alpha=True)
        cmap = np.array(
            [
                (
                    0,
                    0,
                    0,
                    0,
                ),
                color,
            ]
        )
        ctable = np.round(cmap * 255).astype(np.uint8)
        mesh.add_overlay(
            scalars=np.ones(len(data)),
            colormap=ctable,
            rng=[0, 1],
            opacity=alpha,
            name="surf",
        )
        # Now show our field pattern
        mesh.add_overlay(
            scalars=data,
            colormap=colormap,
            rng=[-vmax, vmax],
            opacity=1.0,
            name="field",
        )

        # And the field lines on top
        if n_contours > 0:
            renderer.contour(
                surface=surf,
                scalars=data,
                contours=n_contours,
                vmin=-vmax,
                vmax=vmax,
                opacity=alpha,
                colormap=colormap_lines,
            )

    if time_label is not None:
        if "%" in time_label:
            time_label %= 1e3 * evoked.times[time_idx]
        renderer.text2d(x_window=0.01, y_window=0.01, text=time_label)
    renderer.set_camera(azimuth=10, elevation=60)
    renderer.show()
    return renderer.scene()


@verbose
def plot_alignment(
    info=None,
    trans=None,
    subject=None,
    subjects_dir=None,
    surfaces="auto",
    coord_frame="auto",
    meg=None,
    eeg="original",
    fwd=None,
    dig=False,
    ecog=True,
    src=None,
    mri_fiducials=False,
    bem=None,
    seeg=True,
    fnirs=True,
    show_axes=False,
    dbs=True,
    fig=None,
    interaction="terrain",
    sensor_colors=None,
    verbose=None,
):
    """Plot head, sensor, and source space alignment in 3D.

    Parameters
    ----------
    %(info)s If None (default), no sensor information will be shown.
    %(trans)s "auto" will load trans from the FreeSurfer directory
        specified by ``subject`` and ``subjects_dir`` parameters.

        .. versionchanged:: 0.19
            Support for 'fsaverage' argument.
    %(subject)s Can be omitted if ``src`` is provided.
    %(subjects_dir)s
    surfaces : str | list | dict
        Surfaces to plot. Supported values:

        * scalp: one of 'head', 'outer_skin' (alias for 'head'),
          'head-dense', or 'seghead' (alias for 'head-dense')
        * skull: 'outer_skull', 'inner_skull', 'brain' (alias for
          'inner_skull')
        * brain: one of 'pial', 'white', 'inflated', or 'brain'
          (alias for 'pial').

        Can be dict to specify alpha values for each surface. Use None
        to specify default value. Specified values must be between 0 and 1.
        for example::

            surfaces=dict(brain=0.4, outer_skull=0.6, head=None)

        Defaults to 'auto', which will look for a head surface and plot
        it if found.

        .. note:: For single layer BEMs it is recommended to use ``'brain'``.
    coord_frame : 'auto' | 'head' | 'meg' | 'mri'
        The coordinate frame to use. If ``'auto'`` (default), chooses ``'mri'``
        if ``trans`` was passed, and ``'head'`` otherwise.

        .. versionchanged:: 1.0
           Defaults to ``'auto'``.
    %(meg)s
    %(eeg)s
    %(fwd)s
    dig : bool | 'fiducials'
        If True, plot the digitization points; 'fiducials' to plot fiducial
        points only.
    %(ecog)s
    src : instance of SourceSpaces | None
        If not None, also plot the source space points.
    mri_fiducials : bool | str | path-like
        Plot MRI fiducials (default False). If ``True``, look for a file with
        the canonical name (``bem/{subject}-fiducials.fif``). If ``str``,
        it can be ``'estimated'`` to use :func:`mne.coreg.get_mni_fiducials`,
        otherwise it should provide the full path to the fiducials file.

        .. versionadded:: 0.22
           Support for ``'estimated'``.
    bem : list of dict | instance of ConductorModel | None
        Can be either the BEM surfaces (list of dict), a BEM solution or a
        sphere model. If None, we first try loading
        ``'$SUBJECTS_DIR/$SUBJECT/bem/$SUBJECT-$SOURCE.fif'``, and then look
        for ``'$SUBJECT*$SOURCE.fif'`` in the same directory. For
        ``'outer_skin'``, the subjects bem and bem/flash folders are searched.
        Defaults to None.
    %(seeg)s
    %(fnirs)s
        .. versionadded:: 0.20
    show_axes : bool
        If True (default False), coordinate frame axis indicators will be
        shown:

        * head in pink.
        * MRI in gray (if ``trans is not None``).
        * MEG in blue (if MEG sensors are present).

        .. versionadded:: 0.16
    %(dbs)s
    fig : Figure3D | None
        PyVista scene in which to plot the alignment.
        If ``None``, creates a new 600x600 pixel figure with black background.

        .. versionadded:: 0.16
    %(interaction_scene)s

        .. versionadded:: 0.16
        .. versionchanged:: 1.0
           Defaults to ``'terrain'``.
    sensor_colors : array-like | None
        Colors to use for the sensor glyphs. Can be list-like of color strings
        (length ``n_sensors``) or array-like of RGB(A) values (shape
        ``(n_sensors, 3)`` or ``(n_sensors, 4)``). ``None`` (the default) uses
        the default sensor colors for the :func:`~mne.viz.plot_alignment` GUI.
    %(verbose)s

    Returns
    -------
    fig : instance of Figure3D
        The figure.

    See Also
    --------
    mne.viz.plot_bem

    Notes
    -----
    This function serves the purpose of checking the validity of the many
    different steps of source reconstruction:

    - Transform matrix (keywords ``trans``, ``meg`` and ``mri_fiducials``),
    - BEM surfaces (keywords ``bem`` and ``surfaces``),
    - sphere conductor model (keywords ``bem`` and ``surfaces``) and
    - source space (keywords ``surfaces`` and ``src``).

    .. versionadded:: 0.15
    """
    # Update the backend
    from .backends.renderer import _get_renderer

    meg, eeg, fnirs, warn_meg = _handle_sensor_types(meg, eeg, fnirs)
    _check_option("interaction", interaction, ["trackball", "terrain"])

    info = create_info(1, 1000.0, "misc") if info is None else info
    _validate_type(info, "info")

    # Handle surfaces:
    if surfaces == "auto" and trans is None:
        surfaces = list()  # if no `trans` can't plot mri surfaces
    if isinstance(surfaces, str):
        surfaces = [surfaces]
    if isinstance(surfaces, dict):
        user_alpha = surfaces.copy()
        for key, val in user_alpha.items():
            _validate_type(key, "str", f"surfaces key {repr(key)}")
            _validate_type(val, (None, "numeric"), f"surfaces[{repr(key)}]")
            if val is not None:
                user_alpha[key] = float(val)
                if not 0 <= user_alpha[key] <= 1:
                    raise ValueError(
                        f"surfaces[{repr(key)}] ({val}) must be" " between 0 and 1"
                    )
    else:
        user_alpha = {}
    surfaces = list(surfaces)
    for si, s in enumerate(surfaces):
        _validate_type(s, "str", f"surfaces[{si}]")

    bem = _ensure_bem_surfaces(bem, extra_allow=(ConductorModel, None))
    assert isinstance(bem, ConductorModel) or bem is None

    _check_option("coord_frame", coord_frame, ["head", "meg", "mri", "auto"])
    if coord_frame == "auto":
        coord_frame = "head" if trans is None else "mri"

    if src is not None:
        src = _ensure_src(src)
        src_subject = src._subject
        subject = src_subject if subject is None else subject
        if src_subject is not None and subject != src_subject:
            raise ValueError(
                f'subject ("{subject}") did not match the '
                f'subject name in src ("{src_subject}")'
            )
    # configure transforms
    if isinstance(trans, str) and trans == "auto":
        subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
        trans = _find_trans(subject, subjects_dir)
    trans, trans_type = _get_trans(trans, fro="head", to="mri")

    picks = pick_types(
        info,
        meg=("sensors" in meg),
        ref_meg=("ref" in meg),
        eeg=(len(eeg) > 0),
        ecog=ecog,
        seeg=seeg,
        dbs=dbs,
        fnirs=(len(fnirs) > 0),
    )
    if trans_type == "identity":
        # Some stuff is natively in head coords, others in MRI coords
        msg = (
            "A head<->mri transformation matrix (trans) is required "
            f"to plot %s in {coord_frame} coordinates, "
            "`trans=None` is not allowed"
        )
        if fwd is not None:
            fwd_frame = _frame_to_str[fwd["coord_frame"]]
            if fwd_frame != coord_frame:
                raise ValueError(msg % f"a {fwd_frame}-coordinate forward solution")
        if src is not None:
            src_frame = _frame_to_str[src[0]["coord_frame"]]
            if src_frame != coord_frame:
                raise ValueError(msg % f"a {src_frame}-coordinate source space")
        if mri_fiducials is not False and coord_frame != "mri":
            raise ValueError(msg % "mri fiducials")
        # only enforce needing `trans` if there are channels in "head"/"device"
        if picks.size and coord_frame == "mri":
            raise ValueError(msg % "sensors")
        # if only plotting sphere model no trans needed
        if bem is not None:
            if not bem["is_sphere"]:
                if coord_frame != "mri":
                    raise ValueError(msg % "a BEM")
            elif surfaces not in (["brain"], []):  # can only plot these
                raise ValueError(msg % (", ".join(surfaces) + " surfaces"))
        elif len(surfaces) > 0 and coord_frame != "mri":
            raise ValueError(msg % (", ".join(surfaces) + " surfaces"))
        trans = Transform("head", "mri")  # not used so just use identity
    # get transforms
    head_mri_t = _get_trans(trans, "head", "mri")[0]
    to_cf_t = _get_transforms_to_coord_frame(info, head_mri_t, coord_frame=coord_frame)

    # Surfaces:
    # both the head and helmet will be in MRI coordinates after this
    surfs = dict()

    # Brain surface:
    brain = sorted(set(surfaces) & set(["brain", "pial", "white", "inflated"]))
    if len(brain) > 1:
        raise ValueError(f"Only one brain surface can be plotted, got {brain}.")
    brain = brain[0] if brain else False
    if brain is not False:
        surfaces.pop(surfaces.index(brain))
        if bem is not None and bem["is_sphere"] and brain == "brain":
            surfs["lh"] = _bem_find_surface(bem, "brain")
        else:
            brain = "pial" if brain == "brain" else brain
            subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
            for hemi in ["lh", "rh"]:
                brain_fname = subjects_dir / subject / "surf" / f"{hemi}.{brain}"
                if not brain_fname.is_file():
                    raise RuntimeError(
                        f"No brain surface found for subject {subject}, "
                        f"expected {brain_fname} to exist"
                    )
                surfs[hemi] = _read_mri_surface(brain_fname)
            subjects_dir = str(subjects_dir)

    # Head surface:
    head_keys = ("auto", "head", "outer_skin", "head-dense", "seghead")
    head = [s for s in surfaces if s in head_keys]
    if len(head) > 1:
        raise ValueError("Can only supply one head-like surface name, " f"got {head}")
    head = head[0] if head else False
    if head is not False:
        surfaces.pop(surfaces.index(head))
    elif "projected" in eeg:
        raise ValueError(
            "A head surface is required to project EEG, "
            '"head", "outer_skin", "head-dense" or "seghead" '
            'must be in surfaces or surfaces must be "auto"'
        )

    # Skull surface:
    skulls = [s for s in surfaces if s in ("outer_skull", "inner_skull")]
    for skull_name in skulls:
        surfaces.pop(surfaces.index(skull_name))
        skull = _get_skull_surface(
            skull_name.split("_")[0], subject, subjects_dir, bem=bem
        )
        skull["name"] = skull_name  # set name for alpha
        surfs[skull_name] = skull

    # we've looked through all of them, raise if some remain
    if len(surfaces) > 0:
        raise ValueError(f"Unknown surface type{_pl(surfaces)}: {surfaces}")

    # set colors and alphas
    defaults = DEFAULTS["coreg"]
    no_deep = not (dbs or seeg) or pick_types(info, dbs=True, seeg=True).size == 0
    max_alpha = 1.0 if no_deep else 0.75
    hemi_val = 0.5
    if src is None or (brain and any(s["type"] == "surf" for s in src)):
        hemi_val = max_alpha
    alpha_range = np.linspace(max_alpha / 2.0, 0, 5)[: len(skulls) + 1]
    if src is None and brain is False and len(skulls) == 0 and not show_axes:
        head_alpha = max_alpha
    else:
        head_alpha = alpha_range[0]
    alphas = dict(lh=hemi_val, rh=hemi_val)
    colors = dict(lh=(0.5,) * 3, rh=(0.5,) * 3)
    for idx, name in enumerate(skulls):
        alphas[name] = alpha_range[idx + 1]
        colors[name] = (0.95 - idx * 0.2, 0.85, 0.95 - idx * 0.2)
    if brain is not False and brain in user_alpha:
        alphas["lh"] = alphas["rh"] = user_alpha.pop(brain)
    # replace default alphas with specified user_alpha
    for k, v in user_alpha.items():
        if v is not None:
            alphas[k] = v
        if k in head_keys and v is not None:
            head_alpha = v
    fid_colors = tuple(defaults[f"{key}_color"] for key in ("lpa", "nasion", "rpa"))

    # initialize figure
    renderer = _get_renderer(
        fig,
        name=f"Sensor alignment: {subject}",
        bgcolor=(0.5, 0.5, 0.5),
        size=(800, 800),
    )
    renderer.set_interaction(interaction)

    # plot head
    _, _, head_surf = _plot_head_surface(
        renderer,
        head,
        subject,
        subjects_dir,
        bem,
        coord_frame,
        to_cf_t,
        alpha=head_alpha,
    )

    # plot helmet
    if "helmet" in meg and pick_types(info, meg=True).size > 0:
        _, _, src_surf = _plot_helmet(renderer, info, to_cf_t, head_mri_t, coord_frame)

    # plot surfaces
    if brain and "lh" not in surfs:  # one layer sphere
        assert bem["coord_frame"] == FIFF.FIFFV_COORD_HEAD
        center = bem["r0"].copy()
        center = apply_trans(to_cf_t["head"], center)
        renderer.sphere(center, scale=0.01, color=colors["lh"], opacity=alphas["lh"])
    if show_axes:
        _plot_axes(renderer, info, to_cf_t, head_mri_t)

    # plot points
    _check_option("dig", dig, (True, False, "fiducials"))
    if dig:
        if dig is True:
            _plot_hpi_coils(renderer, info, to_cf_t)
            _plot_head_shape_points(renderer, info, to_cf_t)
        _plot_head_fiducials(renderer, info, to_cf_t, fid_colors)

    if mri_fiducials:
        _plot_mri_fiducials(
            renderer, mri_fiducials, subjects_dir, subject, to_cf_t, fid_colors
        )

    for key, surf in surfs.items():
        # Surfs can sometimes be in head coords (e.g., if coming from sphere)
        assert isinstance(surf, dict), f"{key}: {type(surf)}"
        surf = transform_surface_to(
            surf, coord_frame, [to_cf_t["mri"], to_cf_t["head"]], copy=True
        )
        renderer.surface(
            surface=surf,
            color=colors[key],
            opacity=alphas[key],
            backface_culling=(key != "helmet"),
        )

    # plot sensors (NB snapshot_brain_montage relies on the last thing being
    # plotted being the sensors, so we need to do this after the surfaces)
    if picks.size > 0:
        _plot_sensors(
            renderer,
            info,
            to_cf_t,
            picks,
            meg,
            eeg,
            fnirs,
            warn_meg,
            head_surf,
            "m",
            sensor_colors=sensor_colors,
        )

    if src is not None:
        atlas_ids, colors = read_freesurfer_lut()
        for ss in src:
            src_rr = ss["rr"][ss["inuse"].astype(bool)]
            src_nn = ss["nn"][ss["inuse"].astype(bool)]

            # update coordinate frame
            src_trans = to_cf_t[_frame_to_str[src[0]["coord_frame"]]]
            src_rr = apply_trans(src_trans, src_rr)
            src_nn = apply_trans(src_trans, src_nn, move=False)

            # volume sources
            if ss["type"] == "vol":
                seg_name = ss.get("seg_name", None)
                if seg_name is not None and seg_name in colors:
                    color = colors[seg_name][:3]
                    color = tuple(i / 256.0 for i in color)
                else:
                    color = (1.0, 1.0, 0.0)

            # surface and discrete sources
            else:
                color = (1.0, 1.0, 0.0)

            if len(src_rr) > 0:
                renderer.quiver3d(
                    x=src_rr[:, 0],
                    y=src_rr[:, 1],
                    z=src_rr[:, 2],
                    u=src_nn[:, 0],
                    v=src_nn[:, 1],
                    w=src_nn[:, 2],
                    color=color,
                    mode="cylinder",
                    scale=3e-3,
                    opacity=0.75,
                    glyph_height=0.25,
                    glyph_center=(0.0, 0.0, 0.0),
                    glyph_resolution=20,
                    backface_culling=True,
                )

    if fwd is not None:
        _plot_forward(renderer, fwd, to_cf_t[_frame_to_str[fwd["coord_frame"]]])

    renderer.set_camera(
        azimuth=90, elevation=90, distance=0.6, focalpoint=(0.0, 0.0, 0.0)
    )
    renderer.show()
    return renderer.scene()


def _handle_sensor_types(meg, eeg, fnirs):
    """Handle plotting inputs for sensors types."""
    if eeg is True:
        eeg = ["original"]
    elif eeg is False:
        eeg = list()

    warn_meg = meg is not None  # only warn if the value is explicit
    if meg is True:
        meg = ["helmet", "sensors", "ref"]
    elif meg is None:
        meg = ["helmet", "sensors"]
    elif meg is False:
        meg = list()

    if fnirs is True:
        fnirs = ["pairs"]
    elif fnirs is False:
        fnirs = list()

    if isinstance(meg, str):
        meg = [meg]
    if isinstance(eeg, str):
        eeg = [eeg]
    if isinstance(fnirs, str):
        fnirs = [fnirs]

    for kind, var in zip(("eeg", "meg", "fnirs"), (eeg, meg, fnirs)):
        if not isinstance(var, (list, tuple)) or not all(
            isinstance(x, str) for x in var
        ):
            raise TypeError(f"{kind} must be list or tuple of str, got {type(kind)}")
    for xi, x in enumerate(meg):
        _check_option(f"meg[{xi}]", x, ("helmet", "sensors", "ref"))
    for xi, x in enumerate(eeg):
        _check_option(f"eeg[{xi}]", x, ("original", "projected"))
    for xi, x in enumerate(fnirs):
        _check_option(f"fnirs[{xi}]", x, ("channels", "pairs", "sources", "detectors"))
    return meg, eeg, fnirs, warn_meg


@verbose
def _ch_pos_in_coord_frame(info, to_cf_t, warn_meg=True, verbose=None):
    """Transform positions from head/device/mri to a coordinate frame."""
    from ..forward import _create_meg_coils

    chs = dict(ch_pos=dict(), sources=dict(), detectors=dict())
    unknown_chs = list()  # prepare for chs with unknown coordinate frame
    type_counts = dict()
    for idx in range(info["nchan"]):
        ch_type = channel_type(info, idx)
        if ch_type in type_counts:
            type_counts[ch_type] += 1
        else:
            type_counts[ch_type] = 1
        type_slices = dict(ch_pos=slice(0, 3))
        if ch_type in _FNIRS_CH_TYPES_SPLIT:
            # add sensors and detectors too for fNIRS
            type_slices.update(sources=slice(3, 6), detectors=slice(6, 9))
        for type_name, type_slice in type_slices.items():
            if ch_type in _MEG_CH_TYPES_SPLIT + ("ref_meg",):
                coil_trans = _loc_to_coil_trans(info["chs"][idx]["loc"])
                # Here we prefer accurate geometry in case we need to
                # ConvexHull the coil, we want true 3D geometry (and not, for
                # example, a straight line / 1D geometry)
                this_coil = [info["chs"][idx]]
                try:
                    coil = _create_meg_coils(this_coil, acc="accurate")[0]
                except RuntimeError:  # we don't have an accurate one
                    coil = _create_meg_coils(this_coil, acc="normal")[0]
                # store verts as ch_coord
                ch_coord, triangles = _sensor_shape(coil)
                ch_coord = apply_trans(coil_trans, ch_coord)
                if len(ch_coord) == 0 and warn_meg:
                    warn(f"MEG sensor {info.ch_names[idx]} not found.")
            else:
                ch_coord = info["chs"][idx]["loc"][type_slice]
            ch_coord_frame = info["chs"][idx]["coord_frame"]
            if ch_coord_frame not in (
                FIFF.FIFFV_COORD_UNKNOWN,
                FIFF.FIFFV_COORD_DEVICE,
                FIFF.FIFFV_COORD_HEAD,
                FIFF.FIFFV_COORD_MRI,
            ):
                raise RuntimeError(
                    f"Channel {info.ch_names[idx]} has coordinate frame "
                    f'{ch_coord_frame}, must be "meg", "head" or "mri".'
                )
            # set unknown as head first
            if ch_coord_frame == FIFF.FIFFV_COORD_UNKNOWN:
                unknown_chs.append(info.ch_names[idx])
                ch_coord_frame = FIFF.FIFFV_COORD_HEAD
            ch_coord = apply_trans(to_cf_t[_frame_to_str[ch_coord_frame]], ch_coord)
            if ch_type in _MEG_CH_TYPES_SPLIT + ("ref_meg",):
                chs[type_name][info.ch_names[idx]] = (ch_coord, triangles)
            else:
                chs[type_name][info.ch_names[idx]] = ch_coord
    if unknown_chs:
        warn(
            f'Got coordinate frame "unknown" for {unknown_chs}, assuming '
            '"head" coordinates.'
        )
    logger.info(
        "Channel types::\t"
        + ", ".join([f"{ch_type}: {count}" for ch_type, count in type_counts.items()])
    )
    return chs["ch_pos"], chs["sources"], chs["detectors"]


def _plot_head_surface(
    renderer, head, subject, subjects_dir, bem, coord_frame, to_cf_t, alpha, color=None
):
    """Render a head surface in a 3D scene."""
    color = DEFAULTS["coreg"]["head_color"] if color is None else color
    actor = None
    src_surf = dst_surf = None
    if head is not False:
        src_surf = _get_head_surface(head, subject, subjects_dir, bem=bem)
        src_surf = transform_surface_to(
            src_surf, coord_frame, [to_cf_t["mri"], to_cf_t["head"]], copy=True
        )
        actor, dst_surf = renderer.surface(
            surface=src_surf, color=color, opacity=alpha, backface_culling=False
        )
    return actor, dst_surf, src_surf


def _plot_helmet(
    renderer, info, to_cf_t, head_mri_t, coord_frame, alpha=0.25, color=None
):
    color = DEFAULTS["coreg"]["helmet_color"] if color is None else color
    src_surf = get_meg_helmet_surf(info, head_mri_t)
    assert src_surf["coord_frame"] == FIFF.FIFFV_COORD_MRI
    src_surf = transform_surface_to(
        src_surf, coord_frame, [to_cf_t["mri"], to_cf_t["head"]], copy=True
    )
    actor, dst_surf = renderer.surface(
        surface=src_surf, color=color, opacity=alpha, backface_culling=False
    )
    return actor, dst_surf, src_surf


def _plot_axes(renderer, info, to_cf_t, head_mri_t):
    """Render different axes a 3D scene."""
    axes = [(to_cf_t["head"], (0.9, 0.3, 0.3))]  # always show head
    if not np.allclose(head_mri_t["trans"], np.eye(4)):  # Show MRI
        axes.append((to_cf_t["mri"], (0.6, 0.6, 0.6)))
    if pick_types(info, meg=True).size > 0:  # Show MEG
        axes.append((to_cf_t["meg"], (0.0, 0.6, 0.6)))
    actors = list()
    for ax in axes:
        x, y, z = np.tile(ax[0]["trans"][:3, 3], 3).reshape((3, 3)).T
        u, v, w = ax[0]["trans"][:3, :3]
        actor, _ = renderer.sphere(
            center=np.column_stack((x[0], y[0], z[0])), color=ax[1], scale=3e-3
        )
        actors.append(actor)
        actor, _ = renderer.quiver3d(
            x=x,
            y=y,
            z=z,
            u=u,
            v=v,
            w=w,
            mode="arrow",
            scale=2e-2,
            color=ax[1],
            scale_mode="scalar",
            resolution=20,
            scalars=[0.33, 0.66, 1.0],
        )
        actors.append(actor)
    return actors


def _plot_head_fiducials(renderer, info, to_cf_t, fid_colors):
    defaults = DEFAULTS["coreg"]
    car_loc = _fiducial_coords(info["dig"], FIFF.FIFFV_COORD_HEAD)
    car_loc = apply_trans(to_cf_t["head"], car_loc)
    if len(car_loc) == 0:
        warn("Digitization points not found. Cannot plot digitization.")
    actors = list()
    for color, data in zip(fid_colors, car_loc):
        actor, _ = renderer.sphere(
            center=data,
            color=color,
            scale=defaults["dig_fid_scale"],
            opacity=defaults["dig_fid_opacity"],
            backface_culling=True,
        )
        actors.append(actor)
    return actors


def _plot_mri_fiducials(
    renderer, mri_fiducials, subjects_dir, subject, to_cf_t, fid_colors
):
    from ..coreg import get_mni_fiducials

    defaults = DEFAULTS["coreg"]
    if mri_fiducials is True:
        subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
        if subject is None:
            raise ValueError(
                "Subject needs to be specified to "
                "automatically find the fiducials file."
            )
        mri_fiducials = subjects_dir / subject / "bem" / (subject + "-fiducials.fif")
    if isinstance(mri_fiducials, str) and mri_fiducials == "estimated":
        mri_fiducials = get_mni_fiducials(subject, subjects_dir)
    elif isinstance(mri_fiducials, (str, Path, os.PathLike)):
        mri_fiducials, cf = read_fiducials(mri_fiducials)
        if cf != FIFF.FIFFV_COORD_MRI:
            raise ValueError("Fiducials are not in MRI space")
    if isinstance(mri_fiducials, np.ndarray):
        fid_loc = mri_fiducials
    else:
        fid_loc = _fiducial_coords(mri_fiducials, FIFF.FIFFV_COORD_MRI)
    fid_loc = apply_trans(to_cf_t["mri"], fid_loc)
    transform = np.eye(4)
    transform[:3, :3] = to_cf_t["mri"]["trans"][:3, :3] * defaults["mri_fid_scale"]
    # rotate around Z axis 45 deg first
    transform = transform @ rotation(0, 0, np.pi / 4)
    actors = list()
    for color, data in zip(fid_colors, fid_loc):
        actor, _ = renderer.quiver3d(
            x=data[0],
            y=data[1],
            z=data[2],
            u=1.0,
            v=0.0,
            w=0.0,
            color=color,
            mode="oct",
            scale=1.0,
            opacity=defaults["mri_fid_opacity"],
            backface_culling=True,
            solid_transform=transform,
        )
        actors.append(actor)
    return actors


def _plot_hpi_coils(
    renderer,
    info,
    to_cf_t,
    opacity=0.5,
    scale=None,
    orient_glyphs=False,
    scale_by_distance=False,
    surf=None,
    check_inside=None,
    nearest=None,
):
    defaults = DEFAULTS["coreg"]
    scale = defaults["hpi_scale"] if scale is None else scale
    hpi_loc = np.array(
        [
            d["r"]
            for d in (info["dig"] or [])
            if (
                d["kind"] == FIFF.FIFFV_POINT_HPI
                and d["coord_frame"] == FIFF.FIFFV_COORD_HEAD
            )
        ]
    )
    hpi_loc = apply_trans(to_cf_t["head"], hpi_loc)
    actor, _ = _plot_glyphs(
        renderer=renderer,
        loc=hpi_loc,
        color=defaults["hpi_color"],
        scale=scale,
        opacity=opacity,
        orient_glyphs=orient_glyphs,
        scale_by_distance=scale_by_distance,
        surf=surf,
        backface_culling=True,
        check_inside=check_inside,
        nearest=nearest,
    )
    return actor


def _get_nearest(nearest, check_inside, project_to_trans, proj_rr):
    idx = nearest.query(proj_rr)[1]
    proj_pts = apply_trans(project_to_trans, nearest.data[idx])
    proj_nn = apply_trans(project_to_trans, check_inside.surf["nn"][idx], move=False)
    return proj_pts, proj_nn


def _orient_glyphs(
    pts,
    surf,
    project_to_surface=False,
    mark_inside=False,
    check_inside=None,
    nearest=None,
):
    if check_inside is None:
        check_inside = _CheckInside(surf, mode="pyvista")
    if nearest is None:
        nearest = _DistanceQuery(surf["rr"])
    project_to_trans = np.eye(4)
    inv_trans = np.linalg.inv(project_to_trans)
    proj_rr = apply_trans(inv_trans, pts)
    proj_pts, proj_nn = _get_nearest(nearest, check_inside, project_to_trans, proj_rr)
    vec = pts - proj_pts  # point to the surface
    nn = proj_nn
    if mark_inside and not project_to_surface:
        scalars = (~check_inside(proj_rr)).astype(int)
    else:
        scalars = np.ones(len(pts))
    dist = np.linalg.norm(vec, axis=-1, keepdims=True)
    vectors = (250 * dist + 1) * nn
    return scalars, vectors, proj_pts


def _plot_glyphs(
    renderer,
    loc,
    color,
    scale,
    opacity=1,
    mode="cylinder",
    orient_glyphs=False,
    scale_by_distance=False,
    project_points=False,
    mark_inside=False,
    surf=None,
    backface_culling=False,
    check_inside=None,
    nearest=None,
):
    if surf is not None and len(loc) > 0:
        defaults = DEFAULTS["coreg"]
        scalars, vectors, proj_pts = _orient_glyphs(
            loc, surf, project_points, mark_inside, check_inside, nearest
        )
        if mark_inside:
            from matplotlib.colors import ListedColormap

            color = np.append(color, 1)
            colormap = ListedColormap(
                np.array(
                    [
                        (
                            0,
                            0,
                            0,
                            1,
                        ),
                        color,
                    ]
                )
            )
            color = None
            clim = [0, 1]
        else:
            scalars = None
            colormap = None
            clim = None
        mode = "cylinder" if orient_glyphs else "sphere"
        scale_mode = "vector" if scale_by_distance else "none"
        x, y, z = proj_pts.T if project_points else loc.T
        u, v, w = vectors.T
        return renderer.quiver3d(
            x,
            y,
            z,
            u,
            v,
            w,
            color=color,
            scale=scale,
            mode=mode,
            glyph_height=defaults["eegp_height"],
            glyph_center=(0.0, -defaults["eegp_height"], 0),
            resolution=16,
            glyph_resolution=16,
            glyph_radius=None,
            opacity=opacity,
            scale_mode=scale_mode,
            scalars=scalars,
            colormap=colormap,
            clim=clim,
        )
    else:
        return renderer.sphere(
            center=loc,
            color=color,
            scale=scale,
            opacity=opacity,
            backface_culling=backface_culling,
        )


@verbose
def _plot_head_shape_points(
    renderer,
    info,
    to_cf_t,
    opacity=0.25,
    orient_glyphs=False,
    scale_by_distance=False,
    mark_inside=False,
    surf=None,
    mask=None,
    check_inside=None,
    nearest=None,
    verbose=False,
):
    defaults = DEFAULTS["coreg"]
    ext_loc = np.array(
        [
            d["r"]
            for d in (info["dig"] or [])
            if (
                d["kind"] == FIFF.FIFFV_POINT_EXTRA
                and d["coord_frame"] == FIFF.FIFFV_COORD_HEAD
            )
        ]
    )
    ext_loc = apply_trans(to_cf_t["head"], ext_loc)
    ext_loc = ext_loc[mask] if mask is not None else ext_loc
    actor, _ = _plot_glyphs(
        renderer=renderer,
        loc=ext_loc,
        color=defaults["extra_color"],
        scale=defaults["extra_scale"],
        opacity=opacity,
        orient_glyphs=orient_glyphs,
        scale_by_distance=scale_by_distance,
        mark_inside=mark_inside,
        surf=surf,
        backface_culling=True,
        check_inside=check_inside,
        nearest=nearest,
    )
    return actor


def _plot_forward(renderer, fwd, fwd_trans, fwd_scale=1, scale=1.5e-3, alpha=1):
    from ..forward import Forward

    _validate_type(fwd, [Forward])
    n_dipoles = fwd["source_rr"].shape[0]
    fwd_rr = fwd["source_rr"]
    if fwd["source_ori"] == FIFF.FIFFV_MNE_FIXED_ORI:
        fwd_nn = fwd["source_nn"].reshape(-1, 1, 3)
    else:
        fwd_nn = fwd["source_nn"].reshape(-1, 3, 3)
    # update coordinate frame
    fwd_rr = apply_trans(fwd_trans, fwd_rr) * fwd_scale
    fwd_nn = apply_trans(fwd_trans, fwd_nn, move=False)
    red = (1.0, 0.0, 0.0)
    green = (0.0, 1.0, 0.0)
    blue = (0.0, 0.0, 1.0)
    actors = list()
    for ori, color in zip(range(fwd_nn.shape[1]), (red, green, blue)):
        actor, _ = renderer.quiver3d(
            *fwd_rr.T,
            *fwd_nn[:, ori].T,
            color=color,
            mode="arrow",
            scale_mode="scalar",
            scalars=np.ones(n_dipoles),
            scale=scale,
            opacity=alpha,
        )
        actors.append(actor)
    return actors


def _plot_sensors(
    renderer,
    info,
    to_cf_t,
    picks,
    meg,
    eeg,
    fnirs,
    warn_meg,
    head_surf,
    units,
    sensor_opacity=0.8,
    orient_glyphs=False,
    scale_by_distance=False,
    project_points=False,
    surf=None,
    check_inside=None,
    nearest=None,
    sensor_colors=None,
):
    """Render sensors in a 3D scene."""
    defaults = DEFAULTS["coreg"]
    ch_pos, sources, detectors = _ch_pos_in_coord_frame(
        pick_info(info, picks), to_cf_t=to_cf_t, warn_meg=warn_meg
    )

    actors = dict(
        meg=list(),
        ref_meg=list(),
        eeg=list(),
        fnirs=list(),
        ecog=list(),
        seeg=list(),
        dbs=list(),
    )
    locs = dict(
        eeg=list(),
        fnirs=list(),
        ecog=list(),
        seeg=list(),
        source=list(),
        detector=list(),
    )
    scalar = 1 if units == "m" else 1e3
    for ch_name, ch_coord in ch_pos.items():
        ch_type = channel_type(info, info.ch_names.index(ch_name))
        # for default picking
        if ch_type in _FNIRS_CH_TYPES_SPLIT:
            ch_type = "fnirs"
        elif ch_type in _MEG_CH_TYPES_SPLIT:
            ch_type = "meg"
        # only plot sensor locations if channels/original in selection
        plot_sensors = (ch_type != "fnirs" or "channels" in fnirs) and (
            ch_type != "eeg" or "original" in eeg
        )
        color = defaults[ch_type + "_color"]
        # plot sensors
        if isinstance(ch_coord, tuple):  # is meg, plot coil
            verts, triangles = ch_coord
            actor, _ = renderer.surface(
                surface=dict(rr=verts * scalar, tris=triangles),
                color=color,
                opacity=0.25,
                backface_culling=True,
            )
            actors[ch_type].append(actor)
        else:
            if plot_sensors:
                locs[ch_type].append(ch_coord)
        if ch_name in sources and "sources" in fnirs:
            locs["source"].append(sources[ch_name])
        if ch_name in detectors and "detectors" in fnirs:
            locs["detector"].append(detectors[ch_name])
        if ch_name in sources and ch_name in detectors and "pairs" in fnirs:
            actor, _ = renderer.tube(  # array of origin and dest points
                origin=sources[ch_name][np.newaxis] * scalar,
                destination=detectors[ch_name][np.newaxis] * scalar,
                radius=0.001 * scalar,
            )
            actors[ch_type].append(actor)

    # add sensors
    for sensor_type in locs.keys():
        if len(locs[sensor_type]) > 0:
            sens_loc = np.array(locs[sensor_type])
            sens_loc = sens_loc[~np.isnan(sens_loc).any(axis=1)]
            scale = defaults[sensor_type + "_scale"]
            if sensor_colors is None:
                color = defaults[sensor_type + "_color"]
                actor, _ = _plot_glyphs(
                    renderer=renderer,
                    loc=sens_loc * scalar,
                    color=color,
                    scale=scale * scalar,
                    opacity=sensor_opacity,
                    orient_glyphs=orient_glyphs,
                    scale_by_distance=scale_by_distance,
                    project_points=project_points,
                    surf=surf,
                    check_inside=check_inside,
                    nearest=nearest,
                )
                if sensor_type in ("source", "detector"):
                    sensor_type = "fnirs"
                actors[sensor_type].append(actor)
            else:
                actor_list = []
                for idx_sen in range(sens_loc.shape[0]):
                    sensor_colors = np.asarray(sensor_colors)
                    if (
                        sensor_colors.ndim not in (1, 2)
                        or sensor_colors.shape[0] != sens_loc.shape[0]
                    ):
                        raise ValueError(
                            "sensor_colors should either be None or be "
                            "array-like with shape (n_sensors,) or "
                            "(n_sensors, 3) or (n_sensors, 4). Got shape "
                            f"{sensor_colors.shape}."
                        )
                    color = sensor_colors[idx_sen]

                    actor, _ = _plot_glyphs(
                        renderer=renderer,
                        loc=(sens_loc * scalar)[idx_sen, :],
                        color=color,
                        scale=scale * scalar,
                        opacity=sensor_opacity,
                        orient_glyphs=orient_glyphs,
                        scale_by_distance=scale_by_distance,
                        project_points=project_points,
                        surf=surf,
                        check_inside=check_inside,
                        nearest=nearest,
                    )
                    actor_list.append(actor)
                if sensor_type in ("source", "detector"):
                    sensor_type = "fnirs"
                actors[sensor_type].append(actor_list)

    # add projected eeg
    eeg_indices = pick_types(info, eeg=True)
    if eeg_indices.size > 0 and "projected" in eeg:
        logger.info("Projecting sensors to the head surface")
        eeg_loc = np.array([ch_pos[info.ch_names[idx]] for idx in eeg_indices])
        eeg_loc = eeg_loc[~np.isnan(eeg_loc).any(axis=1)]
        eegp_loc, eegp_nn = _project_onto_surface(
            eeg_loc, head_surf, project_rrs=True, return_nn=True
        )[2:4]
        del eeg_loc
        eegp_loc *= scalar
        scale = defaults["eegp_scale"] * scalar
        actor, _ = renderer.quiver3d(
            x=eegp_loc[:, 0],
            y=eegp_loc[:, 1],
            z=eegp_loc[:, 2],
            u=eegp_nn[:, 0],
            v=eegp_nn[:, 1],
            w=eegp_nn[:, 2],
            color=defaults["eegp_color"],
            mode="cylinder",
            scale=scale,
            opacity=0.6,
            glyph_height=defaults["eegp_height"],
            glyph_center=(0.0, -defaults["eegp_height"] / 2.0, 0),
            glyph_resolution=20,
            backface_culling=True,
        )
        actors["eeg"].append(actor)

    return actors


def _make_tris_fan(n_vert):
    """Make tris given a number of vertices of a circle-like obj."""
    tris = np.zeros((n_vert - 2, 3), int)
    tris[:, 2] = np.arange(2, n_vert)
    tris[:, 1] = tris[:, 2] - 1
    return tris


def _sensor_shape(coil):
    """Get the sensor shape vertices."""
    from scipy.spatial import ConvexHull, Delaunay

    try:
        from scipy.spatial import QhullError
    except ImportError:  # scipy < 1.8
        from scipy.spatial.qhull import QhullError
    id_ = coil["type"] & 0xFFFF
    pad = True
    # Square figure eight
    if id_ in (
        FIFF.FIFFV_COIL_NM_122,
        FIFF.FIFFV_COIL_VV_PLANAR_W,
        FIFF.FIFFV_COIL_VV_PLANAR_T1,
        FIFF.FIFFV_COIL_VV_PLANAR_T2,
    ):
        # wound by right hand rule such that +x side is "up" (+z)
        long_side = coil["size"]  # length of long side (meters)
        offset = 0.0025  # offset of the center portion of planar grad coil
        rrs = np.array(
            [
                [offset, -long_side / 2.0],
                [long_side / 2.0, -long_side / 2.0],
                [long_side / 2.0, long_side / 2.0],
                [offset, long_side / 2.0],
                [-offset, -long_side / 2.0],
                [-long_side / 2.0, -long_side / 2.0],
                [-long_side / 2.0, long_side / 2.0],
                [-offset, long_side / 2.0],
            ]
        )
        tris = np.concatenate(
            (_make_tris_fan(4), _make_tris_fan(4)[:, ::-1] + 4), axis=0
        )
    # Square
    elif id_ in (
        FIFF.FIFFV_COIL_POINT_MAGNETOMETER,
        FIFF.FIFFV_COIL_VV_MAG_T1,
        FIFF.FIFFV_COIL_VV_MAG_T2,
        FIFF.FIFFV_COIL_VV_MAG_T3,
        FIFF.FIFFV_COIL_KIT_REF_MAG,
    ):
        # square magnetometer (potentially point-type)
        size = 0.001 if id_ == 2000 else (coil["size"] / 2.0)
        rrs = np.array([[-1.0, 1.0], [1.0, 1.0], [1.0, -1.0], [-1.0, -1.0]]) * size
        tris = _make_tris_fan(4)
    # Circle
    elif id_ in (
        FIFF.FIFFV_COIL_MAGNES_MAG,
        FIFF.FIFFV_COIL_MAGNES_REF_MAG,
        FIFF.FIFFV_COIL_CTF_REF_MAG,
        FIFF.FIFFV_COIL_BABY_MAG,
        FIFF.FIFFV_COIL_BABY_REF_MAG,
        FIFF.FIFFV_COIL_ARTEMIS123_REF_MAG,
    ):
        n_pts = 15  # number of points for circle
        circle = np.exp(2j * np.pi * np.arange(n_pts) / float(n_pts))
        circle = np.concatenate(([0.0], circle))
        circle *= coil["size"] / 2.0  # radius of coil
        rrs = np.array([circle.real, circle.imag]).T
        tris = _make_tris_fan(n_pts + 1)
    # Circle
    elif id_ in (
        FIFF.FIFFV_COIL_MAGNES_GRAD,
        FIFF.FIFFV_COIL_CTF_GRAD,
        FIFF.FIFFV_COIL_CTF_REF_GRAD,
        FIFF.FIFFV_COIL_CTF_OFFDIAG_REF_GRAD,
        FIFF.FIFFV_COIL_MAGNES_REF_GRAD,
        FIFF.FIFFV_COIL_MAGNES_OFFDIAG_REF_GRAD,
        FIFF.FIFFV_COIL_KIT_GRAD,
        FIFF.FIFFV_COIL_BABY_GRAD,
        FIFF.FIFFV_COIL_ARTEMIS123_GRAD,
        FIFF.FIFFV_COIL_ARTEMIS123_REF_GRAD,
    ):
        # round coil 1st order (off-diagonal) gradiometer
        baseline = coil["base"] if id_ in (5004, 4005) else 0.0
        n_pts = 16  # number of points for circle
        # This time, go all the way around circle to close it fully
        circle = np.exp(2j * np.pi * np.arange(-1, n_pts) / float(n_pts - 1))
        circle[0] = 0  # center pt for triangulation
        circle *= coil["size"] / 2.0
        rrs = np.array(
            [  # first, second coil
                np.concatenate(
                    [circle.real + baseline / 2.0, circle.real - baseline / 2.0]
                ),
                np.concatenate([circle.imag, -circle.imag]),
            ]
        ).T
        tris = np.concatenate(
            [_make_tris_fan(n_pts + 1), _make_tris_fan(n_pts + 1) + n_pts + 1]
        )
    else:
        # 3D convex hull (will fail for 2D geometry)
        rrs = coil["rmag_orig"].copy()
        try:
            tris = _reorder_ccw(rrs, ConvexHull(rrs).simplices)
        except QhullError:  # 2D geometry likely
            logger.debug("Falling back to planar geometry")
            u, _, _ = np.linalg.svd(rrs.T, full_matrices=False)
            u[:, 2] = 0
            rr_rot = rrs @ u
            tris = Delaunay(rr_rot[:, :2]).simplices
            tris = np.concatenate((tris, tris[:, ::-1]))
        pad = False

    # Go from (x,y) -> (x,y,z)
    if pad:
        rrs = np.pad(rrs, ((0, 0), (0, 1)), mode="constant")
    assert rrs.ndim == 2 and rrs.shape[1] == 3
    return rrs, tris


def _process_clim(clim, colormap, transparent, data=0.0, allow_pos_lims=True):
    """Convert colormap/clim options to dict.

    This fills in any "auto" entries properly such that round-trip
    calling gives the same results.
    """
    # Based on type of limits specified, get cmap control points
    from matplotlib.colors import Colormap

    _validate_type(colormap, (str, Colormap), "colormap")
    data = np.asarray(data)
    if isinstance(colormap, str):
        if colormap == "auto":
            if clim == "auto":
                if allow_pos_lims and (data < 0).any():
                    colormap = "mne"
                else:
                    colormap = "hot"
            else:
                if "lims" in clim:
                    colormap = "hot"
                else:  # 'pos_lims' in clim
                    colormap = "mne"
        colormap = _get_cmap(colormap)
    assert isinstance(colormap, Colormap)
    diverging_maps = [
        "PiYG",
        "PRGn",
        "BrBG",
        "PuOr",
        "RdGy",
        "RdBu",
        "RdYlBu",
        "RdYlGn",
        "Spectral",
        "coolwarm",
        "bwr",
        "seismic",
    ]
    diverging_maps += [d + "_r" for d in diverging_maps]
    diverging_maps += ["mne", "mne_analyze"]
    if clim == "auto":
        # this is merely a heuristic!
        if allow_pos_lims and colormap.name in diverging_maps:
            key = "pos_lims"
        else:
            key = "lims"
        clim = {"kind": "percent", key: [96, 97.5, 99.95]}
    if not isinstance(clim, dict):
        raise ValueError('"clim" must be "auto" or dict, got %s' % (clim,))

    if ("lims" in clim) + ("pos_lims" in clim) != 1:
        raise ValueError(
            "Exactly one of lims and pos_lims must be specified "
            "in clim, got %s" % (clim,)
        )
    if "pos_lims" in clim and not allow_pos_lims:
        raise ValueError('Cannot use "pos_lims" for clim, use "lims" ' "instead")
    diverging = "pos_lims" in clim
    ctrl_pts = np.array(clim["pos_lims" if diverging else "lims"], float)
    ctrl_pts = np.array(ctrl_pts, float)
    if ctrl_pts.shape != (3,):
        raise ValueError("clim has shape %s, it must be (3,)" % (ctrl_pts.shape,))
    if (np.diff(ctrl_pts) < 0).any():
        raise ValueError(
            "colormap limits must be monotonically " "increasing, got %s" % (ctrl_pts,)
        )
    clim_kind = clim.get("kind", "percent")
    _check_option("clim['kind']", clim_kind, ["value", "values", "percent"])
    if clim_kind == "percent":
        perc_data = np.abs(data) if diverging else data
        ctrl_pts = np.percentile(perc_data, ctrl_pts)
        logger.info("Using control points %s" % (ctrl_pts,))
    assert len(ctrl_pts) == 3
    clim = dict(kind="value")
    clim["pos_lims" if diverging else "lims"] = ctrl_pts
    mapdata = dict(clim=clim, colormap=colormap, transparent=transparent)
    return mapdata


def _separate_map(mapdata):
    """Help plotters that cannot handle limit equality."""
    diverging = "pos_lims" in mapdata["clim"]
    key = "pos_lims" if diverging else "lims"
    ctrl_pts = np.array(mapdata["clim"][key])
    assert ctrl_pts.shape == (3,)
    if len(set(ctrl_pts)) == 1:  # three points match
        if ctrl_pts[0] == 0:  # all are zero
            warn("All data were zero")
            ctrl_pts = np.arange(3, dtype=float)
        else:
            ctrl_pts *= [0.0, 0.5, 1]  # all nonzero pts == max
    elif len(set(ctrl_pts)) == 2:  # two points match
        # if points one and two are identical, add a tiny bit to the
        # control point two; if points two and three are identical,
        # subtract a tiny bit from point two.
        bump = 1e-5 if ctrl_pts[0] == ctrl_pts[1] else -1e-5
        ctrl_pts[1] = ctrl_pts[0] + bump * (ctrl_pts[2] - ctrl_pts[0])
    mapdata["clim"][key] = ctrl_pts


def _linearize_map(mapdata):
    from matplotlib.colors import ListedColormap

    diverging = "pos_lims" in mapdata["clim"]
    scale_pts = mapdata["clim"]["pos_lims" if diverging else "lims"]
    if diverging:
        lims = [-scale_pts[2], scale_pts[2]]
        ctrl_norm = (
            np.concatenate(
                [-scale_pts[::-1] / scale_pts[2], [0], scale_pts / scale_pts[2]]
            )
            / 2
            + 0.5
        )
        linear_norm = [0, 0.25, 0.5, 0.5, 0.5, 0.75, 1]
        trans_norm = [1, 1, 0, 0, 0, 1, 1]
    else:
        lims = [scale_pts[0], scale_pts[2]]
        range_ = scale_pts[2] - scale_pts[0]
        mid = (scale_pts[1] - scale_pts[0]) / range_ if range_ > 0 else 0.5
        ctrl_norm = [0, mid, 1]
        linear_norm = [0, 0.5, 1]
        trans_norm = [0, 1, 1]
    # do the piecewise linear transformation
    interp_to = np.linspace(0, 1, 256)
    colormap = np.array(
        mapdata["colormap"](np.interp(interp_to, ctrl_norm, linear_norm))
    )
    if mapdata["transparent"]:
        colormap[:, 3] = np.interp(interp_to, ctrl_norm, trans_norm)
    lims = np.array([lims[0], np.mean(lims), lims[1]])
    colormap = ListedColormap(colormap)
    return colormap, lims


def _get_map_ticks(mapdata):
    diverging = "pos_lims" in mapdata["clim"]
    ticks = mapdata["clim"]["pos_lims" if diverging else "lims"]
    delta = 1e-2 * (ticks[2] - ticks[0])
    if ticks[1] <= ticks[0] + delta:  # Only two worth showing
        ticks = ticks[::2]
    if ticks[1] <= ticks[0] + delta:  # Actually only one
        ticks = ticks[::2]
    if diverging:
        idx = int(ticks[0] == 0)
        ticks = list(-np.array(ticks[idx:])[::-1]) + [0] + list(ticks[idx:])
    return np.array(ticks)


def _handle_time(time_label, time_unit, times):
    """Handle time label string and units."""
    _validate_type(time_label, (None, str, "callable"), "time_label")
    if time_label == "auto":
        if times is not None and len(times) > 1:
            if time_unit == "s":
                time_label = "time=%0.3fs"
            elif time_unit == "ms":
                time_label = "time=%0.1fms"
        else:
            time_label = None
    # convert to callable
    if isinstance(time_label, str):
        time_label_fmt = time_label

        def time_label(x):
            try:
                return time_label_fmt % x
            except Exception:
                return time_label  # in case it's static

    assert time_label is None or callable(time_label)
    if times is not None:
        _, times = _check_time_unit(time_unit, times)
    return time_label, times


def _key_pressed_slider(event, params):
    """Handle key presses for time_viewer slider."""
    step = 1
    if event.key.startswith("ctrl"):
        step = 5
        event.key = event.key.split("+")[-1]
    if event.key not in ["left", "right"]:
        return
    time_viewer = event.canvas.figure
    value = time_viewer.slider.val
    times = params["stc"].times
    if params["time_unit"] == "ms":
        times = times * 1000.0
    time_idx = np.argmin(np.abs(times - value))
    if event.key == "left":
        time_idx = np.max((0, time_idx - step))
    elif event.key == "right":
        time_idx = np.min((len(times) - 1, time_idx + step))
    this_time = times[time_idx]
    time_viewer.slider.set_val(this_time)


def _smooth_plot(this_time, params):
    """Smooth source estimate data and plot with mpl."""
    from ..morph import _hemi_morph

    ax = params["ax"]
    stc = params["stc"]
    ax.clear()
    times = stc.times
    scaler = 1000.0 if params["time_unit"] == "ms" else 1.0
    if this_time is None:
        time_idx = 0
    else:
        time_idx = np.argmin(np.abs(times - this_time / scaler))

    if params["hemi_idx"] == 0:
        data = stc.data[: len(stc.vertices[0]), time_idx : time_idx + 1]
    else:
        data = stc.data[len(stc.vertices[0]) :, time_idx : time_idx + 1]

    morph = _hemi_morph(
        params["tris"],
        params["inuse"],
        params["vertices"],
        params["smoothing_steps"],
        maps=None,
        warn=True,
    )
    array_plot = morph @ data

    range_ = params["scale_pts"][2] - params["scale_pts"][0]
    colors = (array_plot - params["scale_pts"][0]) / range_

    faces = params["faces"]
    greymap = params["greymap"]
    cmap = params["cmap"]
    polyc = ax.plot_trisurf(
        *params["coords"].T, triangles=faces, antialiased=False, vmin=0, vmax=1
    )
    color_ave = np.mean(colors[faces], axis=1).flatten()
    curv_ave = np.mean(params["curv"][faces], axis=1).flatten()
    colors = cmap(color_ave)
    # alpha blend
    colors[:, :3] *= colors[:, [3]]
    colors[:, :3] += greymap(curv_ave)[:, :3] * (1.0 - colors[:, [3]])
    colors[:, 3] = 1.0
    polyc.set_facecolor(colors)
    if params["time_label"] is not None:
        ax.set_title(
            params["time_label"](
                times[time_idx] * scaler,
            ),
            color="w",
        )
    _set_aspect_equal(ax)
    ax.axis("off")
    ax.set(xlim=[-80, 80], ylim=(-80, 80), zlim=[-80, 80])
    ax.figure.canvas.draw()


def _plot_mpl_stc(
    stc,
    subject=None,
    surface="inflated",
    hemi="lh",
    colormap="auto",
    time_label="auto",
    smoothing_steps=10,
    subjects_dir=None,
    views="lat",
    clim="auto",
    figure=None,
    initial_time=None,
    time_unit="s",
    background="black",
    spacing="oct6",
    time_viewer=False,
    colorbar=True,
    transparent=True,
):
    """Plot source estimate using mpl."""
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    from matplotlib.widgets import Slider
    import nibabel as nib
    from scipy import stats
    from ..morph import _get_subject_sphere_tris

    if hemi not in ["lh", "rh"]:
        raise ValueError(
            "hemi must be 'lh' or 'rh' when using matplotlib. " "Got %s." % hemi
        )
    lh_kwargs = {
        "lat": {"elev": 0, "azim": 180},
        "med": {"elev": 0, "azim": 0},
        "ros": {"elev": 0, "azim": 90},
        "cau": {"elev": 0, "azim": -90},
        "dor": {"elev": 90, "azim": -90},
        "ven": {"elev": -90, "azim": -90},
        "fro": {"elev": 0, "azim": 106.739},
        "par": {"elev": 30, "azim": -120},
    }
    rh_kwargs = {
        "lat": {"elev": 0, "azim": 0},
        "med": {"elev": 0, "azim": 180},
        "ros": {"elev": 0, "azim": 90},
        "cau": {"elev": 0, "azim": -90},
        "dor": {"elev": 90, "azim": -90},
        "ven": {"elev": -90, "azim": -90},
        "fro": {"elev": 16.739, "azim": 60},
        "par": {"elev": 30, "azim": -60},
    }
    time_viewer = False if time_viewer == "auto" else time_viewer
    kwargs = dict(lh=lh_kwargs, rh=rh_kwargs)
    views = "lat" if views == "auto" else views
    _check_option("views", views, sorted(lh_kwargs.keys()))
    mapdata = _process_clim(clim, colormap, transparent, stc.data)
    _separate_map(mapdata)
    colormap, scale_pts = _linearize_map(mapdata)
    del transparent, mapdata

    time_label, times = _handle_time(time_label, time_unit, stc.times)
    fig = plt.figure(figsize=(6, 6)) if figure is None else figure
    try:
        ax = Axes3D(fig, auto_add_to_figure=False)
    except Exception:  # old mpl
        ax = Axes3D(fig)
    else:
        fig.add_axes(ax)
    hemi_idx = 0 if hemi == "lh" else 1
    surf = subjects_dir / subject / "surf" / f"{hemi}.{surface}"
    if spacing == "all":
        coords, faces = nib.freesurfer.read_geometry(surf)
        inuse = slice(None)
    else:
        stype, sval, ico_surf, src_type_str = _check_spacing(spacing)
        surf = _create_surf_spacing(surf, hemi, subject, stype, ico_surf, subjects_dir)
        inuse = surf["vertno"]
        faces = surf["use_tris"]
        coords = surf["rr"][inuse]
        shape = faces.shape
        faces = stats.rankdata(faces, "dense").reshape(shape) - 1
        faces = np.round(faces).astype(int)  # should really be int-like anyway
    del surf
    vertices = stc.vertices[hemi_idx]
    n_verts = len(vertices)
    tris = _get_subject_sphere_tris(subject, subjects_dir)[hemi_idx]
    cmap = _get_cmap(colormap)
    greymap = _get_cmap("Greys")

    curv = nib.freesurfer.read_morph_data(
        subjects_dir / subject / "surf" / f"{hemi}.curv"
    )[inuse]
    curv = np.clip(np.array(curv > 0, np.int64), 0.33, 0.66)
    params = dict(
        ax=ax,
        stc=stc,
        coords=coords,
        faces=faces,
        hemi_idx=hemi_idx,
        vertices=vertices,
        tris=tris,
        smoothing_steps=smoothing_steps,
        n_verts=n_verts,
        inuse=inuse,
        cmap=cmap,
        curv=curv,
        scale_pts=scale_pts,
        greymap=greymap,
        time_label=time_label,
        time_unit=time_unit,
    )
    _smooth_plot(initial_time, params)

    ax.view_init(**kwargs[hemi][views])

    try:
        ax.set_facecolor(background)
    except AttributeError:
        ax.set_axis_bgcolor(background)

    if time_viewer:
        time_viewer = figure_nobar(figsize=(4.5, 0.25))
        fig.time_viewer = time_viewer
        ax_time = plt.axes()
        if initial_time is None:
            initial_time = 0
        slider = Slider(
            ax=ax_time,
            label="Time",
            valmin=times[0],
            valmax=times[-1],
            valinit=initial_time,
        )
        time_viewer.slider = slider
        callback_slider = partial(_smooth_plot, params=params)
        slider.on_changed(callback_slider)
        callback_key = partial(_key_pressed_slider, params=params)
        time_viewer.canvas.mpl_connect("key_press_event", callback_key)

        time_viewer.subplots_adjust(left=0.12, bottom=0.05, right=0.75, top=0.95)
    fig.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.0)

    # add colorbar
    from mpl_toolkits.axes_grid1.inset_locator import inset_axes

    sm = plt.cm.ScalarMappable(
        cmap=cmap, norm=plt.Normalize(scale_pts[0], scale_pts[2])
    )
    cax = inset_axes(ax, width="80%", height="5%", loc=8, borderpad=3.0)
    plt.setp(plt.getp(cax, "xticklabels"), color="w")
    sm.set_array(np.linspace(scale_pts[0], scale_pts[2], 256))
    if colorbar:
        cb = plt.colorbar(sm, cax=cax, orientation="horizontal")
        cb_yticks = plt.getp(cax, "yticklabels")
        plt.setp(cb_yticks, color="w")
        cax.tick_params(labelsize=16)
        cb.ax.set_facecolor("0.5")
        cax.set(xlim=(scale_pts[0], scale_pts[2]))
    plt_show(True)
    return fig


def link_brains(brains, time=True, camera=False, colorbar=True, picking=False):
    """Plot multiple SourceEstimate objects with PyVista.

    Parameters
    ----------
    brains : list, tuple or np.ndarray
        The collection of brains to plot.
    time : bool
        If True, link the time controllers. Defaults to True.
    camera : bool
        If True, link the camera controls. Defaults to False.
    colorbar : bool
        If True, link the colorbar controllers. Defaults to True.
    picking : bool
        If True, link the vertices picked with the mouse. Defaults to False.
    """
    from .backends.renderer import _get_3d_backend

    if _get_3d_backend() != "pyvistaqt":
        raise NotImplementedError(
            "Expected 3d backend is pyvistaqt but"
            " {} was given.".format(_get_3d_backend())
        )
    from ._brain import Brain, _LinkViewer

    if not isinstance(brains, Iterable):
        brains = [brains]
    if len(brains) == 0:
        raise ValueError("The collection of brains is empty.")
    for brain in brains:
        if not isinstance(brain, Brain):
            raise TypeError(
                "Expected type is Brain but" " {} was given.".format(type(brain))
            )
        # enable time viewer if necessary
        brain.setup_time_viewer()
    subjects = [brain._subject for brain in brains]
    if subjects.count(subjects[0]) != len(subjects):
        raise RuntimeError("Cannot link brains from different subjects.")

    # link brains properties
    _LinkViewer(
        brains=brains,
        time=time,
        camera=camera,
        colorbar=colorbar,
        picking=picking,
    )


def _check_volume(stc, src, surface, backend_name):
    from ..source_estimate import _BaseSurfaceSourceEstimate, _BaseMixedSourceEstimate

    if isinstance(stc, _BaseSurfaceSourceEstimate):
        return False
    else:
        _validate_type(
            src,
            SourceSpaces,
            "src",
            "src when stc is a mixed or volume source estimate",
        )
        if isinstance(stc, _BaseMixedSourceEstimate):
            # When showing subvolumes, surfaces that preserve geometry must
            # be used (i.e., no inflated)
            _check_option(
                "surface",
                surface,
                ("white", "pial"),
                extra="when plotting a mixed source estimate",
            )
        return True


@verbose
def plot_source_estimates(
    stc,
    subject=None,
    surface="inflated",
    hemi="lh",
    colormap="auto",
    time_label="auto",
    smoothing_steps=10,
    transparent=True,
    alpha=1.0,
    time_viewer="auto",
    subjects_dir=None,
    figure=None,
    views="auto",
    colorbar=True,
    clim="auto",
    cortex="classic",
    size=800,
    background="black",
    foreground=None,
    initial_time=None,
    time_unit="s",
    backend="auto",
    spacing="oct6",
    title=None,
    show_traces="auto",
    src=None,
    volume_options=1.0,
    view_layout="vertical",
    add_data_kwargs=None,
    brain_kwargs=None,
    verbose=None,
):
    """Plot SourceEstimate.

    Parameters
    ----------
    stc : SourceEstimate
        The source estimates to plot.
    %(subject_none)s
        If ``None``, ``stc.subject`` will be used.
    surface : str
        The type of surface (inflated, white etc.).
    hemi : str
        Hemisphere id (ie ``'lh'``, ``'rh'``, ``'both'``, or ``'split'``). In
        the case of ``'both'``, both hemispheres are shown in the same window.
        In the case of ``'split'`` hemispheres are displayed side-by-side
        in different viewing panes.
    %(colormap)s
        The default ('auto') uses ``'hot'`` for one-sided data and
        'mne' for two-sided data.
    %(time_label)s
    smoothing_steps : int
        The amount of smoothing.
    %(transparent)s
    alpha : float
        Alpha value to apply globally to the overlay. Has no effect with mpl
        backend.
    time_viewer : bool | str
        Display time viewer GUI. Can also be 'auto', which will mean True
        for the PyVista backend and False otherwise.

        .. versionchanged:: 0.20.0
           "auto" mode added.
    %(subjects_dir)s
    figure : instance of Figure3D | instance of matplotlib.figure.Figure | list | int | None
        If None, a new figure will be created. If multiple views or a
        split view is requested, this must be a list of the appropriate
        length. If int is provided it will be used to identify the PyVista
        figure by it's id or create a new figure with the given id. If an
        instance of matplotlib figure, mpl backend is used for plotting.
    %(views)s

        When plotting a standard SourceEstimate (not volume, mixed, or vector)
        and using the PyVista backend, ``views='flat'`` is also supported to
        plot cortex as a flatmap.

        Using multiple views (list) is not supported by the matplotlib backend.

        .. versionchanged:: 0.21.0
           Support for flatmaps.
    colorbar : bool
        If True, display colorbar on scene.
    %(clim)s
    cortex : str | tuple
        Specifies how binarized curvature values are rendered.
        Either the name of a preset Brain cortex colorscheme (one of
        ``'classic'``, ``'bone'``, ``'low_contrast'``, or ``'high_contrast'``),
        or the name of a colormap, or a tuple with values
        ``(colormap, min, max, reverse)`` to fully specify the curvature
        colors. Has no effect with the matplotlib backend.
    size : float or tuple of float
        The size of the window, in pixels. can be one number to specify
        a square window, or the (width, height) of a rectangular window.
        Has no effect with mpl backend.
    background : matplotlib color
        Color of the background of the display window.
    foreground : matplotlib color | None
        Color of the foreground of the display window. Has no effect with mpl
        backend. None will choose white or black based on the background color.
    initial_time : float | None
        The time to display on the plot initially. ``None`` to display the
        first time sample (default).
    time_unit : ``'s'`` | ``'ms'``
        Whether time is represented in seconds ("s", default) or
        milliseconds ("ms").
    backend : ``'auto'`` | ``'pyvistaqt'`` | ``'matplotlib'``
        Which backend to use. If ``'auto'`` (default), tries to plot with
        pyvistaqt, but resorts to matplotlib if no 3d backend is available.

        .. versionadded:: 0.15.0
    spacing : str
        Only affects the matplotlib backend.
        The spacing to use for the source space. Can be ``'ico#'`` for a
        recursively subdivided icosahedron, ``'oct#'`` for a recursively
        subdivided octahedron, or ``'all'`` for all points. In general, you can
        speed up the plotting by selecting a sparser source space.
        Defaults  to 'oct6'.

        .. versionadded:: 0.15.0
    title : str | None
        Title for the figure. If None, the subject name will be used.

        .. versionadded:: 0.17.0
    %(show_traces)s
    %(src_volume_options)s
    %(view_layout)s
    %(add_data_kwargs)s
    %(brain_kwargs)s
    %(verbose)s

    Returns
    -------
    figure : instance of mne.viz.Brain | matplotlib.figure.Figure
        An instance of :class:`mne.viz.Brain` or matplotlib figure.

    Notes
    -----
    Flatmaps are available by default for ``fsaverage`` but not for other
    subjects reconstructed by FreeSurfer. We recommend using
    :func:`mne.compute_source_morph` to morph source estimates to ``fsaverage``
    for flatmap plotting. If you want to construct your own flatmap for a given
    subject, these links might help:

    - https://surfer.nmr.mgh.harvard.edu/fswiki/FreeSurferOccipitalFlattenedPatch
    - https://openwetware.org/wiki/Beauchamp:FreeSurfer
    """  # noqa: E501
    from .backends.renderer import _get_3d_backend, use_3d_backend
    from ..source_estimate import _BaseSourceEstimate, _check_stc_src

    _check_stc_src(stc, src)
    _validate_type(stc, _BaseSourceEstimate, "stc", "source estimate")
    subjects_dir = get_subjects_dir(subjects_dir=subjects_dir, raise_error=True)
    subject = _check_subject(stc.subject, subject)
    _check_option("backend", backend, ["auto", "matplotlib", "pyvistaqt", "notebook"])
    plot_mpl = backend == "matplotlib"
    if not plot_mpl:
        if backend == "auto":
            try:
                backend = _get_3d_backend()
            except (ImportError, ModuleNotFoundError):
                warn("No 3D backend found. Resorting to matplotlib 3d.")
                plot_mpl = True
    kwargs = dict(
        subject=subject,
        surface=surface,
        hemi=hemi,
        colormap=colormap,
        time_label=time_label,
        smoothing_steps=smoothing_steps,
        subjects_dir=subjects_dir,
        views=views,
        clim=clim,
        figure=figure,
        initial_time=initial_time,
        time_unit=time_unit,
        background=background,
        time_viewer=time_viewer,
        colorbar=colorbar,
        transparent=transparent,
    )
    if plot_mpl:
        return _plot_mpl_stc(stc, spacing=spacing, **kwargs)
    else:
        with use_3d_backend(backend):
            return _plot_stc(
                stc,
                overlay_alpha=alpha,
                brain_alpha=alpha,
                vector_alpha=alpha,
                cortex=cortex,
                foreground=foreground,
                size=size,
                scale_factor=None,
                show_traces=show_traces,
                src=src,
                volume_options=volume_options,
                view_layout=view_layout,
                add_data_kwargs=add_data_kwargs,
                brain_kwargs=brain_kwargs,
                **kwargs,
            )


def _plot_stc(
    stc,
    subject,
    surface,
    hemi,
    colormap,
    time_label,
    smoothing_steps,
    subjects_dir,
    views,
    clim,
    figure,
    initial_time,
    time_unit,
    background,
    time_viewer,
    colorbar,
    transparent,
    brain_alpha,
    overlay_alpha,
    vector_alpha,
    cortex,
    foreground,
    size,
    scale_factor,
    show_traces,
    src,
    volume_options,
    view_layout,
    add_data_kwargs,
    brain_kwargs,
):
    from .backends.renderer import _get_3d_backend, get_brain_class
    from ..source_estimate import _BaseVolSourceEstimate

    vec = stc._data_ndim == 3
    subjects_dir = str(get_subjects_dir(subjects_dir=subjects_dir, raise_error=True))
    subject = _check_subject(stc.subject, subject)

    backend = _get_3d_backend()
    del _get_3d_backend
    Brain = get_brain_class()
    views = _check_views(surface, views, hemi, stc, backend)
    _check_option("hemi", hemi, ["lh", "rh", "split", "both"])
    _check_option("view_layout", view_layout, ("vertical", "horizontal"))
    time_label, times = _handle_time(time_label, time_unit, stc.times)
    show_traces, time_viewer = _check_st_tv(show_traces, time_viewer, times)

    # convert control points to locations in colormap
    use = stc.magnitude().data if vec else stc.data
    mapdata = _process_clim(clim, colormap, transparent, use, allow_pos_lims=not vec)

    volume = _check_volume(stc, src, surface, backend)

    # XXX we should not need to do this for PyVista, the plotter should be
    # smart enough to do this separation in the cmap-to-ctab conversion
    _separate_map(mapdata)
    colormap = mapdata["colormap"]
    diverging = "pos_lims" in mapdata["clim"]
    scale_pts = mapdata["clim"]["pos_lims" if diverging else "lims"]
    transparent = mapdata["transparent"]
    del mapdata

    if hemi in ["both", "split"]:
        hemis = ["lh", "rh"]
    else:
        hemis = [hemi]

    if overlay_alpha is None:
        overlay_alpha = brain_alpha
    if overlay_alpha == 0:
        smoothing_steps = 1  # Disable smoothing to save time.

    title = subject if len(hemis) > 1 else "%s - %s" % (subject, hemis[0])
    kwargs = {
        "subject": subject,
        "hemi": hemi,
        "surf": surface,
        "title": title,
        "cortex": cortex,
        "size": size,
        "background": background,
        "foreground": foreground,
        "figure": figure,
        "subjects_dir": subjects_dir,
        "views": views,
        "alpha": brain_alpha,
    }
    if brain_kwargs is not None:
        kwargs.update(brain_kwargs)
    kwargs["show"] = False
    kwargs["view_layout"] = view_layout
    with warnings.catch_warnings(record=True):  # traits warnings
        brain = Brain(**kwargs)
    del kwargs

    if scale_factor is None:
        # Configure the glyphs scale directly
        width = np.mean(
            [
                np.ptp(brain.geo[hemi].coords[:, 1])
                for hemi in hemis
                if hemi in brain.geo
            ]
        )
        scale_factor = 0.025 * width / scale_pts[-1]

    if transparent is None:
        transparent = True
    center = 0.0 if diverging else None
    kwargs = {
        "array": stc,
        "colormap": colormap,
        "smoothing_steps": smoothing_steps,
        "time": times,
        "time_label": time_label,
        "alpha": overlay_alpha,
        "colorbar": colorbar,
        "vector_alpha": vector_alpha,
        "scale_factor": scale_factor,
        "initial_time": initial_time,
        "transparent": transparent,
        "center": center,
        "fmin": scale_pts[0],
        "fmid": scale_pts[1],
        "fmax": scale_pts[2],
        "clim": clim,
        "src": src,
        "volume_options": volume_options,
        "verbose": None,
    }
    if add_data_kwargs is not None:
        kwargs.update(add_data_kwargs)
    for hemi in hemis:
        if isinstance(stc, _BaseVolSourceEstimate):  # no surf data
            break
        vertices = stc.vertices[0 if hemi == "lh" else 1]
        if len(vertices) == 0:  # no surf data for the given hemi
            continue  # no data
        use_kwargs = kwargs.copy()
        use_kwargs.update(hemi=hemi)
        with warnings.catch_warnings(record=True):  # traits warnings
            brain.add_data(**use_kwargs)

    if volume:
        use_kwargs = kwargs.copy()
        use_kwargs.update(hemi="vol")
        brain.add_data(**use_kwargs)
    del kwargs

    if time_viewer:
        brain.setup_time_viewer(time_viewer=time_viewer, show_traces=show_traces)
    else:
        brain.show()

    return brain


def _check_st_tv(show_traces, time_viewer, times):
    # time_viewer and show_traces
    _check_option("time_viewer", time_viewer, (True, False, "auto"))
    _validate_type(show_traces, (str, bool, "numeric"), "show_traces")
    if isinstance(show_traces, str):
        _check_option(
            "show_traces",
            show_traces,
            ("auto", "separate", "vertex", "label"),
            extra="when a string",
        )
    if time_viewer == "auto":
        time_viewer = True
    if show_traces == "auto":
        show_traces = time_viewer and times is not None and len(times) > 1
    if show_traces and not time_viewer:
        raise ValueError("show_traces cannot be used when time_viewer=False")
    return show_traces, time_viewer


def _glass_brain_crosshairs(params, x, y, z):
    for ax, a, b in (
        (params["ax_y"], x, z),
        (params["ax_x"], y, z),
        (params["ax_z"], x, y),
    ):
        ax.axvline(a, color="0.75")
        ax.axhline(b, color="0.75")


def _cut_coords_to_ijk(cut_coords, img):
    ijk = apply_trans(np.linalg.inv(img.affine), cut_coords)
    ijk = np.clip(np.round(ijk).astype(int), 0, np.array(img.shape[:3]) - 1)
    return ijk


def _ijk_to_cut_coords(ijk, img):
    return apply_trans(img.affine, ijk)


def _load_subject_mri(mri, stc, subject, subjects_dir, name):
    import nibabel as nib
    from nibabel.spatialimages import SpatialImage

    _validate_type(mri, ("path-like", SpatialImage), name)
    if isinstance(mri, str):
        subject = _check_subject(stc.subject, subject)
        mri = nib.load(_check_mri(mri, subject, subjects_dir))
    return mri


@verbose
def plot_volume_source_estimates(
    stc,
    src,
    subject=None,
    subjects_dir=None,
    mode="stat_map",
    bg_img="T1.mgz",
    colorbar=True,
    colormap="auto",
    clim="auto",
    transparent=None,
    show=True,
    initial_time=None,
    initial_pos=None,
    verbose=None,
):
    """Plot Nutmeg style volumetric source estimates using nilearn.

    Parameters
    ----------
    stc : VectorSourceEstimate
        The vector source estimate to plot.
    src : instance of SourceSpaces | instance of SourceMorph
        The source space. Can also be a SourceMorph to morph the STC to
        a new subject (see Examples).

        .. versionchanged:: 0.18
           Support for :class:`~nibabel.spatialimages.SpatialImage`.
    %(subject_none)s
        If ``None``, ``stc.subject`` will be used.
    %(subjects_dir)s
    mode : str
        The plotting mode to use. Either 'stat_map' (default) or 'glass_brain'.
        For "glass_brain", activation absolute values are displayed
        after being transformed to a standard MNI brain.
    bg_img : instance of SpatialImage | str
        The background image used in the nilearn plotting function.
        Can also be a string to use the ``bg_img`` file in the subject's
        MRI directory (default is ``'T1.mgz'``).
        Not used in "glass brain" plotting.
    colorbar : bool, optional
        If True, display a colorbar on the right of the plots.
    %(colormap)s
    %(clim)s
    %(transparent)s
    show : bool
        Show figures if True. Defaults to True.
    initial_time : float | None
        The initial time to plot. Can be None (default) to use the time point
        with the maximal absolute value activation across all voxels
        or the ``initial_pos`` voxel (if ``initial_pos is None`` or not,
        respectively).

        .. versionadded:: 0.19
    initial_pos : ndarray, shape (3,) | None
        The initial position to use (in m). Can be None (default) to use the
        voxel with the maximum absolute value activation across all time points
        or at ``initial_time`` (if ``initial_time is None`` or not,
        respectively).

        .. versionadded:: 0.19
    %(verbose)s

    Returns
    -------
    fig : instance of Figure
        The figure.

    Notes
    -----
    Click on any of the anatomical slices to explore the time series.
    Clicking on any time point will bring up the corresponding anatomical map.

    The left and right arrow keys can be used to navigate in time.
    To move in time by larger steps, use shift+left and shift+right.

    In ``'glass_brain'`` mode, values are transformed to the standard MNI
    brain using the FreeSurfer Talairach transformation
    ``$SUBJECTS_DIR/$SUBJECT/mri/transforms/talairach.xfm``.

    .. versionadded:: 0.17

    .. versionchanged:: 0.19
       MRI volumes are automatically transformed to MNI space in
       ``'glass_brain'`` mode.

    Examples
    --------
    Passing a :class:`mne.SourceMorph` as the ``src``
    parameter can be useful for plotting in a different subject's space
    (here, a ``'sample'`` STC in ``'fsaverage'``'s space)::

    >>> morph = mne.compute_source_morph(src_sample, subject_to='fsaverage')  # doctest: +SKIP
    >>> fig = stc_vol_sample.plot(morph)  # doctest: +SKIP
    """  # noqa: E501
    from matplotlib import pyplot as plt, colors
    import nibabel as nib
    from ..source_estimate import VolSourceEstimate
    from ..morph import SourceMorph

    if not check_version("nilearn", "0.4"):
        raise RuntimeError("This function requires nilearn >= 0.4")

    from nilearn.plotting import plot_stat_map, plot_glass_brain
    from nilearn.image import index_img

    _check_option("mode", mode, ("stat_map", "glass_brain"))
    plot_func = dict(stat_map=plot_stat_map, glass_brain=plot_glass_brain)[mode]
    _validate_type(stc, VolSourceEstimate, "stc")
    if isinstance(src, SourceMorph):
        img = src.apply(stc, "nifti1", mri_resolution=False, mri_space=False)
        stc = src.apply(stc, mri_resolution=False, mri_space=False)
        kind, src_subject = "morph.subject_to", src.subject_to
    else:
        src = _ensure_src(src, kind="volume", extra=" or SourceMorph")
        img = stc.as_volume(src, mri_resolution=False)
        kind, src_subject = "src subject", src._subject
    del src
    _print_coord_trans(
        Transform("mri_voxel", "ras", img.affine),
        prefix="Image affine ",
        units="mm",
        level="debug",
    )
    subject = _check_subject(src_subject, subject, first_kind=kind)
    stc_ijk = np.array(np.unravel_index(stc.vertices[0], img.shape[:3], order="F")).T
    assert stc_ijk.shape == (len(stc.vertices[0]), 3)
    del kind

    # XXX this assumes zooms are uniform, should probably mult by zooms...
    dist_to_verts = _DistanceQuery(stc_ijk, allow_kdtree=True)

    def _cut_coords_to_idx(cut_coords, img):
        """Convert voxel coordinates to index in stc.data."""
        ijk = _cut_coords_to_ijk(cut_coords, img)
        del cut_coords
        logger.debug("    Affine remapped cut coords to [%d, %d, %d] idx" % tuple(ijk))
        dist, loc_idx = dist_to_verts.query(ijk[np.newaxis])
        dist, loc_idx = dist[0], loc_idx[0]
        logger.debug(
            "    Using vertex %d at a distance of %d voxels"
            % (stc.vertices[0][loc_idx], dist)
        )
        return loc_idx

    ax_name = dict(x="X (sagittal)", y="Y (coronal)", z="Z (axial)")

    def _click_to_cut_coords(event, params):
        """Get voxel coordinates from mouse click."""
        if event.inaxes is params["ax_x"]:
            ax = "x"
            x = params["ax_z"].lines[0].get_xdata()[0]
            y, z = event.xdata, event.ydata
        elif event.inaxes is params["ax_y"]:
            ax = "y"
            y = params["ax_x"].lines[0].get_xdata()[0]
            x, z = event.xdata, event.ydata
        elif event.inaxes is params["ax_z"]:
            ax = "z"
            x, y = event.xdata, event.ydata
            z = params["ax_x"].lines[1].get_ydata()[0]
        else:
            logger.debug("    Click outside axes")
            return None
        cut_coords = np.array((x, y, z))
        logger.debug("")

        if params["mode"] == "glass_brain":  # find idx for MIP
            # Figure out what XYZ in world coordinates is in our voxel data
            codes = "".join(nib.aff2axcodes(params["img_idx"].affine))
            assert len(codes) == 3
            # We don't care about directionality, just which is which dim
            codes = codes.replace("L", "R").replace("P", "A").replace("I", "S")
            idx = codes.index(dict(x="R", y="A", z="S")[ax])
            img_data = np.abs(_get_img_fdata(params["img_idx"]))
            ijk = _cut_coords_to_ijk(cut_coords, params["img_idx"])
            if idx == 0:
                ijk[0] = np.argmax(img_data[:, ijk[1], ijk[2]])
                logger.debug("    MIP: i = %d idx" % (ijk[0],))
            elif idx == 1:
                ijk[1] = np.argmax(img_data[ijk[0], :, ijk[2]])
                logger.debug("    MIP: j = %d idx" % (ijk[1],))
            else:
                ijk[2] = np.argmax(img_data[ijk[0], ijk[1], :])
                logger.debug("    MIP: k = %d idx" % (ijk[2],))
            cut_coords = _ijk_to_cut_coords(ijk, params["img_idx"])

        logger.debug(
            "    Cut coords for %s: (%0.1f, %0.1f, %0.1f) mm"
            % ((ax_name[ax],) + tuple(cut_coords))
        )
        return cut_coords

    def _press(event, params):
        """Manage keypress on the plot."""
        pos = params["lx"].get_xdata()
        idx = params["stc"].time_as_index(pos)[0]
        if event.key == "left":
            idx = max(0, idx - 2)
        elif event.key == "shift+left":
            idx = max(0, idx - 10)
        elif event.key == "right":
            idx = min(params["stc"].shape[1] - 1, idx + 2)
        elif event.key == "shift+right":
            idx = min(params["stc"].shape[1] - 1, idx + 10)
        _update_timeslice(idx, params)
        params["fig"].canvas.draw()

    def _update_timeslice(idx, params):
        params["lx"].set_xdata([idx / params["stc"].sfreq + params["stc"].tmin])
        ax_x, ax_y, ax_z = params["ax_x"], params["ax_y"], params["ax_z"]
        plot_map_callback = params["plot_func"]
        # Crosshairs are the first thing plotted in stat_map, and the last
        # in glass_brain
        idxs = [0, 0, 1] if mode == "stat_map" else [-2, -2, -1]
        cut_coords = (
            ax_y.lines[idxs[0]].get_xdata()[0],
            ax_x.lines[idxs[1]].get_xdata()[0],
            ax_x.lines[idxs[2]].get_ydata()[0],
        )
        ax_x.clear()
        ax_y.clear()
        ax_z.clear()
        params.update({"img_idx": index_img(img, idx)})
        params.update({"title": "Activation (t=%.3f s.)" % params["stc"].times[idx]})
        plot_map_callback(params["img_idx"], title="", cut_coords=cut_coords)

    def _update_vertlabel(loc_idx):
        vert_legend.get_texts()[0].set_text(f"{stc.vertices[0][loc_idx]}")

    @verbose_dec
    def _onclick(event, params, verbose=None):
        """Manage clicks on the plot."""
        ax_x, ax_y, ax_z = params["ax_x"], params["ax_y"], params["ax_z"]
        plot_map_callback = params["plot_func"]
        if event.inaxes is params["ax_time"]:
            idx = params["stc"].time_as_index(event.xdata, use_rounding=True)[0]
            _update_timeslice(idx, params)

        cut_coords = _click_to_cut_coords(event, params)
        if cut_coords is None:
            return  # not in any axes

        ax_x.clear()
        ax_y.clear()
        ax_z.clear()
        plot_map_callback(params["img_idx"], title="", cut_coords=cut_coords)
        loc_idx = _cut_coords_to_idx(cut_coords, params["img_idx"])
        ydata = stc.data[loc_idx]
        if loc_idx is not None:
            ax_time.lines[0].set_ydata(ydata)
        else:
            ax_time.lines[0].set_ydata([0.0])
        _update_vertlabel(loc_idx)
        params["fig"].canvas.draw()

    if mode == "glass_brain":
        subject = _check_subject(stc.subject, subject)
        ras_mni_t = read_ras_mni_t(subject, subjects_dir)
        if not np.allclose(ras_mni_t["trans"], np.eye(4)):
            _print_coord_trans(ras_mni_t, prefix="Transforming subject ", units="mm")
            logger.info("")
            # To get from voxel coords to world coords (i.e., define affine)
            # we would apply img.affine, then also apply ras_mni_t, which
            # transforms from the subject's RAS to MNI RAS. So we left-multiply
            # these.
            img = nib.Nifti1Image(img.dataobj, np.dot(ras_mni_t["trans"], img.affine))
        bg_img = None  # not used
    else:  # stat_map
        if bg_img is None:
            bg_img = "T1.mgz"
        bg_img = _load_subject_mri(bg_img, stc, subject, subjects_dir, "bg_img")

    if initial_time is None:
        time_sl = slice(0, None)
    else:
        initial_time = float(initial_time)
        logger.info("Fixing initial time: %s s" % (initial_time,))
        initial_time = np.argmin(np.abs(stc.times - initial_time))
        time_sl = slice(initial_time, initial_time + 1)
    if initial_pos is None:  # find max pos and (maybe) time
        loc_idx, time_idx = np.unravel_index(
            np.abs(stc.data[:, time_sl]).argmax(), stc.data[:, time_sl].shape
        )
        time_idx += time_sl.start
    else:  # position specified
        initial_pos = np.array(initial_pos, float)
        if initial_pos.shape != (3,):
            raise ValueError(
                "initial_pos must be float ndarray with shape "
                "(3,), got shape %s" % (initial_pos.shape,)
            )
        initial_pos *= 1000
        logger.info("Fixing initial position: %s mm" % (initial_pos.tolist(),))
        loc_idx = _cut_coords_to_idx(initial_pos, img)
        if initial_time is not None:  # time also specified
            time_idx = time_sl.start
        else:  # find the max
            time_idx = np.argmax(np.abs(stc.data[loc_idx]))
    img_idx = index_img(img, time_idx)
    assert img_idx.shape == img.shape[:3]
    del initial_time, initial_pos
    ijk = stc_ijk[loc_idx]
    cut_coords = _ijk_to_cut_coords(ijk, img_idx)
    np.testing.assert_allclose(_cut_coords_to_ijk(cut_coords, img_idx), ijk)
    logger.info(
        "Showing: t = %0.3f s, (%0.1f, %0.1f, %0.1f) mm, "
        "[%d, %d, %d] vox, %d vertex"
        % (
            (stc.times[time_idx],)
            + tuple(cut_coords)
            + tuple(ijk)
            + (stc.vertices[0][loc_idx],)
        )
    )
    del ijk

    # Plot initial figure
    fig, (axes, ax_time) = plt.subplots(2)
    axes.set(xticks=[], yticks=[])
    marker = "o" if len(stc.times) == 1 else None
    ydata = stc.data[loc_idx]
    h = ax_time.plot(stc.times, ydata, color="k", marker=marker)[0]
    if len(stc.times) > 1:
        ax_time.set(xlim=stc.times[[0, -1]])
    ax_time.set(xlabel="Time (s)", ylabel="Activation")
    vert_legend = ax_time.legend([h], [""], title="Vertex")
    _update_vertlabel(loc_idx)
    lx = ax_time.axvline(stc.times[time_idx], color="g")
    fig.tight_layout()

    allow_pos_lims = mode != "glass_brain"
    mapdata = _process_clim(clim, colormap, transparent, stc.data, allow_pos_lims)
    _separate_map(mapdata)
    diverging = "pos_lims" in mapdata["clim"]
    ticks = _get_map_ticks(mapdata)
    colormap, scale_pts = _linearize_map(mapdata)
    del mapdata

    ylim = [min((scale_pts[0], ydata.min())), max((scale_pts[-1], ydata.max()))]
    ylim = np.array(ylim) + np.array([-1, 1]) * 0.05 * np.diff(ylim)[0]
    dup_neg = False
    if stc.data.min() < 0:
        ax_time.axhline(0.0, color="0.5", ls="-", lw=0.5, zorder=2)
        dup_neg = not diverging  # glass brain with signed data
    yticks = list(ticks)
    if dup_neg:
        yticks += [0] + list(-np.array(ticks))
    yticks = np.unique(yticks)
    ax_time.set(yticks=yticks)
    ax_time.set(ylim=ylim)
    del yticks

    if not diverging:  # set eq above iff one-sided
        # there is a bug in nilearn where this messes w/transparency
        # Need to double the colormap
        if (scale_pts < 0).any():
            # XXX We should fix this, but it's hard to get nilearn to
            # use arbitrary bounds :(
            # Should get them to support non-mirrored colorbars, or
            # at least a proper `vmin` for one-sided things.
            # Hopefully this is a sufficiently rare use case!
            raise ValueError(
                "Negative colormap limits for sequential "
                'control points clim["lims"] not supported '
                "currently, consider shifting or flipping the "
                "sign of your data for visualization purposes"
            )
        # due to nilearn plotting weirdness, extend this to go
        # -scale_pts[2]->scale_pts[2] instead of scale_pts[0]->scale_pts[2]
        colormap = _get_cmap(colormap)
        colormap = colormap(
            np.interp(np.linspace(-1, 1, 256), scale_pts / scale_pts[2], [0, 0.5, 1])
        )
        colormap = colors.ListedColormap(colormap)
    vmax = scale_pts[-1]

    # black_bg = True is needed because of some matplotlib
    # peculiarity. See: https://stackoverflow.com/a/34730204
    # Otherwise, event.inaxes does not work for ax_x and ax_z
    plot_kwargs = dict(
        threshold=None,
        axes=axes,
        resampling_interpolation="nearest",
        vmax=vmax,
        figure=fig,
        colorbar=colorbar,
        bg_img=bg_img,
        cmap=colormap,
        black_bg=True,
        symmetric_cbar=True,
    )

    def plot_and_correct(*args, **kwargs):
        axes.clear()
        if params.get("fig_anat") is not None and plot_kwargs["colorbar"]:
            params["fig_anat"]._cbar.ax.clear()
        with warnings.catch_warnings(record=True):  # nilearn bug; ax recreated
            warnings.simplefilter("ignore", DeprecationWarning)
            params["fig_anat"] = partial(plot_func, **plot_kwargs)(*args, **kwargs)
        params["fig_anat"]._cbar.outline.set_visible(False)
        for key in "xyz":
            params.update({"ax_" + key: params["fig_anat"].axes[key].ax})
        # Fix nilearn bug w/cbar background being white
        if plot_kwargs["colorbar"]:
            params["fig_anat"]._cbar.ax.set_facecolor("0.5")
            # adjust one-sided colorbars
            if not diverging:
                _crop_colorbar(params["fig_anat"]._cbar, *scale_pts[[0, -1]])
            params["fig_anat"]._cbar.set_ticks(params["cbar_ticks"])
        if mode == "glass_brain":
            _glass_brain_crosshairs(params, *kwargs["cut_coords"])

    params = dict(
        stc=stc,
        ax_time=ax_time,
        plot_func=plot_and_correct,
        img_idx=img_idx,
        fig=fig,
        lx=lx,
        mode=mode,
        cbar_ticks=ticks,
    )

    plot_and_correct(stat_map_img=params["img_idx"], title="", cut_coords=cut_coords)

    if show:
        plt.show()
    fig.canvas.mpl_connect(
        "button_press_event", partial(_onclick, params=params, verbose=verbose)
    )
    fig.canvas.mpl_connect("key_press_event", partial(_press, params=params))

    return fig


def _check_views(surf, views, hemi, stc=None, backend=None):
    from ._brain.view import views_dicts
    from ..source_estimate import SourceEstimate

    _validate_type(views, (list, tuple, str), "views")
    views = [views] if isinstance(views, str) else list(views)
    if surf == "flat":
        _check_option("views", views, (["auto"], ["flat"]))
        views = ["flat"]
    elif len(views) == 1 and views[0] == "auto":
        views = ["lateral"]
    if views == ["flat"]:
        if stc is not None:
            _validate_type(
                stc, SourceEstimate, "stc", "SourceEstimate when a flatmap is used"
            )
        if backend is not None:
            if backend not in ("pyvistaqt", "notebook"):
                raise RuntimeError(
                    "The PyVista 3D backend must be used to " "plot a flatmap"
                )
    if (views == ["flat"]) ^ (surf == "flat"):  # exactly only one of the two
        raise ValueError(
            'surface="flat" must be used with views="flat", got '
            f"surface={repr(surf)} and views={repr(views)}"
        )
    _check_option("hemi", hemi, ("split", "both", "lh", "rh", "vol", None))
    use_hemi = "lh" if hemi == "split" or hemi is None else hemi
    for vi, v in enumerate(views):
        _check_option(f"views[{vi}]", v, sorted(views_dicts[use_hemi]))
    return views


@verbose
def plot_vector_source_estimates(
    stc,
    subject=None,
    hemi="lh",
    colormap="hot",
    time_label="auto",
    smoothing_steps=10,
    transparent=None,
    brain_alpha=0.4,
    overlay_alpha=None,
    vector_alpha=1.0,
    scale_factor=None,
    time_viewer="auto",
    subjects_dir=None,
    figure=None,
    views="lateral",
    colorbar=True,
    clim="auto",
    cortex="classic",
    size=800,
    background="black",
    foreground=None,
    initial_time=None,
    time_unit="s",
    show_traces="auto",
    src=None,
    volume_options=1.0,
    view_layout="vertical",
    add_data_kwargs=None,
    brain_kwargs=None,
    verbose=None,
):
    """Plot VectorSourceEstimate with PyVista.

    A "glass brain" is drawn and all dipoles defined in the source estimate
    are shown using arrows, depicting the direction and magnitude of the
    current moment at the dipole. Additionally, an overlay is plotted on top of
    the cortex with the magnitude of the current.

    Parameters
    ----------
    stc : VectorSourceEstimate | MixedVectorSourceEstimate
        The vector source estimate to plot.
    %(subject_none)s
        If ``None``, ``stc.subject`` will be used.
    hemi : str, 'lh' | 'rh' | 'split' | 'both'
        The hemisphere to display.
    %(colormap)s
        This should be a sequential colormap.
    %(time_label)s
    smoothing_steps : int
        The amount of smoothing.
    %(transparent)s
    brain_alpha : float
        Alpha value to apply globally to the surface meshes. Defaults to 0.4.
    overlay_alpha : float
        Alpha value to apply globally to the overlay. Defaults to
        ``brain_alpha``.
    vector_alpha : float
        Alpha value to apply globally to the vector glyphs. Defaults to 1.
    scale_factor : float | None
        Scaling factor for the vector glyphs. By default, an attempt is made to
        automatically determine a sane value.
    time_viewer : bool | str
        Display time viewer GUI. Can be "auto", which is True for the PyVista
        backend and False otherwise.

        .. versionchanged:: 0.20
           Added "auto" option and default.
    subjects_dir : str
        The path to the freesurfer subjects reconstructions.
        It corresponds to Freesurfer environment variable SUBJECTS_DIR.
    figure : instance of Figure3D | list | int | None
        If None, a new figure will be created. If multiple views or a
        split view is requested, this must be a list of the appropriate
        length. If int is provided it will be used to identify the PyVista
        figure by it's id or create a new figure with the given id.
    %(views)s
    colorbar : bool
        If True, display colorbar on scene.
    %(clim_onesided)s
    cortex : str or tuple
        Specifies how binarized curvature values are rendered.
        either the name of a preset Brain cortex colorscheme (one of
        'classic', 'bone', 'low_contrast', or 'high_contrast'), or the
        name of a colormap, or a tuple with values (colormap, min,
        max, reverse) to fully specify the curvature colors.
    size : float or tuple of float
        The size of the window, in pixels. can be one number to specify
        a square window, or the (width, height) of a rectangular window.
    background : matplotlib color
        Color of the background of the display window.
    foreground : matplotlib color | None
        Color of the foreground of the display window.
        None will choose black or white based on the background color.
    initial_time : float | None
        The time to display on the plot initially. ``None`` to display the
        first time sample (default).
    time_unit : 's' | 'ms'
        Whether time is represented in seconds ("s", default) or
        milliseconds ("ms").
    %(show_traces)s
    %(src_volume_options)s
    %(view_layout)s
    %(add_data_kwargs)s
    %(brain_kwargs)s
    %(verbose)s

    Returns
    -------
    brain : mne.viz.Brain
        A instance of :class:`mne.viz.Brain`.

    Notes
    -----
    .. versionadded:: 0.15

    If the current magnitude overlay is not desired, set ``overlay_alpha=0``
    and ``smoothing_steps=1``.
    """
    from ..source_estimate import _BaseVectorSourceEstimate

    _validate_type(stc, _BaseVectorSourceEstimate, "stc", "vector source estimate")
    return _plot_stc(
        stc,
        subject=subject,
        surface="white",
        hemi=hemi,
        colormap=colormap,
        time_label=time_label,
        smoothing_steps=smoothing_steps,
        subjects_dir=subjects_dir,
        views=views,
        clim=clim,
        figure=figure,
        initial_time=initial_time,
        time_unit=time_unit,
        background=background,
        time_viewer=time_viewer,
        colorbar=colorbar,
        transparent=transparent,
        brain_alpha=brain_alpha,
        overlay_alpha=overlay_alpha,
        vector_alpha=vector_alpha,
        cortex=cortex,
        foreground=foreground,
        size=size,
        scale_factor=scale_factor,
        show_traces=show_traces,
        src=src,
        volume_options=volume_options,
        view_layout=view_layout,
        add_data_kwargs=add_data_kwargs,
        brain_kwargs=brain_kwargs,
    )


@verbose
def plot_sparse_source_estimates(
    src,
    stcs,
    colors=None,
    linewidth=2,
    fontsize=18,
    bgcolor=(0.05, 0, 0.1),
    opacity=0.2,
    brain_color=(0.7,) * 3,
    show=True,
    high_resolution=False,
    fig_name=None,
    fig_number=None,
    labels=None,
    modes=("cone", "sphere"),
    scale_factors=(1, 0.6),
    verbose=None,
    **kwargs,
):
    """Plot source estimates obtained with sparse solver.

    Active dipoles are represented in a "Glass" brain.
    If the same source is active in multiple source estimates it is
    displayed with a sphere otherwise with a cone in 3D.

    Parameters
    ----------
    src : dict
        The source space.
    stcs : instance of SourceEstimate or list of instances of SourceEstimate
        The source estimates.
    colors : list
        List of colors.
    linewidth : int
        Line width in 2D plot.
    fontsize : int
        Font size.
    bgcolor : tuple of length 3
        Background color in 3D.
    opacity : float in [0, 1]
        Opacity of brain mesh.
    brain_color : tuple of length 3
        Brain color.
    show : bool
        Show figures if True.
    high_resolution : bool
        If True, plot on the original (non-downsampled) cortical mesh.
    fig_name : str
        PyVista figure name.
    fig_number : int
        Matplotlib figure number.
    labels : ndarray or list of ndarray
        Labels to show sources in clusters. Sources with the same
        label and the waveforms within each cluster are presented in
        the same color. labels should be a list of ndarrays when
        stcs is a list ie. one label for each stc.
    modes : list
        Should be a list, with each entry being ``'cone'`` or ``'sphere'``
        to specify how the dipoles should be shown.
        The pivot for the glyphs in ``'cone'`` mode is always the tail
        whereas the pivot in ``'sphere'`` mode is the center.
    scale_factors : list
        List of floating point scale factors for the markers.
    %(verbose)s
    **kwargs : kwargs
        Keyword arguments to pass to renderer.mesh.

    Returns
    -------
    surface : instance of Figure3D
        The 3D figure containing the triangular mesh surface.
    """
    import matplotlib.pyplot as plt

    # Update the backend
    from .backends.renderer import _get_renderer

    linestyles = [
        ("solid", "solid"),  # noqa: E241
        ("dashed", "dashed"),  # noqa: E241
        ("dotted", "dotted"),  # noqa: E241
        ("dashdot", "dashdot"),  # noqa: E241
        ("loosely dotted", (0, (1, 10))),  # noqa: E241
        ("dotted", (0, (1, 1))),  # noqa: E241
        ("densely dotted", (0, (1, 1))),  # noqa: E241
        ("loosely dashed", (0, (5, 10))),  # noqa: E241
        ("dashed", (0, (5, 5))),  # noqa: E241
        ("densely dashed", (0, (5, 1))),  # noqa: E241
        ("loosely dashdotted", (0, (3, 10, 1, 10))),  # noqa: E241
        ("dashdotted", (0, (3, 5, 1, 5))),  # noqa: E241
        ("densely dashdotted", (0, (3, 1, 1, 1))),  # noqa: E241
        ("dashdotdotted", (0, (3, 5, 1, 5, 1, 5))),  # noqa: E241
        ("loosely dashdotdotted", (0, (3, 10, 1, 10, 1, 10))),  # noqa: E241
        ("densely dashdotdotted", (0, (3, 1, 1, 1, 1, 1))),  # noqa: E241
    ]

    known_modes = ["cone", "sphere"]
    if not isinstance(modes, (list, tuple)) or not all(
        mode in known_modes for mode in modes
    ):
        raise ValueError("mode must be a list containing only " '"cone" or "sphere"')
    if not isinstance(stcs, list):
        stcs = [stcs]
    if labels is not None and not isinstance(labels, list):
        labels = [labels]

    if colors is None:
        colors = _get_color_list()

    linestyles = cycle(linestyles)
    linestyles = [next(linestyles)[1] for _ in range(len(stcs))]

    # Show 3D
    lh_points = src[0]["rr"]
    rh_points = src[1]["rr"]
    points = np.r_[lh_points, rh_points]

    lh_normals = src[0]["nn"]
    rh_normals = src[1]["nn"]
    normals = np.r_[lh_normals, rh_normals]

    if high_resolution:
        use_lh_faces = src[0]["tris"]
        use_rh_faces = src[1]["tris"]
    else:
        use_lh_faces = src[0]["use_tris"]
        use_rh_faces = src[1]["use_tris"]

    use_faces = np.r_[use_lh_faces, lh_points.shape[0] + use_rh_faces]

    points *= 170

    vertnos = [np.r_[stc.lh_vertno, lh_points.shape[0] + stc.rh_vertno] for stc in stcs]
    unique_vertnos = np.unique(np.concatenate(vertnos).ravel())

    renderer = _get_renderer(bgcolor=bgcolor, size=(600, 600), name=fig_name)
    renderer.mesh(
        x=points[:, 0],
        y=points[:, 1],
        z=points[:, 2],
        triangles=use_faces,
        color=brain_color,
        opacity=opacity,
        backface_culling=True,
        shading=True,
        normals=normals,
        **kwargs,
    )

    # Show time courses
    fig = plt.figure(fig_number)
    fig.clf()
    ax = fig.add_subplot(111)

    colors = cycle(colors)

    logger.info("Total number of active sources: %d" % len(unique_vertnos))

    if labels is not None:
        colors = [
            next(colors) for _ in range(np.unique(np.concatenate(labels).ravel()).size)
        ]

    for idx, v in enumerate(unique_vertnos):
        # get indices of stcs it belongs to
        ind = [k for k, vertno in enumerate(vertnos) if v in vertno]
        is_common = len(ind) > 1

        if labels is None:
            c = next(colors)
        else:
            # if vertex is in different stcs than take label from first one
            c = colors[labels[ind[0]][vertnos[ind[0]] == v]]

        mode = modes[1] if is_common else modes[0]
        scale_factor = scale_factors[1] if is_common else scale_factors[0]

        if isinstance(scale_factor, (np.ndarray, list, tuple)) and len(
            unique_vertnos
        ) == len(scale_factor):
            scale_factor = scale_factor[idx]

        x, y, z = points[v]
        nx, ny, nz = normals[v]
        renderer.quiver3d(
            x=x,
            y=y,
            z=z,
            u=nx,
            v=ny,
            w=nz,
            color=_to_rgb(c),
            mode=mode,
            scale=scale_factor,
        )

        for k in ind:
            vertno = vertnos[k]
            mask = vertno == v
            assert np.sum(mask) == 1
            linestyle = linestyles[k]
            ax.plot(
                1e3 * stcs[k].times,
                1e9 * stcs[k].data[mask].ravel(),
                c=c,
                linewidth=linewidth,
                linestyle=linestyle,
            )

    ax.set_xlabel("Time (ms)", fontsize=18)
    ax.set_ylabel("Source amplitude (nAm)", fontsize=18)

    if fig_name is not None:
        ax.set_title(fig_name)
    plt_show(show)

    renderer.show()
    return renderer.scene()


@verbose
def plot_dipole_locations(
    dipoles,
    trans=None,
    subject=None,
    subjects_dir=None,
    mode="orthoview",
    coord_frame="mri",
    idx="gof",
    show_all=True,
    ax=None,
    block=False,
    show=True,
    scale=None,
    color=None,
    *,
    highlight_color="r",
    fig=None,
    title=None,
    head_source="seghead",
    surf="pial",
    width=None,
    verbose=None,
):
    """Plot dipole locations.

    If mode is set to 'arrow' or 'sphere', only the location of the first
    time point of each dipole is shown else use the show_all parameter.

    Parameters
    ----------
    dipoles : list of instances of Dipole | Dipole
        The dipoles to plot.
    trans : dict | None
        The mri to head trans.
        Can be None with mode set to '3d'.
    subject : str | None
        The FreeSurfer subject name (will be used to set the FreeSurfer
        environment variable ``SUBJECT``).
        Can be ``None`` with mode set to ``'3d'``.
    %(subjects_dir)s
    mode : str
        Can be:

        ``'arrow'`` or ``'sphere'``
            Plot in 3D mode using PyVista with the given glyph type.
        ``'orthoview'``
            Plot in matplotlib ``Axes3D`` using matplotlib with MRI slices
            shown on the sides of a cube, with the dipole(s) shown as arrows
            extending outward from a dot (i.e., the arrows pivot on the tail).
        ``'outlines'``
            Plot in matplotlib ``Axes`` using a quiver of arrows for the
            dipoles in three axes (axial, coronal, and sagittal views),
            with the arrow pivoting in the middle of the arrow.

        .. versionchanged:: 1.1
           Added support for ``'outlines'``.
    coord_frame : str
        Coordinate frame to use: 'head' or 'mri'. Can also be 'mri_rotated'
        when mode equals ``'outlines'``. Defaults to 'mri'.

        .. versionadded:: 0.14.0
        .. versionchanged:: 1.1
           Added support for ``'mri_rotated'``.
    idx : int | 'gof' | 'amplitude'
        Index of the initially plotted dipole. Can also be 'gof' to plot the
        dipole with highest goodness of fit value or 'amplitude' to plot the
        dipole with the highest amplitude. The dipoles can also be browsed
        through using up/down arrow keys or mouse scroll. Defaults to 'gof'.
        Only used if mode equals 'orthoview'.

        .. versionadded:: 0.14.0
    show_all : bool
        Whether to always plot all the dipoles. If ``True`` (default), the
        active dipole is plotted as a red dot and its location determines the
        shown MRI slices. The non-active dipoles are plotted as small blue
        dots. If ``False``, only the active dipole is plotted.
        Only used if ``mode='orthoview'``.

        .. versionadded:: 0.14.0
    ax : instance of matplotlib Axes3D | list of matplotlib Axes | None
        Axes to plot into. If None (default), axes will be created.
        If mode equals ``'orthoview'``, must be a single ``Axes3D``.
        If mode equals ``'outlines'``, must be a list of three ``Axes``.

        .. versionadded:: 0.14.0
    block : bool
        Whether to halt program execution until the figure is closed. Defaults
        to False.
        Only used if mode equals 'orthoview'.

        .. versionadded:: 0.14.0
    show : bool
        Show figure if True. Defaults to True.
        Only used if mode equals 'orthoview'.
    scale : float
        The scale (size in meters) of the dipoles if ``mode`` is not
        ``'orthoview'``. The default is 0.03 when mode is ``'outlines'`` and
        0.005 otherwise.
    color : tuple
        The color of the dipoles.
        The default (None) will use ``'y'`` if mode is ``'orthoview'`` and
        ``show_all`` is True, else 'r'. Can also be a list of colors to use
        when mode is ``'outlines'``.

        .. versionchanged:: 0.19.0
           Color is now passed in orthoview mode.
    highlight_color : color
        The highlight color. Only used in orthoview mode with
        ``show_all=True``.

        .. versionadded:: 0.19.0
    fig : instance of Figure3D | None
        3D figure in which to plot the alignment.
        If ``None``, creates a new 600x600 pixel figure with black background.
        Only used when mode is ``'arrow'`` or ``'sphere'``.

        .. versionadded:: 0.19.0
    title : str | None
        The title of the figure if ``mode='orthoview'`` (ignored for all other
        modes). If ``None``, dipole number and its properties (amplitude,
        orientation etc.) will be shown. Defaults to ``None``.

        .. versionadded:: 0.21.0
    %(head_source)s
        Only used when mode equals ``'outlines'``.

        .. versionadded:: 1.1
    surf : str | None
        Brain surface to show outlines for, can be ``'white'``, ``'pial'``, or
        ``None``. Only used when mode is ``'outlines'``.

        .. versionadded:: 1.1
    width : float | None
        Width of the matplotlib quiver arrow, see
        :meth:`matplotlib:matplotlib.axes.Axes.quiver`. If None (default),
        when mode is ``'outlines'`` 0.015 will be used, and when mode is
        ``'orthoview'`` the matplotlib default is used.
    %(verbose)s

    Returns
    -------
    fig : instance of Figure3D or matplotlib.figure.Figure
        The PyVista figure or matplotlib Figure.

    Notes
    -----
    .. versionadded:: 0.9.0
    """
    _validate_type(mode, str, "mode")
    _validate_type(coord_frame, str, "coord_frame")
    _check_option("mode", mode, ("orthoview", "outlines", "arrow", "sphere"))
    if mode in ("orthoview", "outlines"):
        subjects_dir = str(get_subjects_dir(subjects_dir, raise_error=True))
    kwargs = dict(
        trans=trans,
        subject=subject,
        subjects_dir=subjects_dir,
        coord_frame=coord_frame,
        ax=ax,
        block=block,
        show=show,
        color=color,
        title=title,
        width=width,
    )
    dipoles = _check_concat_dipoles(dipoles)
    if mode == "orthoview":
        fig = _plot_dipole_mri_orthoview(
            dipoles,
            idx=idx,
            show_all=show_all,
            highlight_color=highlight_color,
            **kwargs,
        )
    elif mode == "outlines":
        fig = _plot_dipole_mri_outlines(
            dipoles, head_source=head_source, surf=surf, scale=scale, **kwargs
        )
    else:
        assert mode in ("arrow", "sphere"), mode
        fig = _plot_dipole_3d(
            dipoles,
            trans=trans,
            coord_frame=coord_frame,
            color=color,
            fig=fig,
            scale=scale,
            mode=mode,
        )

    return fig


def snapshot_brain_montage(fig, montage, hide_sensors=True):
    """Take a snapshot of a PyVista Scene and project channels onto 2d coords.

    Note that this will take the raw values for 3d coordinates of each channel,
    without applying any transforms. If brain images are flipped up/dn upon
    using `~matplotlib.pyplot.imshow`, check your matplotlib backend as this
    behavior changes.

    Parameters
    ----------
    fig : instance of Figure3D
        The figure on which you've plotted electrodes using
        :func:`mne.viz.plot_alignment`.
    montage : instance of DigMontage or Info | dict
        The digital montage for the electrodes plotted in the scene. If
        :class:`~mne.Info`, channel positions will be pulled from the ``loc``
        field of ``chs``. dict should have ch:xyz mappings.
    hide_sensors : bool
        Whether to remove the spheres in the scene before taking a snapshot.
        The sensors will always be shown in the final figure. If you want an
        image of just the brain, use :class:`mne.viz.Brain` instead.

    Returns
    -------
    xy : array, shape (n_channels, 2)
        The 2d location of each channel on the image of the current scene view.
    im : array, shape (m, n, 3)
        The screenshot of the current scene view.
    """
    from ..channels import DigMontage
    from .. import Info

    # Update the backend
    from .backends.renderer import _get_renderer

    if fig is None:
        raise ValueError("The figure must have a scene")
    if isinstance(montage, DigMontage):
        chs = montage._get_ch_pos()
        ch_names, xyz = zip(*[(ich, ixyz) for ich, ixyz in chs.items()])
    elif isinstance(montage, Info):
        xyz = [ich["loc"][:3] for ich in montage["chs"]]
        ch_names = [ich["ch_name"] for ich in montage["chs"]]
    elif isinstance(montage, dict):
        if not all(len(ii) == 3 for ii in montage.values()):
            raise ValueError("All electrode positions must be length 3")
        ch_names, xyz = zip(*[(ich, ixyz) for ich, ixyz in montage.items()])
    else:
        raise TypeError(
            "montage must be an instance of `DigMontage`, `Info`," " or `dict`"
        )

    # initialize figure
    renderer = _get_renderer(fig, show=True)

    xyz = np.vstack(xyz)
    proj = renderer.project(xyz=xyz, ch_names=ch_names)
    if hide_sensors is True:
        proj.visible(False)

    im = renderer.screenshot()
    proj.visible(True)
    return proj.xy, im


def _plot_dipole_mri_orthoview(
    dipole,
    trans,
    subject,
    subjects_dir=None,
    coord_frame="head",
    idx="gof",
    show_all=True,
    ax=None,
    block=False,
    show=True,
    color=None,
    highlight_color="r",
    title=None,
    width=None,
):
    """Plot dipoles on top of MRI slices in 3-D."""
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    _import_nibabel("plotting MRI slices")

    _check_option("coord_frame", coord_frame, ["head", "mri"])

    if idx == "gof":
        idx = np.argmax(dipole.gof)
    elif idx == "amplitude":
        idx = np.argmax(np.abs(dipole.amplitude))
    else:
        idx = _ensure_int(idx, "idx", 'an int or one of ["gof", "amplitude"]')

    vox, ori, pos, data = _get_dipole_loc(
        dipole, trans, subject, subjects_dir, coord_frame
    )

    dims = len(data)  # Symmetric size assumed.
    dd = dims // 2
    if ax is None:
        fig, ax = plt.subplots(1, subplot_kw=dict(projection="3d"))
    else:
        _validate_type(ax, Axes3D, "ax", "Axes3D", extra='when mode is "orthoview"')
        fig = ax.get_figure()

    gridx, gridy = np.meshgrid(
        np.linspace(-dd, dd, dims), np.linspace(-dd, dd, dims), indexing="ij"
    )
    params = {
        "ax": ax,
        "data": data,
        "idx": idx,
        "dipole": dipole,
        "vox": vox,
        "gridx": gridx,
        "gridy": gridy,
        "ori": ori,
        "coord_frame": coord_frame,
        "show_all": show_all,
        "pos": pos,
        "color": color,
        "highlight_color": highlight_color,
        "title": title,
        "width": width,
    }
    _plot_dipole(**params)
    ax.view_init(elev=30, azim=-140)

    callback_func = partial(_dipole_changed, params=params)
    fig.canvas.mpl_connect("scroll_event", callback_func)
    fig.canvas.mpl_connect("key_press_event", callback_func)

    plt_show(show, block=block)
    return fig


RAS_AFFINE = np.eye(4)
RAS_AFFINE[:3, 3] = [-128] * 3
RAS_SHAPE = (256, 256, 256)


def _get_dipole_loc(dipole, trans, subject, subjects_dir, coord_frame):
    """Get the dipole locations and orientations."""
    import nibabel as nib
    from nibabel.processing import resample_from_to

    _check_option("coord_frame", coord_frame, ["head", "mri"])

    subjects_dir = str(get_subjects_dir(subjects_dir=subjects_dir, raise_error=True))
    t1_fname = op.join(subjects_dir, subject, "mri", "T1.mgz")
    t1 = nib.load(t1_fname)
    # Do everything in mm here to make life slightly easier
    vox_ras_t, _, mri_ras_t, _, _ = _read_mri_info(t1_fname, units="mm")
    head_mri_t = _get_trans(trans, fro="head", to="mri")[0].copy()
    head_mri_t["trans"][:3, 3] *= 1000  # m→mm
    del trans
    pos = dipole.pos * 1e3  # m→mm
    ori = dipole.ori
    # Figure out how to always resample to an identity, 256x256x256 RAS:
    #
    # 1. Resample to head or MRI surface RAS (the conditional), but also
    # 2. Resample to what will work for the standard 1mm** RAS_AFFINE (resamp)
    #
    # We could do this with two resample_from_to calls, but it's cleaner,
    # faster, and we get fewer boundary artifacts if we do it in one shot.
    # So first olve usamp s.t. ``upsamp @ vox_ras_t == RAS_AFFINE`` (2):
    upsamp = np.linalg.solve(vox_ras_t["trans"].T, RAS_AFFINE.T).T
    # Now figure out how we would resample from RAS to head or MRI coords:
    if coord_frame == "head":
        dest_ras_t = combine_transforms(head_mri_t, mri_ras_t, "head", "ras")["trans"]
    else:
        pos = apply_trans(head_mri_t, pos)
        ori = apply_trans(head_mri_t, dipole.ori, move=False)
        dest_ras_t = mri_ras_t["trans"]
    # The order here is wacky because we need `resample_from_to` to operate
    # in a reverse order
    affine = np.dot(np.dot(dest_ras_t, upsamp), vox_ras_t["trans"])
    t1 = resample_from_to(t1, (RAS_SHAPE, affine), order=0)
    # Now we could do:
    #
    #    t1 = SpatialImage(t1.dataobj, RAS_AFFINE)
    #
    # And t1 would be in our destination (mri or head) space. But we don't
    # need to construct the image -- let's just get our voxel coords and data:
    vox = apply_trans(np.linalg.inv(RAS_AFFINE), pos)
    t1_data = _get_img_fdata(t1)
    return vox, ori, pos, t1_data


def _plot_dipole(
    ax,
    data,
    vox,
    idx,
    dipole,
    gridx,
    gridy,
    ori,
    coord_frame,
    show_all,
    pos,
    color,
    highlight_color,
    title,
    width,
):
    """Plot dipoles."""
    import matplotlib.pyplot as plt

    xidx, yidx, zidx = np.round(vox[idx]).astype(int)
    xslice = data[xidx]
    yslice = data[:, yidx]
    zslice = data[:, :, zidx]

    ori = ori[idx]
    if color is None:
        color = "y" if show_all else "r"
    color = np.array(_to_rgb(color, alpha=True))
    highlight_color = np.array(
        _to_rgb(highlight_color, name="highlight_color", alpha=True)
    )
    if show_all:
        colors = np.repeat(color[np.newaxis], len(vox), axis=0)
        colors[idx] = highlight_color
        size = np.repeat(5, len(vox))
        size[idx] = 20
        visible = np.arange(len(vox))
    else:
        colors = color
        size = 20
        visible = idx

    offset = np.min(gridx)
    xyz = pos
    ax.scatter(
        xs=xyz[visible, 0],
        ys=xyz[visible, 1],
        zs=xyz[visible, 2],
        zorder=2,
        s=size,
        facecolor=colors,
    )
    xx = np.linspace(offset, xyz[idx, 0], xidx)
    yy = np.linspace(offset, xyz[idx, 1], yidx)
    zz = np.linspace(offset, xyz[idx, 2], zidx)
    ax.plot(
        xx,
        np.repeat(xyz[idx, 1], len(xx)),
        zs=xyz[idx, 2],
        zorder=1,
        linestyle="-",
        color=highlight_color,
    )
    ax.plot(
        np.repeat(xyz[idx, 0], len(yy)),
        yy,
        zs=xyz[idx, 2],
        zorder=1,
        linestyle="-",
        color=highlight_color,
    )
    ax.plot(
        np.repeat(xyz[idx, 0], len(zz)),
        np.repeat(xyz[idx, 1], len(zz)),
        zs=zz,
        zorder=1,
        linestyle="-",
        color=highlight_color,
    )
    q_kwargs = dict(length=50, color=highlight_color, pivot="tail")
    if width is not None:
        q_kwargs["width"] = width
    ax.quiver(xyz[idx, 0], xyz[idx, 1], xyz[idx, 2], ori[0], ori[1], ori[2], **q_kwargs)
    dims = np.array([(len(data) / -2.0), (len(data) / 2.0)])
    ax.set(xlim=-dims, ylim=-dims, zlim=dims)

    # Plot slices
    ax.contourf(
        xslice, gridx, gridy, offset=offset, zdir="x", cmap="gray", zorder=0, alpha=0.5
    )
    ax.contourf(
        gridx, yslice, gridy, offset=offset, zdir="y", cmap="gray", zorder=0, alpha=0.5
    )
    ax.contourf(
        gridx, gridy, zslice, offset=offset, zdir="z", cmap="gray", zorder=0, alpha=0.5
    )

    # Plot orientations
    args = np.array([list(xyz[idx]) + list(ori)] * 3)
    for ii in range(3):
        args[ii, [ii, ii + 3]] = [offset + 0.5, 0]  # half a mm inward  (z ord)
    ax.quiver(*args.T, alpha=0.75, **q_kwargs)

    # These are the only two options
    coord_frame_name = "Head" if coord_frame == "head" else "MRI"

    if title is None:
        title = "Dipole #%s / %s @ %.3fs, GOF: %.1f%%, %.1fnAm\n%s: " % (
            idx + 1,
            len(dipole.times),
            dipole.times[idx],
            dipole.gof[idx],
            dipole.amplitude[idx] * 1e9,
            coord_frame_name,
        ) + "(%0.1f, %0.1f, %0.1f) mm" % tuple(xyz[idx])

    ax.get_figure().suptitle(title)

    ax.set_xlabel("x")
    ax.set_ylabel("y")
    ax.set_zlabel("z")

    plt.draw()


def _dipole_changed(event, params):
    """Handle dipole plotter scroll/key event."""
    if event.key is not None:
        if event.key == "up":
            params["idx"] += 1
        elif event.key == "down":
            params["idx"] -= 1
        else:  # some other key
            return
    elif event.step > 0:  # scroll event
        params["idx"] += 1
    else:
        params["idx"] -= 1
    params["idx"] = min(max(0, params["idx"]), len(params["dipole"].pos) - 1)
    params["ax"].clear()
    _plot_dipole(**params)


@fill_doc
def plot_brain_colorbar(
    ax,
    clim,
    colormap="auto",
    transparent=True,
    orientation="vertical",
    label="Activation",
    bgcolor="0.5",
):
    """Plot a colorbar that corresponds to a brain activation map.

    Parameters
    ----------
    ax : instance of Axes
        The Axes to plot into.
    %(clim)s
    %(colormap)s
    %(transparent)s
    orientation : str
        Orientation of the colorbar, can be "vertical" or "horizontal".
    label : str
        The colorbar label.
    bgcolor : color
        The color behind the colorbar (for alpha blending).

    Returns
    -------
    cbar : instance of ColorbarBase
        The colorbar.

    Notes
    -----
    .. versionadded:: 0.19
    """
    from matplotlib.colorbar import ColorbarBase
    from matplotlib.colors import Normalize

    mapdata = _process_clim(clim, colormap, transparent)
    ticks = _get_map_ticks(mapdata)
    colormap, lims = _linearize_map(mapdata)
    del mapdata
    norm = Normalize(vmin=lims[0], vmax=lims[2])
    cbar = ColorbarBase(
        ax, cmap=colormap, norm=norm, ticks=ticks, label=label, orientation=orientation
    )
    # make the colorbar background match the brain color
    cbar.ax.set(facecolor=bgcolor)
    # remove the colorbar frame except for the line containing the ticks
    cbar.outline.set_visible(False)
    cbar.ax.set_frame_on(True)
    for key in ("left", "top", "bottom" if orientation == "vertical" else "right"):
        ax.spines[key].set_visible(False)
    return cbar


@dataclass()
class _3d_Options:
    antialias: Optional[bool]
    depth_peeling: Optional[bool]
    smooth_shading: Optional[bool]
    multi_samples: Optional[int]


_3d_options = _3d_Options(
    antialias=None,
    depth_peeling=None,
    smooth_shading=None,
    multi_samples=None,
)
_3d_default = _3d_Options(
    antialias="true",
    depth_peeling="true",
    smooth_shading="true",
    multi_samples="4",
)


def set_3d_options(
    antialias=None, depth_peeling=None, smooth_shading=None, *, multi_samples=None
):
    """Set 3D rendering options.

    Parameters
    ----------
    antialias : bool | None
        If bool, whether to enable or disable full-screen anti-aliasing.
        False is useful when renderers have problems (such as software
        MESA renderers). If None, use the default setting. This option
        can also be controlled using an environment variable, e.g.,
        ``MNE_3D_OPTION_ANTIALIAS=false``.
    depth_peeling : bool | None
        If bool, whether to enable or disable accurate transparency.
        False is useful when renderers have problems (for instance
        while X forwarding on remote servers). If None, use the default
        setting. This option can also be controlled using an environment
        variable, e.g., ``MNE_3D_OPTION_DEPTH_PEELING=false``.
    smooth_shading : bool | None
        If bool, whether to enable or disable smooth color transitions
        between polygons. False is useful on certain configurations
        where this type of shading is not supported or for performance
        reasons. This option can also be controlled using an environment
        variable, e.g., ``MNE_3D_OPTION_SMOOTH_SHADING=false``.
    multi_samples : int
        Number of multi-samples. Should be 1 for MESA for volumetric rendering
        to work properly.

        .. versionadded:: 1.1

    Notes
    -----
    .. versionadded:: 0.21.0
    """
    if antialias is not None:
        _3d_options.antialias = bool(antialias)
    if depth_peeling is not None:
        _3d_options.depth_peeling = bool(depth_peeling)
    if smooth_shading is not None:
        _3d_options.smooth_shading = bool(smooth_shading)
    if multi_samples is not None:
        _3d_options.multi_samples = int(multi_samples)


def _get_3d_option(key):
    _validate_type(key, "str", "key")
    opt = getattr(_3d_options, key)
    if opt is None:  # parse get_config (and defaults)
        default_value = getattr(_3d_default, key)
        opt = get_config(f"MNE_3D_OPTION_{key.upper()}", default_value)
        if key == "multi_samples":
            opt = int(opt)
        else:
            opt = opt.lower() == "true"
    return opt