[MRG+1] BUG: projection vectors (mne-tools#5641)

doc/whats_new.rst

@@ -89,6 +89,7 @@ Changelog
 
 Bug
 ~~~
+- Fix :func:`mne.io.Raw.plot_proj_topomap` by `Joan Massich`_
 
 - Fix reading edf file annotations by `Joan Massich`_
 

examples/io/plot_read_proj.py

@@ -0,0 +1,91 @@
+
+"""
+==============================================
+Read and visualize projections (SSP and other)
+==============================================
+
+This example shows how to read and visualize Signal Subspace Projectors (SSP)
+vector. Such projections are sometimes referred to as PCA projections.
+"""
+
+# Author: Joan Massich <[email protected]>
+#
+# License: BSD (3-clause)
+
+import matplotlib.pyplot as plt
+
+import mne
+from mne import read_proj
+from mne.io import read_raw_fif
+
+from mne.datasets import sample
+
+print(__doc__)
+
+data_path = sample.data_path()
+
+subjects_dir = data_path + '/subjects'
+fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+ecg_fname = data_path + '/MEG/sample/sample_audvis_ecg-proj.fif'
+
+###############################################################################
+# Load the FIF file and display the projections present in the file. Here the
+# projections are added to the file during the acquisition and are obtained
+# from empty room recordings.
+raw = read_raw_fif(fname)
+empty_room_proj = raw.info['projs']
+
+# Display the projections stored in `info['projs']` from the raw object
+raw.plot_projs_topomap()
+
+###############################################################################
+# Display the projections one by one
+fig, axes = plt.subplots(1, len(empty_room_proj))
+for proj, ax in zip(empty_room_proj, axes):
+    proj.plot_topomap(axes=ax)
+
+###############################################################################
+# 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)
+
+###############################################################################
+# As shown in the tutorial on how to
+# :ref:`sphx_glr_auto_tutorials_plot_visualize_raw.py`
+# the ECG projections can be loaded from a file and added to the raw object
+
+# read the projections
+ecg_projs = read_proj(ecg_fname)
+
+# 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)

mne/viz/tests/test_utils.py

@@ -9,7 +9,8 @@
 import pytest
 
 from mne.viz.utils import (compare_fiff, _fake_click, _compute_scalings,
-                           _validate_if_list_of_axes, _get_color_list)
+                           _validate_if_list_of_axes, _get_color_list,
+                           _setup_vmin_vmax)
 from mne.viz import ClickableImage, add_background_image, mne_analyze_colormap
 from mne.utils import run_tests_if_main
 from mne.io import read_raw_fif
@@ -26,6 +27,13 @@
 ev_fname = op.join(base_dir, 'test_raw-eve.fif')
 
 
+def test_setup_vmin_vmax_warns():
+    """Test that _setup_vmin_vmax warns properly."""
+    expected_msg = r'\(min=0.0, max=1\) range.*minimum of data is -1'
+    with pytest.warns(UserWarning, match=expected_msg):
+        _setup_vmin_vmax(data=[-1, 0], vmin=None, vmax=None, norm=True)
+
+
 def test_get_color_list():
     """Test getting a colormap from rcParams."""
     colors = _get_color_list()

mne/viz/topomap.py

@@ -145,13 +145,36 @@ def _add_colorbar(ax, im, cmap, side="right", pad=.05, title=None,
     divider = make_axes_locatable(ax)
     cax = divider.append_axes(side, size=size, pad=pad)
     cbar = plt.colorbar(im, cax=cax, cmap=cmap, format=format)
-    if cmap[1]:
+    if cmap is not None and cmap[1]:
         ax.CB = DraggableColorbar(cbar, im)
     if title is not None:
         cax.set_title(title, y=1.05, fontsize=10)
     return cbar, cax
 
 
+def _eliminate_zeros(proj):
+    """Remove grad or mag data if only contains 0s (gh 5641)."""
+    GRAD_ENDING = ('2', '3')
+    MAG_ENDING = '1'
+
+    proj = copy.deepcopy(proj)
+    proj['data']['data'] = np.atleast_2d(proj['data']['data'])
+
+    for ending in (GRAD_ENDING, MAG_ENDING):
+        names = proj['data']['col_names']
+        idx = [i for i, name in enumerate(names) if name.endswith(ending)]
+
+        # if all 0, remove the 0s an their labels
+        if not proj['data']['data'][0][idx].any():
+            new_col_names = np.delete(np.array(names), idx).tolist()
+            new_data = np.delete(np.array(proj['data']['data'][0]), idx)
+            proj['data']['col_names'] = new_col_names
+            proj['data']['data'] = np.array([new_data])
+
+    proj['data']['ncol'] = len(proj['data']['col_names'])
+    return proj
+
+
 def plot_projs_topomap(projs, layout=None, cmap=None, sensors=True,
                        colorbar=False, res=64, size=1, show=True,
                        outlines='head', contours=6, image_interp='bilinear',
@@ -231,6 +254,10 @@ def plot_projs_topomap(projs, layout=None, cmap=None, sensors=True,
                                    _pair_grad_sensors_ch_names_neuromag122,
                                    Layout, _merge_grad_data)
     from ..channels import _get_ch_type
+
+    is_layout_parameter_none = layout is None
+    is_info_parameter_none = info is None
+
     if info is not None:
         if not isinstance(info, Info):
             raise TypeError('info must be an instance of Info, got %s'
@@ -254,7 +281,6 @@ def plot_projs_topomap(projs, layout=None, cmap=None, sensors=True,
     nrows = math.floor(math.sqrt(n_projs))
     ncols = math.ceil(n_projs / nrows)
 
-    cmap = _setup_cmap(cmap)
     if axes is None:
         plt.figure()
         axes = list()
@@ -269,6 +295,7 @@ def plot_projs_topomap(projs, layout=None, cmap=None, sensors=True,
         title = proj['desc']
         title = '\n'.join(title[ii:ii + 22] for ii in range(0, len(title), 22))
         axes[proj_idx].set_title(title, fontsize=10)
+        proj = _eliminate_zeros(proj)  # gh 5641
         ch_names = _clean_names(proj['data']['col_names'],
                                 remove_whitespace=True)
         data = proj['data']['data'].ravel()
@@ -310,12 +337,30 @@ def plot_projs_topomap(projs, layout=None, cmap=None, sensors=True,
                     data = _merge_grad_data(data[grad_pairs]).ravel()
                 break
             if len(idx) == 0:
-                raise RuntimeError('Cannot find a proper layout for '
-                                   'projection %s, consider explicitly '
-                                   'passing a Layout or Info as the layout '
-                                   'parameter.' % proj['desc'])
+                if ch_names[0].startswith('EEG'):
+                    msg = ('Cannot find a proper layout for projection {0}.'
+                           ' The proper layout of an EEG topomap cannot be'
+                           ' inferred from the data. '.format(proj['desc']))
+                    if is_layout_parameter_none and is_info_parameter_none:
+                        msg += (' For EEG data, valid `layout` or `info` is'
+                                ' required. None was provided, please consider'
+                                ' passing one of them.')
+                    elif not is_layout_parameter_none:
+                        msg += (' A `layout` was provided but could not be'
+                                ' used for display. Please review the `layout`'
+                                ' parameter.')
+                    else:  # layout is none, but we have info
+                        msg += (' The `info` parameter was provided but could'
+                                ' not be for display. Please review the `info`'
+                                ' parameter.')
+                    raise RuntimeError(msg)
+                else:
+                    raise RuntimeError('Cannot find a proper layout for '
+                                       'projection %s, consider explicitly '
+                                       'passing a Layout or Info as the layout'
+                                       ' parameter.' % proj['desc'])
 
-        im = plot_topomap(data, pos[:, :2], vmax=None, cmap=cmap[0],
+        im = plot_topomap(data, pos[:, :2], vmax=None, cmap=cmap,
                           sensors=sensors, res=res, axes=axes[proj_idx],
                           outlines=outlines, contours=contours,
                           image_interp=image_interp, show=False)[0]

mne/viz/utils.py

@@ -46,25 +46,42 @@
 
 
 def _setup_vmin_vmax(data, vmin, vmax, norm=False):
-    """Handle vmin and vmax parameters."""
+    """Handle vmin and vmax parameters for visualizing topomaps.
+
+    For the normal use-case (when `vmin` and `vmax` are None), the parameter
+    `norm` drives the computation. When norm=False, data is supposed to come
+    from a mag and the output tuple (vmin, vmax) is symmetric range
+    (-x, x) where x is the max(abs(data)). When norm=False (aka data is the L2
+    norm of a gradiometer pair) the output tuple corresponds to (0, x).
+
+    Otherwise, vmin and vmax are callables that drive the operation.
+    """
+    should_warn = False
     if vmax is None and vmin is None:
         vmax = np.abs(data).max()
-        if norm:
-            vmin = 0.
-        else:
-            vmin = -vmax
+        vmin = 0. if norm else -vmax
+        if vmin == 0 and np.min(data) < 0:
+            should_warn = True
+
     else:
         if callable(vmin):
             vmin = vmin(data)
         elif vmin is None:
-            if norm:
-                vmin = 0.
-            else:
-                vmin = np.min(data)
+            vmin = 0. if norm else np.min(data)
+            if vmin == 0 and np.min(data) < 0:
+                should_warn = True
+
         if callable(vmax):
             vmax = vmax(data)
         elif vmax is None:
             vmax = np.max(data)
+
+    if should_warn:
+        warn_msg = ("_setup_vmin_vmax output a (min={vmin}, max={vmax})"
+                    " range whereas the minimum of data is {data_min}")
+        warn_val = {'vmin': vmin, 'vmax': vmax, 'data_min': np.min(data)}
+        warn(warn_msg.format(**warn_val), UserWarning)
+
     return vmin, vmax
 
 

tutorials/plot_artifacts_correction_ssp.py

@@ -2,6 +2,12 @@
 Artifact Correction with SSP
 ============================
 
+This tutorial explains how to estimate Signal Subspace Projectors (SSP)
+for correction of ECG and EOG artifacts.
+
+See :ref:`sphx_glr_auto_examples_io_plot_read_proj.py` for how to read
+and visualize already present SSP projection vectors.
+
 """
 import numpy as np
 

tutorials/plot_visualize_raw.py

@@ -105,6 +105,15 @@
 raw.add_proj(projs)
 raw.plot_projs_topomap()
 
+###############################################################################
+# Note that the projections in ``raw.info['projs']`` can be visualized using
+# :meth:`raw.plot_projs_topomap <mne.io.Raw.plot_projs_topomap>` or calling
+# :func:`proj.plot_topomap <mne.Projection.plot_topomap>`
+#
+# more examples can be found in
+# :ref:`sphx_glr_auto_examples_io_plot_read_proj.py`
+projs[0].plot_topomap()
+
 ###############################################################################
 # The first three projectors that we see are the SSP vectors from empty room
 # measurements to compensate for the noise. The fourth one is the average EEG