[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` to 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

@@ -28,7 +28,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 +402,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 +422,47 @@ 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))
+
+
+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 +673,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.')
+    subjects_dir = get_subjects_dir(subjects_dir, raise_error=True)
+    fname = _check_fname(op.join(subjects_dir, subject, 'bem',
+                                 surf + '_skull.surf'),
+                         overwrite='read', must_exist=True,
+                         name=f'{surf} skull surface')
+    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(

mne/surface.py

@@ -26,7 +26,7 @@
 from .io.pick import pick_types
 from .parallel import parallel_func
 from .transforms import (transform_surface_to, _pol_to_cart, _cart_to_sph,
-                         _get_trans, apply_trans, Transform,
+                         _get_trans, apply_trans, Transform, _frame_to_str,
                          apply_volume_registration)
 from .utils import (logger, verbose, get_subjects_dir, warn, _check_fname,
                     _check_option, _ensure_int, _TempDir, run_subprocess,
@@ -1736,27 +1736,11 @@ def _marching_cubes(image, level, smooth=0):
     return out
 
 
-def _get_surface_RAS_volumes(base_image, subject_from, subjects_dir):
-    """Check if the image is in Freesurfer surface RAS space."""
-    import nibabel as nib
-    if not isinstance(base_image, nib.spatialimages.SpatialImage):
-        base_image = nib.load(base_image)
-    base_image = nib.MGHImage(
-        _get_img_fdata(base_image).astype(np.float32), base_image.affine)
-    fs_t1 = nib.load(op.join(subjects_dir, subject_from, 'mri', 'T1.mgz'))
-    if not np.allclose(base_image.affine, fs_t1.affine, atol=1e-6):
-        raise RuntimeError('The `base_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 base_image, fs_t1
-
-
 def _warn_missing_chs(montage, dig_image, after_warp):
     """Warn that channels are missing."""
     # ensure that each electrode contact was marked in at least one voxel
     missing = set(np.arange(1, len(montage.ch_names) + 1)).difference(
-        set(np.unique(_get_img_fdata(dig_image))))
+        set(np.unique(np.asanyarray(dig_image.dataobj))))
     missing_ch = [montage.ch_names[idx - 1] for idx in missing]
     if missing_ch:
         warn('Channels ' + ', '.join(missing_ch) + ' were not assigned '
@@ -1822,7 +1806,7 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
     """
     _require_version('nibabel', 'SDR morph', '2.1.0')
     _require_version('dipy', 'SDR morph', '0.10.1')
-    from .channels import DigMontage
+    from .channels import DigMontage, make_dig_montage
     from ._freesurfer import _check_subject_dir
     import nibabel as nib
 
@@ -1840,28 +1824,40 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
     _check_subject_dir(subject_to, subjects_dir)
 
     # load image and make sure it's in surface RAS
-    image, fs_t1 = _get_surface_RAS_volumes(
-        base_image, subject_from, subjects_dir)
+    if not isinstance(base_image, nib.spatialimages.SpatialImage):
+        base_image = nib.load(base_image)
+    fs_from_img = nib.load(
+        op.join(subjects_dir, subject_from, 'mri', 'brain.mgz'))
+    if not np.allclose(base_image.affine, fs_from_img.affine, atol=1e-6):
+        raise RuntimeError('The `base_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')
 
     # get montage channel coordinates
     ch_dict = montage.get_positions()
     if ch_dict['coord_frame'] != 'mri':
+        bad_coord_frames = np.unique([d['coord_frame'] for d in montage.dig])
+        bad_coord_frames = ', '.join([
+            _frame_to_str[cf] if cf in _frame_to_str else str(cf)
+            for cf in bad_coord_frames])
         raise RuntimeError('Coordinate frame not supported, expected '
-                           '"mri", got ' + str(ch_dict['coord_frame']))
-    ch_coords = np.array(list(ch_dict['ch_pos'].values()))
+                           f'"mri", got {bad_coord_frames}')
+    ch_names = list(ch_dict['ch_pos'].keys())
+    ch_coords = np.array([ch_dict['ch_pos'][name] for name in ch_names])
 
     # convert to freesurfer voxel space
     ch_coords = apply_trans(
-        np.linalg.inv(fs_t1.header.get_vox2ras_tkr()), ch_coords * 1000)
+        np.linalg.inv(fs_from_img.header.get_vox2ras_tkr()), ch_coords * 1000)
 
     # take channel coordinates and use the image to transform them
     # into a volume where all the voxels over a threshold nearby
     # are labeled with an index
-    image_data = _get_img_fdata(image)
+    image_data = _get_img_fdata(base_image)
     if use_min:
         image_data *= -1
     thresh = np.quantile(image_data, thresh)
-    image_from = np.zeros(image.shape, dtype=int)
+    image_from = np.zeros(base_image.shape, dtype=int)
     for i, ch_coord in enumerate(ch_coords):
         # this looks up to a voxel away, it may be marked imperfectly
         volume = _voxel_neighbors(ch_coord, image_data, thresh=thresh,
@@ -1880,7 +1876,7 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
                 image_from[voxel] = i + 1
 
     # apply the mapping
-    image_from = nib.Nifti1Image(image_from, fs_t1.affine)
+    image_from = nib.Nifti1Image(image_from, fs_from_img.affine)
     _warn_missing_chs(montage, image_from, after_warp=False)
 
     template_brain = nib.load(
@@ -1893,25 +1889,20 @@ def warp_montage_volume(montage, base_image, reg_affine, sdr_morph,
     _warn_missing_chs(montage, image_to, after_warp=True)
 
     # recover the contact positions as the center of mass
-    warped_data = _get_img_fdata(image_to)
+    warped_data = np.asanyarray(image_to.dataobj)
     for val, ch_coord in enumerate(ch_coords, 1):
-        ch_coord[:] = np.array(
+        ch_coord[:] = np.asanyarray(
             np.where(warped_data == val), float).mean(axis=1)
 
-    # convert back to surface RAS
+    # convert back to surface RAS of the template
+    fs_to_img = nib.load(
+        op.join(subjects_dir, subject_to, 'mri', 'brain.mgz'))
     ch_coords = apply_trans(
-        fs_t1.header.get_vox2ras_tkr(), ch_coords) / 1000
-
-    # copy before modifying
-    montage_warped = montage.copy()
-
-    # modify montage to be returned
-    idx = 0
-    for dig_point in montage_warped.dig:
-        if dig_point['kind'] == FIFF.FIFFV_POINT_EEG:
-            assert dig_point['ident'] == idx + 1  # 1 indexed
-            dig_point['r'] = ch_coords[idx]
-            idx += 1
+        fs_to_img.header.get_vox2ras_tkr(), ch_coords) / 1000
+
+    # make warped montage
+    montage_warped = make_dig_montage(
+        dict(zip(ch_names, ch_coords)), coord_frame='mri')
     return montage_warped, image_from, image_to