From 6afcfb12c968b61804f3c4f4f5e9852b96adb64c Mon Sep 17 00:00:00 2001 From: payam Date: Fri, 23 Oct 2026 14:28:33 +0200 Subject: [PATCH 1/7] add subcortical surf data --- mne/viz/_brain/_brain.py | 55 +++++++++++++++++++++++++++++----- mne/viz/_brain/_time_viewer.py | 5 ++-- 2 files changed, 51 insertions(+), 9 deletions(-) diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py index 760066ab1bf..4c8d61ce031 100644 --- a/mne/viz/_brain/_brain.py +++ b/mne/viz/_brain/_brain.py @@ -1167,7 +1167,7 @@ def _add_volume_data(self, hemi, src, volume_options): from ..backends._pyvista import _hide_testing_actor _validate_type(src, SourceSpaces, "src") - _check_option("src.kind", src.kind, ("volume",)) + _check_option("src.kind", src.kind, ("volume", "subcortical_surf")) _validate_type(volume_options, (dict, "numeric", None), "volume_options") assert hemi == "vol" if not isinstance(volume_options, dict): @@ -1203,6 +1203,9 @@ def _add_volume_data(self, hemi, src, volume_options): if alpha is None: alpha = 0.4 if self._data[hemi]["array"].ndim == 3 else 1.0 alpha = np.clip(float(alpha), 0.0, 1.0) + if src.kind == "subcortical_surf": + self._data[hemi]["alpha"] = alpha + return self._add_subcortical_surf_data(hemi, src), None resolution = volume_options["resolution"] surface_alpha = volume_options["surface_alpha"] if surface_alpha is None: @@ -1295,6 +1298,40 @@ def _add_volume_data(self, hemi, src, volume_options): return actor_pos, actor_neg + def _add_subcortical_surf_data(self, hemi, src): + from ...morph import _hemi_morph + + hemi_data = self._data[hemi] + mesh = hemi_data.get("subcortical_mesh") + if mesh is not None: # already created in another view + self._renderer.plotter.add_actor(mesh._actor, render=False) + return mesh._actor + # all meshes go in one LayeredMesh, so offset each one's vertex numbers + offsets = np.cumsum([0] + [s["np"] for s in src])[:-1] + rr = np.concatenate([s["rr"] for s in src]) + if self._units == "mm": + rr = rr * 1000 + tris = np.concatenate([s["tris"] + off for s, off in zip(src, offsets)]) + vertices = np.concatenate([s["vertno"] + off for s, off in zip(src, offsets)]) + if len(vertices) != len(hemi_data["vertices"]): + raise ValueError( + f"The stc has {len(hemi_data['vertices'])} subcortical vertices but " + f"src has {len(vertices)}, they must match" + ) + mesh = LayeredMesh( + self._renderer, + vertices=rr, + triangles=tris, + normals=np.concatenate([s["nn"] for s in src]), + ) + mesh._map() + # the data overlay itself is added by _update_current_time_idx + hemi_data["smooth_mat"] = _hemi_morph( + tris, np.arange(len(rr)), vertices, "nearest", maps=None, warn=False + ) + hemi_data["subcortical_mesh"] = mesh + return mesh._actor + def add_label( self, label, @@ -2957,12 +2994,13 @@ def _update_colormap_range(self, fmin=None, fmid=None, fmax=None, alpha=None): for hemi in ["lh", "rh", "vol"]: hemi_data = self._data.get(hemi) if hemi_data is not None: - if hemi in self.layered_meshes: - mesh = self.layered_meshes[hemi] + mesh = self.layered_meshes.get(hemi, hemi_data.get("subcortical_mesh")) + if mesh is not None: mesh.update_overlay( name=self._active_data_key, colormap=self._data["ctable"], - opacity=alpha, + # like volumes, subcortical meshes use volume_options alpha + opacity=alpha if hemi != "vol" else hemi_data["alpha"], rng=rng, ) self._renderer._set_colormap_range( @@ -3138,8 +3176,8 @@ def _update_current_time_idx(self, time_idx): # grid.point_data["values"][21334] = values.max() # update the mesh scalar values - if hemi in self.layered_meshes: - mesh = self.layered_meshes[hemi] + mesh = self.layered_meshes.get(hemi, hemi_data.get("subcortical_mesh")) + if mesh is not None: mesh.smooth_mat = hemi_data.get("smooth_mat") key_rng = [ ( @@ -3149,7 +3187,10 @@ def _update_current_time_idx(self, time_idx): ), key_data["fmax"], ] - if data_key in mesh._overlays: + if data_key in mesh._overlays and hemi not in self.layered_meshes: + # a subcortical mesh holds only this key's overlay + mesh.update_overlay(data_key, scalars=act_data) + elif data_key in mesh._overlays: # Stage without recomposing; a single update below handles # all overlays in O(N) instead of O(N²). mesh.update_overlay(data_key, scalars=act_data, update=False) diff --git a/mne/viz/_brain/_time_viewer.py b/mne/viz/_brain/_time_viewer.py index 05f6969c2f0..817a0875ab6 100644 --- a/mne/viz/_brain/_time_viewer.py +++ b/mne/viz/_brain/_time_viewer.py @@ -753,8 +753,7 @@ def _update_act_data_smooth(self): act_data = np.linalg.norm(act_data, axis=1) smooth_mat = hemi_data.get("smooth_mat") vertices = hemi_data["vertices"] - if hemi == "vol": - assert smooth_mat is None + if hemi == "vol" and smooth_mat is None: # not subcortical meshes smooth_mat = csr_array( (np.ones(len(vertices)), (vertices, np.arange(len(vertices)))) ) @@ -1067,6 +1066,8 @@ def _on_pick(self, vtk_picker, event): hemi = "vol" if self.act_data_smooth[hemi][0] is None: # no data to add for hemi return + if hemi == "vol" and "grid" not in self._data[hemi]: + return # picking on subcortical surfaces is not supported yet pos = np.array(vtk_picker.GetPickPosition()) if hemi == "vol": # VTK will give us the point closest to the viewer in the vol. From c23a3176b0fbac06e1817a069f3cfc0486564cbe Mon Sep 17 00:00:00 2001 From: payam Date: Fri, 23 Oct 2026 14:28:58 +0200 Subject: [PATCH 2/7] add test --- mne/viz/tests/test_3d.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index df767f83897..491f84cf183 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -24,6 +24,7 @@ read_evokeds, read_forward_solution, read_trans, + setup_subcortical_source_space, setup_volume_source_space, use_coil_def, ) @@ -1556,6 +1557,27 @@ def test_mixed_sources_plot_surface(renderer_interactive): brain.close() del brain + # cortex + subcortical surface source space + src += setup_subcortical_source_space( + "sample", label="Left-Amygdala", aseg="aseg", subjects_dir=subjects_dir + ) + assert src.kind == "mixed" + vertno = [s["vertno"] for s in src] + data = rng.standard_normal((sum(len(v) for v in vertno), T)) + stc = MixedSourceEstimate(data, vertno, 0, 1, subject="sample") + brain = stc.plot( + src=src, surface="white", subjects_dir=subjects_dir, colorbar=False, alpha=0.1 + ) + mesh = brain._data["vol"]["subcortical_mesh"] + assert mesh._polydata.n_points == src[2]["np"] + # the cortex alpha does not apply, as with volumes + assert mesh._overlays[brain._active_data_key]._opacity == 1.0 + colors = mesh._current_colors.copy() + brain.set_time_point(1) + assert not np.array_equal(colors, mesh._current_colors) + brain.close() + del brain + @testing.requires_testing_data @pytest.mark.slowtest From c732d0efa5d9b4631755f170daf2d62155727b20 Mon Sep 17 00:00:00 2001 From: payam Date: Fri, 23 Oct 2026 14:29:24 +0200 Subject: [PATCH 3/7] small explanation as warning --- mne/source_space/_source_space.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/source_space/_source_space.py b/mne/source_space/_source_space.py index 133dc37975d..1985c548367 100644 --- a/mne/source_space/_source_space.py +++ b/mne/source_space/_source_space.py @@ -2154,9 +2154,9 @@ def setup_subcortical_source_space( created by this function are not (yet) compatible with morphing (:class:`~mne.SourceMorph`), :func:`mne.extract_label_time_course` does not yet know how to select vertices within a subcortical-surface - label, and plotting support is limited (for example, the - :meth:`~mne.MixedSourceEstimate.plot` method does not yet support - these source spaces). + label, and plotting support is limited (for example, + :meth:`~mne.MixedSourceEstimate.plot` shows these source spaces, but + picking time courses on them is not supported). Parameters ---------- From bd1b241f302ad8a4867d32435a630f57782ecb6f Mon Sep 17 00:00:00 2001 From: payam Date: Fri, 23 Oct 2026 14:29:43 +0200 Subject: [PATCH 4/7] add it to the example --- examples/visualization/brain.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/examples/visualization/brain.py b/examples/visualization/brain.py index 49719eb927e..ce3444e4fdf 100644 --- a/examples/visualization/brain.py +++ b/examples/visualization/brain.py @@ -269,3 +269,33 @@ def gaussian_patch(coords, center, sigma=15.0): smoothing_steps=5, ) brain.show_view(azimuth=180, elevation=70, distance=380, focalpoint=(0, 10, 20)) + +# %% +# Plot data on subcortical structures +# ----------------------------------- +# +# :func:`mne.setup_subcortical_source_space` builds surface source spaces for +# subcortical structures. Their data is added with ``hemi="vol"`` together with +# the source space. Because the meshes are in real anatomical coordinates, use a +# non-inflated cortical surface such as ``surf="white"``. Here we show a synthetic +# gradient along the anterior-posterior axis of the left amygdala and hippocampus. +# :meth:`mne.MixedSourceEstimate.plot` shows them in the same way. + +src = mne.setup_subcortical_source_space( + "sample", label=["Left-Hippocampus", "Left-Amygdala"], subjects_dir=subjects_dir +) +rr = np.concatenate([s["rr"][s["vertno"]] for s in src]) +data = (rr[:, 1] - rr[:, 1].min()) / np.ptp(rr[:, 1]) +brain = mne.viz.Brain( + "sample", subjects_dir=subjects_dir, hemi="lh", surf="white", **brain_kwargs +) +brain.add_data( + data, + hemi="vol", + vertices=np.concatenate([s["vertno"] for s in src]), + src=src, + fmin=0, + fmax=1, + colormap="viridis", +) +brain.show_view(azimuth=150, elevation=90, distance=400, focalpoint=(0, -10, 10)) From 2960816c9b3d94810c88237dc559425837d813bc Mon Sep 17 00:00:00 2001 From: payam Date: Fri, 23 Oct 2026 14:31:47 +0200 Subject: [PATCH 5/7] changelog --- doc/changes/dev/14450.newfeature.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 doc/changes/dev/14450.newfeature.rst diff --git a/doc/changes/dev/14450.newfeature.rst b/doc/changes/dev/14450.newfeature.rst new file mode 100644 index 00000000000..1164e5e9ef6 --- /dev/null +++ b/doc/changes/dev/14450.newfeature.rst @@ -0,0 +1 @@ +:meth:`mne.MixedSourceEstimate.plot` now shows source spaces created by :func:`mne.setup_subcortical_source_space` as surface meshes alongside the cortex, by `Payam Sadeghi-Shabestari`_. From 57e7224438884fd31c8777406f652af8c2b8b983 Mon Sep 17 00:00:00 2001 From: payam Date: Sat, 24 Oct 2026 14:11:53 +0200 Subject: [PATCH 6/7] copy marching cubes output so it does not keep VTK arrays alive --- mne/surface.py | 5 ++--- mne/tests/test_surface.py | 1 + 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/mne/surface.py b/mne/surface.py index 7c9bee6a3ff..833aa59d348 100644 --- a/mne/surface.py +++ b/mne/surface.py @@ -1960,9 +1960,8 @@ def _marching_cubes(image, level, smooth=0, fill_hole_size=None, use_flying_edge polydata = geometry.GetOutput() rr = vtk_to_numpy(polydata.GetPoints().GetData()) tris = vtk_to_numpy(polydata.GetPolys().GetConnectivityArray()).reshape(-1, 3) - rr = np.ascontiguousarray(rr) - tris = np.ascontiguousarray(tris) - out.append((rr, tris)) + # copy, as views on VTK memory would keep the VTK arrays alive + out.append((rr.copy(), tris.copy())) return out diff --git a/mne/tests/test_surface.py b/mne/tests/test_surface.py index 7db6a19bdfe..27104503297 100644 --- a/mne/tests/test_surface.py +++ b/mne/tests/test_surface.py @@ -286,6 +286,7 @@ def test_marching_cubes(dtype, value, smooth, order): out = _marching_cubes(data, level, smooth=smooth) assert len(out) == 1 verts, triangles = out[0] + assert verts.flags.owndata and triangles.flags.owndata # not views on VTK data # verts and faces are rather large so use checksum rtol = 1e-2 if smooth else 1e-9 assert_allclose(verts.sum(axis=0), [14700, 14700, 14700], rtol=rtol) From fcda7c1e5aab366e0da4f3219b1888f0476b5f86 Mon Sep 17 00:00:00 2001 From: payam Date: Tue, 27 Oct 2026 15:16:56 +0100 Subject: [PATCH 7/7] fix coordinate and add a base layer --- mne/viz/_brain/_brain.py | 59 ++++++++++++++++++++++++++++++++-------- mne/viz/tests/test_3d.py | 14 ++++++++-- 2 files changed, 59 insertions(+), 14 deletions(-) diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py index 4c8d61ce031..ca05fe4af94 100644 --- a/mne/viz/_brain/_brain.py +++ b/mne/viz/_brain/_brain.py @@ -12,6 +12,7 @@ import numpy as np from scipy.interpolate import interp1d +from ..._fiff.constants import FIFF from ..._fiff.meas_info import Info from ..._fiff.pick import pick_types from ..._freesurfer import ( @@ -34,6 +35,7 @@ _get_trans, _get_transforms_to_coord_frame, apply_trans, + fit_matched_points, ) from ...utils import ( Bunch, @@ -1309,8 +1311,29 @@ def _add_subcortical_surf_data(self, hemi, src): # all meshes go in one LayeredMesh, so offset each one's vertex numbers offsets = np.cumsum([0] + [s["np"] for s in src])[:-1] rr = np.concatenate([s["rr"] for s in src]) - if self._units == "mm": - rr = rr * 1000 + nn = np.concatenate([s["nn"] for s in src]) + scale = 1000.0 if self._units == "mm" else 1.0 + if src[0]["coord_frame"] == FIFF.FIFFV_COORD_HEAD: # e.g., fwd["src"] + for hemi, s in zip(("lh", "rh"), self._data["src"][:2]): + if hemi in self.geo and s["np"] == len(self.geo[hemi].coords): + try: + head_mri_t = fit_matched_points( + s["rr"], self.geo[hemi].coords / scale, tol=1e-4 + ) + break + except RuntimeError: # not the same surface + pass + else: + raise ValueError( + "The subcortical source spaces are in head coordinates, and " + "the head to MRI transform could not be determined from the " + "cortical source spaces. Transform them to MRI coordinates " + "first, e.g. with mne.transform_surface_to(s, 'mri', trans) " + "for each source space." + ) + rr = apply_trans(head_mri_t, rr) + nn = apply_trans(head_mri_t, nn, move=False) + rr = rr * scale tris = np.concatenate([s["tris"] + off for s, off in zip(src, offsets)]) vertices = np.concatenate([s["vertno"] + off for s, off in zip(src, offsets)]) if len(vertices) != len(hemi_data["vertices"]): @@ -1318,20 +1341,31 @@ def _add_subcortical_surf_data(self, hemi, src): f"The stc has {len(hemi_data['vertices'])} subcortical vertices but " f"src has {len(vertices)}, they must match" ) - mesh = LayeredMesh( - self._renderer, - vertices=rr, - triangles=tris, - normals=np.concatenate([s["nn"] for s in src]), - ) + mesh = LayeredMesh(self._renderer, vertices=rr, triangles=tris, normals=nn) mesh._map() - # the data overlay itself is added by _update_current_time_idx + # gray base layer like the cortex has, so the mesh stays visible where the + # data is transparent; the data overlay is added by _update_current_time_idx + mesh.add_overlay( + scalars=np.zeros(len(rr)), + colormap=np.array([self._brain_color]) * 255, + rng=[0, 1], + opacity=self._alpha, + name="curv", + ) hemi_data["smooth_mat"] = _hemi_morph( tris, np.arange(len(rr)), vertices, "nearest", maps=None, warn=False ) hemi_data["subcortical_mesh"] = mesh return mesh._actor + @property + def _subcortical_meshes(self): + return [ + key_data["vol"]["subcortical_mesh"] + for key_data in self._all_data.values() + if "subcortical_mesh" in key_data.get("vol", {}) + ] + def add_label( self, label, @@ -2755,8 +2789,9 @@ def set_cortex_alpha(self, alpha): The opacity of the cortical surface, between 0 and 1. """ self._alpha = float(alpha) - for hemi in self._hemis: - self.layered_meshes[hemi].update_overlay("curv", opacity=self._alpha) + meshes = [self.layered_meshes[hemi] for hemi in self._hemis] + for mesh in meshes + self._subcortical_meshes: + mesh.update_overlay("curv", opacity=self._alpha) self._renderer._update() def set_cortex_colormap(self, cortex): @@ -2786,6 +2821,8 @@ def set_cortex_colormap(self, cortex): colormap=geo_kwargs["colormap"], rng=[geo_kwargs["vmin"], geo_kwargs["vmax"]], ) + for mesh in self._subcortical_meshes: + mesh.update_overlay("curv", colormap=np.array([self._brain_color]) * 255) self._renderer._update() def set_surf(self, surf): diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index 491f84cf183..20800246c15 100644 --- a/mne/viz/tests/test_3d.py +++ b/mne/viz/tests/test_3d.py @@ -1565,13 +1565,21 @@ def test_mixed_sources_plot_surface(renderer_interactive): vertno = [s["vertno"] for s in src] data = rng.standard_normal((sum(len(v) for v in vertno), T)) stc = MixedSourceEstimate(data, vertno, 0, 1, subject="sample") + src_head = src.copy()._transform_to("head", read_trans(trans_fname)) brain = stc.plot( - src=src, surface="white", subjects_dir=subjects_dir, colorbar=False, alpha=0.1 + src=src_head, + surface="white", + subjects_dir=subjects_dir, + colorbar=False, + alpha=0.1, ) mesh = brain._data["vol"]["subcortical_mesh"] - assert mesh._polydata.n_points == src[2]["np"] - # the cortex alpha does not apply, as with volumes + assert_allclose(mesh._polydata.points, src[2]["rr"] * 1000, atol=1e-2) + # gray base layer with the cortex alpha, data layer with the volume alpha + assert mesh._overlays["curv"]._opacity == 0.1 assert mesh._overlays[brain._active_data_key]._opacity == 1.0 + brain.set_cortex_alpha(0.5) + assert mesh._overlays["curv"]._opacity == 0.5 colors = mesh._current_colors.copy() brain.set_time_point(1) assert not np.array_equal(colors, mesh._current_colors)