[ENH, MRG] Add mne.viz.Brain.plot_sensors and refactor mne.viz.plot_alignment

doc/_includes/forward.rst

@@ -690,7 +690,7 @@ EEG forward solution in the sphere model
    dataset tutorial <plt_brainstorm_phantom_ctf_eeg_sphere_geometry>`,
    :ref:`Brainstorm Elekta phantom dataset tutorial
    <plt_brainstorm_phantom_elekta_eeg_sphere_geometry>`, and
-   :ref:`plot_source_alignment_without_mri`.
+   :ref:`tut-source-alignment-without-mri`.
 
 When the sphere model is employed, the computation of the EEG solution can be
 substantially accelerated by using approximation methods described by Mosher

doc/changes/0.21.inc

@@ -42,7 +42,7 @@ Enhancements
 
 - Speed up reading of annotations in EDF+ files **by new contributor** |Jeroen Van Der Donckt|_
 
-- Add head to mri and mri to voxel space transform details to :ref:`plot_source_alignment` tutorial, by `Alex Rockhill`_
+- Add head to mri and mri to voxel space transform details to :ref:`tut-source-alignment` tutorial, by `Alex Rockhill`_
 
 - Improve memory efficiency of :func:`mne.concatenate_epochs` by `Eric Larson`_
 

doc/changes/latest.inc

@@ -133,6 +133,10 @@ Enhancements
 
 - Add :func:`mne.channels.read_dig_localite` to read Localite electrode location files (:gh:`9658` by `Clemens Brunner`_)
 
+- Add :meth:`mne.viz.Brain.add_sensors` to plot sensor locations (:gh:`9585` by `Alex Rockhill`_)
+
+- Add :func:`mne.coreg.estimate_head_mri_t` estimate the head->mri transform from fsaverage fiducials (:gh:`9585` by `Alex Rockhill`_)
+
 Bugs
 ~~~~
 - Fix bug with :meth:`mne.Epochs.crop` and :meth:`mne.Evoked.crop` when ``include_tmax=False``, where the last sample was always cut off, even when ``tmax > epo.times[-1]`` (:gh:`9378` **by new contributor** |Jan Sosulski|_)
@@ -187,6 +191,10 @@ Bugs
 
 - Fix bug where setting of a montage with fNIRS data got set to "unknown" coordinate frame when it should have been in "head" (:gh:`9630` by `Alex Rockhill`_)
 
+- 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`_)
+
+- Raise error when no ``trans`` is provided to :func:`mne.viz.plot_alignment` when required instead of assuming identitiy head->mri transform (:gh:`9585` 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`_)

doc/mri.rst

@@ -15,6 +15,7 @@ Step by step instructions for using :func:`gui.coregistration`:
    :toctree: generated/
 
    coreg.get_mni_fiducials
+   coreg.estimate_head_mri_t
    get_montage_volume_labels
    gui.coregistration
    gui.fiducials

doc/overview/cookbook.rst

@@ -314,7 +314,7 @@ for each experimental session, *i.e.*, each time when new head
 digitization data are employed.
 
 The corregistration is stored in ``-trans.fif`` file. If is present,
-you can follow :ref:`plot_source_alignment` to validate its correctness.
+you can follow :ref:`tut-source-alignment` to validate its correctness.
 If the ``-trans.fif`` is not present or the alignment is not correct
 you need to use :func:`mne.gui.coregistration` (or its convenient command line
 equivalent :ref:`mne coreg`) to generate it.

mne/_freesurfer.py

@@ -10,13 +10,16 @@
 from gzip import GzipFile
 
 from .bem import _bem_find_surface, read_bem_surfaces
+from .io import _loc_to_coil_trans
 from .io.constants import FIFF
 from .io.meas_info import read_fiducials
+from .io.pick import channel_type, _FNIRS_CH_TYPES_SPLIT, _MEG_CH_TYPES_SPLIT
 from .transforms import (apply_trans, invert_transform, combine_transforms,
-                         _ensure_trans, read_ras_mni_t, Transform)
+                         _ensure_trans, read_ras_mni_t, Transform, _get_trans,
+                         _frame_to_str)
 from .surface import read_surface, _read_mri_surface
 from .utils import (verbose, _validate_type, _check_fname, _check_option,
-                    get_subjects_dir, _require_version, logger)
+                    get_subjects_dir, _require_version, logger, warn)
 
 
 def _check_subject_dir(subject, subjects_dir):
@@ -28,7 +31,7 @@ def _check_subject_dir(subject, subjects_dir):
             raise ValueError('Freesurfer recon-all subject folder '
                              'is incorrect or improperly formatted, '
                              f'got {op.join(subjects_dir, subject)}')
-    return subjects_dir
+    return op.join(subjects_dir, subject)
 
 
 def _get_aseg(aseg, subject, subjects_dir):
@@ -402,7 +405,7 @@ def get_mni_fiducials(subject, subjects_dir=None, verbose=None):
 
     For more details about the coordinate systems and transformations involved,
     see https://surfer.nmr.mgh.harvard.edu/fswiki/CoordinateSystems and
-    :ref:`plot_source_alignment`.
+    :ref:`tut-source-alignment`.
     """
     # Eventually we might want to allow using the MNI Talairach with-skull
     # transformation rather than the standard brain-based MNI Talaranch
@@ -422,6 +425,120 @@ def get_mni_fiducials(subject, subjects_dir=None, verbose=None):
     return fids
 
 
+@verbose
+def estimate_head_mri_t(subject, subjects_dir=None, verbose=None):
+    """Estimate the head->mri transform from fsaverage fiducials.
+
+    A subject's fiducials can be estimated given a Freesurfer ``recon-all``
+    by transforming ``fsaverage`` fiducials using the inverse Talairach
+    transform, see :func:`mne.coreg.get_mni_fiducials`.
+
+    Parameters
+    ----------
+    %(subject)s
+    %(subjects_dir)s
+    %(verbose)s
+
+    Returns
+    -------
+    %(trans_not_none)s
+    """
+    from .channels.montage import make_dig_montage, compute_native_head_t
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    lpa, nasion, rpa = get_mni_fiducials(subject, subjects_dir)
+    montage = make_dig_montage(lpa=lpa['r'], nasion=nasion['r'], rpa=rpa['r'],
+                               coord_frame='mri')
+    return invert_transform(compute_native_head_t(montage))
+
+
+@verbose
+def _get_transforms_to_coord_frame(info, trans, coord_frame='mri',
+                                   verbose=None):
+    """Get the transforms to a coordinate frame from device, head and mri."""
+    head_mri_t = _get_trans(trans, 'head', 'mri')[0]
+    dev_head_t = _get_trans(info['dev_head_t'], 'meg', 'head')[0]
+    mri_dev_t = invert_transform(combine_transforms(
+        dev_head_t, head_mri_t, 'meg', 'mri'))
+    to_cf_t = dict(
+        meg=_ensure_trans([dev_head_t, mri_dev_t, Transform('meg', 'meg')],
+                          fro='meg', to=coord_frame),
+        head=_ensure_trans([dev_head_t, head_mri_t, Transform('head', 'head')],
+                           fro='head', to=coord_frame),
+        mri=_ensure_trans([head_mri_t, mri_dev_t, Transform('mri', 'mri')],
+                          fro='mri', to=coord_frame))
+    return to_cf_t
+
+
+@verbose
+def _ch_pos_in_coord_frame(info, to_cf_t, warn_meg=True, verbose=None):
+    """Transform positions from head/device/mri to a coordinate frame."""
+    from .forward import _create_meg_coils
+    from .viz._3d import _sensor_shape
+    chs = dict(ch_pos=dict(), sources=dict(), detectors=dict())
+    unknown_chs = list()  # prepare for chs with unknown coordinate frame
+    type_counts = dict()
+    for idx in range(info['nchan']):
+        ch_type = channel_type(info, idx)
+        if ch_type in type_counts:
+            type_counts[ch_type] += 1
+        else:
+            type_counts[ch_type] = 1
+        type_slices = dict(ch_pos=slice(0, 3))
+        if ch_type in _FNIRS_CH_TYPES_SPLIT:
+            # add sensors and detectors too for fNIRS
+            type_slices.update(sources=slice(3, 6), detectors=slice(6, 9))
+        for type_name, type_slice in type_slices.items():
+            if ch_type in _MEG_CH_TYPES_SPLIT + ('ref_meg',):
+                coil_trans = _loc_to_coil_trans(info['chs'][idx]['loc'])
+                coil = _create_meg_coils([info['chs'][idx]], acc='normal')[0]
+                # store verts as ch_coord
+                ch_coord, triangles = _sensor_shape(coil)
+                ch_coord = apply_trans(coil_trans, ch_coord)
+                if len(ch_coord) == 0 and warn_meg:
+                    warn(f'MEG sensor {info.ch_names[idx]} not found.')
+            else:
+                ch_coord = info['chs'][idx]['loc'][type_slice]
+            ch_coord_frame = info['chs'][idx]['coord_frame']
+            if ch_coord_frame not in (FIFF.FIFFV_COORD_UNKNOWN,
+                                      FIFF.FIFFV_COORD_DEVICE,
+                                      FIFF.FIFFV_COORD_HEAD,
+                                      FIFF.FIFFV_COORD_MRI):
+                raise RuntimeError(
+                    f'Channel {info.ch_names[idx]} has coordinate frame '
+                    f'{ch_coord_frame}, must be "meg", "head" or "mri".')
+            # set unknown as head first
+            if ch_coord_frame == FIFF.FIFFV_COORD_UNKNOWN:
+                unknown_chs.append(info.ch_names[idx])
+                ch_coord_frame = FIFF.FIFFV_COORD_HEAD
+            ch_coord = apply_trans(
+                to_cf_t[_frame_to_str[ch_coord_frame]], ch_coord)
+            if ch_type in _MEG_CH_TYPES_SPLIT + ('ref_meg',):
+                chs[type_name][info.ch_names[idx]] = (ch_coord, triangles)
+            else:
+                chs[type_name][info.ch_names[idx]] = ch_coord
+    if unknown_chs:
+        warn(f'Got coordinate frame "unknown" for {unknown_chs}, assuming '
+             '"head" coordinates.')
+    logger.info('Channel types::\t' + ', '.join(
+        [f'{ch_type}: {count}' for ch_type, count in type_counts.items()]))
+    return chs['ch_pos'], chs['sources'], chs['detectors']
+
+
+def _ensure_image_in_surface_RAS(image, subject, subjects_dir):
+    """Check if the image is in Freesurfer surface RAS space."""
+    import nibabel as nib
+    if not isinstance(image, nib.spatialimages.SpatialImage):
+        image = nib.load(image)
+    image = nib.MGHImage(image.dataobj.astype(np.float32), image.affine)
+    fs_img = nib.load(op.join(subjects_dir, subject, 'mri', 'brain.mgz'))
+    if not np.allclose(image.affine, fs_img.affine, atol=1e-6):
+        raise RuntimeError('The `image` is not aligned to Freesurfer '
+                           'surface RAS space. This space is required as '
+                           'it is the space where the anatomical '
+                           'segmentation and reconstructed surfaces are')
+    return image  # returns MGH image for header
+
+
 @verbose
 def read_talxfm(subject, subjects_dir=None, verbose=None):
     """Compute MRI-to-MNI transform from FreeSurfer talairach.xfm file.
@@ -632,3 +749,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/montage.py

@@ -392,7 +392,7 @@ def add_estimated_fiducials(self, subject, subjects_dir=None,
 
         See Also
         --------
-        :ref:`plot_source_alignment`
+        :ref:`tut-source-alignment`
 
         Notes
         -----

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/coreg.py

@@ -22,7 +22,8 @@
 from .io.constants import FIFF
 # keep get_mni_fiducials for backward compat (no burden to keep in this
 # namespace, too)
-from ._freesurfer import _read_mri_info, get_mni_fiducials  # noqa: F401
+from ._freesurfer import (_read_mri_info, get_mni_fiducials,  # noqa: F401
+                          estimate_head_mri_t)  # noqa: F401
 from .label import read_label, Label
 from .source_space import (add_source_space_distances, read_source_spaces,  # noqa: E501,F401
                            write_source_spaces)

mne/defaults.py

@@ -77,6 +77,8 @@
         head_color=(0.988, 0.89, 0.74),
         hpi_color=(1., 0., 1.),
         extra_color=(1., 1., 1.),
+        meg_color=(0., 0.25, 0.5), ref_meg_color=(0.5, 0.5, 0.5),
+        helmet_color=(0.0, 0.0, 0.6),
         eeg_color=(1., 0.596, 0.588), eegp_color=(0.839, 0.15, 0.16),
         ecog_color=(1., 1., 1.),
         dbs_color=(0.82, 0.455, 0.659),

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(