1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
|
# -*- coding: utf-8 -*-
"""
.. _ex-mixed-norm-inverse:
================================================================
Compute sparse inverse solution with mixed norm: MxNE and irMxNE
================================================================
Runs an (ir)MxNE (L1/L2 :footcite:`GramfortEtAl2012` or L0.5/L2
:footcite:`StrohmeierEtAl2014` mixed norm) inverse solver.
L0.5/L2 is done with irMxNE which allows for sparser source estimates with less
amplitude bias due to the non-convexity of the L0.5/L2 mixed norm penalty.
"""
# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
# Daniel Strohmeier <daniel.strohmeier@tu-ilmenau.de>
#
# License: BSD-3-Clause
# %%
# sphinx_gallery_thumbnail_number = 2
import numpy as np
import mne
from mne.datasets import sample
from mne.inverse_sparse import mixed_norm, make_stc_from_dipoles
from mne.minimum_norm import make_inverse_operator, apply_inverse
from mne.viz import (plot_sparse_source_estimates,
plot_dipole_locations, plot_dipole_amplitudes)
print(__doc__)
data_path = sample.data_path()
meg_path = data_path / 'MEG' / 'sample'
fwd_fname = meg_path / 'sample_audvis-meg-eeg-oct-6-fwd.fif'
ave_fname = meg_path / 'sample_audvis-ave.fif'
cov_fname = meg_path / 'sample_audvis-shrunk-cov.fif'
subjects_dir = data_path / 'subjects'
# Read noise covariance matrix
cov = mne.read_cov(cov_fname)
# Handling average file
condition = 'Left Auditory'
evoked = mne.read_evokeds(ave_fname, condition=condition, baseline=(None, 0))
evoked.crop(tmin=0, tmax=0.3)
# Handling forward solution
forward = mne.read_forward_solution(fwd_fname)
# %%
# Run solver with SURE criterion :footcite:`DeledalleEtAl2014`
alpha = "sure" # regularization parameter between 0 and 100 or SURE criterion
loose, depth = 0.9, 0.9 # loose orientation & depth weighting
n_mxne_iter = 10 # if > 1 use L0.5/L2 reweighted mixed norm solver
# if n_mxne_iter > 1 dSPM weighting can be avoided.
# Compute dSPM solution to be used as weights in MxNE
inverse_operator = make_inverse_operator(evoked.info, forward, cov,
depth=depth, fixed=True,
use_cps=True)
stc_dspm = apply_inverse(evoked, inverse_operator, lambda2=1. / 9.,
method='dSPM')
# Compute (ir)MxNE inverse solution with dipole output
dipoles, residual = mixed_norm(
evoked, forward, cov, alpha, loose=loose, depth=depth, maxit=3000,
tol=1e-4, active_set_size=10, debias=False, weights=stc_dspm,
weights_min=8., n_mxne_iter=n_mxne_iter, return_residual=True,
return_as_dipoles=True, verbose=True, random_state=0,
# for this dataset we know we should use a high alpha, so avoid some
# of the slower (lower) alpha values
sure_alpha_grid=np.linspace(100, 40, 10),
)
t = 0.083
tidx = evoked.time_as_index(t)
for di, dip in enumerate(dipoles, 1):
print(f'Dipole #{di} GOF at {1000 * t:0.1f} ms: '
f'{float(dip.gof[tidx]):0.1f}%')
# %%
# Plot dipole activations
plot_dipole_amplitudes(dipoles)
# Plot dipole location of the strongest dipole with MRI slices
idx = np.argmax([np.max(np.abs(dip.amplitude)) for dip in dipoles])
plot_dipole_locations(dipoles[idx], forward['mri_head_t'], 'sample',
subjects_dir=subjects_dir, mode='orthoview',
idx='amplitude')
# Plot dipole locations of all dipoles with MRI slices
for dip in dipoles:
plot_dipole_locations(dip, forward['mri_head_t'], 'sample',
subjects_dir=subjects_dir, mode='orthoview',
idx='amplitude')
# %%
# Plot residual
ylim = dict(eeg=[-10, 10], grad=[-400, 400], mag=[-600, 600])
evoked.pick_types(meg=True, eeg=True, exclude='bads')
evoked.plot(ylim=ylim, proj=True, time_unit='s')
residual.pick_types(meg=True, eeg=True, exclude='bads')
residual.plot(ylim=ylim, proj=True, time_unit='s')
# %%
# Generate stc from dipoles
stc = make_stc_from_dipoles(dipoles, forward['src'])
# %%
# View in 2D and 3D ("glass" brain like 3D plot)
solver = "MxNE" if n_mxne_iter == 1 else "irMxNE"
plot_sparse_source_estimates(forward['src'], stc, bgcolor=(1, 1, 1),
fig_name="%s (cond %s)" % (solver, condition),
opacity=0.1)
# %%
# Morph onto fsaverage brain and view
morph = mne.compute_source_morph(stc, subject_from='sample',
subject_to='fsaverage', spacing=None,
sparse=True, subjects_dir=subjects_dir)
stc_fsaverage = morph.apply(stc)
src_fsaverage_fname = (
subjects_dir / 'fsaverage' / 'bem' / 'fsaverage-ico-5-src.fif')
src_fsaverage = mne.read_source_spaces(src_fsaverage_fname)
plot_sparse_source_estimates(src_fsaverage, stc_fsaverage, bgcolor=(1, 1, 1),
fig_name="Morphed %s (cond %s)" % (solver,
condition), opacity=0.1)
# %%
# References
# ----------
# .. footbibliography::
|
