Note
Click here to download the full example code
Compute LCMV inverse solution in volume source space¶
Compute LCMV beamformer on an auditory evoked dataset in a volume source space,
and show activation on fsaverage
.
# Author: Alexandre Gramfort <alexandre.gramfort@inria.fr>
#
# License: BSD (3-clause)
import mne
from mne.datasets import sample
from mne.beamformer import make_lcmv, apply_lcmv
print(__doc__)
Data preprocessing:
data_path = sample.data_path()
subjects_dir = data_path + '/subjects'
raw_fname = data_path + '/MEG/sample/sample_audvis_filt-0-40_raw.fif'
fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-vol-7-fwd.fif'
# Get epochs
event_id, tmin, tmax = [1, 2], -0.2, 0.5
# Read forward model
forward = mne.read_forward_solution(fname_fwd)
# Setup for reading the raw data
raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw.info['bads'] = ['MEG 2443', 'EEG 053'] # 2 bads channels
events = mne.find_events(raw)
# Pick the channels of interest
raw.pick(['meg', 'eog'])
# Read epochs
proj = False # already applied
epochs = mne.Epochs(raw, events, event_id, tmin, tmax,
baseline=(None, 0), preload=True, proj=proj,
reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))
evoked = epochs.average()
# Visualize sensor space data
evoked.plot_joint()
Out:
Reading forward solution from /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis-meg-vol-7-fwd.fif...
Reading a source space...
[done]
1 source spaces read
Desired named matrix (kind = 3523) not available
Read MEG forward solution (3757 sources, 306 channels, free orientations)
Source spaces transformed to the forward solution coordinate frame
Opening raw data file /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_filt-0-40_raw.fif...
Read a total of 4 projection items:
PCA-v1 (1 x 102) idle
PCA-v2 (1 x 102) idle
PCA-v3 (1 x 102) idle
Average EEG reference (1 x 60) idle
Range : 6450 ... 48149 = 42.956 ... 320.665 secs
Ready.
Current compensation grade : 0
Reading 0 ... 41699 = 0.000 ... 277.709 secs...
319 events found
Event IDs: [ 1 2 3 4 5 32]
145 matching events found
Applying baseline correction (mode: mean)
Not setting metadata
Created an SSP operator (subspace dimension = 3)
Loading data for 145 events and 106 original time points ...
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on MAG : ['MEG 1411']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on MAG : ['MEG 1711']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on MAG : ['MEG 1711']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
Rejecting epoch based on EOG : ['EOG 061']
29 bad epochs dropped
Compute covariance matrices, fit and apply spatial filter.
# Read regularized noise covariance and compute regularized data covariance
noise_cov = mne.compute_covariance(epochs, tmin=tmin, tmax=0, method='shrunk',
rank=None)
data_cov = mne.compute_covariance(epochs, tmin=0.04, tmax=0.15,
method='shrunk', rank=None)
# Compute weights of free orientation (vector) beamformer with weight
# normalization (neural activity index, NAI). Providing a noise covariance
# matrix enables whitening of the data and forward solution. Source orientation
# is optimized by setting pick_ori to 'max-power'.
# weight_norm can also be set to 'unit-noise-gain'. Source orientation can also
# be 'normal' (but only when using a surface-based source space) or None,
# which computes a vector beamfomer. Note, however, that not all combinations
# of orientation selection and weight normalization are implemented yet.
filters = make_lcmv(evoked.info, forward, data_cov, reg=0.05,
noise_cov=noise_cov, pick_ori='max-power',
weight_norm='nai', rank=None)
print(filters)
# You can save these with:
# filters.save('filters-lcmv.h5')
# Apply this spatial filter to the evoked data.
stc = apply_lcmv(evoked, filters, max_ori_out='signed')
Out:
Computing data rank from raw with rank=None
Using tolerance 4.7e-09 (2.2e-16 eps * 305 dim * 6.9e+04 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Created an SSP operator (subspace dimension = 3)
Setting small MEG eigenvalues to zero (without PCA)
Reducing data rank from 305 -> 302
Estimating covariance using SHRUNK
Done.
Number of samples used : 3596
[done]
Computing data rank from raw with rank=None
Using tolerance 4e-09 (2.2e-16 eps * 305 dim * 5.8e+04 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Created an SSP operator (subspace dimension = 3)
Setting small MEG eigenvalues to zero (without PCA)
Reducing data rank from 305 -> 302
Estimating covariance using SHRUNK
Done.
Number of samples used : 2088
[done]
Computing data rank from covariance with rank=None
Using tolerance 1e-12 (2.2e-16 eps * 305 dim * 15 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Computing data rank from covariance with rank=None
Using tolerance 3.8e-13 (2.2e-16 eps * 305 dim * 5.6 max singular value)
Estimated rank (mag + grad): 302
MEG: rank 302 computed from 305 data channels with 3 projectors
Making LCMV beamformer with rank {'meg': 302}
Computing inverse operator with 305 channels.
305 out of 306 channels remain after picking
Selected 305 channels
Whitening the forward solution.
Created an SSP operator (subspace dimension = 3)
Computing data rank from covariance with rank={'meg': 302}
Setting small MEG eigenvalues to zero (without PCA)
Creating the source covariance matrix
Adjusting source covariance matrix.
Computing beamformer filters for 3757 sources
Filter computation complete
<Beamformer | LCMV, subject unknown, 3757 vert, 305 ch, max-power ori, nai norm, rank 302>
Plot source space activity:
# You can save result in stc files with:
# stc.save('lcmv-vol')
lims = [0.3, 0.6, 0.9]
stc.plot(
src=forward['src'], subject='sample', subjects_dir=subjects_dir,
clim=dict(kind='value', pos_lims=lims), mode='stat_map',
initial_time=0.1, verbose=True)

Out:
Fixing initial time: 0.1 sec
Showing: t = 0.100 s, (-61.3, -19.0, -6.3) mm, [2, 9, 10] vox, 5651 vertex
Now let’s plot this on a glass brain, which will automatically transform the data to MNI Talairach space:
stc.plot(
src=forward['src'], subject='sample', subjects_dir=subjects_dir,
mode='glass_brain', clim=dict(kind='value', lims=lims),
initial_time=0.1, verbose=True)

Out:
Transforming subject RAS (non-zero origin) -> MNI Talairach
1.022485 -0.008449 -0.036217 5.60 mm
0.071071 0.914866 0.406098 -19.82 mm
0.008756 -0.433700 1.028119 -1.55 mm
0.000000 0.000000 0.000000 1.00
Fixing initial time: 0.1 sec
Showing: t = 0.100 s, (-56.7, -44.1, -0.3) mm, [2, 9, 10] vox, 5651 vertex
Finally let’s get another view, this time plotting again a 'stat_map'
style but using volumetric morphing to get data to fsaverage space,
which we can get by passing a mne.SourceMorph
as the src
argument to mne.VolSourceEstimate.plot
. To save a bit of speed when
applying the morph, we will crop the STC:
morph = mne.compute_source_morph(
forward['src'], 'sample', 'fsaverage', subjects_dir=subjects_dir,
zooms=7, verbose=True)
stc.copy().crop(0.05, 0.18).plot(
src=morph, subject='fsaverage', subjects_dir=subjects_dir,
mode='stat_map', clim=dict(kind='value', pos_lims=lims),
initial_time=0.1, verbose=True)

Out:
volume source space inferred...
loading /home/circleci/mne_data/MNE-sample-data/subjects/sample/mri/brain.mgz as "from" volume
loading /home/circleci/mne_data/MNE-sample-data/subjects/fsaverage/mri/brain.mgz as "to" volume
Computing nonlinear Symmetric Diffeomorphic Registration...
Optimizing translation:
Optimizing level 2 [max iter: 100]
Optimizing level 1 [max iter: 100]
Optimizing level 0 [max iter: 10]
Optimizing rigid-body:
Optimizing level 2 [max iter: 100]
Optimizing level 1 [max iter: 100]
Optimizing level 0 [max iter: 10]
Optimizing full affine:
Optimizing level 2 [max iter: 100]
Optimizing level 1 [max iter: 100]
Optimizing level 0 [max iter: 10]
Optimizing SDR:
Creating scale space from the moving image. Levels: 3. Sigma factor: 0.200000.
Creating scale space from the static image. Levels: 3. Sigma factor: 0.200000.
Optimizing level 2
Optimizing level 1
Optimizing level 0
[done]
Fixing initial time: 0.1 sec
Showing: t = 0.100 s, (-54.0, -37.0, 2.0) mm, [26, 18, 13] vox, 18489 vertex
Total running time of the script: ( 0 minutes 34.552 seconds)
Estimated memory usage: 721 MB