diff --git a/mne/viz/_3d.py b/mne/viz/_3d.py
index 4af33fa3005..5b7ced94f2b 100644
--- a/mne/viz/_3d.py
+++ b/mne/viz/_3d.py
@@ -425,7 +425,8 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
                    meg=None, eeg='original', fwd=None,
                    dig=False, ecog=True, src=None, mri_fiducials=False,
                    bem=None, seeg=True, fnirs=True, show_axes=False, dbs=True,
-                   fig=None, interaction='trackball', verbose=None):
+                   fig=None, interaction='trackball', silhouette=None,
+                   verbose=None):
     """Plot head, sensor, and source space alignment in 3D.
 
     Parameters
@@ -522,6 +523,10 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
         .. versionadded:: 0.16
     dbs : bool
         If True (default), show DBS (deep brain stimulation) electrodes.
+    silhouette : dict | None
+       As a dict, it contains the ``color``, ``linewidth`` and ``alpha``opacity
+       of the brain's silhouette to display, otherwise it is None.
+       Defaults to None.
     fig : mayavi.mlab.Figure | None
         Mayavi Scene in which to plot the alignment.
         If ``None``, creates a new 600x600 pixel figure with black background.
@@ -1032,7 +1037,8 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
     for key, surf in surfs.items():
         renderer.surface(surface=surf, color=colors[key],
                          opacity=alphas[key],
-                         backface_culling=(key != 'helmet'))
+                         backface_culling=(key != 'helmet'),
+                         silhouette=silhouette)
     if brain and 'lh' not in surfs:  # one layer sphere
         assert bem['coord_frame'] == FIFF.FIFFV_COORD_HEAD
         center = bem['r0'].copy()
diff --git a/mne/viz/backends/_pyvista.py b/mne/viz/backends/_pyvista.py
index 1b340d2a974..8f595b91bde 100644
--- a/mne/viz/backends/_pyvista.py
+++ b/mne/viz/backends/_pyvista.py
@@ -302,7 +302,7 @@ def polydata(self, mesh, color=None, opacity=1.0, normals=None,
                  backface_culling=False, scalars=None, colormap=None,
                  vmin=None, vmax=None, interpolate_before_map=True,
                  representation='surface', line_width=1.,
-                 polygon_offset=None, **kwargs):
+                 polygon_offset=None, silhouette=None, **kwargs):
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", category=FutureWarning)
             rgba = False
@@ -344,6 +344,25 @@ def polydata(self, mesh, color=None, opacity=1.0, normals=None,
                 mapper.SetRelativeCoincidentTopologyPolygonOffsetParameters(
                     polygon_offset, polygon_offset)
 
+            if silhouette is not None:
+                silhouette_filter = vtk.vtkPolyDataSilhouette()
+                silhouette_filter.SetInputData(mesh.decimate(0.9))
+                silhouette_filter.SetCamera(
+                    self.plotter.renderer.GetActiveCamera())
+                silhouette_filter.SetEnableFeatureAngle(0)
+                silhouette_mapper = vtk.vtkPolyDataMapper()
+                silhouette_mapper.SetInputConnection(
+                    silhouette_filter.GetOutputPort())
+                _, prop = self.plotter.add_actor(
+                    silhouette_mapper, reset_camera=False, name=None,
+                    culling=False, pickable=False)
+                if "color" in silhouette:
+                    prop.SetColor(silhouette["color"])
+                if "alpha" in silhouette:
+                    prop.SetOpacity(silhouette["alpha"])
+                if "linewidth" in silhouette:
+                    prop.SetLineWidth(silhouette["linewidth"])
+
             return actor, mesh
 
     def mesh(self, x, y, z, triangles, color, opacity=1.0, shading=False,
@@ -408,7 +427,8 @@ def contour(self, surface, scalars, contours, width=1.0, opacity=1.0,
     def surface(self, surface, color=None, opacity=1.0,
                 vmin=None, vmax=None, colormap=None,
                 normalized_colormap=False, scalars=None,
-                backface_culling=False, polygon_offset=None):
+                backface_culling=False, polygon_offset=None,
+                silhouette=None):
         with warnings.catch_warnings():
             warnings.filterwarnings("ignore", category=FutureWarning)
             normals = surface.get('nn', None)
@@ -430,6 +450,7 @@ def surface(self, surface, color=None, opacity=1.0,
             vmin=vmin,
             vmax=vmax,
             polygon_offset=polygon_offset,
+            silhouette=silhouette,
         )
 
     def sphere(self, center, color, scale, opacity=1.0,
diff --git a/tutorials/misc/plot_seeg.py b/tutorials/misc/plot_seeg.py
index b7370f80693..832a380cc35 100644
--- a/tutorials/misc/plot_seeg.py
+++ b/tutorials/misc/plot_seeg.py
@@ -118,9 +118,15 @@
 ###############################################################################
 # Let's check to make sure everything is aligned.
 
-fig = mne.viz.plot_alignment(raw.info, trans, 'fsaverage',
-                             subjects_dir=subjects_dir, show_axes=True,
-                             surfaces=["pial", "head"])
+fig = mne.viz.create_3d_figure(size=(800, 600))
+mne.viz.plot_alignment(raw.info, trans, 'fsaverage',
+                       subjects_dir=subjects_dir, show_axes=True,
+                       surfaces=["pial", "head"],
+                       silhouette=dict(
+                           color=(0, 0, 0),
+                           alpha=1.0,
+                           linewidth=4),
+                       fig=fig)
 
 ###############################################################################
 # Next, we'll get the raw data and plot its amplitude over time.