WIP: Fix topomap plotting for EEG

[larsoner]

Closes #3987
Closes #4880
Closes #5190
Closes #5472
Closes #6304

Todo:

• ensure our montages actually use the same head sizes (done by dig refactoring)
• Make the plot_sensors be in physical coordinates
• make all code assume a 0.095 m head radius
• check that all examples in gist are fixed
• deprecate params: transform in montage, outlines, ...
• get tests to pass again
• check plot_topomap, plot_sensors, montage.plot examples via circle full

cc @jona-sassenhagen @sappelhoff @cbrnr you might be interested in this plan.

After this, due to using physical units for the head, #5471 should be simpler.

[cbrnr]

Sounds like a good plan. What about #5269 - should we compute the locations ourselves so we have exactly one location per label? Otherwise, we should make sure that whenever we talk about a specific location, e.g. C3, its location should be the same in all our montages.

[larsoner]

I think that's an orthogonal issue, let's discuss that over in #5269

[cbrnr]

The question whether to compute the locations or whether to use a "master montage" from which all montages are derived is orthogonal, but I'd still add the point "ensure all our montages actually use the same underlying locations" to your to do list.

[larsoner]

That seems like an issue that can and should be tackled separately, this PR will be annoying enough to review and test as is

[cbrnr]

OK, then let's not forget about this and discuss in #5269.

[larsoner]

Too risky to rush to finish this for 0.17, better to make it the first thing to work on in 0.18 to make sure we have months to catch and fix bugs

[mmagnuski]

I've editted the TODO list to remove the fix that is included in #5754.

[mmagnuski]

make all code assume a 8cm head radius (or whatever we use in Montage)

Regarding this second point : maybe this should be considered during Montage refactor: to have explicit units in the Montage and maybe also head radius? If there are no head digitization points then head could be prepresented in Montage as a sphere with specified radius.

[agramfort]

as a reminder the idea of a refactor is to read from montage file a dig object that has all necessary information (coordinate system, units etc.) so yes it's part of the puzzle

[mmagnuski]

regarding:

ensure our montages actually use the same head sizes

I have checked all the montages shipped with mne, and:

• all HydroCel montages have average channel distance from the center between 8. - 10. so it is safe to assume these files use cm units.
• EGI_256 has radius of 1.00
• all biosemi and easycap montages have radius of 85. so they seem to use mm
• mgh60 and mgh70 and all standard_... have radius of 0.1 so they probably use m

Here are all the results:

Montage: EGI_256
n channels: 256
average distance from (0, 0, 0): 1.00
distance std: 0.00

Montage: GSN-HydroCel-128
n channels: 131
average distance from (0, 0, 0): 8.72
distance std: 0.58

Montage: GSN-HydroCel-129
n channels: 132
average distance from (0, 0, 0): 8.72
distance std: 0.58

Montage: GSN-HydroCel-256
n channels: 259
average distance from (0, 0, 0): 9.75
distance std: 1.01

Montage: GSN-HydroCel-257
n channels: 260
average distance from (0, 0, 0): 9.75
distance std: 1.01

Montage: GSN-HydroCel-32
n channels: 36
average distance from (0, 0, 0): 8.75
distance std: 0.66

Montage: GSN-HydroCel-64_1.0
n channels: 67
average distance from (0, 0, 0): 9.79
distance std: 0.73

Montage: GSN-HydroCel-65_1.0
n channels: 68
average distance from (0, 0, 0): 9.80
distance std: 0.72

Montage: biosemi128
n channels: 131
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: biosemi16
n channels: 19
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: biosemi160
n channels: 163
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: biosemi256
n channels: 259
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: biosemi32
n channels: 35
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: biosemi64
n channels: 67
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: easycap-M1
n channels: 74
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: easycap-M10
n channels: 61
average distance from (0, 0, 0): 85.00
distance std: 0.00

Montage: mgh60
n channels: 63
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: mgh70
n channels: 73
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_1005
n channels: 346
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_1020
n channels: 97
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_alphabetic
n channels: 68
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_postfixed
n channels: 103
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_prefixed
n channels: 77
average distance from (0, 0, 0): 0.10
distance std: 0.01

Montage: standard_primed
n channels: 103
average distance from (0, 0, 0): 0.10
distance std: 0.01

[mmagnuski]

@larsoner

make all code assume a 8cm head radius

I was recently reading through the topomap code and my thoughts went back to this issue.
Currently topomap is plotted in such a way (at least for EEG) that channel positions from cartesian are turned to spherical, and then the spherical radius is not used at all. This means for example that montages that assume non-spherical head (ellipsoidal for example like EGI) could be plotted a bit better (more on this in the next post, later today). But it got me thinking:

• the radius shoud not matter for outlines='head' anyway (because it is not used currently), only the head center. In this scenario it is enough to not scale the channels to the whole head and it should look ok. However it seems that the essence of outlines='head' is to scale channels to the whole head space, so not scaling would make it almost identical to outlines='skirt' (see next point). This feature (stretching/scaling) leads to strange topos both for channels that cover just part of the head and channels that extend below head circumference: EGI for example has some channels on the cheeks and due to this all channels squeezed further into the head circle with outlines='head', which does not look particulary good or realistic.
• however the outlines='skirt' option can be fixed relatively easily. We can just assume, as I think is done in eeglab, the the head outline represnets the circumference of the head. This would mean that the radius of the head outline should be np.pi / 2 and channel positions projected to 2d should not be rescaled (before rescaling their distance from the center of the topomap corresponds to their polar angle in 3d).
• it is hard to say what would be the right fix for outlines='head'. One option is to have the head outline be at least np.pi / 2 but wider if some channels extend beyond that.

All this relates to how channels are projected to 2d for EEG currently. I did not check how or whether it differes for MEG. It also makes the most sense if generic montage is used, not digitized channel positions. In the case of digitized positions the head shape would not be spherical so the current way channels are projected to 2d may look weird. But anyway, for the proposed fix to work for digitization, only the head center would have to be defined.
So, to summarize, I think the real problem is what to do in outlines='head' option. Even if after Montage/Digitization refactor specific head radius is assumed for topomap every time, it would still not be clear how the behavior should change for outlines='head' (at least if we don't want to make it identical to outlines='skirt'). My vote would be to fix the outlines='skirt' option anyway now, not waiting for Montage refactor and we would have at least one easy solution for this issue.
I can paste some code for the proposed fix later today.

[mmagnuski]

(btw - I meant that radius is not used in to_sphere path, I didn't look at the behavior when to_sphere is not True)

[larsoner]

But anyway, for the proposed fix to work for digitization, only the head center would have to be defined.

Yes for plotting perhaps, but we still want a way of controlling it e.g. for source localization / coreg with surrogate MRI (like fsaverage). Standardizing around some standard default (80cm?) and method for changing it should only help these issues and people's use of these objects, even if it ends up being not necessary to fix some of them (e.g., plotting in some modes).

[sappelhoff]

I just got triggered by the 80 :-) and wanted to cross ref to this PR and the comment that @massich made: #6600 (review)

[mmagnuski]

@larsoner
Oh, yes, I meant this for topomap plotting only. I'm not saying Montage refactor is not needed!

[mmagnuski]

So, to make what I described more tangible (and actionable perhaps) here is a sketch of the proposed solution for outlines='skirt'. Notice that the radius could be scaled to whatever value we prefer, but the general idea is to keep correspondence between head outline (graphical representation in the topomap) and head circumferece (approximately at widest z axis level):

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.channels import read_montage
from mne.channels.layout import _find_topomap_coords


# modified outlines construction
# ------------------------------
def vars_for_skirt(info):
    '''Create outlines dictionary for outlines='skirt' mode with head outline
    corresponding to head circumference.'''
    radius = np.pi / 2
    ll = np.linspace(0, 2 * np.pi, 101)
    head_x = np.cos(ll) * radius
    head_y = np.sin(ll) * radius
    nose_x = np.array([0.18, 0, -0.18]) * radius
    nose_y = np.array([radius - .004, radius * 1.15, radius - .004])
    ear_x = np.array([.497, .510, .518, .5299, .5419, .54, .547,
                      .532, .510, .489]) * (radius / 0.5)
    ear_y = np.array([.0555, .0775, .0783, .0746, .0555, -.0055, -.0932,
                      -.1313, -.1384, -.1199]) * (radius / 0.5)
    
    outlines_dict = dict(head=(head_x, head_y), nose=(nose_x, nose_y),
                         ear_left=(ear_x, ear_y), ear_right=(-ear_x, ear_y))
    outlines_dict['autoshrink'] = False
    
    mask_outlines = np.c_[head_x, head_y]
    picks = range(len(info['ch_names']))
    pos = _find_topomap_coords(info, picks=picks)
    head_pos = dict(center=(0., 0.))
    
    mask_scale = max(1.25, np.linalg.norm(pos, axis=1).max() / radius)
    mask_outlines *= mask_scale
    outlines_dict['clip_radius'] = (mask_scale * radius,) * 2    
    outlines_dict['mask_pos'] = (mask_outlines[:, 0], mask_outlines[:, 1])
    return pos, outlines_dict, head_pos


# create data
# -----------
ch_names = [f'E{idx}' for idx in range(1, 65)] + ['Cz']
mntg = read_montage('GSN-HydroCel-65_1.0', ch_names=ch_names)
info = mne.create_info(ch_names, sfreq=250., ch_types='eeg', montage=mntg)

picks = [0, 1, 2, 4, 5, 7, 8, 9, 10, 55, 56, 57, 58, 59, 60, 61]
info_picked = mne.pick_info(info, sel=picks)

activity = [0.26, -0.02, -0.16, 0.2, -0.13, -0.15, 0.1,  0.03, 0.18, 0.12,
            0.38, 0.2, 0.39, 0.33, 0.48, 0.46, 0.65, 0.84, 0.74, 0.33, 0.61,
            0.46, 1.02, 0.02, -0.02, -0.06, -0.49, 0.08, -0.07, -0.31, -0.09,
            0.01, -0.12, 0.13, -0.25, 0.16, 0.01, -0.01, 0.02, 0.04, 0.49,
            -0.29, 0.04, -0.24, -0.13, 0., 0.18, 0.19, 0.28, 0.11, 0.26, 0.28,
            0.05, -0.06, -0.04, -0.11, -0.11, -0.03, -0.24, 0.01, 0.24, 0.43,
            0.44, 0.76, 0.61]
activity = np.array(activity)

# figure
# ------
fig, ax = plt.subplots(figsize=(10, 10), ncols=2, nrows=2)
topo_kwargs = dict(sensors='kx', outlines='skirt', vmin=-1.5, vmax=1.5)

ax[0, 0].set_title('current, all channels', fontsize=14)
mne.viz.plot_topomap(activity, info, axes=ax[0, 0], **topo_kwargs)

ax[0, 1].set_title('current, selected channels', fontsize=14)
mne.viz.plot_topomap(activity[picks], info_picked, axes=ax[0, 1],
                     **topo_kwargs)


# proposed fix: use circumference as head radius when outlines='skirt'
topo_kwargs2 = topo_kwargs.copy()

ax[1, 0].set_title('"fix", all channels', fontsize=14)
pos, outlines, head_pos = vars_for_skirt(info)
topo_kwargs2['outlines'] = outlines
mne.viz.plot_topomap(activity, pos, axes=ax[1, 0], head_pos=head_pos,
                     **topo_kwargs2)

ax[1, 1].set_title('"fix", selected channels', fontsize=14)
pos, outlines, head_pos = vars_for_skirt(info_picked)
topo_kwargs2['outlines'] = outlines
mne.viz.plot_topomap(activity[picks], pos, axes=ax[1, 1],
                     head_pos=head_pos, **topo_kwargs2)

This is the resulting figure:

The improvement for second column is obvoius but the differences in channel placement for the first column are also an improvement. The two most frontal channels (the two closest to the nose) are just above the eyes on real human heads (as the lower row suggests), not above the forhead (as the upper row seems to suggest):

(also the middle of the cap is better aligned with the head and the cheek channels are further outside head outlines, as they are usually pretty low - often below the nose tip if one is not smiling :) )

[mmagnuski]

@larsoner
The same as in the code above can be done in the case of outlines='head' but then the radius would be defined as max(np.linalg.norm(pos, axis=1).max(), np.pi / 2) (maybe plus 5% or some small value) so that all channels would be within radius. WDYT?
Another thing: do you know how often to_sphere=False option was used in mne topomap plotting?

[agramfort]

to move forward here I would suggest the following: - open a PR with an example that shows all the things that are broken or could be a bit improved (plot topos with different options, caps etc.) - change the mne behavior and report the updated example we can look at to see that things are actually better now.

[mmagnuski]

@agramfort Ok, I will!

[larsoner]

I'll close this but others feel free to take over!

[larsoner]

the radius shoud not matter for outlines='head' anyway (because it is not used currently), only the head center. In this scenario it is enough to not scale the channels to the whole head and it should look ok. However it seems that the essence of outlines='head' is to scale channels to the whole head space, so not scaling would make it almost identical to outlines='skirt'

Based on this and how many bugs there are, I think we should get rid of this outlines stuff. It just complicates things without much benefit once we fix the behavior.

So I think we should build off of what you did @mmagnuski, specifically:

1. Make the one parameter head_radius or radius. We make this default to 0.095 as in make_standard_montage, make 'auto' an option. (It should not be the default because the fit will behave badly when very few electrodes are present.)
2. Plot everything in head coords, with the head circle representing the head at z=0 (i.e., representing a circle with radius 0.095 in physical space, where the LPA/Nasion/RPA are on the same plane as the circle)

I have made a gist containing what I think are all of the bad behaviors:

https://gist.github.com/larsoner/f30e6b3b2f4e7c7979abbd367ea7ee54

My plan of action is to take @mmagnuski's code and:

1. Swap it in for our current code
2. Ensure it does points 1 and 2 above
3. Ensure it fixes all of these problems in the gist
4. Fix tests
5. Properly deprecate no-longer-needed arguments
6. Run a circle full build and look at examples

Let me know if anyone sees problems with this plan.

[agramfort]

sounds like a great plan !

…

[mmagnuski]

@larsoner
Hi there, I have been MIA, but I am back. :)
I can work on this starting this weekend if that is ok with all of you.

[mmagnuski]

I am all in favor of throwing out outlines='head' and having outlines='skirt' as the default.
And your plan @larsoner, is great, so I don't want to take all that fun away from you if you already invested some time and motivation in it. :)

[larsoner]

@mmagnuski I think I have it almost working after a lot of refactoring. Will push here shortly then you can try it this weekend

[larsoner]

Okay @mmagnuski that's all I have time for, feel free to try it. If there are obvious bugs feel free to push commits to fix them.

[larsoner]

Argh push force then reopen is problematic, will open a new PR