diff --git a/doc/changes/dev/14016.newfeature.rst b/doc/changes/dev/14016.newfeature.rst new file mode 100644 index 00000000000..51bbbeb8fc0 --- /dev/null +++ b/doc/changes/dev/14016.newfeature.rst @@ -0,0 +1 @@ +:class:`mne.viz.Brain`'s ``silhouette`` option now accepts a spacing string (e.g. ``"ico5"``, now the default) to decimate using the same vertex-picking as :func:`mne.setup_source_space` instead of a generic mesh decimation. :func:`mne.viz.plot_alignment` and :class:`mne.viz.Brain` should also be faster when plotting many sensors (e.g., whole-head MEG systems), by `Eric Larson`_. diff --git a/doc/sphinxext/mne_doc_utils.py b/doc/sphinxext/mne_doc_utils.py index 8e7415b3375..9575ef72cf1 100644 --- a/doc/sphinxext/mne_doc_utils.py +++ b/doc/sphinxext/mne_doc_utils.py @@ -7,6 +7,7 @@ import functools import gc import os +import sys import time import warnings from pathlib import Path @@ -14,14 +15,11 @@ import numpy as np import pyvista import sphinx.util.logging +from refleak.testing import assert_no_instances from sphinx.errors import ExtensionError import mne -from mne.utils import ( - _assert_no_instances, - _get_extra_data_path, - sizeof_fmt, -) +from mne.utils import Bunch, _get_extra_data_path, sizeof_fmt from mne.viz import Brain sphinx_logger = sphinx.util.logging.getLogger("mne") @@ -177,6 +175,18 @@ def reset_modules(gallery_conf, fname, when): neo.io.stimfitio.STFIO_ERR = None except Exception: pass + # IPython.core.completer does `import __main__` at import time, permanently + # capturing whatever sys.modules['__main__'] was at that moment. Since SG + # temporarily swaps sys.modules['__main__'] to each example's throwaway + # namespace while it runs, if that import happens to fire during one of + # those windows, it pins that example's globals (and anything reachable + # from them) alive for the rest of the process. + try: + import IPython.core.completer + + IPython.core.completer.__main__ = sys.modules["__main__"] + except Exception: + pass gc.collect() # Agg does not call close_event so let's clean up on our own :( @@ -192,19 +202,21 @@ def reset_modules(gallery_conf, fname, when): # to just test MNEQtBrowser skips = os.getenv("MNE_SKIP_INSTANCE_ASSERTIONS", "").lower() prefix = "" + request = Bunch() # just give it something to say "we have done GC already" + request.node = Bunch() if skips not in ("true", "1", "all"): prefix = "Clean " skips = skips.split(",") if "brain" not in skips: - _assert_no_instances(Brain, when) # calls gc.collect() + assert_no_instances(Brain, when=when, request=request) if Plotter is not None and "plotter" not in skips: - _assert_no_instances(Plotter, when) + assert_no_instances(Plotter, when=when, request=request) if BackgroundPlotter is not None and "backgroundplotter" not in skips: - _assert_no_instances(BackgroundPlotter, when) + assert_no_instances(BackgroundPlotter, when=when, request=request) if vtkPolyData is not None and "vtkpolydata" not in skips: - _assert_no_instances(vtkPolyData, when) + assert_no_instances(vtkPolyData, when=when, request=request) if "_renderer" not in skips: - _assert_no_instances(_Renderer, when) + assert_no_instances(_Renderer, when=when, request=request) if MNEQtBrowser is not None and "mneqtbrowser" not in skips: # Ensure any manual fig.close() events get properly handled from mne_qt_browser._pg_figure import QApplication @@ -213,7 +225,7 @@ def reset_modules(gallery_conf, fname, when): if inst is not None: for _ in range(2): inst.processEvents() - _assert_no_instances(MNEQtBrowser, when) + assert_no_instances(MNEQtBrowser, when=when, request=request) # This will overwrite some Sphinx printing but it's useful # for memory timestamps if os.getenv("SG_STAMP_STARTS", "").lower() == "true": diff --git a/doc/sphinxext/related_software.py b/doc/sphinxext/related_software.py index 1fe0c62ebe2..bc3d2a34cc3 100644 --- a/doc/sphinxext/related_software.py +++ b/doc/sphinxext/related_software.py @@ -29,10 +29,10 @@ from sphinx.errors import ExtensionError from sphinx.util.display import status_iterator -# If a package is in MNE-Installers (preferred method), no need to add it here. -# But still add it to doc/sphinxext/related_software.txt! +# 1. If a package is in MNE-Installers (preferred method), no need to add it here. +# But still add it to doc/sphinxext/related_software.txt! -# If it's available on PyPI, add it to this set: +# 2. If it's available on PyPI, add it to this set: PYPI_PACKAGES = { "cross-domain-saliency-maps", "meggie", @@ -41,7 +41,7 @@ "zuna", } -# If it's not available on PyPI, add it to this dict: +# 3. If it's not available on PyPI, add it to this dict: MANUAL_PACKAGES = { # TODO: These packages are not pip-installable as of 2025/11/19, so we have to # manually populate them -- should open issues on their package repos. diff --git a/mne/conftest.py b/mne/conftest.py index b6da95824e3..b42f4123ae1 100644 --- a/mne/conftest.py +++ b/mne/conftest.py @@ -22,6 +22,7 @@ import pytest from packaging.version import Version from pytest import StashKey, register_assert_rewrite +from refleak.testing import assert_no_instances, gc_collect_once # Any `assert` statements in our testing functions should be verbose versions register_assert_rewrite("mne.utils._testing") @@ -38,7 +39,6 @@ from mne.stats import cluster_level from mne.utils import ( Bunch, - _assert_no_instances, _check_qt_version, _chmod_rw_R, _pl, @@ -628,10 +628,10 @@ def triaxial_evoked(triaxial_raw): @pytest.fixture -def garbage_collect(): +def garbage_collect(request): """Garbage collect on exit.""" yield - gc.collect() + gc_collect_once(request) @pytest.fixture @@ -687,7 +687,9 @@ def pg_backend(request, garbage_collect): mne_qt_browser._browser_instances.clear() if not _test_passed(request): return - _assert_no_instances(MNEQtBrowser, f"Closure of {request.node.name}") + assert_no_instances( + MNEQtBrowser, f"Closure of {request.node.name}", request=request + ) @pytest.fixture( @@ -1081,8 +1083,13 @@ def brain_gc(request): close_func() if not _test_passed(request): return - _assert_no_instances(Brain, "after") - # Check VTK + gc_collect_once(request) + # Brain._instances is a WeakSet populated only when MNE_3D_BACKEND_TESTING + # is set (see Brain.__init__), so use it instead of a slow gc.get_objects() + # scan of the whole process to check for lingering Brain instances. + assert_no_instances(Brain, "after", request=request, objs=list(Brain._instances)) + # Check VTK -- these aren't individually tracked, so we still need a full + # heap scan here. objs = gc.get_objects() bad = list() for o in objs: diff --git a/mne/gui/tests/test_coreg.py b/mne/gui/tests/test_coreg.py index feb28ed770f..13ef01ca827 100644 --- a/mne/gui/tests/test_coreg.py +++ b/mne/gui/tests/test_coreg.py @@ -263,7 +263,10 @@ def test_coreg_gui_pyvista_basic(tmp_path, monkeypatch, renderer_interactive_pyv coreg._redraw(verbose="debug") assert "Drawing meg sensors" in log.getvalue() assert coreg._actors["helmet"] is not None - assert len(coreg._actors["sensors"]) == 306 + # coil surfaces sharing a shape are GPU-instanced into a single + # mesh/actor; this Neuromag sample dataset has 2 distinct MEG coil + # shapes (magnetometer + planar gradiometer) + assert len(coreg._actors["sensors"]) == 2 assert coreg._orient_glyphs assert coreg._scale_by_distance assert coreg._mark_inside diff --git a/mne/surface.py b/mne/surface.py index c3075bf7d9d..9153e080202 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1269,6 +1269,20 @@ def _tessellate_sphere(mylevel): return rr, tris +def _decimate_surface_ico_oct(subject, subjects_dir, hemi, surf, spacing): + from .source_space._source_space import _check_spacing + + stype, _, ico_surf, _ = _check_spacing(spacing, verbose=False) + subjects_dir = Path(subjects_dir) + surf_fname = subjects_dir / subject / "surf" / f"{hemi}.{surf}" + dec = _create_surf_spacing(surf_fname, hemi, subject, stype, ico_surf, subjects_dir) + vertno, use_tris = dec["vertno"], dec["use_tris"] + lut = np.zeros(dec["np"], int) + lut[vertno] = np.arange(len(vertno)) + tris = lut[use_tris] + return vertno, tris + + def _create_surf_spacing(surf, hemi, subject, stype, ico_surf, subjects_dir): """Load a surf and use the subdivided icosahedron to get points.""" # Based on load_source_space_surf_spacing() in load_source_space.c diff --git a/mne/utils/config.py b/mne/utils/config.py index b46a8d5c9c5..8b57fd837cf 100644 --- a/mne/utils/config.py +++ b/mne/utils/config.py @@ -863,6 +863,7 @@ def sys_info( "pytest-qt", "pytest-rerunfailures", "pytest-timeout", + "refleak", "codespell", "ipython", "mypy", diff --git a/mne/utils/misc.py b/mne/utils/misc.py index 38c14cbceb0..bf8478094b0 100644 --- a/mne/utils/misc.py +++ b/mne/utils/misc.py @@ -5,7 +5,6 @@ # Copyright the MNE-Python contributors. import fnmatch -import gc import hashlib import inspect import os @@ -353,80 +352,13 @@ def _file_like(obj): return all(callable(getattr(obj, name, None)) for name in ("read", "seek")) -def _fullname(obj, *, referent=None): - klass = obj.__class__ - module = klass.__module__ - name = klass.__qualname__ - if module != "builtins": - name = f"{module}.{name}" - if referent is not None: - if isinstance(obj, list | tuple): - for ii, item in enumerate(obj): - if item is referent: - name += f"[{ii}]" - break - elif isinstance(obj, dict): - for key, value in obj.items(): - if key is referent: - name += "-key" - break - if value is referent: - name += f"[{key!r}]" - break - return name - - +# Low-effort backward compat wrapper in case other MNE libraries use this function def _assert_no_instances(cls, when=""): + from refleak.testing import assert_no_instances + __tracebackhide__ = True - n = 0 - ref = list() - gc.collect() - objs = gc.get_objects() - for obj in objs: # e.g., vtkPolyData, Brain, Plotter, etc. - try: - check = isinstance(obj, cls) - except Exception: # such as a weakref - check = False - if check: - if cls.__name__ == "Brain": - ref.append(f"Brain._cleaned = {getattr(obj, '_cleaned', None)}") - rr = gc.get_referrers(obj) - count = 0 - for r in rr: # e.g., list, dict, etc. that holds the reference to obj - if ( - r is not objs - and r is not globals() - and r is not locals() - and not inspect.isframe(r) - ): - name = _fullname(r, referent=obj) - if isinstance(r, list | dict | tuple): - rep = f"len={len(r)}" - r_ = gc.get_referrers(r) - types = list() - for x in r_: - types.append(_fullname(x, referent=r)) - types = " / ".join(sorted(types)) - rep += f" | {len(r_)} referrers: {types}" - del r_ - else: - rep = "repr=" - rep += repr(r)[:100].replace("\n", " ") - # If it's a __closure__, get more information - if rep.startswith(" 0 - del obj - del objs - gc.collect() - assert n == 0, f"\n{n} {cls.__name__} @ {when}:\n" + "\n".join(ref) + + assert_no_instances(cls, when=when) def _resource_path(submodule, filename): diff --git a/mne/utils/tests/test_misc.py b/mne/utils/tests/test_misc.py index c3bc63dea76..64f00512664 100644 --- a/mne/utils/tests/test_misc.py +++ b/mne/utils/tests/test_misc.py @@ -10,7 +10,25 @@ import pytest import mne -from mne.utils import _clean_names, catch_logging, run_subprocess, sizeof_fmt +from mne.utils import ( + _assert_no_instances, + _clean_names, + catch_logging, + run_subprocess, + sizeof_fmt, +) + + +def test_assert_no_instances(): + """Test that our wrapper works.""" + + class _Foo: + pass + + _holder = {"key": _Foo()} + + with pytest.raises(AssertionError, match="after closing"): + _assert_no_instances(_Foo, "after closing") def test_sizeof_fmt(): diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py index 2c5e081b8f0..6b153b1d341 100644 --- a/mne/viz/_3d.py +++ b/mne/viz/_3d.py @@ -51,6 +51,7 @@ _angle_between_quats, _ensure_trans, _find_trans, + _find_vector_rotation, _frame_to_str, _get_trans, _get_transforms_to_coord_frame, @@ -1103,7 +1104,8 @@ def _ch_pos_in_coord_frame(info, to_cf_t, warn_meg=True, verbose=None): # 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",): + is_meg = ch_type in _MEG_CH_TYPES_SPLIT + ("ref_meg",) + if is_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 @@ -1117,11 +1119,15 @@ def _ch_pos_in_coord_frame(info, to_cf_t, warn_meg=True, verbose=None): coil = _create_meg_coils(this_coil, acc="normal", coilset=coilset)[ 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: + # Keep the local (template) geometry and the per-channel + # transform separate (rather than baking into absolute + # vertex positions) so that channels sharing a coil shape + # can share one template/actor when rendered. + local_rr, triangles, extra_z = _sensor_shape(coil) + if len(local_rr) == 0 and warn_meg: warn(f"MEG sensor {info.ch_names[idx]} not found.") + coil_trans = coil_trans.copy() + coil_trans[:3, 3] += coil_trans[:3, :3] @ [0, 0, extra_z] else: ch_coord = info["chs"][idx]["loc"][type_slice] ch_coord_frame = info["chs"][idx]["coord_frame"] @@ -1139,10 +1145,21 @@ def _ch_pos_in_coord_frame(info, to_cf_t, warn_meg=True, verbose=None): 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) + if is_meg: + frame_trans = to_cf_t[_frame_to_str[ch_coord_frame]]["trans"] + transform = frame_trans @ coil_trans + position = transform[:3, 3] + quat = rot_to_quat(transform[:3, :3]) + shape_key = (local_rr.round(10).tobytes(), triangles.tobytes()) + chs[type_name][info.ch_names[idx]] = ( + local_rr, + triangles, + position, + quat, + shape_key, + ) else: + ch_coord = apply_trans(to_cf_t[_frame_to_str[ch_coord_frame]], ch_coord) chs[type_name][info.ch_names[idx]] = ch_coord if unknown_chs: warn( @@ -1227,7 +1244,6 @@ def _plot_axes(renderer, info, to_cf_t, head_mri_t): scale=2e-2, color=ax[1], scale_mode="scalar", - resolution=20, scalars=[0.33, 0.66, 1.0], ) actors.append(actor) @@ -1327,11 +1343,11 @@ def _plot_hpi_coils( ] ) hpi_loc = apply_trans(to_cf_t["head"], hpi_loc) - actor, _ = _plot_glyphs( + return _plot_glyphs( renderer=renderer, loc=hpi_loc, - color=defaults["hpi_color"], - scale=scale, + colors=defaults["hpi_color"], + scales=scale, opacity=opacity, orient_glyphs=orient_glyphs, scale_by_distance=scale_by_distance, @@ -1340,7 +1356,6 @@ def _plot_hpi_coils( check_inside=check_inside, nearest=nearest, ) - return actor def _get_nearest(nearest, check_inside, project_to_trans, proj_rr): @@ -1379,68 +1394,80 @@ def _orient_glyphs( def _plot_glyphs( renderer, loc, - color, - scale, - opacity=1, - mode="cylinder", + colors, + scales, + opacity, + *, orient_glyphs=False, scale_by_distance=False, project_points=False, mark_inside=False, surf=None, + orient_nn=None, + cylinder_geom=None, backface_culling=False, check_inside=None, nearest=None, ): - from matplotlib.colors import ListedColormap, to_rgba + """Plot glyphs (sensors, HPI coils, head-shape points) as one instanced actor. + + ``loc`` is already in render units. ``colors`` (any Matplotlib color spec or + an ``(n, 4)`` / ``(1, 4)`` RGBA array) and ``scales`` may be per-instance or + a single value to broadcast; a single GPU-instanced actor is produced + regardless. When ``surf`` is set (coreg), glyphs are oriented along, scaled + by distance to, projected onto, and/or marked inside the surface. Callers + that have already projected/oriented the points (e.g. projected EEG) can + instead pass per-instance orientation vectors as ``orient_nn`` and, if the + default coreg cylinder is not wanted, a ``cylinder_geom`` dict of + ``_cylinder_geom`` keyword arguments. + """ + from matplotlib.colors import to_rgba, to_rgba_array - _validate_type(mark_inside, bool, "mark_inside") - if surf is not None and len(loc) > 0: - defaults = DEFAULTS["coreg"] - scalars, vectors, proj_pts = _orient_glyphs( + defaults = DEFAULTS["coreg"] + n = len(loc) + if n == 0: + return None + colors = np.array(np.broadcast_to(to_rgba_array(colors), (n, 4)), float) + colors[:, 3] *= opacity + scales = np.broadcast_to(np.asarray(scales, float).reshape(-1), (n,)) + positions = loc + vectors = orient_nn # explicit per-instance orientation, if given + if surf is not None: + scalars, surf_vectors, proj_pts = _orient_glyphs( loc, surf, project_points, mark_inside, check_inside, nearest ) - if mark_inside: - colormap = ListedColormap([to_rgba("darkslategray"), to_rgba(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), + if project_points: + positions = proj_pts + if scale_by_distance: + scales = scales * np.linalg.norm(surf_vectors, axis=1) + if mark_inside: # recolor points that fall inside the surface + colors[scalars < 0.5, :3] = to_rgba("darkslategray")[:3] + if orient_glyphs: # point cylinders along the surface normal + vectors = surf_vectors + kind, template_kw = "sphere", dict() + quats = np.zeros((n, 3)) # identity: unrotated (sphere is orientation-free) + if vectors is not None: # orient cylinders along the given direction + kind = "cylinder" + template_kw = cylinder_geom or dict( + height=defaults["eegp_height"], + center=(0.0, -defaults["eegp_height"], 0.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, ) + x_axis = np.array([1.0, 0.0, 0.0]) + nn = vectors / np.linalg.norm(vectors, axis=1, keepdims=True) + rots = np.array([_find_vector_rotation(x_axis, this_nn) for this_nn in nn]) + quats = rot_to_quat(rots) + rr, tris = renderer._glyph_template(kind, **template_kw) + actor, _ = renderer.instanced_mesh( + rr=rr, + tris=tris, + positions=positions, + quats=quats, + colors=colors, + scales=scales, + backface_culling=backface_culling, + ) + return actor @verbose @@ -1471,11 +1498,11 @@ def _plot_head_shape_points( ) 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( + return _plot_glyphs( renderer=renderer, loc=ext_loc, - color=defaults["extra_color"], - scale=defaults["extra_scale"], + colors=defaults["extra_color"], + scales=defaults["extra_scale"], opacity=opacity, orient_glyphs=orient_glyphs, scale_by_distance=scale_by_distance, @@ -1485,7 +1512,6 @@ def _plot_head_shape_points( check_inside=check_inside, nearest=nearest, ) - return actor def _plot_forward(renderer, fwd, fwd_trans, fwd_scale=1, scale=1.5e-3, alpha=1): @@ -1572,7 +1598,14 @@ def _plot_sensors_3d( plot_sensors = False # plot sensors if isinstance(ch_coord, tuple): # is meg, plot coil - ch_coord = dict(rr=ch_coord[0] * unit_scalar, tris=ch_coord[1]) + local_rr, triangles, position, quat, shape_key = ch_coord + ch_coord = dict( + rr=local_rr * unit_scalar, + tris=triangles, + position=position * unit_scalar, + quat=quat, + shape_key=shape_key, + ) if plot_sensors: if ch_type == "eeg": if "original" in eeg: @@ -1650,113 +1683,73 @@ def _plot_sensors_3d( if isinstance(sens_loc[0], dict): # meg coil if len(colors) == 1: colors = [colors[0]] * len(sens_loc) - for surface, color in zip(sens_loc, colors): - actor, _ = renderer.surface( - surface=surface, - color=color[:3], - opacity=this_alpha * color[3], + colors = np.array(colors, float) + colors[:, 3] *= this_alpha + # Group coils by shape: one GPU-instanced actor per shape + groups = defaultdict(list) + for si, surface in enumerate(sens_loc): + if len(surface["rr"]) == 0: + continue # coil geometry not found (already warned above) + groups[surface["shape_key"]].append(si) + for idxs in groups.values(): + template = sens_loc[idxs[0]] + positions = np.array([sens_loc[i]["position"] for i in idxs]) + quats = np.array([sens_loc[i]["quat"] for i in idxs]) + actor, _ = renderer.instanced_mesh( + rr=template["rr"], + tris=template["tris"], + positions=positions, + quats=quats, + colors=colors[idxs], backface_culling=False, # visible from all sides ) actors[ch_type].append(actor) else: + # One GPU-instanced actor regardless of how many distinct + # colors/scales are requested (broadcasting handles 1-vs-N). sens_loc = np.array(sens_loc, float) mask = ~np.isnan(sens_loc).any(axis=1) - if ch_type == "eegp": # special case, need to project + loc = sens_loc[mask] + these_colors = colors[mask] if len(colors) == len(mask) else colors + these_scales = scales[mask] if len(scales) == len(mask) else scales + orient_nn = cylinder_geom = None + backface_culling = False + actor_key = ch_type + if ch_type == "eegp": + # Projected EEG differs: reproject onto the scalp, orient along + # its normals, and use a single color/scale + flat disc glyph. logger.info("Projecting sensors to the head surface") - eegp_loc, eegp_nn = _project_onto_surface( - sens_loc[mask], head_surf, project_rrs=True, return_nn=True + loc, orient_nn = _project_onto_surface( + loc, head_surf, project_rrs=True, return_nn=True )[2:4] - eegp_loc *= unit_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=colors[0], # TODO: Maybe eventually support multiple - mode="cylinder", - scale=scales[0] * unit_scalar, # TODO: Also someday maybe multiple - opacity=sensor_alpha[ch_type], - glyph_height=defaults["eegp_height"], - glyph_center=(0.0, -defaults["eegp_height"] / 2.0, 0), - glyph_resolution=20, - backface_culling=True, + these_colors = colors[0] # TODO: someday maybe support multiple + these_scales = scales[0] * unit_scalar # TODO: someday maybe multiple + height = defaults["eegp_height"] + cylinder_geom = dict( + radius=0.15, + height=height, + center=(0.0, -height / 2.0, 0.0), + resolution=20, ) - actors["eeg"].append(actor) - elif len(colors) == 1 and len(scales) == 1: - # Single color mode (one actor) - actor, _ = _plot_glyphs( - renderer=renderer, - loc=sens_loc[mask] * unit_scalar, - color=colors[0, :3], - scale=scales[0], - opacity=this_alpha * colors[0, 3], - orient_glyphs=orient_glyphs, - scale_by_distance=scale_by_distance, - project_points=project_points, - surf=surf, - check_inside=check_inside, - nearest=nearest, - ) - actors[ch_type].append(actor) - elif len(colors) == len(sens_loc) and len(scales) == 1: - # Multi-color single scale mode (multiple actors) - for loc, color, usable in zip(sens_loc, colors, mask): - if not usable: - continue - actor, _ = _plot_glyphs( - renderer=renderer, - loc=loc * unit_scalar, - color=color[:3], - scale=scales[0], - opacity=this_alpha * color[3], - orient_glyphs=orient_glyphs, - scale_by_distance=scale_by_distance, - project_points=project_points, - surf=surf, - check_inside=check_inside, - nearest=nearest, - ) - actors[ch_type].append(actor) - elif len(colors) == 1 and len(scales) == len(sens_loc): - # Multi-scale single color mode (multiple actors) - for loc, scale, usable in zip(sens_loc, scales, mask): - if not usable: - continue - actor, _ = _plot_glyphs( - renderer=renderer, - loc=loc * unit_scalar, - color=colors[0, :3], - scale=scale, - opacity=this_alpha * colors[0, 3], - orient_glyphs=orient_glyphs, - scale_by_distance=scale_by_distance, - project_points=project_points, - surf=surf, - check_inside=check_inside, - nearest=nearest, - ) - actors[ch_type].append(actor) - else: - # Multi-color multi-scale mode (multiple actors) - for loc, color, scale, usable in zip(sens_loc, colors, scales, mask): - if not usable: - continue - actor, _ = _plot_glyphs( - renderer=renderer, - loc=loc * unit_scalar, - color=color[:3], - scale=scale, - opacity=this_alpha * color[3], - orient_glyphs=orient_glyphs, - scale_by_distance=scale_by_distance, - project_points=project_points, - surf=surf, - check_inside=check_inside, - nearest=nearest, - ) - actors[ch_type].append(actor) + backface_culling = True + actor_key = "eeg" + actor = _plot_glyphs( + renderer=renderer, + loc=loc * unit_scalar, + colors=these_colors, + scales=these_scales, + opacity=this_alpha, + orient_glyphs=orient_glyphs, + scale_by_distance=scale_by_distance, + project_points=project_points, + surf=surf, + orient_nn=orient_nn, + cylinder_geom=cylinder_geom, + backface_culling=backface_culling, + check_inside=check_inside, + nearest=nearest, + ) + actors[actor_key].append(actor) actors = dict(actors) # get rid of defaultdict @@ -1778,7 +1771,8 @@ def _sensor_shape(coil): except ImportError: # scipy < 1.8 from scipy.spatial.qhull import QhullError id_ = coil["type"] & 0xFFFF - z_value = 0 + add_z_coord = True # almost all geometry is planar + extra_z = 0.0 # Square figure eight if id_ in ( FIFF.FIFFV_COIL_NM_122, @@ -1805,7 +1799,7 @@ def _sensor_shape(coil): (_make_tris_fan(4), _make_tris_fan(4)[:, ::-1] + 4), axis=0 ) # Offset for visibility (using heuristic for sanely named Neuromag coils) - z_value = 0.001 * (1 + coil["chname"].endswith("2")) + extra_z = 0.001 * (1 + coil["chname"].endswith("2")) # Square elif id_ in ( FIFF.FIFFV_COIL_POINT_MAGNETOMETER, @@ -1876,13 +1870,13 @@ def _sensor_shape(coil): rr_rot = rrs @ u tris = Delaunay(rr_rot[:, :2]).simplices tris = np.concatenate((tris, tris[:, ::-1])) - z_value = None + add_z_coord = False # Go from (x,y) -> (x,y,z) - if z_value is not None: - rrs = np.pad(rrs, ((0, 0), (0, 1)), mode="constant", constant_values=z_value) + if add_z_coord: + rrs = np.pad(rrs, ((0, 0), (0, 1)), mode="constant", constant_values=0.0) assert rrs.ndim == 2 and rrs.shape[1] == 3 - return rrs, tris + return rrs, tris, extra_z def _process_clim(clim, colormap, transparent, data=0.0, allow_pos_lims=True): diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py index b66bdf5c436..4d64eacdd5f 100644 --- a/mne/viz/_brain/_brain.py +++ b/mne/viz/_brain/_brain.py @@ -8,6 +8,7 @@ import time import traceback import warnings +import weakref from functools import partial from io import BytesIO @@ -29,7 +30,12 @@ ) from ...defaults import DEFAULTS, _handle_default from ...fixes import _reshape_view -from ...surface import _marching_cubes, _mesh_borders, mesh_edges +from ...surface import ( + _decimate_surface_ico_oct, + _marching_cubes, + _mesh_borders, + mesh_edges, +) from ...transforms import ( Transform, _frame_to_str, @@ -164,9 +170,15 @@ class Brain: %(view_layout)s silhouette : dict | bool As a dict, it contains the ``color``, ``linewidth``, ``alpha`` opacity - and ``decimate`` (level of decimation between 0 and 1 or None) of the - brain's silhouette to display. If True, the default values are used - and if False, no silhouette will be displayed. Defaults to False. + and ``decimate`` of the brain's silhouette to display. ``decimate`` can be + a level of decimation between 0 and 1, None for no decimation, or a spacing + string (e.g. ``"ico4"``, ``"oct6"``) to instead pick vertices the same way + :func:`mne.setup_source_space` does, which is faster and better preserves + the shape than plain decimation. If True, ``"ico5"`` is used and + if False, no silhouette will be displayed. Defaults to False. + + .. versionchanged:: 1.13 + The default ``decimate`` value changed from ``0.9`` to ``"ico5"``. %(theme_3d)s show : bool Display the window as soon as it is ready. Defaults to True. @@ -275,6 +287,11 @@ class Brain: +-------------------------------------+--------------+---------------+ """ + # Tracks live instances when MNE_3D_BACKEND_TESTING is set (see + # __init__), so that tests can check for lingering instances without + # having to scan gc.get_objects() for the whole process. + _instances = weakref.WeakSet() + def __init__( self, subject, @@ -298,6 +315,7 @@ def __init__( theme=None, show=True, ): + from ..backends import renderer as _renderer_mod from ..backends.renderer import _get_renderer, backend _validate_type(subject, str, "subject") @@ -346,6 +364,10 @@ def __init__( self.time_viewer = False self._hash = time.time_ns() + # This is only used in testing! Needs self._hash to already be set, + # since Brain.__hash__ (used when adding to the WeakSet) depends on it. + if _renderer_mod.MNE_3D_BACKEND_TESTING: + Brain._instances.add(self) self._hemi = hemi self._units = units self._alpha = float(alpha) @@ -366,7 +388,7 @@ def __init__( "color": self._bg_color, "line_width": 2, "alpha": alpha, - "decimate": 0.9, + "decimate": "ico5", } _validate_type(silhouette, (dict, bool), "silhouette") if isinstance(silhouette, dict): @@ -403,7 +425,6 @@ def __init__( self._renderer._window_set_theme(theme) self.plotter = self._renderer.plotter self.widgets = dict() - self._setup_canonical_rotation() # plot hemis @@ -450,12 +471,30 @@ def __init__( self._renderer.plotter.add_actor(actor, render=False) if self.silhouette: mesh = self.layered_meshes[h] + decimate = self._silhouette["decimate"] + if isinstance(decimate, str): + import pyvista as pv + + vertno, tris = _decimate_surface_ico_oct( + self._subject, + self._subjects_dir, + h, + self.geo[h].surf, + decimate, + ) + sil_mesh = pv.PolyData( + self.geo[h].coords[vertno], + np.c_[np.full(len(tris), 3), tris], + ) + decimate = None # already decimated + else: + sil_mesh = mesh._polydata self._renderer._silhouette( - mesh=mesh._polydata, + mesh=sil_mesh, color=self._silhouette["color"], line_width=self._silhouette["line_width"], alpha=self._silhouette["alpha"], - decimate=self._silhouette["decimate"], + decimate=decimate, ) self._set_camera(**views_dicts[h][v]) @@ -1426,7 +1465,13 @@ def _add_vertex_glyph(self, hemi, mesh, vertex_id, update=True): def _remove_vertex_glyph(self, *, hemi, vertex_id, render=True): _ensure_int(vertex_id) assert isinstance(hemi, str), f"got {type(hemi)} for {hemi=}" - spheres = self._picked_points.pop((hemi, vertex_id)) + # When linked via _LinkViewer, removing a vertex on one brain cascades + # to all linked brains, so by the time a given brain's own loop (e.g. + # in clear_glyphs) reaches this (hemi, vertex_id) it may already be + # gone; just no-op in that case. + spheres = self._picked_points.pop((hemi, vertex_id), None) + if spheres is None: + return color, line = spheres[0]["color"], spheres[0]["line"] line.remove() self.mpl_canvas.update_plot() @@ -3023,7 +3068,12 @@ def add_annotation( ) if hover: - obs = self.plotter.AddObserver("MouseMoveEvent", self._on_annotation_hover) + + @_auto_weakref + def on_annotation_hover(iren, event): + self._on_annotation_hover(iren, event) + + obs = self.plotter.AddObserver("MouseMoveEvent", on_annotation_hover) for hemi in hemis: caption = self._create_caption() self._annots[hemi][-1].update(caption=caption, obs=obs) diff --git a/mne/viz/_brain/tests/test_brain.py b/mne/viz/_brain/tests/test_brain.py index d50495b14ab..05ded04fbb6 100644 --- a/mne/viz/_brain/tests/test_brain.py +++ b/mne/viz/_brain/tests/test_brain.py @@ -638,7 +638,7 @@ def test_add_sensors_scales(renderer_interactive_pyvistaqt): title=title, cortex=cortex, units="m", - silhouette=dict(decimate=0.95), + silhouette=dict(decimate="oct6"), **kwargs, ) @@ -957,7 +957,7 @@ def test_brain_time_viewer(renderer_interactive_pyvistaqt, pixel_ratio, brain_gc brain = _create_testing_brain( hemi="both", show_traces=False, - brain_kwargs=dict(silhouette=dict(decimate=0.95)), + brain_kwargs=dict(silhouette=dict(decimate="ico5")), ) # test sub routines when show_traces=False brain._on_pick(None, None) diff --git a/mne/viz/backends/_abstract.py b/mne/viz/backends/_abstract.py index 7c9d97167aa..9a4d764eb7e 100644 --- a/mne/viz/backends/_abstract.py +++ b/mne/viz/backends/_abstract.py @@ -97,7 +97,7 @@ def set_interaction(self, interaction): @classmethod @abstractmethod - def legend(self, labels, border=False, size=0.1, face="triangle", loc="upper left"): + def legend(self, labels, size=0.1, face="triangle", loc="upper left"): """Add a legend to the scene. Parameters @@ -106,9 +106,6 @@ def legend(self, labels, border=False, size=0.1, face="triangle", loc="upper lef Each entry must contain two strings, (label, color), where ``label`` is the name of the item to add, and ``color`` is the color of the label to add. - border : bool - Controls if there will be a border around the legend. - The default is False. size : float The size of the entire figure window. loc : str @@ -144,7 +141,6 @@ def mesh( representation="surface", line_width=1.0, normals=None, - polygon_offset=None, name=None, **kwargs, ): @@ -190,8 +186,6 @@ def mesh( The width of the lines when representation='wireframe'. normals : array, shape (n_vertices, 3) The array containing the normal of each vertex. - polygon_offset : float - If not None, the factor used to resolve coincident topology. name : str | None The name of the mesh. kwargs : args @@ -266,7 +260,6 @@ def surface( normalized_colormap=False, scalars=None, backface_culling=False, - polygon_offset=None, *, name=None, ): @@ -294,8 +287,6 @@ def surface( The scalar valued associated to the vertices. backface_culling : bool If True, enable backface culling on the surface. - polygon_offset : float - If not None, the factor used to resolve coincident topology. name : str | None Name of the surface. """ @@ -407,7 +398,6 @@ def quiver3d( color, scale, mode, - resolution=8, glyph_height=None, glyph_center=None, glyph_resolution=None, @@ -447,10 +437,6 @@ def quiver3d( The given value specifies the maximum glyph size in drawing units. mode : 'arrow', 'cone' or 'cylinder' The type of the quiver. - resolution : int - The resolution of the glyph created. Depending on the type of - glyph, it represents the number of divisions in its geometric - representation. glyph_height : float The height of the glyph used with the quiver. glyph_center : tuple @@ -485,6 +471,73 @@ def quiver3d( """ pass + @classmethod + @abstractmethod + def instanced_mesh( + self, + rr, + tris, + positions, + quats, + colors, + scales=None, + opacity=1.0, + backface_culling=False, + name=None, + ): + """Add one mesh GPU-instanced at multiple positions/orientations. + + Unlike :meth:`quiver3d` or :meth:`sphere`, which bake a merged, + static glyph mesh, this draws one copy of a single template mesh at + each of ``n_instances`` locations using per-instance position, + orientation, and color that can be changed later without rebuilding + the geometry (see ``mesh`` in "Returns" below). + + Parameters + ---------- + rr : array, shape (n_vertices, 3) + The vertices of the template mesh, in the local (object) frame + shared by all instances. + tris : array, shape (n_tris, 3) + The triangles of the template mesh. + positions : array, shape (n_instances, 3) + The position of each instance. + quats : array, shape (n_instances, 3) + The orientation of each instance as a unit quaternion, given as + only the ``(x, y, z)`` vector part (``w`` omitted, recoverable + as ``sqrt(max(1 - x**2 - y**2 - z**2, 0))``) -- the same + convention used by :func:`~mne.transforms.rot_to_quat` / + :func:`~mne.transforms.quat_to_rot`. + colors : array, shape (n_instances, 4) + The per-instance RGBA color (float values between 0 and 1) to + use for each instance. Per-instance alpha (the fourth column) + is respected. + scales : array, shape (n_instances,) | None + The per-instance isotropic scale factor applied to the template + geometry. If ``None`` (the default), no per-instance scaling is + applied and the size baked into ``rr`` is used as-is (e.g. for + MEG coils). + opacity : float + A uniform opacity multiplier applied on top of the per-instance + alpha values in ``colors``. + backface_culling : bool + If True, enable backface culling on the mesh. + name : str | None + Name of the mesh. + + Returns + ------- + actor : + The actor in the scene. + mesh : + The point cloud (one point per instance) backing the glyph + mapper. Its per-instance color and orientation arrays can be + mutated in place (followed by a render update) to recolor or + re-orient individual instances live, without rebuilding the + actor or its geometry. + """ + pass + @classmethod @abstractmethod def text2d( diff --git a/mne/viz/backends/_pyvista.py b/mne/viz/backends/_pyvista.py index ca61d601cf7..d020f468e5a 100644 --- a/mne/viz/backends/_pyvista.py +++ b/mne/viz/backends/_pyvista.py @@ -17,45 +17,20 @@ import numpy as np import pyvista -from pyvista import Line, Plotter, PolyData, UnstructuredGrid, close_all +from pyvista import Line, Plotter, PolyData, close_all +from pyvista.plotting.plotter import _ALL_PLOTTERS from pyvistaqt import BackgroundPlotter - -from ...fixes import _compare_version -from ...surface import _vtk_smooth -from ...transforms import _cart_to_sph, _sph_to_cart, apply_trans -from ...utils import ( - _check_option, - _require_version, - _validate_type, - warn, -) -from ._abstract import Figure3D, _AbstractRenderer -from ._utils import ( - ALLOWED_QUIVER_MODES, - _alpha_blend_background, - _get_colormap_from_array, - _init_mne_qtapp, -) - -try: - from pyvista.plotting.plotter import _ALL_PLOTTERS -except Exception: # PV < 0.40 - from pyvista.plotting.plotting import _ALL_PLOTTERS - from vtkmodules.util.numpy_support import numpy_to_vtk from vtkmodules.vtkCommonCore import VTK_UNSIGNED_CHAR, vtkCommand, vtkLookupTable -from vtkmodules.vtkCommonDataModel import VTK_VERTEX, vtkPiecewiseFunction +from vtkmodules.vtkCommonDataModel import vtkPiecewiseFunction from vtkmodules.vtkCommonTransforms import vtkTransform from vtkmodules.vtkFiltersCore import vtkGlyph3D from vtkmodules.vtkFiltersGeneral import vtkMarchingContourFilter from vtkmodules.vtkFiltersHybrid import vtkPolyDataSilhouette from vtkmodules.vtkFiltersSources import ( vtkArrowSource, - vtkConeSource, vtkCylinderSource, vtkGlyphSource2D, - vtkPlatonicSolidSource, - vtkSphereSource, ) from vtkmodules.vtkImagingCore import vtkImageReslice from vtkmodules.vtkRenderingCore import ( @@ -64,12 +39,30 @@ vtkColorTransferFunction, vtkCoordinate, vtkDataSetMapper, + vtkGlyph3DMapper, vtkMapper, vtkPolyDataMapper, vtkVolume, ) from vtkmodules.vtkRenderingVolumeOpenGL2 import vtkSmartVolumeMapper +from ...fixes import _compare_version +from ...surface import _vtk_smooth +from ...transforms import _cart_to_sph, _sph_to_cart, apply_trans +from ...utils import ( + _check_option, + _require_version, + _validate_type, + warn, +) +from ._abstract import Figure3D, _AbstractRenderer +from ._utils import ( + ALLOWED_QUIVER_MODES, + _alpha_blend_background, + _get_colormap_from_array, + _init_mne_qtapp, +) + try: from vtkmodules.vtkFiltersGeneral import vtkTransformFilter except ImportError: # TODO VERSION VTK 9.7+ @@ -314,22 +307,18 @@ def scene(self): def update_lighting(self): # Inspired from Mayavi's version of Raymond Maple 3-lights illumination + # below and centered, left and above, right and above + az_el_in = ((0, -45, 0.7), (-60, 30, 0.7), (60, 30, 0.7)) for renderer in self._all_renderers: - lights = list(renderer.GetLights()) - headlight = lights.pop(0) - headlight.SetSwitch(False) - # below and centered, left and above, right and above - az_el_in = ((0, -45, 0.7), (-60, 30, 0.7), (60, 30, 0.7)) - for li, light in enumerate(lights): - if li < len(az_el_in): - light.SetSwitch(True) - light.SetPosition(_to_pos(*az_el_in[li][:2])) - light.SetIntensity(az_el_in[li][2]) - else: - light.SetSwitch(False) - light.SetPosition(_to_pos(0.0, 0.0)) - light.SetIntensity(0.0) - light.SetColor(1.0, 1.0, 1.0) + renderer.remove_all_lights() + for azimuth, elevation, intensity in az_el_in: + light = pyvista.Light( + position=_to_pos(azimuth, elevation), + color="white", + light_type="camera light", + intensity=intensity, + ) + renderer.add_light(light) def set_interaction(self, interaction): if not hasattr(self.plotter, "iren") or self.plotter.iren is None: @@ -346,7 +335,7 @@ def set_interaction(self, interaction): kwargs["mouse_wheel_zooms"] = True getattr(self.plotter, f"enable_{interaction}_style")(**kwargs) - def legend(self, labels, border=False, size=0.1, face="triangle", loc="upper left"): + def legend(self, labels, size=0.1, face="triangle", loc="upper left"): return self.plotter.add_legend(labels, size=(size, size), face=face, loc=loc) def polydata( @@ -363,7 +352,6 @@ def polydata( interpolate_before_map=True, representation="surface", line_width=1.0, - polygon_offset=None, *, name=None, **kwargs, @@ -413,13 +401,6 @@ def polydata( **kwargs, ) - if polygon_offset is not None: - mapper = actor.GetMapper() - mapper.SetResolveCoincidentTopologyToPolygonOffset() - mapper.SetRelativeCoincidentTopologyPolygonOffsetParameters( - polygon_offset, polygon_offset - ) - return actor, mesh def mesh( @@ -440,7 +421,6 @@ def mesh( representation="surface", line_width=1.0, normals=None, - polygon_offset=None, name=None, **kwargs, ): @@ -460,7 +440,6 @@ def mesh( interpolate_before_map=interpolate_before_map, representation=representation, line_width=line_width, - polygon_offset=polygon_offset, name=name, **kwargs, ) @@ -516,7 +495,6 @@ def surface( normalized_colormap=False, scalars=None, backface_culling=False, - polygon_offset=None, *, name=None, ): @@ -538,7 +516,6 @@ def surface( colormap=colormap, vmin=vmin, vmax=vmax, - polygon_offset=polygon_offset, name=name, ) @@ -552,21 +529,17 @@ def sphere( backface_culling=False, radius=None, ): - from vtkmodules.vtkFiltersSources import vtkSphereSource - factor = 1.0 if radius is not None else scale center = np.array(center, dtype=float) if len(center) == 0: return None, None _check_option("center.ndim", center.ndim, (1, 2)) _check_option("center.shape[-1]", center.shape[-1], (3,)) - sphere = vtkSphereSource() - sphere.SetThetaResolution(resolution) - sphere.SetPhiResolution(resolution) - if radius is not None: - sphere.SetRadius(radius) - sphere.Update() - geom = sphere.GetOutput() + geom = pyvista.Sphere( + radius=0.5 if radius is None else radius, + theta_resolution=resolution, + phi_resolution=resolution, + ) mesh = PolyData(center) glyph = mesh.glyph(orient=False, scale=False, factor=factor, geom=geom) actor = _add_mesh( @@ -628,7 +601,6 @@ def quiver3d( color, scale, mode, - resolution=8, *, glyph_height=None, glyph_center=None, @@ -648,12 +620,9 @@ def quiver3d( _check_option("scale_mode", scale_mode, list(scale_map)) factor = scale vectors = np.c_[u, v, w] - points = np.vstack(np.c_[x, y, z]) + points = np.vstack(np.c_[x, y, z]).astype(float) n_points = len(points) - cell_type = np.full(n_points, VTK_VERTEX) - cells = np.c_[np.full(n_points, 1), range(n_points)] - args = (cells, cell_type, points) - grid = UnstructuredGrid(*args) + grid = PolyData(points) if scalars is None: scalars = np.ones((n_points,)) mesh_scalars = None @@ -667,45 +636,25 @@ def quiver3d( alg = _glyph(grid, orient="vec", scalars="scalars", factor=factor) mesh = pyvista.wrap(alg.GetOutput()) else: - tr = None if mode == "cone": - glyph = vtkConeSource() - glyph.SetCenter(0.5, 0, 0) - if glyph_radius is not None: - glyph.SetRadius(glyph_radius) + geom = pyvista.Cone(center=(0.5, 0, 0), radius=glyph_radius) elif mode == "cylinder": - glyph = vtkCylinderSource() - if glyph_radius is not None: - glyph.SetRadius(glyph_radius) - elif mode == "oct": - glyph = vtkPlatonicSolidSource() - glyph.SetSolidTypeToOctahedron() - else: - assert mode == "sphere", mode # guaranteed above - glyph = vtkSphereSource() - if mode == "cylinder": - if glyph_height is not None: - glyph.SetHeight(glyph_height) - if glyph_center is not None: - glyph.SetCenter(glyph_center) - if glyph_resolution is not None: - glyph.SetResolution(glyph_resolution) - tr = vtkTransform() - tr.RotateWXYZ(90, 0, 0, 1) + geom = _cylinder_geom( + radius=glyph_radius, + height=glyph_height, + center=glyph_center, + resolution=glyph_resolution, + ) elif mode == "oct": + geom = pyvista.PlatonicSolid(kind="octahedron") if solid_transform is not None: assert solid_transform.shape == (4, 4) - tr = vtkTransform() - tr.SetMatrix(solid_transform.astype(np.float64).ravel()) - if tr is not None: - # fix orientation - glyph.Update() - trp = vtkTransformFilter() - trp.SetInputData(glyph.GetOutput()) - trp.SetTransform(tr) - glyph = trp - glyph.Update() - geom = glyph.GetOutput() + geom = geom.transform( + solid_transform.astype(np.float64), inplace=True + ) + else: + assert mode == "sphere", mode # guaranteed above + geom = pyvista.Sphere(theta_resolution=8, phi_resolution=8) mesh = grid.glyph( orient="vec", scale=scale_map[scale_mode], @@ -725,6 +674,78 @@ def quiver3d( ) return actor, mesh + # quiver3d (above) and instanced_mesh (below) split along principled lines: + # quiver3d bakes a static, merged glyph mesh with direction-vector + # orientation and scalar/colormap coloring (arrows and friends), while + # instanced_mesh GPU-instances one template with per-instance, updatable + # quaternion orientation and RGBA coloring (sensors, MEG coils). + def instanced_mesh( + self, + rr, + tris, + positions, + quats, + colors, + scales=None, + opacity=1.0, + backface_culling=False, + *, + name=None, + ): + faces = np.c_[np.full(len(tris), 3), tris] + geom = PolyData(np.asarray(rr, float), faces) + _compute_normals(geom) + + cloud = PolyData(np.asarray(positions, float).copy()) + cloud.point_data["orientation"] = _quat_to_vtk_wxyz(np.asarray(quats, float)) + cloud.point_data["colors"] = (np.asarray(colors, float) * 255).astype(np.uint8) + + mapper = vtkGlyph3DMapper() + mapper.SetInputData(cloud) + mapper.SetSourceData(geom) + mapper.SetOrientationArray("orientation") + mapper.SetOrientationModeToQuaternion() + if scales is None: + # size is baked into rr (e.g. MEG coils); no per-instance scaling + mapper.ScalingOff() + else: + cloud.point_data["scale"] = np.asarray(scales, float) + mapper.SetScaleArray("scale") + mapper.SetScaleModeToScaleByMagnitude() + mapper.SetScaleFactor(1.0) + mapper.ScalingOn() + mapper.SetScalarModeToUsePointFieldData() + mapper.SelectColorArray("colors") + mapper.SetColorModeToDirectScalars() + + actor = self._actor(mapper) + prop = actor.GetProperty() + prop.SetOpacity(opacity) + prop.SetBackfaceCulling(backface_culling) + if self.smooth_shading and "Normals" in geom.point_data: + prop.SetInterpolationToPhong() + self.plotter.add_actor( + actor, name=name, render=False, reset_camera=False, pickable=True + ) + return actor, cloud + + def _glyph_template(self, kind, **kwargs): + """Return (rr, tris) for a standard template mesh for instanced_mesh. + + ``kind`` is ``"sphere"`` (unit-diameter, i.e. radius 0.5) or + ``"cylinder"`` (see ``_cylinder_geom`` for ``**kwargs``). The template + is oriented along +x so per-instance quaternions can point it anywhere. + """ + if kind == "sphere": + geom = pyvista.Sphere(radius=0.5, theta_resolution=8, phi_resolution=8) + else: + assert kind == "cylinder", kind + geom = pyvista.wrap(_cylinder_geom(**kwargs)) + geom = geom.triangulate() + rr = np.asarray(geom.points, float) + tris = np.asarray(geom.faces).reshape(-1, 4)[:, 1:] + return rr, tris + def text2d( self, x_window, @@ -761,7 +782,7 @@ def text2d( return actor def text3d(self, x, y, z, text, scale, color="white"): - kwargs = dict( + actor = self.plotter.add_point_labels( points=np.array([x, y, z]).astype(float), labels=[text], point_size=scale, @@ -769,10 +790,8 @@ def text3d(self, x, y, z, text, scale, color="white"): font_family=self.font_family, name=text, shape_opacity=0, + always_visible=True, ) - if "always_visible" in signature(self.plotter.add_point_labels).parameters: - kwargs["always_visible"] = True - actor = self.plotter.add_point_labels(**kwargs) _hide_testing_actor(actor) return actor @@ -946,15 +965,9 @@ def _set_volume_range(self, volume, ctable, alpha, scalar_bar, rng): scalar_bar.SetLookupTable(lut) def _sphere(self, center, color, radius): - from vtkmodules.vtkFiltersSources import vtkSphereSource - - sphere = vtkSphereSource() - sphere.SetThetaResolution(8) - sphere.SetPhiResolution(8) - sphere.SetRadius(radius) - sphere.SetCenter(center) - sphere.Update() - mesh = pyvista.wrap(sphere.GetOutput()) + mesh = pyvista.Sphere( + radius=radius, center=center, theta_resolution=8, phi_resolution=8 + ) actor = _add_mesh(self.plotter, mesh=mesh, color=color) return actor, mesh @@ -1084,6 +1097,13 @@ def _compute_normals(mesh): ) +def _quat_to_vtk_wxyz(quat): + """Convert MNE's (..., 3) unit quaternions to VTK-order (w, x, y, z).""" + assert quat.ndim >= 2 and quat.shape[-1] == 3, quat.shape + w = np.sqrt(np.clip(1.0 - (quat * quat).sum(-1), 0.0, 1.0)) + return np.concatenate([w[..., np.newaxis], quat], axis=-1) + + def _add_mesh(plotter, **kwargs): """Patch PyVista add_mesh.""" mesh = kwargs.get("mesh") @@ -1268,7 +1288,33 @@ def _process_events(plotter): def _add_camera_callback(camera, callback): - camera.AddObserver(vtkCommand.ModifiedEvent, callback) + return camera.AddObserver(vtkCommand.ModifiedEvent, callback) + + +def _cylinder_geom(radius=None, height=None, center=None, resolution=None): + """Build a cylinder vtkPolyData with its axis along +x. + + vtkCylinderSource's axis is along y, so we rotate 90 degrees about z to + match the arrow/cone convention of pointing along x (the axis that + vtkGlyph3D/vtkGlyph3DMapper orient toward the per-instance vector). + """ + source = vtkCylinderSource() + if radius is not None: + source.SetRadius(radius) + if height is not None: + source.SetHeight(height) + if center is not None: + source.SetCenter(center) + if resolution is not None: + source.SetResolution(resolution) + source.Update() + tr = vtkTransform() + tr.RotateWXYZ(90, 0, 0, 1) + trp = vtkTransformFilter() + trp.SetInputData(source.GetOutput()) + trp.SetTransform(tr) + trp.Update() + return trp.GetOutput() def _arrow_glyph(grid, factor): @@ -1339,14 +1385,7 @@ def _glyph( @contextmanager def _disabled_depth_peeling(): - try: - from pyvista import global_theme - except Exception: # workaround for older PyVista - from pyvista import rcParams - - depth_peeling = rcParams["depth_peeling"] - else: - depth_peeling = global_theme.depth_peeling + depth_peeling = pyvista.global_theme.depth_peeling depth_peeling_enabled = depth_peeling["enabled"] depth_peeling["enabled"] = False try: diff --git a/mne/viz/backends/tests/test_renderer.py b/mne/viz/backends/tests/test_renderer.py index 399aced2cfd..661eacf9a2d 100644 --- a/mne/viz/backends/tests/test_renderer.py +++ b/mne/viz/backends/tests/test_renderer.py @@ -10,6 +10,7 @@ import pytest from matplotlib.font_manager import findfont +from mne.transforms import rot_to_quat from mne.utils import run_subprocess from mne.viz import Figure3D, get_3d_backend, set_3d_backend from mne.viz.backends._utils import ALLOWED_QUIVER_MODES @@ -150,6 +151,22 @@ def test_3d_backend(renderer): with pytest.raises(ValueError, match="Invalid value"): rend.quiver3d(mode="foo", **kwargs) + # use instanced_mesh + inst_positions = np.array([[0.0, 0.0, 0.0], [tet_size, 0.0, 0.0]]) + inst_quats = np.array([rot_to_quat(np.eye(3)), rot_to_quat(np.eye(3))]) + inst_colors = np.array([[1.0, 0.0, 0.0, 1.0], [0.0, 1.0, 0.0, 1.0]]) + _, inst_cloud = rend.instanced_mesh( + rr=sph_center * sph_scale, + tris=tet_indices, + positions=inst_positions, + quats=inst_quats, + colors=inst_colors, + ) + # colors can be updated in place (e.g. for future sensor + # highlighting/hover) without rebuilding the actor or its geometry + inst_cloud.point_data["colors"][0] = [0, 0, 255, 255] + inst_cloud.Modified() + # use tube rend.tube(origin=np.array([[0, 0, 0]]), destination=np.array([[0, 1, 0]])) _, tube = rend.tube( diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 154313a1eeb..4e6cfce4504 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -485,13 +485,27 @@ def test_plot_alignment_meg(renderer, system): sensor_colors=sensor_colors, ) assert isinstance(fig, Figure3D) - # count the number of objects: should be n_meg_ch + 1 (helmet) + 1 (head) - use_info = pick_info( - this_info, - pick_types(this_info, meg=True, eeg=False, ref_meg="ref" in meg, exclude=()), - ) - n_actors = use_info["nchan"] + 2 - _assert_n_actors(fig, renderer, n_actors) + # one actor per distinct coil shape + helmet + head + # Neuromag has 2 shapes, CTF and BTi expose only mags, KIT's mag + ref mag are same + n_shapes = 1 if system in ("CTF", "BTi") else 2 + _assert_n_actors(fig, renderer, n_shapes + 2) + + if system == "Neuromag": + # Passing fully-random per-channel colors should not blow up the + # actor count (proves per-channel MEG coil coloring no longer + # requires one actor per unique color). + rng = np.random.default_rng(0) + n_meg = len(pick_types(this_info, meg=True)) + fig2 = plot_alignment( + this_info, + read_trans(trans_fname), + subject="sample", + subjects_dir=subjects_dir, + meg=meg, + eeg=False, + sensor_colors=dict(meg=rng.random((n_meg, 4))), + ) + _assert_n_actors(fig2, renderer, n_shapes + 2) @testing.requires_testing_data diff --git a/mne/viz/tests/test_raw.py b/mne/viz/tests/test_raw.py index dab3deb5165..b4adce79f47 100644 --- a/mne/viz/tests/test_raw.py +++ b/mne/viz/tests/test_raw.py @@ -13,6 +13,7 @@ import pytest from matplotlib import backend_bases, rc_context from numpy.testing import assert_allclose, assert_array_equal +from refleak.testing import assert_no_instances import mne from mne import Annotations, create_info, pick_types @@ -21,7 +22,6 @@ from mne.datasets import testing from mne.io import RawArray from mne.utils import ( - _assert_no_instances, _dt_to_stamp, _record_warnings, catch_logging, @@ -1357,7 +1357,7 @@ def test_plotting_memory_garbage_collection(raw, pg_backend): from mne_qt_browser._pg_figure import MNEQtBrowser assert len(mne_qt_browser._browser_instances) == 0 - _assert_no_instances(MNEQtBrowser, "after closing") + assert_no_instances(MNEQtBrowser, "after closing") def test_plotting_scalebars(browser_backend, qtbot): diff --git a/pyproject.toml b/pyproject.toml index 64e0a55fc5c..f84e22f8ced 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ doc = [ "pyvistaqt >= 0.11", # released 2023-06-30, will become 0.12 on 2028-07-01 "pyxdf", "pyzmq != 24.0.0", + "refleak", "scikit-learn >= 1.5", # released 2024-05-21, will become 1.6 on 2026-12-09 "seaborn >= 0.5, != 0.11.2", "selenium >= 4.27.1", @@ -53,6 +54,7 @@ test = [ "pytest-qt >= 4.3", "pytest-rerunfailures", "pytest-timeout >= 2.2", + "refleak", "ruff >= 0.1", "twine", "vulture", diff --git a/tools/pylock.ci-old.toml b/tools/pylock.ci-old.toml index 8980f0beb51..ce3d975e03e 100644 --- a/tools/pylock.ci-old.toml +++ b/tools/pylock.ci-old.toml @@ -1,5 +1,5 @@ # This file was autogenerated by uv via the following command: -# uv pip compile pyproject.toml --python 3.10 --python-platform x86_64-unknown-linux-gnu --group test --group lockfile_extras --resolution lowest-direct --format pylock.toml --output-file tools/pylock.ci-old.toml +# uv pip compile pyproject.toml --python 3.10 --python-platform x86_64-unknown-linux-gnu --group test --group lockfile_extras --resolution lowest-direct --format pylock.toml --output-file tools/pylock.ci-old.toml --no-cache lock-version = "1.0" created-by = "uv" requires-python = ">=3.10" @@ -428,6 +428,12 @@ version = "0.37.0" sdist = { url = "https://files.pythonhosted.org/packages/22/f5/df4e9027acead3ecc63e50fe1e36aca1523e1719559c499951bb4b53188f/referencing-0.37.0.tar.gz", upload-time = 2025-10-13T15:30:48Z, size = 78036, hashes = { sha256 = "44aefc3142c5b842538163acb373e24cce6632bd54bdb01b21ad5863489f50d8" } } wheels = [{ url = "https://files.pythonhosted.org/packages/2c/58/ca301544e1fa93ed4f80d724bf5b194f6e4b945841c5bfd555878eea9fcb/referencing-0.37.0-py3-none-any.whl", upload-time = 2025-10-13T15:30:47Z, size = 26766, hashes = { sha256 = "381329a9f99628c9069361716891d34ad94af76e461dcb0335825aecc7692231" } }] +[[packages]] +name = "refleak" +version = "0.1.1" +sdist = { url = "https://files.pythonhosted.org/packages/48/3d/b1fe0b6bd9fccf9666c40beafff4e042c5f83eaa54d5408ccc5559e782be/refleak-0.1.1.tar.gz", upload-time = 2026-07-04T00:57:10Z, size = 17584, hashes = { sha256 = "05c7ec81497c6f0dd80028fb933136a7d9b055a978f74b17f8bc787e1602b764" } } +wheels = [{ url = "https://files.pythonhosted.org/packages/50/fe/b26d9cad0eb21bbbaf345b9543eabc7c28f093ed022a4df9b3027df54bcd/refleak-0.1.1-py3-none-any.whl", upload-time = 2026-07-04T00:57:09Z, size = 11937, hashes = { sha256 = "aafe8c1141b59df7c10eaa9bdf509ddefcf00624be644118adb678c4af30c38f" } }] + [[packages]] name = "requests" version = "2.32.5" diff --git a/tutorials/inverse/50_beamformer_lcmv.py b/tutorials/inverse/50_beamformer_lcmv.py index bd7228fa840..0add1434fd2 100644 --- a/tutorials/inverse/50_beamformer_lcmv.py +++ b/tutorials/inverse/50_beamformer_lcmv.py @@ -23,11 +23,10 @@ # %% # Introduction to beamformers # --------------------------- -# A beamformer is a spatial filter that reconstructs source activity by -# scanning through a grid of pre-defined source points and estimating activity -# at each of those source points independently. A set of weights is -# constructed for each defined source location which defines the contribution -# of each sensor to this source. +# A beamformer is a spatial filter that reconstructs source activity by scanning through +# a grid of pre-defined source points and estimating activity at each of those source +# points independently. A set of weights is constructed for each defined source location +# which defines the contribution of each sensor to this source. # # Beamformers are often used for their focal reconstructions and their ability # to reconstruct deeper sources. They can also suppress external noise sources. diff --git a/tutorials/preprocessing/25_background_filtering.py b/tutorials/preprocessing/25_background_filtering.py index 537da6a22b2..bbb00ad18cc 100644 --- a/tutorials/preprocessing/25_background_filtering.py +++ b/tutorials/preprocessing/25_background_filtering.py @@ -150,7 +150,6 @@ from scipy import signal import mne -from mne.fixes import minimum_phase from mne.time_frequency.tfr import morlet from mne.viz import plot_filter, plot_ideal_filter @@ -323,11 +322,9 @@ # # We can construct a minimum-phase filter from our existing linear-phase # filter, and note that the falloff is not as steep. Here we do this with function -# ``mne.fixes.minimum_phase()`` to avoid a SciPy bug; once SciPy 1.14.0 is released you -# could directly use # :func:`scipy.signal.minimum_phase(..., half=False) `. -h_min = minimum_phase(h, half=False) +h_min = signal.minimum_phase(h, half=False) plot_filter(h_min, sfreq, freq, gain, "Minimum-phase", **kwargs) # %% diff --git a/tutorials/preprocessing/70_fnirs_processing.py b/tutorials/preprocessing/70_fnirs_processing.py index afafa9c7bdd..efae1fe04ed 100644 --- a/tutorials/preprocessing/70_fnirs_processing.py +++ b/tutorials/preprocessing/70_fnirs_processing.py @@ -18,8 +18,6 @@ # %% -from itertools import compress - import matplotlib.pyplot as plt import numpy as np @@ -118,7 +116,7 @@ # In this example we will mark all channels with a SCI less than 0.5 as bad # (this dataset is quite clean, so no channels are marked as bad). -raw_od.info["bads"] = list(compress(raw_od.ch_names, sci < 0.5)) +raw_od.info["bads"] = [raw_od.ch_names[idx] for idx in np.where(sci < 0.5)[0]] # %% # At this stage it is appropriate to inspect your data