Deep source modeling (phase 2)

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`_.

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))

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
     ----------

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
 
 

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)

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)

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.

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