[med-svn] [python-mne] 210/353: ENH : add implementation of LCMV
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:24:59 UTC 2015
This is an automated email from the git hooks/post-receive script.
yoh pushed a commit to tag 0.4
in repository python-mne.
commit 2923459327ef3dbfda0f17e5f0788c38e3e570b7
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Sat Jun 18 13:27:19 2011 -0400
ENH : add implementation of LCMV
---
doc/source/whats_new.rst | 2 +
examples/inverse/plot_lcmv_beamformer.py | 73 +++++++++++++++++++++
mne/beamformer/__init__.py | 4 ++
mne/beamformer/_lcmv.py | 105 +++++++++++++++++++++++++++++++
mne/beamformer/tests/__init__.py | 0
mne/beamformer/tests/test_lcmv.py | 66 +++++++++++++++++++
mne/cov.py | 44 +++++++++++++
7 files changed, 294 insertions(+)
diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst
index be59acc..55157f6 100644
--- a/doc/source/whats_new.rst
+++ b/doc/source/whats_new.rst
@@ -7,6 +7,8 @@ Current
Changelog
~~~~~~~~~
+ - LCMV Beamformer by `Alex Gramfort`_.
+
- Add support for reading named channel selections by `Martin Luessi`_.
- Add Raw.filter method to more easily band pass data by `Alex Gramfort`_.
diff --git a/examples/inverse/plot_lcmv_beamformer.py b/examples/inverse/plot_lcmv_beamformer.py
new file mode 100644
index 0000000..39af9c6
--- /dev/null
+++ b/examples/inverse/plot_lcmv_beamformer.py
@@ -0,0 +1,73 @@
+"""
+======================================
+Compute LCMV beamformer on evoked data
+======================================
+
+Compute LCMV beamformer solution on evoked dataset
+and stores the solution in stc files for visualisation.
+
+"""
+
+# Author: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+print __doc__
+
+import pylab as pl
+import numpy as np
+
+import mne
+from mne.datasets import sample
+from mne.fiff import Raw, pick_types
+from mne.beamformer import lcmv
+
+data_path = sample.data_path('..')
+raw_fname = data_path + '/MEG/sample/sample_audvis_raw.fif'
+event_fname = data_path + '/MEG/sample/sample_audvis_raw-eve.fif'
+fname_fwd = data_path + '/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif'
+fname_cov = data_path + '/MEG/sample/sample_audvis-cov.fif'
+label_name = 'Aud-lh'
+fname_label = data_path + '/MEG/sample/labels/%s.label' % label_name
+
+###############################################################################
+# Get epochs
+event_id, tmin, tmax = 1, -0.2, 0.5
+
+# Setup for reading the raw data
+raw = Raw(raw_fname)
+raw.info['bads'] = ['MEG 2443', 'EEG 053'] # 2 bads channels
+events = mne.read_events(event_fname)
+
+# Set up pick list: EEG + MEG - bad channels (modify to your needs)
+left_temporal_channels = mne.read_selection('Left-temporal')
+picks = pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,
+ exclude=raw.info['bads'], selection=left_temporal_channels)
+
+# Read epochs
+epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
+ picks=picks, baseline=(None, 0), preload=True,
+ reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))
+evoked = epochs.average()
+
+forward = mne.read_forward_solution(fname_fwd)
+
+noise_cov = mne.read_cov(fname_cov)
+noise_cov = mne.cov.regularize(noise_cov, evoked.info,
+ mag=0.05, grad=0.05, eeg=0.1, proj=True)
+
+data_cov = mne.compute_covariance(epochs, tmin=0.04, tmax=0.15)
+stc = lcmv(evoked, forward, noise_cov, data_cov, reg=0.01)
+
+# Save result in stc files
+stc.save('lcmv')
+
+###############################################################################
+# View activation time-series
+data, times, _ = mne.label_time_courses(fname_label, "lcmv-lh.stc")
+pl.close('all')
+pl.plot(1e3 * times, np.mean(data, axis=0))
+pl.xlabel('time (ms)')
+pl.ylabel('LCMV value')
+pl.title('LCMV in %s' % label_name)
+pl.show()
diff --git a/mne/beamformer/__init__.py b/mne/beamformer/__init__.py
new file mode 100644
index 0000000..7c2df4c
--- /dev/null
+++ b/mne/beamformer/__init__.py
@@ -0,0 +1,4 @@
+"""Artifacts finding/correction related functions
+"""
+
+from ._lcmv import lcmv
diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py
new file mode 100644
index 0000000..40b8afa
--- /dev/null
+++ b/mne/beamformer/_lcmv.py
@@ -0,0 +1,105 @@
+"""Compute Linearly constrained minimum variance (LCMV) beamformer.
+"""
+
+# Authors: Alexandre Gramfort <gramfort at nmr.mgh.harvard.edu>
+#
+# License: BSD (3-clause)
+
+import numpy as np
+from scipy import linalg
+
+from ..fiff.constants import FIFF
+from ..fiff.proj import make_projector
+from ..fiff.pick import pick_types, pick_channels_forward
+from ..minimum_norm.inverse import _make_stc, _get_vertno, combine_xyz
+from ..cov import compute_whitener
+
+
+def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01):
+ """Linearly Constrained Minimum Variance (LCMV) beamformer.
+
+ Compute LCMV on evoked data starting from
+ a forward operator.
+
+ Parameters
+ ----------
+ evoked : Evoked
+ Evoked data to invert
+ forward : dict
+ Forward operator
+
+ Returns
+ -------
+ stc : dict
+ Source time courses
+
+ Notes
+ -----
+ The original reference is:
+ Van Veen et al. Localization of brain electrical activity via linearly
+ constrained minimum variance spatial filtering.
+ Biomedical Engineering (1997) vol. 44 (9) pp. 867--880
+ """
+ is_free_ori = forward['source_ori'] == FIFF.FIFFV_MNE_FREE_ORI
+
+ picks = pick_types(evoked.info, meg=True, eeg=True,
+ exclude=evoked.info['bads'])
+ ch_names = [evoked.ch_names[k] for k in picks]
+
+ forward = pick_channels_forward(forward, include=ch_names)
+
+ M = evoked.data
+ G = forward['sol']['data']
+
+ # Handle SSPs
+ proj, ncomp, _ = make_projector(evoked.info['projs'], evoked.ch_names)
+ M = np.dot(proj, M)
+ G = np.dot(proj, G)
+
+ # Handle whitening + data covariance
+ W, _ = compute_whitener(noise_cov, evoked.info, picks)
+
+ # whiten data and leadfield
+ M = np.dot(W, M)
+ G = np.dot(W, G)
+
+ # Apply SSPs + whitener to data covariance
+ Cm = data_cov['data']
+ Cm = np.dot(proj, np.dot(Cm, proj.T))
+ Cm = np.dot(W, np.dot(Cm, W.T))
+
+ # Cm += reg * np.trace(Cm) / len(Cm) * np.eye(len(Cm))
+ Cm_inv = linalg.pinv(Cm, reg)
+
+ # Compute spatial filters
+ W = np.dot(G.T, Cm_inv)
+ n_orient = 3 if is_free_ori else 1
+ n_sources = G.shape[1] // n_orient
+ for k in range(n_sources):
+ Wk = W[n_orient * k: n_orient * k + n_orient]
+ Gk = G[:, n_orient * k: n_orient * k + n_orient]
+ Ck = np.dot(Wk, Gk)
+ Wk[:] = np.dot(linalg.pinv(Ck, 0.01), Wk)
+
+ # noise normalization
+ noise_norm = np.sum(W ** 2, axis=1)
+ if is_free_ori:
+ noise_norm = np.sum(np.reshape(noise_norm, (-1, 3)), axis=1)
+
+ noise_norm = np.sqrt(noise_norm)
+
+ sol = np.dot(W, M)
+
+ if is_free_ori:
+ print 'combining the current components...',
+ sol = combine_xyz(sol)
+
+ sol /= noise_norm[:, None]
+
+ tstep = 1.0 / evoked.info['sfreq']
+ tmin = float(evoked.first) / evoked.info['sfreq']
+ vertno = _get_vertno(forward['src'])
+ stc = _make_stc(sol, tmin, tstep, vertno)
+ print '[done]'
+
+ return stc
diff --git a/mne/beamformer/tests/__init__.py b/mne/beamformer/tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/mne/beamformer/tests/test_lcmv.py b/mne/beamformer/tests/test_lcmv.py
new file mode 100644
index 0000000..99e03e7
--- /dev/null
+++ b/mne/beamformer/tests/test_lcmv.py
@@ -0,0 +1,66 @@
+import os.path as op
+
+from nose.tools import assert_true
+import numpy as np
+# from numpy.testing import assert_array_almost_equal, assert_equal
+
+import mne
+from mne.datasets import sample
+from mne.beamformer import lcmv
+
+
+examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples')
+data_path = sample.data_path(examples_folder)
+fname_data = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis-ave.fif')
+fname_raw = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis_raw.fif')
+fname_cov = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis-cov.fif')
+fname_fwd = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis-meg-oct-6-fwd.fif')
+fname_event = op.join(data_path, 'MEG', 'sample',
+ 'sample_audvis_raw-eve.fif')
+label = 'Aud-lh'
+fname_label = op.join(data_path, 'MEG', 'sample', 'labels', '%s.label' % label)
+
+label = mne.read_label(fname_label)
+noise_cov = mne.read_cov(fname_cov)
+raw = mne.fiff.Raw(fname_raw)
+forward = mne.read_forward_solution(fname_fwd)
+events = mne.read_events(fname_event)
+
+
+def test_lcmv():
+ """Test LCMV
+ """
+ event_id, tmin, tmax = 1, -0.2, 0.2
+
+ # Setup for reading the raw data
+ raw.info['bads'] = ['MEG 2443', 'EEG 053'] # 2 bads channels
+
+ # Set up pick list: EEG + MEG - bad channels (modify to your needs)
+ left_temporal_channels = mne.read_selection('Left-temporal')
+ picks = mne.fiff.pick_types(raw.info, meg=True, eeg=False, stim=True, eog=True,
+ exclude=raw.info['bads'], selection=left_temporal_channels)
+
+ # Read epochs
+ epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
+ picks=picks, baseline=(None, 0), preload=True,
+ reject=dict(grad=4000e-13, mag=4e-12, eog=150e-6))
+ evoked = epochs.average()
+
+ noise_cov = mne.read_cov(fname_cov)
+ noise_cov = mne.cov.regularize(noise_cov, evoked.info,
+ mag=0.05, grad=0.05, eeg=0.1, proj=True)
+
+ data_cov = mne.compute_covariance(epochs, tmin=0.04, tmax=0.15)
+ stc = lcmv(evoked, forward, noise_cov, data_cov, reg=0.01)
+
+ stc_pow = np.sum(stc.data, axis=1)
+ idx = np.argmax(stc_pow)
+ max_stc = stc.data[idx]
+ tmax = stc.times[np.argmax(max_stc)]
+
+ assert_true(0.09 < tmax < 0.1)
+ assert_true(2. < np.max(max_stc) < 3.)
diff --git a/mne/cov.py b/mne/cov.py
index d740240..c63a831 100644
--- a/mne/cov.py
+++ b/mne/cov.py
@@ -518,3 +518,47 @@ def regularize(cov, info, mag=0.1, grad=0.1, eeg=0.1, exclude=None,
cov['data'] = C
return cov
+
+
+def compute_whitener(noise_cov, info, picks=None):
+ """Compute whitening matrix
+
+ Parameters
+ ----------
+ noise_cov : Covariance
+ The noise covariance
+ info : dict
+ The measurement info
+ picks : array of int | None
+ The channels indices to include. If None the data
+ channels in info, except bad channels, are used.
+
+ Returns
+ -------
+ W : 2d array
+ The whitening matrix
+ ch_names : list
+ The channel names
+ """
+ if picks is None:
+ picks = pick_types(info, meg=True, eeg=True, exclude=info['bads'])
+
+ ch_names = [info['chs'][k]['ch_name'] for k in picks]
+
+ noise_cov = copy.deepcopy(noise_cov)
+ noise_cov = prepare_noise_cov(noise_cov, info, ch_names)
+ n_chan = len(ch_names)
+
+ W = np.zeros((n_chan, n_chan), dtype=np.float)
+ #
+ # Omit the zeroes due to projection
+ #
+ eig = noise_cov['eig']
+ nzero = (eig > 0)
+ W[nzero, nzero] = 1.0 / np.sqrt(eig[nzero])
+ #
+ # Rows of eigvec are the eigenvectors
+ #
+ W = np.dot(W, noise_cov['eigvec'])
+ W = np.dot(noise_cov['eigvec'].T, W)
+ return W, ch_names
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/python-mne.git
More information about the debian-med-commit
mailing list