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`_. 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)) 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 ---------- 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) diff --git a/mne/viz/_brain/_brain.py b/mne/viz/_brain/_brain.py index 46567fb21c0..9869c6e16a1 100644 --- a/mne/viz/_brain/_brain.py +++ b/mne/viz/_brain/_brain.py @@ -11,10 +11,11 @@ import numpy as np from scipy.interpolate import interp1d +from ..._fiff.constants import FIFF from ..._freesurfer import _estimate_talxfm_rigid from ...defaults import _handle_default from ...surface import _decimate_surface_ico_oct, _mesh_borders, mesh_edges -from ...transforms import apply_trans +from ...transforms import apply_trans, fit_matched_points from ...utils import ( Bunch, _auto_weakref, @@ -1141,7 +1142,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): @@ -1177,6 +1178,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: @@ -1269,6 +1273,72 @@ 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]) + 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"]): + 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=nn) + mesh._map() + # 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, @@ -1976,8 +2046,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): @@ -2007,6 +2078,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): @@ -2215,12 +2288,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( @@ -2396,8 +2470,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 = [ ( @@ -2407,7 +2481,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. diff --git a/mne/viz/tests/test_3d.py b/mne/viz/tests/test_3d.py index df767f83897..20800246c15 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,35 @@ 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") + src_head = src.copy()._transform_to("head", read_trans(trans_fname)) + brain = stc.plot( + src=src_head, + surface="white", + subjects_dir=subjects_dir, + colorbar=False, + alpha=0.1, + ) + mesh = brain._data["vol"]["subcortical_mesh"] + 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) + brain.close() + del brain + @testing.requires_testing_data @pytest.mark.slowtest