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 <[email protected]> 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 <[email protected]> +# +# 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 <[email protected]> +# +# 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
