ENH: Add upsampling for MEG helmet surface (#13179)

Co-authored-by: autofix-ci[bot] <114827586+autofix-ci[bot]@users.noreply.github.com>
Co-authored-by: Daniel McCloy <dan@mccloy.info>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

Author: 4 people

doc/changes/devel/13179.bugfix.rst

@@ -0,0 +1 @@
+Fix bug with :func:`mne.viz.plot_alignment` where ``eeg="projected"`` was not plotted, by `Eric Larson`_.

doc/changes/devel/13179.newfeature.rst

@@ -0,0 +1 @@
+Add ``upsampling`` option to :func:`mne.make_field_map` to allow upsampling MEG helmet surfaces for plotting, by `Eric Larson`_.

examples/visualization/eeg_on_scalp.py

@@ -5,7 +5,7 @@
 Plotting EEG sensors on the scalp
 =================================
 
-In this example, digitized EEG sensor locations are shown on the scalp.
+In this example, digitized EEG sensor locations are shown on the scalp surface.
 """
 # Author: Eric Larson <larson.eric.d@gmail.com>
 #

examples/visualization/mne_helmet.py

@@ -35,9 +35,10 @@
     ch_type="meg",
     subject="sample",
     subjects_dir=subjects_dir,
+    upsampling=2,
 )
 time = 0.083
-fig = mne.viz.create_3d_figure((256, 256))
+fig = mne.viz.create_3d_figure((512, 512), bgcolor="w", title="MNE helmet")
 mne.viz.plot_alignment(
     evoked.info,
     subject="sample",
@@ -50,13 +51,18 @@
     coord_frame="mri",
 )
 evoked.plot_field(
-    maps, time=time, fig=fig, time_label=None, vmax=5e-13, time_viewer=False
+    maps,
+    time=time,
+    fig=fig,
+    time_label=None,
+    vmax=5e-13,
+    time_viewer=False,
 )
 mne.viz.set_3d_view(
     fig,
     azimuth=40,
     elevation=87,
     focalpoint=(0.0, -0.01, 0.04),
-    roll=-25,
+    roll=-105,
     distance=0.55,
 )

mne/forward/_field_interpolation.py

@@ -446,6 +446,7 @@ def make_field_map(
     origin=(0.0, 0.0, 0.04),
     n_jobs=None,
     *,
+    upsampling=1,
     head_source=("bem", "head"),
     verbose=None,
 ):
@@ -483,6 +484,9 @@ def make_field_map(
 
         .. versionadded:: 0.11
     %(n_jobs)s
+    %(helmet_upsampling)s
+
+        .. versionadded:: 1.10
     %(head_source)s
 
         .. versionadded:: 1.1
@@ -527,7 +531,7 @@ def make_field_map(
     surfs = []
     for this_type in types:
         if this_type == "meg" and meg_surf == "helmet":
-            surf = get_meg_helmet_surf(info, trans)
+            surf = get_meg_helmet_surf(info, trans, upsampling=upsampling)
         else:
             surf = get_head_surf(subject, source=head_source, subjects_dir=subjects_dir)
         surfs.append(surf)

mne/surface.py

@@ -44,6 +44,7 @@
     _hashable_ndarray,
     _import_nibabel,
     _pl,
+    _soft_import,
     _TempDir,
     _validate_type,
     fill_doc,
@@ -173,7 +174,7 @@ def _get_head_surface(subject, source, subjects_dir, on_defects, raise_error=Tru
 
 
 @verbose
-def get_meg_helmet_surf(info, trans=None, *, verbose=None):
+def get_meg_helmet_surf(info, trans=None, *, upsampling=1, verbose=None):
     """Load the MEG helmet associated with the MEG sensors.
 
     Parameters
@@ -183,6 +184,7 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
         The head<->MRI transformation, usually obtained using
         read_trans(). Can be None, in which case the surface will
         be in head coordinates instead of MRI coordinates.
+    %(helmet_upsampling)s
     %(verbose)s
 
     Returns
@@ -198,7 +200,10 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
     from .bem import _fit_sphere, read_bem_surfaces
     from .channels.channels import _get_meg_system
 
+    _validate_type(upsampling, "int", "upsampling")
+
     system, have_helmet = _get_meg_system(info)
+    incomplete = False
     if have_helmet:
         logger.info(f"Getting helmet for system {system}")
         fname = _helmet_path / f"{system}.fif.gz"
@@ -231,8 +236,24 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
         # remove the frontal point we added from the simplices
         tris = tris[(tris != len(sph) - 1).all(-1)]
         tris = _reorder_ccw(rr, tris)
-
         surf = dict(rr=rr, tris=tris)
+        incomplete = True
+    if upsampling > 1:
+        # Use VTK (could also use Butterfly but Loop is smoother)
+        pv = _soft_import("pyvista", "upsample a mesh")
+        factor = 4 ** (upsampling - 1)
+        rr, tris = surf["rr"], surf["tris"]
+        logger.info(
+            f"Upsampling from {len(rr)} to {len(rr) * factor} vertices ({upsampling=})"
+        )
+        tris = np.c_[np.full(len(tris), 3), tris]
+        mesh = pv.PolyData(rr, tris)
+        mesh = mesh.subdivide(upsampling - 1, subfilter="linear")
+        rr, tris = mesh.points, mesh.faces.reshape(-1, 4)[:, 1:]
+        tris = _reorder_ccw(rr, tris)
+        surf = dict(rr=rr, tris=tris)
+        incomplete = True
+    if incomplete:
         complete_surface_info(surf, copy=False, verbose=False)
 
     # Ignore what the file says, it's in device coords and we want MRI coords

mne/utils/docs.py

@@ -2009,6 +2009,13 @@ def _reflow_param_docstring(docstring, has_first_line=True, width=75):
     :func:`mne.get_head_surf` for more information.
 """
 
+docdict["helmet_upsampling"] = """
+upsampling : int
+    The upsampling factor to use for the helmet mesh. The default (1) does no
+    upsampling. Larger integers lead to more densely sampled helmet surfaces, and
+    the number of vertices increases as a factor of ``4**(upsampling-1)``.
+"""
+
 docdict["hitachi_fname"] = """
 fname : list | str
     Path(s) to the Hitachi CSV file(s). This should only be a list for

mne/viz/_3d.py

@@ -915,6 +915,7 @@ def plot_alignment(
             surf, coord_frame, [to_cf_t["mri"], to_cf_t["head"]], copy=True
         )
         renderer.surface(
+            name=key,
             surface=surf,
             color=colors[key],
             opacity=alphas[key],
@@ -1021,7 +1022,7 @@ def _handle_sensor_types(meg, eeg, fnirs):
 
     alpha_map = dict(
         meg=dict(sensors="meg", helmet="meg_helmet", ref="ref_meg"),
-        eeg=dict(original="eeg", projected="eeg_projected"),
+        eeg=dict(original="eeg", projected="eegp"),
         fnirs=dict(channels="fnirs", pairs="fnirs_pairs"),
     )
     sensor_alpha = {
@@ -1169,7 +1170,11 @@ def _plot_helmet(
             src_surf, coord_frame, [to_cf_t["mri"], to_cf_t["head"]], copy=True
         )
     actor, dst_surf = renderer.surface(
-        surface=src_surf, color=color, opacity=alpha, backface_culling=False
+        surface=src_surf,
+        color=color,
+        opacity=alpha,
+        backface_culling=False,
+        name="helmet",
     )
     return actor, dst_surf, src_surf
 
@@ -1538,7 +1543,7 @@ def _plot_sensors_3d(
             if not fnirs or "channels" not in fnirs:
                 plot_sensors = False
         elif ch_type == "eeg":
-            if not eeg or "original" not in eeg:
+            if not eeg:
                 plot_sensors = False
         elif ch_type == "meg":
             if not meg or "sensors" not in meg:
@@ -1547,7 +1552,13 @@ def _plot_sensors_3d(
         if isinstance(ch_coord, tuple):  # is meg, plot coil
             ch_coord = dict(rr=ch_coord[0] * unit_scalar, tris=ch_coord[1])
         if plot_sensors:
-            locs[ch_type].append(ch_coord)
+            if ch_type == "eeg":
+                if "original" in eeg:
+                    locs[ch_type].append(ch_coord)
+                if "projected" in eeg:
+                    locs["eegp"].append(ch_coord)
+            else:
+                locs[ch_type].append(ch_coord)
         if ch_name in sources and "sources" in fnirs:
             locs["source"].append(sources[ch_name])
         if ch_name in detectors and "detectors" in fnirs:
@@ -1628,7 +1639,30 @@ def _plot_sensors_3d(
         else:
             sens_loc = np.array(sens_loc, float)
             mask = ~np.isnan(sens_loc).any(axis=1)
-            if len(colors) == 1 and len(scales) == 1:
+            if ch_type == "eegp":  # special case, need to project
+                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
+                )[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,
+                )
+                actors["eeg"].append(actor)
+            elif len(colors) == 1 and len(scales) == 1:
                 # Single color mode (one actor)
                 actor, _ = _plot_glyphs(
                     renderer=renderer,
@@ -1701,29 +1735,7 @@ def _plot_sensors_3d(
                         nearest=nearest,
                     )
                     actors[ch_type].append(actor)
-        if ch_type == "eeg" and "projected" in eeg:
-            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
-            )[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=defaults["eegp_color"],
-                mode="cylinder",
-                scale=defaults["eegp_scale"] * unit_scalar,
-                opacity=sensor_alpha["eeg_projected"],
-                glyph_height=defaults["eegp_height"],
-                glyph_center=(0.0, -defaults["eegp_height"] / 2.0, 0),
-                glyph_resolution=20,
-                backface_culling=True,
-            )
-            actors["eeg"].append(actor)
+
     actors = dict(actors)  # get rid of defaultdict
 
     return actors

mne/viz/backends/_abstract.py

@@ -127,6 +127,7 @@ def mesh(
         triangles,
         color,
         opacity=1.0,
+        *,
         backface_culling=False,
         scalars=None,
         colormap=None,
@@ -137,6 +138,7 @@ def mesh(
         line_width=1.0,
         normals=None,
         polygon_offset=None,
+        name=None,
         **kwargs,
     ):
         """Add a mesh in the scene.
@@ -183,6 +185,8 @@ def mesh(
             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
             The arguments to pass to triangular_mesh
 
@@ -254,6 +258,8 @@ def surface(
         scalars=None,
         backface_culling=False,
         polygon_offset=None,
+        *,
+        name=None,
     ):
         """Add a surface in the scene.
 
@@ -281,6 +287,8 @@ def surface(
             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.
         """
         pass
 

mne/viz/backends/_pyvista.py

@@ -356,6 +356,8 @@ def polydata(
         representation="surface",
         line_width=1.0,
         polygon_offset=None,
+        *,
+        name=None,
         **kwargs,
     ):
         from matplotlib.colors import to_rgba_array
@@ -386,6 +388,7 @@ def polydata(
         actor = _add_mesh(
             plotter=self.plotter,
             mesh=mesh,
+            name=name,
             color=color,
             scalars=scalars,
             edge_color=color,
@@ -430,6 +433,7 @@ def mesh(
         line_width=1.0,
         normals=None,
         polygon_offset=None,
+        name=None,
         **kwargs,
     ):
         vertices = np.c_[x, y, z].astype(float)
@@ -449,6 +453,7 @@ def mesh(
             representation=representation,
             line_width=line_width,
             polygon_offset=polygon_offset,
+            name=name,
             **kwargs,
         )
 
@@ -504,6 +509,8 @@ def surface(
         scalars=None,
         backface_culling=False,
         polygon_offset=None,
+        *,
+        name=None,
     ):
         normals = surface.get("nn", None)
         vertices = np.array(surface["rr"])
@@ -524,6 +531,7 @@ def surface(
             vmin=vmin,
             vmax=vmax,
             polygon_offset=polygon_offset,
+            name=name,
         )
 
     def sphere(
@@ -1019,7 +1027,6 @@ def _silhouette(self, mesh, color=None, line_width=None, alpha=None, decimate=No
         silhouette_mapper.SetInputConnection(silhouette_filter.GetOutputPort())
         actor, prop = self.plotter.add_actor(
             silhouette_mapper,
-            name=None,
             culling=False,
             pickable=False,
             reset_camera=False,