MRG+1, ENH: add border argument for topomap extrapolation

[mmagnuski]

Reference issue

Continuation of topomap extrapolation pull request.

What does this implement/fix?

This adds a border argument to mne.viz.plot_topomap that allows to change the values set to extrapolation points. The default is zero, just as before, but 'mean' allows to extrapolate to mean of channels that are neighbours to given extrapolation point. See the example below:

(I know, I will have to fix the channel positions after the recent changes to plot_topomap channel position handling)

border='mean' is especially useful when plotting log-power (and setting vmin, vmax):

This is a draft at this point, but feel free to comment.

TODOs:

• add tests
• %()s docs
• expose border in other topomap functions?

[agramfort]

this seems useful. Is this option offered by any other software?

[mmagnuski]

@agramfort Unfortunatelly I don't know. Probably some packages do something like border='mean' here by default, but I've never checked.

[mmagnuski]

It just became an issue for me lately when plotting log-transformed power - all the topo values were negative (in the range of -12 or -11) so the default extrapolation border of 0 looked very weird (the topo values could not be discerned).

[agramfort]

all other packages use zero outside interpolation? I would like to make sure MNE users
get plots that they would get with historical packages. Here for mean option use case seems legit

[larsoner]

border='mean' seems like it might make a sensible new default. Our topomaps will already be changing a bit in this cycle because of #7066 so it's not a bad time to do it. Arguably it could be a bug fix...

You could try setting border='auto' in these PRs and doing a commit with [circle full] in the message so we could see how all topomap examples look with the new default. If they all look at least as good or better, then it's a good reason to switch / consider it a bugfix.

[mmagnuski]

@larsoner BTW, I am writing a test for this PR and it seems I do not understand well the way channel positions are handled in plot_topomap currently. I create 5 fake channels and try to mach the sphere radius to their positions:

import numpy as np
import mne

ch_names = list('abcde')
ch_pos = np.array([[0, 0, 1], [1, 0, 0], [-1, 0, 0],
                   [0, -1, 0], [0, 1, 0]])
ch_pos_dict = {name: pos for name, pos in zip(ch_names, ch_pos)}
dig = mne.channels.make_dig_montage(ch_pos_dict, coord_frame='head')

data = np.ones(5) * 5
info = mne.create_info(ch_names, 250, ['eeg'] * 5)
info.set_montage(dig)

ax, _ = mne.viz.plot_topomap(data, info, extrapolate='head',
                             border=0, sphere=1.)

my expectation was that the radius of 1. in sphere and the channel positions extending to 1 on x and y axis would give me 4 channels on the head outline, however I get this:

are channel positions and sphere in different units by default? It doesn't look like cm vs m, and when I do dig.plot() I get a plot with meters on axis labels (for example x (m)).

However, if change the radius to be 0.1:

ax, _ = mne.viz.plot_topomap(data, info, extrapolate='head',
                             border=0, sphere=0.1)

the channel position is better but the exrapolation is weird:

I expected that extrapolate='head' would give me extrapolation to the head contour.

[mmagnuski]

If that is unexpected I can try to fix that in a subsequent PR.

[mmagnuski]

@agramfort - I can check eeglab and fieldtrip later this week. I think some of the figures in eeglab are in dB by default so negative values should be expected - however, I never got a "cosmic frog" topomap image as in my example with log-power and border=0 in mne. This would mean that eeglab does something like border='mean' but I would have to check to be sure.

[larsoner]

my expectation was that the radius of 1. in sphere and the channel positions extending to 1 on x and y axis would give me 4 channels on the head outline

IIRC the points with z=0 should be on the head outline, points with z>0 should be inside, z<0 should be outside (i.e., radius only depends on z) and the x/y should determine the angle. I think this was the intent. Though there might also be some interaction with the head size / radius in there, too -- for example you might need both that z=0 and np.linalg.norm([x, y, z]) == head_radius for the point to be on the head border.

So you might have indeed found a bug. Any unit tests you can write to check these cases would be welcome. But first we should decide on the behavior. To me it makes sense to have the radius be given by some multiplicative factor dependent on z: z=0 gives a factor of 1, z > 0 gives factors < 1, z < 0 gives factors > 1. And then there is also a multiplicative factor like norm(pos) / head_radius. Something like this...

But IIRC when you go to implement, I think all the computations no longer use normalized units, so the head sphere will extend to a radius of, e.g., 0.1 instead of 1.

If this does not make sense or you want to hash it out interactively, I'm on gitter

[mmagnuski]

To me it makes sense to have the radius be given by some multiplicative factor dependent on z: z=0 gives a factor of 1, z > 0 gives factors < 1, z < 0 gives factors > 1. And then there is also a multiplicative factor like norm(pos) / head_radius. Something like this...

Yes, that was my expectation at least, so I think we are on the same page.

But IIRC when you go to implement, I think all the computations no longer use normalized units, so the head sphere will extend to a radius of, e.g., 0.1 instead of 1.

I'm not sure I understand it here - I had my channel positions be on z=0 and x=0, y=1 for example and asked for sphere=1. but the channel was not on the head outline, but deep inside the head (first example plot).

[larsoner]

I'm not sure I understand it here - I had my channel positions be on z=0 and x=0, y=1 for example and asked for sphere=1. but the channel was not on the head outline, but deep inside the head (first example plot).

Yes that sounds like errant behavior. I'm just saying that when you eventually go to write unit tests for this sort of thing, if you use a head radius of something other than 1. you shouldn't expect z=0 points with norm([x, y, z]) == head_radius to be at locations with norm 1, but rather of norm head_radius.

[mmagnuski]

Yes that sounds like errant behavior. I'm just saying that when you eventually go to write unit tests for this sort of thing, if you use a head radius of something other than 1. you shouldn't expect z=0 points with norm([x, y, z]) == head_radius to be at locations with norm 1, but rather of norm head_radius.

Ah, sure. :)

[codecov]

Codecov Report

Merging #7193 into master will increase coverage by <.01%.
The diff coverage is 100%.

@@            Coverage Diff             @@
##           master    #7193      +/-   ##
==========================================
+ Coverage   89.76%   89.77%   +<.01%     
==========================================
  Files         445      445              
  Lines       79874    80013     +139     
  Branches    12773    12792      +19     
==========================================
+ Hits        71702    71832     +130     
- Misses       5374     5375       +1     
- Partials     2798     2806       +8

[mmagnuski]

Hm, I have a strange Travis error in the topomap test. I can't reproduce locally. It seems no contour levels are created, but matplotlib does: self.vmin = np.amin(self.levels). I'm not sure whether it's a problem on our side or matplotlib's.

――――――――――――――――――――――――――― test_plot_topomap_basic ――――――――――――――――――――――――――――
mne/viz/tests/test_topomap.py:214: in test_plot_topomap_basic
    ax, _ = plot_topomap(data, info, extrapolate='head', border='mean', sphere=1)
mne/viz/topomap.py:686: in plot_topomap
    sphere=sphere, border=border)[:2]
mne/viz/topomap.py:837: in _plot_topomap
    linewidths=linewidth / 2.)
../../../miniconda/envs/test/lib/python3.5/site-packages/matplotlib/__init__.py:1898: in inner
    return func(ax, *args, **kwargs)
../../../miniconda/envs/test/lib/python3.5/site-packages/matplotlib/axes/_axes.py:5825: in contour
    contours = mcontour.QuadContourSet(self, *args, **kwargs)
../../../miniconda/envs/test/lib/python3.5/site-packages/matplotlib/contour.py:865: in __init__
    self._process_levels()
../../../miniconda/envs/test/lib/python3.5/site-packages/matplotlib/contour.py:1199: in _process_levels
    self.vmin = np.amin(self.levels)
../../../miniconda/envs/test/lib/python3.5/site-packages/numpy/core/fromnumeric.py:2352: in amin
    out=out, **kwargs)
../../../miniconda/envs/test/lib/python3.5/site-packages/numpy/core/_methods.py:29: in _amin
    return umr_minimum(a, axis, None, out, keepdims)
E   ValueError: zero-size array to reduction operation minimum which has no identity

[larsoner]

Probably a matplotlib bug but we usually try to protect against this sort of thing by checking to see if we are going to get no levels and if so don't even call contour.

Or if you to punt on this and work around it, you can probably make your data out of np.random.RandomState(0).randn(*shape) instead of np.zeros(shape).

[mmagnuski]

I will soon add border to other topo functions (all that allow using extrapolate argument) and run circle full.

[codecov-io]

Codecov Report

Merging #7193 into master will increase coverage by 0.03%.
The diff coverage is 97.77%.

@@            Coverage Diff            @@
##           master   #7193      +/-   ##
=========================================
+ Coverage   89.76%   89.8%   +0.03%     
=========================================
  Files         445     445              
  Lines       79874   80338     +464     
  Branches    12773   12903     +130     
=========================================
+ Hits        71702   72144     +442     
- Misses       5374    5380       +6     
- Partials     2798    2814      +16

[larsoner]

@mmagnuski did you sort out the scaling issues? Curious to know what the problem was if you did / when you do

[mmagnuski]

@larsoner, no, I didn't get to it yet. I will do that in a separate PR just after this one. The joy of debugging awaits! :)

[mmagnuski]

@larsoner @agramfort

I've added/exposed the border argument where extrapolate was already present/exposed.
Because of this I did not add border to:

• plot_topomap of tfr
• plot_psd_topomap of Epochs
• _plot_corrmap (viz.topomap) and corrmap from preprocessing.ica
• plot_csd (viz.topomap) and CrossSpectralDensity.plot()

I have changed ProjMixin.plot_projs_topomap a little more:

• added extrapolate and border (because mne.viz.plot_projs_topomap had extrapolate already)
• changed the order of arguemnts passed to mne.viz.plot_projs_topomap (only inside the function, no API change)
• found out that ch_type argument is unused

I would also prefer to test changing the default border to 'mean' (border='auto') in a subsequent PR.
Let me know if this is ok.

[larsoner]

Yes that sounds like a good plan to me

Ready for review?

[mmagnuski]

Yes :)

[larsoner]

LGTM +1 for merge

[agramfort]

@mmagnuski can you jsut update latest.inc to update what's new page? thx

[mmagnuski]

Ah, right, sure!

[mmagnuski]

@agramfort - done.

[mmagnuski]

thanks @larsoner and @agramfort ! 🚀