MRG, ENH: Fix plot_vector_source_estimates

[GuillaumeFavelier]

This PR closes #7020:

• Scale automatically the glyphs when scale_factor=None
• Fix the bug where the glyphs disappears when brain_alpha is in [0.8, 1]
• Experiment with the type of the glyphs (2D or 3D)

[larsoner]

Can you touch the examples that use vector source estimate plotting so we can see if the auto mode is working for them?

[codecov]

Codecov Report

Merging #7030 into master will increase coverage by <.01%.
The diff coverage is 61.53%.

@@            Coverage Diff             @@
##           master    #7030      +/-   ##
==========================================
+ Coverage   89.71%   89.71%   +<.01%     
==========================================
  Files         434      438       +4     
  Lines       77293    77480     +187     
  Branches    12559    12587      +28     
==========================================
+ Hits        69344    69513     +169     
- Misses       5148     5156       +8     
- Partials     2801     2811      +10

[GuillaumeFavelier]

From what I tried locally, the difference in Y axis is too big:

y = np.array(brain.geo[hemi].coords[:, 1])
y_diff = (np.max(y) - np.min(y)) / 5.  # presumably get fifth of the width 

There is probably something I missed.

[GuillaumeFavelier]

This is what I mean by too big:

Taken from: https://16646-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/source-modeling/plot_mne_dspm_source_localization.html#sphx-glr-auto-tutorials-source-modeling-plot-mne-dspm-source-localization-py

[drammock]

Sometimes that's what my brain feels like though. Are we sure it's not correct? 😆

[GuillaumeFavelier]

By discussing offline with @larsoner, I should check if the units get converted and use the convenient numpy.ptp() function to get the range. I expect the result to be different 😅

[GuillaumeFavelier]

I checked the units and everything is fine. Plus I noticed that the scaling of the glyphs is modified twice since there is time information so I decided to use a simpler approach. This one sets the scale_factor of the glyphs as a post-process. Hopefully, it's better.

[GuillaumeFavelier]

This version looks more reasonable to me:

Taken from: https://16651-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/source-modeling/plot_mne_dspm_source_localization.html#sphx-glr-auto-tutorials-source-modeling-plot-mne-dspm-source-localization-py

[larsoner]

Maybe a bit too big? I think it might look better if you shrink to 0.75 or 0.5 the size -- do either of those look better (or at least as good) to you?

[GuillaumeFavelier]

Locally, by varying the factor on the brain width:

20%	15%	10%

In my opinion, 10% is the best among those and decreases the clutter.

[larsoner]

Agreed 10% is clearest

[drammock]

I also like 10% best

[GuillaumeFavelier]

Don't hesitate to tell me if you want me to modify anything else or test another parameters. Otherwise, this PR is ready for reviews.

[larsoner]

@GuillaumeFavelier don't forget to push the 10% change

[larsoner]

... and also update the BUG section of changes/latest.inc to say that we improved the logic for how scale is automatically computed when plotting vector source estimates

[GuillaumeFavelier]

Any idea how I can fix this?

changes/latest.inc:98: WARNING: py:func reference target not found: plot_vector_source_estimates

[larsoner]

mne.viz.plot_vector_source_estimates

It must match the namespace in python_reference.rst

[GuillaumeFavelier]

Oh I see!

[jasmainak]

change brain_alpha to something like 0.8 or 1.0 and you see that this does not completely fix the issue.

[GuillaumeFavelier]

I didn't notice that this PR was supposed to close #7020, I'll update right now.

[larsoner]

I put that in when I thought it was sufficient. Feel free to remove the label if you want to proceed with multiple PRs, or keep at it in this one. Either way is fine by me

[GuillaumeFavelier]

Thanks @larsoner I decided to fix everything in this PR.

[GuillaumeFavelier]

It might require more reviews as it moves forward though.

[GuillaumeFavelier]

Now I suggest a patch enabling depth peeling when brain_alpha is less than 1.0. This modification did the trick for me. Is it working for you too @jasmainak ?

[larsoner]

That's fine, I'm sure @jasmainak is up for it :)

[GuillaumeFavelier]

For reference on depth peeling: enthought/mayavi#491 (comment)

[GuillaumeFavelier]

Why this is not enabled by default actually? Any ideas? 🤔

[larsoner]

It's slower (requires multiple draw passes) and is not supported on all hardware. So better not to use it unless you need to. For example if we can fix it by setting some opacity, translucency, or depth testing params it would be faster and work on more graphics cards.

[GuillaumeFavelier]

I understand. But the number of parameters to modify to obtain a similar result plus the maintenance of those which is not trivial in my opinion has a certain cost.

And I wouldn't be worried about performance for such a lightweight 3D scene especially considering how modern hardware gets more powerful each year. I did not have to work on a very demanding example so far on MNE but maybe there is.

[GuillaumeFavelier]

So far, this is the result on my configuration by experimenting with 3D arrows:

[larsoner]

Looks (and sounds) reasonable to me, @jasmainak can you try it?

[drammock]

just wondering if it's the right choice to scale both length and thickness, or whether keeping thickness constant and scaling only length would be better. Personally I'm not sure.

[larsoner]

It might not be so easy to have Mayavi (or PyVista) just scale one of the two params. @GuillaumeFavelier do you know?

I'm okay with it either way, agreed we should at least think about it.

[jasmainak]

okay okay, certified tester @jasmainak is coming. Give me a moment to try :)

[jasmainak]

I like the 3D glyphs far far better than the 2D ones and I don't see any noticeable lag on my mac.

However, I think some tiny tweaks would certainly help

• make the default transparency less: something like brain_alpha=0.7
• make the arrows smaller: it's too big right now to make any sense from them. On the somato data, scale_factor=10 looked much better than the default.

[agramfort]

https://16699-1301584-gh.circle-artifacts.com/0/dev/auto_tutorials/source-modeling/plot_mne_dspm_source_localization.html#sphx-glr-auto-tutorials-source-modeling-plot-mne-dspm-source-localization-py LGTM merging. thx @GuillaumeFavelier

[rythorpe]

I found a small bug in this commit when trying to plot both hemispheres of the MNE source estimate in this script. Shown below is my traceback.

Also, as a user I have to say that implementing the 3D glyphs would be much easier to view even though they take up slightly more space. Scaling by both length and thickness makes it immediately clear which vertices contribute the most.

(mne_env) [rthorpe@gpu701 rthorpe]$ python somato_mne_vector_plot.py 
/gpfs_home/rthorpe/envs/mne_env/lib/python3.7/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.
  warnings.warn(msg, category=DeprecationWarning)
Opening raw data file /users/rthorpe/mne_data/MNE-somato-data/sub-01/meg/sub-01_task-somato_meg.fif...
    Range : 237600 ... 506999 =    791.189 ...  1688.266 secs
Ready.
Current compensation grade : 0
Reading 0 ... 269399  =      0.000 ...   897.077 secs...
Filtering raw data in 1 contiguous segment
Setting up band-pass filter from 1 - 40 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 1.00
- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)
- Upper passband edge: 40.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 45.00 Hz)
- Filter length: 991 samples (3.300 sec)

111 events found
Event IDs: [1]
111 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
0 projection items activated
Loading data for 111 events and 121 original time points ...
0 bad epochs dropped
Reading forward solution from /users/rthorpe/mne_data/MNE-somato-data/derivatives/sub-01/sub-01_task-somato-fwd.fif...
    Reading a source space...
    [done]
    Reading a source space...
    [done]
    2 source spaces read
    Desired named matrix (kind = 3523) not available
    Read MEG forward solution (8155 sources, 306 channels, free orientations)
    Source spaces transformed to the forward solution coordinate frame
Computing data rank from raw with rank=None
    Using tolerance 1.6e-08 (2.2e-16 eps * 306 dim * 2.3e+05  max singular value)
    Estimated rank (mag + grad): 306
    MEG: rank 306 computed from 306 data channels with 0 projectors
somato_mne_vector_plot.py:27: RuntimeWarning: Something went wrong in the data-driven estimation of the data rank as it exceeds the theoretical rank from the info (306 > 64). Consider setting rank to "auto" or setting it explicitly as an integer.
  cov = mne.compute_covariance(epochs)
Reducing data rank from 306 -> 306
Estimating covariance using EMPIRICAL
Done.
Number of samples used : 13431
[done]
Converting forward solution to surface orientation
    No patch info available. The standard source space normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Computing inverse operator with 306 channels.
    306 out of 306 channels remain after picking
Selected 306 channels
Creating the depth weighting matrix...
    204 planar channels
    limit = 7615/8155 = 10.004172
    scale = 5.17919e-08 exp = 0.8
Applying loose dipole orientations. Loose value of 0.2.
Whitening the forward solution.
Computing data rank from covariance with rank=None
    Using tolerance 2.1e-12 (2.2e-16 eps * 306 dim * 31  max singular value)
    Estimated rank (mag + grad): 64
    MEG: rank 64 computed from 306 data channels with 0 projectors
    Setting small MEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing SVD of whitened and weighted lead field matrix.
    largest singular value = 2.42763
    scaling factor to adjust the trace = 3.87378e+18
Preparing the inverse operator for use...
    Scaled noise and source covariance from nave = 1 to nave = 111
    Created the regularized inverter
    The projection vectors do not apply to these channels.
    Created the whitener using a noise covariance matrix with rank 64 (242 small eigenvalues omitted)
Applying inverse operator to "1"...
    Picked 306 channels from the data
    Computing inverse...
    Eigenleads need to be weighted ...
    Computing residual...
    Explained  85.2% variance
[done]
QStandardPaths: XDG_RUNTIME_DIR not set, defaulting to '/tmp/runtime-rthorpe'
Using control points [2.84792465e-11 3.74674855e-11 1.22291182e-10]
Traceback (most recent call last):
  File "somato_mne_vector_plot.py", line 41, in <module>
    brain = stc.plot(subjects_dir=subjects_dir,subject='01',hemi='both',views='par',initial_time=0.04)
  File "/gpfs_home/rthorpe/envs/mne_env/lib/python3.7/site-packages/mne/source_estimate.py", line 2117, in plot
    initial_time=initial_time, time_unit=time_unit
  File "/gpfs_home/rthorpe/envs/mne_env/lib/python3.7/site-packages/mne/viz/_3d.py", line 2318, in plot_vector_source_estimates
    layer_id = brain.data['layer_id']
  File "/gpfs_home/rthorpe/envs/mne_env/lib/python3.7/site-packages/surfer/viz.py", line 767, in data
    self._get_one_brain([[]], 'data')
  File "/gpfs_home/rthorpe/envs/mne_env/lib/python3.7/site-packages/surfer/viz.py", line 728, in _get_one_brain
    'or brain.brains.' % name)
ValueError: Cannot access brain.data when more than one view is plotted. Use brain.brain_matrix or brain.brains.

[jasmainak]

I can confirm that I can reproduce the issue!

[GuillaumeFavelier]

I can reproduce as well, thanks for reporting 👍 I opened an issue about it in #7083

Thanks also for the feedback about the 3D glyphs, I'm actually refactoring this part of the code in #7060, I can continue experimenting there and if it turns out to be good enough, we can integrate it I think.