MRG, BUG: Fix M/EEG topomap plotting (mne-tools#7066)

* MAINT: Refactor to route through _find

* BUG: Fix plotting of sensors relative to head

* FIX: Remove problematic workaround

* FIX: Org

* FIX: Comparison

* ENH: Use sphere [circle full]

* FIX: Units [circle full]

* API: Cleaner deprecation [circle full]

* FIX: Debugging cruft [ci skip]

* FIX: Example

* FIX: Eradicate layout [circle full]

* FIX: Standardize [circle full]

* FIX: Fix for Travis

* BUG: Fix projs topomap [circle full]

* FIX: Fix for sphere and draw [circle full]

* FIX: Warning

* revise sensor locations tutorial

* fix crossref

* FIX: Better min/max [circle full]

* FIX: Fix default [circle full]

* FIX: Test

* FIX: Pass info [circle full]

* FIX: Fix animation

* FIX: Example

* FIX: Example

* DOC: Many fixes [circle full]

* FIX: Fix for changed example [skip travis]

doc/changes/latest.inc

@@ -126,6 +126,8 @@ Bug
 
 - Fix handling of repeated events in :class:`mne.Epochs` by `Fahimeh Mamashli`_ and `Alex Gramfort`_
 
+- Fix many bugs with plotting sensors overlaid on a head outline. All plotting is now done in head coordinates and scaled by ``head_radius``, which defaults to 0.095, by `Eric Larson`_
+
 - Fix :func:`mne.io.anonymize_info` to allow shifting dates of service and to match anticipated changes in mne-cpp by `Luke Bloy`_
 
 - Fix reading of cardinals in .htps files (identifier are int not strings) by `Alex Gramfort`_
@@ -155,9 +157,13 @@ API
 
 - :meth:`mne.Epochs.shift_time` and :meth:`mne.Evoked.shift_time` now allow shifting times by arbitrary amounts (previously only by integer multiples of the sampling period), by `Daniel McCloy`_ and `Eric Larson`_.
 
+- The ``head_pos`` argument of :func:`mne.Evoked.plot_topomap` and related functions has been deprecated in favor of ``head_radius``, by `Eric Larson`_.
+
+- The ``layout`` argument to topomap-related functions such as :meth:`mne.Evoked.plot_topomap` and :func:`mne.viz.plot_tfr_topomap` has been deprecated in favor of channel-position based flattening based on the ``info`` and ``sphere`` argument, by `Eric Larson`_.
+
 - The APIs of :meth:`mne.io.Raw.plot_projs_topomap`, :meth:`mne.Epochs.plot_projs_topomap` and :meth:`mne.Evoked.plot_projs_topomap` are now more similar to :func:`mne.viz.plot_projs_topomap` by `Daniel McCloy`_.
 
-- :func:`mne.viz.plot_projs_topomap` now accepts both :class:`~mne.channels.Layout` and :class:`~mne.Info` (previously one or the other), and will ignore :class:`~mne.Info` if :class:`~mne.channels.Layout` is provided `Daniel McCloy`_.
+- The function :func:`mne.setup_volume_source_space` has a ``sphere_units`` argument that defaults to ``'mm'`` in 0.20 but will change to ``'m'`` in 0.21, set it to avoid a warning by `Eric Larson`_.
 
 - :func:`mne.viz.plot_projs_topomap` and the related methods :meth:`mne.io.Raw.plot_projs_topomap`, :meth:`mne.Epochs.plot_projs_topomap` and :meth:`mne.Evoked.plot_projs_topomap` now accept parameter ``vlim`` to control the colormap, with keyword ``'joint'`` computing the colormap jointly across all projectors of a given channel type, by `Daniel McCloy`_.
 

doc/python_reference.rst

@@ -185,6 +185,7 @@ Datasets
    brainstorm.bst_resting.data_path
    brainstorm.bst_raw.data_path
    eegbci.load_data
+   eegbci.standardize
    fetch_aparc_sub_parcellation
    fetch_fsaverage
    fetch_hcp_mmp_parcellation

examples/decoding/plot_decoding_csp_eeg.py

@@ -37,7 +37,7 @@
 from sklearn.model_selection import ShuffleSplit, cross_val_score
 
 from mne import Epochs, pick_types, events_from_annotations
-from mne.channels import read_layout
+from mne.channels import make_standard_montage
 from mne.io import concatenate_raws, read_raw_edf
 from mne.datasets import eegbci
 from mne.decoding import CSP
@@ -56,6 +56,9 @@
 
 raw_fnames = eegbci.load_data(subject, runs)
 raw = concatenate_raws([read_raw_edf(f, preload=True) for f in raw_fnames])
+eegbci.standardize(raw)  # set channel names
+montage = make_standard_montage('standard_1005')
+raw.set_montage(montage)
 
 # strip channel names of "." characters
 raw.rename_channels(lambda x: x.strip('.'))
@@ -102,9 +105,7 @@
 # plot CSP patterns estimated on full data for visualization
 csp.fit_transform(epochs_data, labels)
 
-layout = read_layout('EEG1005')
-csp.plot_patterns(epochs.info, layout=layout, ch_type='eeg',
-                  units='Patterns (AU)', size=1.5)
+csp.plot_patterns(epochs.info, ch_type='eeg', units='Patterns (AU)', size=1.5)
 
 ###############################################################################
 # Look at performance over time

examples/decoding/plot_decoding_spoc_CMC.py

@@ -31,7 +31,6 @@
 from mne import Epochs
 from mne.decoding import SPoC
 from mne.datasets.fieldtrip_cmc import data_path
-from mne.channels import read_layout
 
 from sklearn.pipeline import make_pipeline
 from sklearn.linear_model import Ridge
@@ -87,5 +86,4 @@
 # Plot the contributions to the detected components (i.e., the forward model)
 
 spoc.fit(X, y)
-layout = read_layout('CTF151.lay')
-spoc.plot_patterns(meg_epochs.info, layout=layout)
+spoc.plot_patterns(meg_epochs.info)

examples/decoding/plot_decoding_xdawn_eeg.py

@@ -7,17 +7,6 @@
 spatial Xdawn filters are trained and applied on the signal. Channels are
 concatenated and rescaled to create features vectors that will be fed into
 a logistic regression.
-
-References
-----------
-.. [1] Rivet, B., Souloumiac, A., Attina, V., & Gibert, G. (2009). xDAWN
-       algorithm to enhance evoked potentials: application to brain-computer
-       interface. Biomedical Engineering, IEEE Transactions on, 56(8),
-       2035-2043.
-.. [2] Rivet, B., Cecotti, H., Souloumiac, A., Maby, E., & Mattout, J. (2011,
-       August). Theoretical analysis of xDAWN algorithm: application to an
-       efficient sensor selection in a P300 BCI. In Signal Processing
-       Conference, 2011 19th European (pp. 1382-1386). IEEE.
 """
 # Authors: Alexandre Barachant <[email protected]>
 #
@@ -36,7 +25,6 @@
 from mne.datasets import sample
 from mne.preprocessing import Xdawn
 from mne.decoding import Vectorizer
-from mne.viz import tight_layout
 
 
 print(__doc__)
@@ -48,7 +36,8 @@
 raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
 event_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw-eve.fif'
 tmin, tmax = -0.1, 0.3
-event_id = dict(aud_l=1, aud_r=2, vis_l=3, vis_r=4)
+event_id = {'Auditory/Left': 1, 'Auditory/Right': 2,
+            'Visual/Left': 3, 'Visual/Right': 4}
 n_filter = 3
 
 # Setup for reading the raw data
@@ -92,37 +81,43 @@
 cm_normalized = cm.astype(float) / cm.sum(axis=1)[:, np.newaxis]
 
 # Plot confusion matrix
-plt.imshow(cm_normalized, interpolation='nearest', cmap=plt.cm.Blues)
-plt.title('Normalized Confusion matrix')
-plt.colorbar()
+fig, ax = plt.subplots(1)
+im = ax.imshow(cm_normalized, interpolation='nearest', cmap=plt.cm.Blues)
+ax.set(title='Normalized Confusion matrix')
+fig.colorbar(im)
 tick_marks = np.arange(len(target_names))
 plt.xticks(tick_marks, target_names, rotation=45)
 plt.yticks(tick_marks, target_names)
-tight_layout()
-plt.ylabel('True label')
-plt.xlabel('Predicted label')
-plt.show()
-
+fig.tight_layout()
+ax.set(ylabel='True label', xlabel='Predicted label')
 
 ###############################################################################
 # The ``patterns_`` attribute of a fitted Xdawn instance (here from the last
 # cross-validation fold) can be used for visualization.
 
-fig, axes = plt.subplots(nrows=len(event_id), ncols=n_filter)
+fig, axes = plt.subplots(nrows=len(event_id), ncols=n_filter,
+                         figsize=(n_filter, len(event_id) * 2))
 fitted_xdawn = clf.steps[0][1]
-for i, (cur_class, cur_patterns) in enumerate(fitted_xdawn.patterns_.items()):
-    tmp_info = epochs.info.copy()
-    tmp_info['sfreq'] = 1.
+tmp_info = epochs.info.copy()
+tmp_info['sfreq'] = 1.
+for ii, cur_class in enumerate(sorted(event_id)):
+    cur_patterns = fitted_xdawn.patterns_[cur_class]
     pattern_evoked = EvokedArray(cur_patterns[:n_filter].T, tmp_info, tmin=0)
     pattern_evoked.plot_topomap(
-        times=np.arange(n_filter), ch_type=None,
-        scalings=None,
-        time_format='{} / comp. %01d'.format(cur_class),
-        colorbar=False,
-        head_pos={'center': [0, 0]},
-        show_names=False,
-        axes=axes[i, :],
-        extrapolate='head',
-        show=False)
-fig.subplots_adjust(hspace=0.3)
-plt.show()
+        times=np.arange(n_filter),
+        time_format='Component %d' if ii == 0 else '', colorbar=False,
+        show_names=False, axes=axes[ii], show=False)
+    axes[ii, 0].set(ylabel=cur_class)
+fig.tight_layout(h_pad=1.0, w_pad=1.0, pad=0.1)
+
+###############################################################################
+# References
+# ----------
+# .. [1] Rivet, B., Souloumiac, A., Attina, V., & Gibert, G. (2009). xDAWN
+#        algorithm to enhance evoked potentials: application to brain-computer
+#        interface. Biomedical Engineering, IEEE Transactions on, 56(8),
+#        2035-2043.
+# .. [2] Rivet, B., Cecotti, H., Souloumiac, A., Maby, E., & Mattout, J. (2011,
+#        August). Theoretical analysis of xDAWN algorithm: application to an
+#        efficient sensor selection in a P300 BCI. In Signal Processing
+#        Conference, 2011 19th European (pp. 1382-1386). IEEE.

examples/forward/plot_left_cerebellum_volume_source.py

@@ -33,10 +33,10 @@
 
 # setup a volume source space of the left cerebellum cortex
 volume_label = 'Left-Cerebellum-Cortex'
-sphere = (0, 0, 0, 120)
-lh_cereb = setup_volume_source_space(subject, mri=aseg_fname, sphere=sphere,
-                                     volume_label=volume_label,
-                                     subjects_dir=subjects_dir)
+sphere = (0, 0, 0, 0.12)
+lh_cereb = setup_volume_source_space(
+    subject, mri=aseg_fname, sphere=sphere, volume_label=volume_label,
+    subjects_dir=subjects_dir, sphere_units='m')
 
 # Combine the source spaces
 src = surf + lh_cereb

examples/io/plot_read_proj.py

@@ -40,14 +40,15 @@
 
 ###############################################################################
 # Display the projections one by one
-fig, axes = plt.subplots(1, len(empty_room_proj))
+n_cols = len(empty_room_proj)
+fig, axes = plt.subplots(1, n_cols, figsize=(2 * n_cols, 2))
 for proj, ax in zip(empty_room_proj, axes):
-    proj.plot_topomap(axes=ax)
+    proj.plot_topomap(axes=ax, info=raw.info)
 
 ###############################################################################
 # Use the function in `mne.viz` to display a list of projections
 assert isinstance(empty_room_proj, list)
-mne.viz.plot_projs_topomap(empty_room_proj)
+mne.viz.plot_projs_topomap(empty_room_proj, info=raw.info)
 
 ###############################################################################
 # .. TODO: add this when the tutorial is up: "As shown in the tutorial
@@ -61,32 +62,3 @@
 # add them to raw and plot everything
 raw.add_proj(ecg_projs)
 raw.plot_projs_topomap()
-
-###############################################################################
-# Displaying the projections from a raw object requires no extra information
-# since all the layout information is present in `raw.info`.
-# MNE is able to automatically determine the layout for some magnetometer and
-# gradiometer configurations but not the layout of EEG electrodes.
-#
-# Here we display the `ecg_projs` individually and we provide extra parameters
-# for EEG. (Notice that planar projection refers to the gradiometers and axial
-# refers to magnetometers.)
-#
-# Notice that the conditional is just for illustration purposes. We could
-# `raw.info` in all cases to avoid the guesswork in `plot_topomap` and ensure
-# that the right layout is always found
-fig, axes = plt.subplots(1, len(ecg_projs))
-for proj, ax in zip(ecg_projs, axes):
-    if proj['desc'].startswith('ECG-eeg'):
-        proj.plot_topomap(axes=ax, info=raw.info)
-    else:
-        proj.plot_topomap(axes=ax)
-
-###############################################################################
-# The correct layout or a list of layouts from where to choose can also be
-# provided. Just for illustration purposes, here we generate the
-# `possible_layouts` from the raw object itself, but it can come from somewhere
-# else.
-possible_layouts = [mne.find_layout(raw.info, ch_type=ch_type)
-                    for ch_type in ('grad', 'mag', 'eeg')]
-mne.viz.plot_projs_topomap(ecg_projs, layout=possible_layouts)

examples/visualization/plot_evoked_topomap.py

@@ -19,6 +19,7 @@
 
 import numpy as np
 import matplotlib.pyplot as plt
+
 from mne.datasets import sample
 from mne import read_evokeds
 
@@ -80,12 +81,15 @@
 # Compare this with ``extrapolate='head'`` (second topography below) where
 # extrapolation goes to 0 at the head circle and ``extrapolate='local'`` where
 # extrapolation is performed only within some distance from channels:
+
 extrapolations = ['box', 'head', 'local']
 fig, axes = plt.subplots(figsize=(7.5, 2.5), ncols=3)
 
+# Here we look at EEG channels, and use a custom head sphere to get all the
+# sensors to be well within the drawn head surface
 for ax, extr in zip(axes, extrapolations):
-    evoked.plot_topomap(0.1, ch_type='mag', size=2, extrapolate=extr, axes=ax,
-                        show=False, colorbar=False)
+    evoked.plot_topomap(0.1, ch_type='eeg', size=2, extrapolate=extr, axes=ax,
+                        show=False, colorbar=False, sphere=(0., 0., 0., 0.09))
     ax.set_title(extr, fontsize=14)
 
 ###############################################################################

mne/beamformer/tests/test_lcmv.py

@@ -287,8 +287,9 @@ def test_make_lcmv(tmpdir, reg, proj):
     # selection and rank reduction of the leadfield
     sphere = mne.make_sphere_model(r0=(0., 0., 0.), head_radius=0.080)
     src = mne.setup_volume_source_space(subject=None, pos=15., mri=None,
-                                        sphere=(0.0, 0.0, 0.0, 80.0),
-                                        bem=None, mindist=5.0, exclude=2.0)
+                                        sphere=(0.0, 0.0, 0.0, 0.08),
+                                        bem=None, mindist=5.0, exclude=2.0,
+                                        sphere_units='m')
 
     fwd_sphere = mne.make_forward_solution(evoked.info, trans=None, src=src,
                                            bem=sphere, eeg=False, meg=True)

mne/beamformer/tests/test_rap_music.py

@@ -157,7 +157,8 @@ def test_rap_music_sphere():
     sphere = mne.make_sphere_model(r0=(0., 0., 0.04))
     src = mne.setup_volume_source_space(subject=None, pos=10.,
                                         sphere=(0.0, 0.0, 40, 65.0),
-                                        mindist=5.0, exclude=0.0)
+                                        mindist=5.0, exclude=0.0,
+                                        sphere_units='mm')
     forward = mne.make_forward_solution(evoked.info, trans=None, src=src,
                                         bem=sphere)
 

mne/bem.py

@@ -839,16 +839,15 @@ def fit_sphere_to_headshape(info, dig_kinds='auto', units='m', verbose=None):
     ----------
     info : instance of Info
         Measurement info.
-    dig_kinds : list of str | str
-        Kind of digitization points to use in the fitting. These can be any
-        combination of ('cardinal', 'hpi', 'eeg', 'extra'). Can also
-        be 'auto' (default), which will use only the 'extra' points if
-        enough (more than 10) are available, and if not, uses 'extra' and
-        'eeg' points.
+    %(dig_kinds)s
     units : str
         Can be "m" (default) or "mm".
 
         .. versionadded:: 0.12
+    move_origin : bool
+        If True, allow the origin to vary. Otherwise, fix it at (0, 0, 0).
+
+        .. versionadded:: 0.20
 
     %(verbose)s
 
@@ -1004,8 +1003,7 @@ def constraint(center_rad):
     x_opt = fmin_cobyla(cost_fun, x0, constraint, rhobeg=radius_init,
                         rhoend=radius_init * 1e-6, disp=disp)
 
-    origin = x_opt[:3]
-    radius = x_opt[3]
+    origin, radius = x_opt[:3], x_opt[3]
     return radius, origin
 
 

mne/channels/__init__.py

@@ -3,6 +3,7 @@
 Can be used for setting of sensor locations used for processing and plotting.
 """
 
+from ..defaults import HEAD_SIZE_DEFAULT
 from .layout import (Layout, make_eeg_layout, make_grid_layout, read_layout,
                      find_layout, generate_2d_layout)
 from .montage import (DigMontage,
@@ -11,8 +12,7 @@
                       read_dig_polhemus_isotrak, read_polhemus_fastscan,
                       compute_dev_head_t, make_standard_montage,
                       read_custom_montage, read_dig_hpts,
-                      compute_native_head_t,
-                      )
+                      compute_native_head_t)
 from .channels import (equalize_channels, rename_channels, fix_mag_coil_types,
                        read_ch_connectivity, _get_ch_type,
                        find_ch_connectivity, make_1020_channel_selections)