MRG, ENH: Make use of fiducials (mne-tools#6791)

Add fiducials to all standard templates
Update coil files

Author: larsoner

.travis.yml

@@ -175,8 +175,8 @@ script:
         pip install -e .;
         python mne/tests/test_evoked.py;
       fi;
-    - echo pytest -m "${CONDITION}" ${USE_DIRS}
-    - pytest -m "${CONDITION}" ${USE_DIRS}
+    - echo pytest -m "${CONDITION}" --cov=mne ${USE_DIRS}
+    - pytest -m "${CONDITION}" --cov=mne ${USE_DIRS}
     - if [ "${DEPS}" == "nodata" ]; then
         make pep;
       fi;

azure-pipelines.yml

@@ -76,7 +76,7 @@ jobs:
     displayName: 'Print config'
   - script: python -c "import mne; mne.datasets.testing.data_path(verbose=True)"
     displayName: 'Get test data'
-  - script: pytest -m "not ultraslowtest" mne
+  - script: pytest -m "not ultraslowtest" --cov=mne mne
     displayName: 'Run tests'
   - script: codecov --root %BUILD_REPOSITORY_LOCALPATH% -t %CODECOV_TOKEN%
     displayName: 'Codecov'

doc/python_reference.rst

@@ -190,6 +190,7 @@ Datasets
    fetch_hcp_mmp_parcellation
    hf_sef.data_path
    kiloword.data_path
+   limo.load_data
    misc.data_path
    mtrf.data_path
    multimodal.data_path

examples/visualization/plot_montage.py

@@ -20,55 +20,39 @@
 from mne.datasets import fetch_fsaverage
 from mne.viz import set_3d_title, set_3d_view
 
-import warnings
-
-warnings.simplefilter("ignore")
 
 ###############################################################################
 # Check all montages against a sphere
-#
-
-sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.085)
 
 for current_montage in get_builtin_montages():
-
     montage = mne.channels.make_standard_montage(current_montage)
+    info = mne.create_info(
+        ch_names=montage.ch_names, sfreq=100., ch_types='eeg', montage=montage)
+    sphere = mne.make_sphere_model(r0='auto', head_radius='auto', info=info)
     fig = mne.viz.plot_alignment(
         # Plot options
-        show_axes=True, dig=True, surfaces='head', bem=sphere,
-
-        # Create dummy info
-        info=mne.create_info(
-            ch_names=montage.ch_names,
-            sfreq=1,
-            ch_types='eeg',
-            montage=montage,
-        ),
-    )
+        show_axes=True, dig='fiducials', surfaces='head',
+        bem=sphere, info=info)
     set_3d_view(figure=fig, azimuth=135, elevation=80)
     set_3d_title(figure=fig, title=current_montage)
 
 
 ###############################################################################
 # Check all montages against fsaverage
-#
 
 subjects_dir = op.dirname(fetch_fsaverage())
 
 for current_montage in get_builtin_montages():
     montage = mne.channels.make_standard_montage(current_montage)
+    # Create dummy info
+    info = mne.create_info(
+        ch_names=montage.ch_names, sfreq=100., ch_types='eeg', montage=montage)
     fig = mne.viz.plot_alignment(
         # Plot options
-        show_axes=True, dig=True, surfaces='head', trans=None,
-        subject='fsaverage', subjects_dir=subjects_dir,
-
-        # Create dummy info
-        info=mne.create_info(
-            ch_names=montage.ch_names,
-            sfreq=1,
-            ch_types='eeg',
-            montage=montage,
-        ),
+        show_axes=True, dig='fiducials', surfaces='head', mri_fiducials=True,
+        subject='fsaverage', subjects_dir=subjects_dir, info=info,
+        coord_frame='mri',
+        trans='fsaverage',  # transform from head coords to fsaverage's MRI
     )
     set_3d_view(figure=fig, azimuth=135, elevation=80)
     set_3d_title(figure=fig, title=current_montage)

mne/_digitization/_utils.py

@@ -7,6 +7,7 @@
 #
 # License: BSD (3-clause)
 
+from collections import OrderedDict
 import datetime
 import os.path as op
 import re
@@ -337,7 +338,9 @@ def _make_dig_points(nasion=None, lpa=None, rpa=None, hpi=None,
                         'kind': FIFF.FIFFV_POINT_EXTRA,
                         'coord_frame': coord_frame})
     if dig_ch_pos is not None:
-        keys = sorted(dig_ch_pos.keys())
+        keys = dig_ch_pos.keys()
+        if not isinstance(dig_ch_pos, OrderedDict):
+            keys = sorted(keys)
         try:  # use the last 3 as int if possible (e.g., EEG001->1)
             idents = []
             for key in keys:

mne/channels/_standard_montage_utils.py

@@ -14,11 +14,13 @@
 
 MONTAGE_PATH = op.join(op.dirname(_CHANNELS_INIT_FILE), 'data', 'montages')
 
-HEAD_SIZE_DEFAULT = 0.085  # in [m]
 _str = 'U100'
 
 
-def _egi_256():
+# In standard_1020, T9=LPA, T10=RPA, Nasion is the same as Iz with a
+# sign-flipped Y value
+
+def _egi_256(head_size):
     fname = op.join(MONTAGE_PATH, 'EGI_256.csd')
     # Label, Theta, Phi, Radius, X, Y, Z, off sphere surface
     options = dict(comments='//',
@@ -27,45 +29,66 @@ def _egi_256():
     pos = np.stack([xs, ys, zs], axis=-1)
 
     # Fix pos to match Montage code
-    # pos -= np.mean(pos, axis=0)
-    pos = HEAD_SIZE_DEFAULT * (pos / np.linalg.norm(pos, axis=1).mean())
+    pos *= head_size / np.median(np.linalg.norm(pos, axis=1))
+
+    # For this cap, the Nasion is the frontmost electrode,
+    # LPA/RPA we approximate by putting 75% of the way (toward the front)
+    # between the two electrodes that are halfway down the ear holes
+    nasion = pos[ch_names.index('E31')]
+    lpa = (0.75 * pos[ch_names.index('E67')] +
+           0.25 * pos[ch_names.index('E94')])
+    rpa = (0.75 * pos[ch_names.index('E219')] +
+           0.25 * pos[ch_names.index('E190')])
 
     return make_dig_montage(
         ch_pos=OrderedDict(zip(ch_names, pos)),
-        coord_frame='head',
+        coord_frame='unknown', nasion=nasion, lpa=lpa, rpa=rpa,
     )
 
 
-def _easycap(basename):
+def _easycap(basename, head_size):
     fname = op.join(MONTAGE_PATH, basename)
     options = dict(skip_header=1, dtype=(_str, 'i4', 'i4'))
     ch_names, theta, phi = _safe_np_loadtxt(fname, **options)
 
-    radii = np.full_like(phi, 1)  # XXX: HEAD_SIZE_DEFAULT should work
+    radii = np.full(len(phi), head_size)
     pos = _sph_to_cart(np.stack(
         [radii, np.deg2rad(phi), np.deg2rad(theta)],
         axis=-1,
     ))
-
-    # scale up to realistic head radius (8.5cm == 85mm):
-    pos *= HEAD_SIZE_DEFAULT  # XXXX: this should be done through radii
+    nasion = np.concatenate([[0],  pos[ch_names.index('Fpz'), 1:]])
+    nasion *= head_size / np.linalg.norm(nasion)
+    lpa = np.mean([pos[ch_names.index('FT9')],
+                   pos[ch_names.index('TP9')]], axis=0)
+    lpa *= head_size / np.linalg.norm(lpa)  # on sphere
+    rpa = np.mean([pos[ch_names.index('FT10')],
+                   pos[ch_names.index('TP10')]], axis=0)
+    rpa *= head_size / np.linalg.norm(rpa)
 
     return make_dig_montage(
         ch_pos=OrderedDict(zip(ch_names, pos)),
-        coord_frame='head',
+        coord_frame='unknown', nasion=nasion, lpa=lpa, rpa=rpa,
     )
 
 
-def _hydrocel(basename, fid_names=('FidNz', 'FidT9', 'FidT10')):
+def _hydrocel(basename, head_size):
+    fid_names = ('FidNz', 'FidT9', 'FidT10')
     fname = op.join(MONTAGE_PATH, basename)
     options = dict(dtype=(_str, 'f4', 'f4', 'f4'))
     ch_names, xs, ys, zs = _safe_np_loadtxt(fname, **options)
 
-    pos = np.stack([xs, ys, zs], axis=-1) * 0.01
+    pos = np.stack([xs, ys, zs], axis=-1)
     ch_pos = OrderedDict(zip(ch_names, pos))
-    _ = [ch_pos.pop(n, None) for n in fid_names]
+    nasion, lpa, rpa = [ch_pos.pop(n) for n in fid_names]
+    scale = head_size / np.median(np.linalg.norm(pos, axis=-1))
+    for value in ch_pos.values():
+        value *= scale
+    nasion *= scale
+    lpa *= scale
+    rpa *= scale
 
-    return make_dig_montage(ch_pos=ch_pos, coord_frame='head')
+    return make_dig_montage(ch_pos=ch_pos, coord_frame='unknown',
+                            nasion=nasion, rpa=rpa, lpa=lpa)
 
 
 def _str_names(ch_names):
@@ -79,39 +102,32 @@ def _safe_np_loadtxt(fname, **kwargs):
     return (ch_names,) + others
 
 
-def _biosemi(basename, fid_names=('Nz', 'LPA', 'RPA')):
+def _biosemi(basename, head_size):
+    fid_names = ('Nz', 'LPA', 'RPA')
     fname = op.join(MONTAGE_PATH, basename)
     options = dict(skip_header=1, dtype=(_str, 'i4', 'i4'))
     ch_names, theta, phi = _safe_np_loadtxt(fname, **options)
 
-    radii = np.full_like(phi, 1)  # XXX: HEAD_SIZE_DEFAULT should work
+    radii = np.full(len(phi), head_size)
     pos = _sph_to_cart(np.stack(
         [radii, np.deg2rad(phi), np.deg2rad(theta)],
         axis=-1,
     ))
 
-    # scale up to realistic head radius (8.5cm == 85mm):
-    pos *= HEAD_SIZE_DEFAULT  # XXXX: this should be done through radii
-
     ch_pos = OrderedDict(zip(ch_names, pos))
-    _ = [ch_pos.pop(n, None) for n in fid_names]
+    nasion, lpa, rpa = [ch_pos.pop(n) for n in fid_names]
 
-    return make_dig_montage(ch_pos=ch_pos, coord_frame='head')
+    return make_dig_montage(ch_pos=ch_pos, coord_frame='unknown',
+                            nasion=nasion, lpa=lpa, rpa=rpa)
 
 
-def _mgh_or_standard(basename, fid_names=('Nz', 'LPA', 'RPA')):
+def _mgh_or_standard(basename, head_size):
+    fid_names = ('Nz', 'LPA', 'RPA')
     fname = op.join(MONTAGE_PATH, basename)
 
     ch_names_, pos = [], []
     with open(fname) as fid:
-        # Default units are meters
-        for line in fid:
-            if 'UnitPosition' in line:
-                units = line.split()[1]
-                scale_factor = dict(m=1., mm=1e-3)[units]
-                break
-        else:
-            raise RuntimeError('Could not detect units in file %s' % fname)
+        # Ignore units as we will scale later using the norms anyway
         for line in fid:
             if 'Positions\n' in line:
                 break
@@ -125,12 +141,18 @@ def _mgh_or_standard(basename, fid_names=('Nz', 'LPA', 'RPA')):
                 break
             ch_names_.append(line.strip(' ').strip('\n'))
 
-    pos = np.array(pos) * scale_factor
-
+    pos = np.array(pos)
     ch_pos = OrderedDict(zip(ch_names_, pos))
-    _ = [ch_pos.pop(n, None) for n in fid_names]
-
-    return make_dig_montage(ch_pos=ch_pos, coord_frame='head')
+    nasion, lpa, rpa = [ch_pos.pop(n) for n in fid_names]
+    scale = head_size / np.median(np.linalg.norm(pos, axis=1))
+    for value in ch_pos.values():
+        value *= scale
+    nasion *= scale
+    lpa *= scale
+    rpa *= scale
+
+    return make_dig_montage(ch_pos=ch_pos, coord_frame='unknown',
+                            nasion=nasion, lpa=lpa, rpa=rpa)
 
 
 standard_montage_look_up_table = {

mne/channels/montage.py

@@ -50,6 +50,7 @@
 
 from .channels import DEPRECATED_PARAM
 
+HEAD_SIZE_DEFAULT = 0.095  # in [m]
 
 _BUILT_IN_MONTAGES = [
     'EGI_256',
@@ -840,6 +841,7 @@ def dig_ch_pos(self):
 
     def _get_ch_pos(self):
         pos = [d['r'] for d in _get_dig_eeg(self.dig)]
+        assert len(self.ch_names) == len(pos)
         return dict(zip(self.ch_names, pos))
 
     def _get_dig_names(self):
@@ -1718,7 +1720,7 @@ def compute_dev_head_t(montage):
     return Transform(fro='meg', to='head', trans=trans)
 
 
-def make_standard_montage(kind):
+def make_standard_montage(kind, head_size=HEAD_SIZE_DEFAULT):
     """Read a generic (built-in) montage.
 
     Individualized (digitized) electrode positions should be read in using
@@ -1729,6 +1731,8 @@ def make_standard_montage(kind):
     kind : str
         The name of the montage file without the file extension (e.g.
         kind='easycap-M10' for 'easycap-M10.txt'). See notes for valid kinds.
+    head_size : float
+        The head size (in meters) to use for spherical montages.
 
     Returns
     -------
@@ -1795,5 +1799,4 @@ def make_standard_montage(kind):
         raise ValueError('Could not find the montage %s. Please provide one '
                          'among: %s' % (kind,
                                         standard_montage_look_up_table.keys()))
-    else:
-        return standard_montage_look_up_table[kind]()
+    return standard_montage_look_up_table[kind](head_size=head_size)

mne/channels/tests/test_montage.py

@@ -1034,7 +1034,7 @@ def test_set_montage():
     assert ((orig_pos != new_pos).all())
 
     r0 = _fit_sphere(new_pos)[1]
-    assert_allclose(r0, [0., -0.016, 0.], atol=1e-3)
+    assert_allclose(r0, [0.000775, 0.006881, 0.047398], atol=1e-3)
 
 
 # XXX: this does not check ch_names + it cannot work because of write_dig