refactor plot_alignment to the extent possible

doc/changes/latest.inc

@@ -175,6 +175,8 @@ Bugs
 
 - Fix bug in :func:`mne.viz.set_3d_backend` and :func:`mne.viz.get_3d_backend` where the PyVistaQt-based backend was ambiguously named ``'pyvista'`` instead of ``'pyvistaqt'``; use ``set_3d_backend('pyvistaqt')`` and expect ``'pyvistaqt'`` as the output of :func:`mne.viz.get_3d_backend` instead of ``'pyvista'``, and consider using ``get_3d_backend().startswith('pyvista')`` for example for backward-compatible conditionals (:gh:`9607` by `Guillaume Favelier`_)
 
+- Fix bug where "seeg", "ecog", "dbs" and "fnirs" data had coordinate frame unknown upon loading from a file when it should have been in "head" (:gh:`9580` by `Alex Rockhill`_)
+
 API changes
 ~~~~~~~~~~~
 - In `mne.compute_source_morph`, the ``niter_affine`` and ``niter_sdr`` parameters have been replaced by ``niter`` and ``pipeline`` parameters for more consistent and finer-grained control of registration/warping steps and iteration (:gh:`9505` by `Alex Rockhill`_ and `Eric Larson`_)

mne/_freesurfer.py

@@ -702,3 +702,43 @@ def _get_head_surface(surf, subject, subjects_dir, bem=None, verbose=None):
                 return _read_mri_surface(fname)
     raise IOError('No head surface found for subject '
                   f'{subject} after trying:\n' + '\n'.join(try_fnames))
+
+
+@verbose
+def _get_skull_surface(surf, subject, subjects_dir, bem=None, verbose=None):
+    """Get a skull surface from the Freesurfer subject directory.
+
+    Parameters
+    ----------
+    surf : str
+        The name of the surface 'outer' or 'inner'.
+    %(subject)s
+    %(subjects_dir)s
+    bem : mne.bem.ConductorModel | None
+        The conductor model that stores information about the skull surface.
+    %(verbose)s
+
+    Returns
+    -------
+    skull_surf : dict | None
+        A dictionary with keys 'rr', 'tris', 'ntri', 'use_tris', 'np'
+        and 'coord_frame' that store information for mesh plotting and other
+        useful information about the head surface.
+
+    Notes
+    -----
+    .. versionadded: 0.24
+    """
+    if bem is not None:
+        try:
+            return _bem_find_surface(bem, surf + '_skull')
+        except RuntimeError:
+            logger.info('Could not find the surface for '
+                        'skull in the provided BEM model, '
+                        'looking in the subject directory.')
+    fname = op.join(get_subjects_dir(subjects_dir, raise_error=True),
+                    subject, 'bem', surf + '_skull.surf')
+    if not op.isfile(fname):
+        raise ValueError('Skull surface cannot be found, should be '
+                         f'located at {fname}')
+    return _read_mri_surface(fname)

mne/channels/tests/test_montage.py

@@ -1362,10 +1362,8 @@ def test_montage_head_frame(ch_type):
     # gh-9446
     data = np.random.randn(2, 100)
     info = create_info(['a', 'b'], 512, ch_type)
-    coord_frame = getattr(
-        FIFF, f'FIFFV_COORD_{"HEAD" if ch_type == "eeg" else "UNKNOWN"}')
     for ch in info['chs']:
-        assert ch['coord_frame'] == coord_frame
+        assert ch['coord_frame'] == FIFF.FIFFV_COORD_HEAD
     raw = RawArray(data, info)
     ch_pos = dict(a=[-0.00250136, 0.04913788, 0.05047056],
                   b=[-0.00528394, 0.05066484, 0.05061559])

mne/io/tag.py

@@ -333,6 +333,10 @@ def _read_coord_trans_struct(fid, tag, shape, rlims):
     FIFF.FIFFV_MEG_CH: FIFF.FIFFV_COORD_DEVICE,
     FIFF.FIFFV_REF_MEG_CH: FIFF.FIFFV_COORD_DEVICE,
     FIFF.FIFFV_EEG_CH: FIFF.FIFFV_COORD_HEAD,
+    FIFF.FIFFV_ECOG_CH: FIFF.FIFFV_COORD_HEAD,
+    FIFF.FIFFV_SEEG_CH: FIFF.FIFFV_COORD_HEAD,
+    FIFF.FIFFV_DBS_CH: FIFF.FIFFV_COORD_HEAD,
+    FIFF.FIFFV_FNIRS_CH: FIFF.FIFFV_COORD_HEAD
 }
 
 

mne/io/tests/test_meas_info.py

@@ -317,6 +317,18 @@ def test_io_dig_points(tmpdir):
             read_polhemus_fastscan(dest, on_header_missing='warn')
 
 
+def test_io_coord_frame(tmpdir):
+    """Test round trip for coordinate frame."""
+    fname = tmpdir.join('test.fif')
+    for ch_type in ('eeg', 'seeg', 'ecog', 'dbs', 'hbo', 'hbr'):
+        info = create_info(
+            ch_names=['Test Ch'], sfreq=1000., ch_types=[ch_type])
+        info['chs'][0]['loc'][:3] = [0.05, 0.01, -0.03]
+        write_info(fname, info)
+        info2 = read_info(fname)
+        assert info2['chs'][0]['coord_frame'] == FIFF.FIFFV_COORD_HEAD
+
+
 def test_make_dig_points():
     """Test application of Polhemus HSP to info."""
     extra_points = read_polhemus_fastscan(

mne/viz/_3d.py

@@ -21,7 +21,8 @@
 
 from ..defaults import DEFAULTS
 from ..fixes import _crop_colorbar, _get_img_fdata, _get_args
-from .._freesurfer import _read_mri_info, _check_mri, _get_head_surface
+from .._freesurfer import (_read_mri_info, _check_mri, _get_head_surface,
+                           compute_mri_head_t, _get_skull_surface)
 from ..io import _loc_to_coil_trans
 from ..io.pick import pick_types, _picks_to_idx
 from ..io.constants import FIFF
@@ -30,10 +31,9 @@
                             SourceSpaces, read_freesurfer_lut)
 
 from ..surface import (get_meg_helmet_surf, _read_mri_surface, _DistanceQuery,
-                       transform_surface_to, _project_onto_surface,
-                       _reorder_ccw)
-from ..transforms import (_find_trans, apply_trans, rot_to_quat,
-                          combine_transforms, _get_trans, _ensure_trans,
+                       _project_onto_surface, _reorder_ccw)
+from ..transforms import (apply_trans, rot_to_quat, combine_transforms,
+                          _get_trans, _ensure_trans, transform_surface_to,
                           invert_transform, Transform, rotation,
                           read_ras_mni_t, _print_coord_trans)
 from ..utils import (get_subjects_dir, logger, _check_subject, verbose, warn,
@@ -414,6 +414,12 @@ def plot_evoked_field(evoked, surf_maps, time=None, time_label='t = %0.0f ms',
     return renderer.scene()
 
 
+def _apply_trans_to_points(trans, points):
+    """Apply the transform to change the space."""
+    for ch in points:
+        points[ch] = apply_trans(trans, points)
+
+
 @verbose
 def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
                    surfaces='auto', coord_frame='head',
@@ -656,7 +662,7 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
     if trans == 'auto':
         # let's try to do this in MRI coordinates so they're easy to plot
         subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
-        trans = _find_trans(subject, subjects_dir)
+        trans = compute_mri_head_t(subject, subjects_dir)
     head_mri_t, _ = _get_trans(trans, 'head', 'mri')
     dev_head_t, _ = _get_trans(info['dev_head_t'], 'meg', 'head')
     del trans
@@ -692,23 +698,12 @@ def plot_alignment(info=None, trans=None, subject=None, subjects_dir=None,
             surfs['head'] = head_surf
 
     # Skull:
-    skull = list()
+    skulls = list()
     for name in ('outer_skull', 'inner_skull'):
         if name in surfaces:
             surfaces.pop(surfaces.index(name))
-            if bem is None:
-                fname = op.join(
-                    get_subjects_dir(subjects_dir, raise_error=True),
-                    subject, 'bem', name + '.surf')
-                if not op.isfile(fname):
-                    raise ValueError('bem is None and the the %s file cannot '
-                                     'be found:\n%s' % (name, fname))
-                surf = _read_mri_surface(fname)
-            else:
-                surf = _bem_find_surface(bem, name).copy()
-            surf['name'] = name
-            skull.append(surf)
-    assert all(isinstance(s, dict) for s in skull)
+            skulls.append(_get_skull_surface(
+                name.split('_')[0], subject, subjects_dir, bem=bem))
 
     if mri_fiducials:
         if mri_fiducials is True:

mne/viz/_brain/_brain.py

@@ -35,7 +35,7 @@
 from ...externals.decorator import decorator
 from ..._freesurfer import (vertex_to_mni, read_talxfm, read_freesurfer_lut,
                             _get_head_surface, compute_mri_head_t,
-                            _ch_pos_in_mri)
+                            _ch_pos_in_mri, _get_skull_surface)
 from ...io.pick import (pick_types, channel_type, _picks_to_idx,
                         _FNIRS_CH_TYPES_SPLIT)
 from ...io.meas_info import Info
@@ -2361,6 +2361,45 @@ def add_head(self, dense=True, color=None, alpha=0.5):
 
         self._renderer._update()
 
+    @fill_doc
+    def add_skull(self, outer=True, inner=False, color=None, alpha=0.5):
+        """Add a mesh to render the skull surface.
+
+        Parameters
+        ----------
+        outer : bool
+            Whether render the outer skull.
+        inner : bool
+            Whether to render the inner skull.
+        color : matplotlib-style color | None
+            A list of anything matplotlib accepts: string, RGB, hex, etc.
+            (default: "gray").
+        alpha : float in [0, 1]
+            Alpha level to control opacity.
+
+        Notes
+        -----
+        .. versionadded:: 0.24
+        """
+        from matplotlib.colors import colorConverter
+
+        for surf_name in ['outer'] * outer + ['inner'] * inner:
+            surf = _get_skull_surface(surf_name, self._subject_id,
+                                      self._subjects_dir)
+            verts, triangles = surf['rr'], surf['tris']
+            verts *= 1e3 if self._units == 'mm' else 1
+            if color is None:
+                color = 'gray'
+            color = colorConverter.to_rgba(color, alpha)
+
+            for h in self._hemis:
+                for ri, ci, v in self._iter_views(h):
+                    self._renderer.mesh(
+                        *verts.T, triangles=triangles, color=color,
+                        opacity=alpha, reset_camera=False, render=False)
+
+        self._renderer._update()
+
     @fill_doc
     def add_volume_labels(self, aseg='aparc+aseg', labels=None, colors=None,
                           alpha=0.5, smooth=0.9, legend=None):

mne/viz/_brain/tests/test_brain.py

@@ -315,8 +315,9 @@ def __init__(self):
     brain.add_foci([0], coords_as_verts=True,
                    hemi=hemi, color='blue')
 
-    # add head
+    # add head and skull
     brain.add_head(color='red', alpha=0.1)
+    brain.add_skull(outer=True, inner=True, color='green', alpha=0.1)
 
     # add volume labels
     brain.add_volume_labels(