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 db09465a5d83070d4aa1ef8c0dc8788fbeb8f3d9 Author: Martin Luessi <[email protected]> Date: Mon Jul 30 14:16:40 2012 -0400 ENH: LCMV for raw data --- mne/beamformer/__init__.py | 2 +- mne/beamformer/_lcmv.py | 197 +++++++++++++++++++++++++++++--------- mne/beamformer/tests/test_lcmv.py | 52 +++++++++- mne/cov.py | 11 ++- mne/tests/test_cov.py | 9 +- 5 files changed, 218 insertions(+), 53 deletions(-) diff --git a/mne/beamformer/__init__.py b/mne/beamformer/__init__.py index 7c2df4c..c2b177a 100644 --- a/mne/beamformer/__init__.py +++ b/mne/beamformer/__init__.py @@ -1,4 +1,4 @@ """Artifacts finding/correction related functions """ -from ._lcmv import lcmv +from ._lcmv import lcmv, lcmv_raw diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py index 389cf88..6915cf6 100644 --- a/mne/beamformer/_lcmv.py +++ b/mne/beamformer/_lcmv.py @@ -8,71 +8,81 @@ import numpy as np from scipy import linalg +from ..fiff import Evoked, Raw from ..fiff.constants import FIFF from ..fiff.proj import make_projector -from ..fiff.pick import pick_types, pick_channels_forward +from ..fiff.pick import pick_types, pick_channels_forward, pick_channels_cov 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 Linearly Constrained Minimum Variance (LCMV) beamformer - on evoked data. - - NOTE : This implementation is heavilly tested so please - report any issue or suggestions. +def _lcmv_any(evraw, forward, noise_cov, data_cov, reg, label, start, stop, + picks): + """ LCMV beamformer for evoked or raw data """ - Parameters - ---------- - evoked : Evoked - Evoked data to invert - forward : dict - Forward operator - noise_cov : Covariance - The noise covariance - data_cov : Covariance - The data covariance - reg : float - The regularization for the whitened data covariance. - - 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] + if picks is None: + picks = pick_types(evraw.info, meg=True, eeg=True, + exclude=evraw.info['bads']) + + ch_names = [evraw.ch_names[k] for k in picks] + # restrict forward solution to selected channels forward = pick_channels_forward(forward, include=ch_names) - M = evoked.data - G = forward['sol']['data'] + # get gain matrix (forward operator) + if label is not None: + if forward['src'][0]['type'] != 'surf': + return Exception('Labels are only supported with surface ' + 'source spaces') + + vertno_fwd = _get_vertno(forward['src']) + if label['hemi'] == 'lh': + vertno_sel = np.intersect1d(vertno_fwd[0], label['vertices']) + idx_sel = np.searchsorted(vertno_fwd[0], vertno_sel) + vertno = [vertno_sel, np.empty(0, dtype=vertno_sel.dtype)] + elif label['hemi'] == 'rh': + vertno_sel = np.intersect1d(vertno_fwd[1], label['vertices']) + idx_sel = len(vertno_fwd[0]) + np.searchsorted(vertno_fwd[1], + vertno_sel) + vertno = [np.empty(0, dtype=vertno_sel.dtype), vertno_sel] + else: + raise Exception("Unknown hemisphere type") + + if is_free_ori: + idx_sel_free = np.zeros(3 * len(idx_sel), dtype=idx_sel.dtype) + for i in range(3): + idx_sel_free[i::3] = 3 * idx_sel + i + idx_sel = idx_sel_free + + G = forward['sol']['data'][:, idx_sel] + else: + vertno = _get_vertno(forward['src']) + G = forward['sol']['data'] + + if isinstance(evraw, Raw): + M, times = evraw[picks, start:stop] + elif isinstance(evraw, Evoked): + M = evraw.data[picks, start:stop] + times = evraw.times[start:stop] + else: + raise ValueError('evraw has to be of type Evoked or Raw') # Handle SSPs - proj, ncomp, _ = make_projector(evoked.info['projs'], evoked.ch_names) + proj, ncomp, _ = make_projector(evraw.info['projs'], ch_names) M = np.dot(proj, M) G = np.dot(proj, G) # Handle whitening + data covariance - W, _ = compute_whitener(noise_cov, evoked.info, picks) + W, _ = compute_whitener(noise_cov, evraw.info, picks) # whiten data and leadfield M = np.dot(W, M) G = np.dot(W, G) # Apply SSPs + whitener to data covariance + data_cov = pick_channels_cov(data_cov, include=ch_names) Cm = data_cov['data'] Cm = np.dot(proj, np.dot(Cm, proj.T)) Cm = np.dot(W, np.dot(Cm, W.T)) @@ -105,10 +115,109 @@ def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01): sol /= noise_norm[:, None] - tstep = 1.0 / evoked.info['sfreq'] - tmin = float(evoked.first) / evoked.info['sfreq'] - vertno = _get_vertno(forward['src']) + tstep = 1.0 / evraw.info['sfreq'] + tmin = times[0] stc = _make_stc(sol, tmin, tstep, vertno) print '[done]' return stc + + +def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None, + start=None, stop=None): + """Linearly Constrained Minimum Variance (LCMV) beamformer. + + Compute Linearly Constrained Minimum Variance (LCMV) beamformer + on evoked data. + + NOTE : This implementation is heavilly tested so please + report any issue or suggestions. + + Parameters + ---------- + evoked : Evoked + Evoked data to invert + forward : dict + Forward operator + noise_cov : Covariance + The noise covariance + data_cov : Covariance + The data covariance + reg : float + The regularization for the whitened data covariance. + label : Label + Restricts the LCMV solution to a given label + start : int + Index of first time sample (index not time is seconds) + stop : int + Index of first time sample not to include (index not time is seconds) + + 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 + """ + + stc = _lcmv_any(evoked, forward, noise_cov, data_cov, reg, label, + start, stop, None) + + return stc + + +def lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=None, + start=None, stop=None, picks=None): + """Linearly Constrained Minimum Variance (LCMV) beamformer. + + Compute Linearly Constrained Minimum Variance (LCMV) beamformer + on raw data. + + NOTE : This implementation is heavilly tested so please + report any issue or suggestions. + + Parameters + ---------- + raw : mne.fiff.Raw + Raw data to invert + forward : dict + Forward operator + noise_cov : Covariance + The noise covariance + data_cov : Covariance + The data covariance + reg : float + The regularization for the whitened data covariance. + label : Label + Restricts the LCMV solution to a given label + start : int + Index of first time sample (index not time is seconds) + stop : int + Index of first time sample not to include (index not time is seconds) + picks: aray of int + Channel indices in raw to use for beamforming (if None all channels + are used) + + 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 + """ + + stc = _lcmv_any(raw, forward, noise_cov, data_cov, reg, label, start, stop, + picks) + + return stc + diff --git a/mne/beamformer/tests/test_lcmv.py b/mne/beamformer/tests/test_lcmv.py index 99e03e7..69ddd82 100644 --- a/mne/beamformer/tests/test_lcmv.py +++ b/mne/beamformer/tests/test_lcmv.py @@ -2,11 +2,11 @@ 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 +from numpy.testing import assert_array_almost_equal import mne from mne.datasets import sample -from mne.beamformer import lcmv +from mne.beamformer import lcmv, lcmv_raw examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') @@ -32,7 +32,7 @@ events = mne.read_events(fname_event) def test_lcmv(): - """Test LCMV + """Test LCMV with evoked data """ event_id, tmin, tmax = 1, -0.2, 0.2 @@ -41,8 +41,10 @@ def test_lcmv(): # 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) + 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, @@ -64,3 +66,43 @@ def test_lcmv(): assert_true(0.09 < tmax < 0.1) assert_true(2. < np.max(max_stc) < 3.) + + +def test_lcmv_raw(): + """Test LCMV with raw data + """ + tmin, tmax = 0, 20 + # 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) + + noise_cov = mne.read_cov(fname_cov) + noise_cov = mne.cov.regularize(noise_cov, raw.info, + mag=0.05, grad=0.05, eeg=0.1, proj=True) + + start, stop = raw.time_to_index(tmin, tmax) + + # use only the left-temporal MEG channels for LCMV + picks = mne.fiff.pick_types(raw.info, meg=True, exclude=raw.info['bads'], + selection=left_temporal_channels) + + data_cov = mne.compute_raw_data_covariance(raw, tmin=tmin, tmax=tmax) + + stc = lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=label, + start=start, stop=stop, picks=picks) + + assert_array_almost_equal(np.array([tmin, tmax]), + np.array([stc.times[0], stc.times[-1]]), + decimal=2) + + # make sure we get an stc with vertices only in the lh + vertno = [forward['src'][0]['vertno'], forward['src'][1]['vertno']] + assert_true(len(stc.vertno[0]) == len(np.intersect1d(vertno[0], + label['vertices']))) + assert_true(len(stc.vertno[1]) == 0) + # TODO: test more things diff --git a/mne/cov.py b/mne/cov.py index c63a831..3c2c30c 100644 --- a/mne/cov.py +++ b/mne/cov.py @@ -133,7 +133,7 @@ def read_cov(fname): # Estimate from data def compute_raw_data_covariance(raw, tmin=None, tmax=None, tstep=0.2, - reject=None, flat=None): + reject=None, flat=None, picks=None): """Estimate noise covariance matrix from a continuous segment of raw data It is typically useful to estimate a noise covariance @@ -164,6 +164,9 @@ def compute_raw_data_covariance(raw, tmin=None, tmax=None, tstep=0.2, Rejection parameters based on flatness of signal Valid keys are 'grad' | 'mag' | 'eeg' | 'eog' | 'ecg' If flat is None then no rejection is done. + picks : array of int + Indices of channels to include (if None, all channels + are used) Returns ------- @@ -180,8 +183,12 @@ def compute_raw_data_covariance(raw, tmin=None, tmax=None, tstep=0.2, stop = int(ceil(tmax * sfreq)) step = int(ceil(tstep * raw.info['sfreq'])) + if picks is None: + picks_data = pick_types(raw.info, meg=True, eeg=True, eog=False) + else: + picks_data = picks + picks = pick_types(raw.info, meg=True, eeg=True, eog=True) - picks_data = pick_types(raw.info, meg=True, eeg=True, eog=False) idx = [list(picks).index(k) for k in picks_data] data = 0 diff --git a/mne/tests/test_cov.py b/mne/tests/test_cov.py index 9341330..353208b 100644 --- a/mne/tests/test_cov.py +++ b/mne/tests/test_cov.py @@ -9,7 +9,7 @@ from ..cov import regularize from .. import read_cov, Epochs, merge_events, \ find_events, compute_raw_data_covariance, \ compute_covariance -from ..fiff import Raw, pick_channels_cov +from ..fiff import Raw, pick_channels_cov, pick_channels cov_fname = op.join(op.dirname(__file__), '..', 'fiff', 'tests', 'data', 'test-cov.fif') @@ -56,6 +56,13 @@ def test_cov_estimation_on_raw_segment(): assert_true((linalg.norm(cov.data - cov_read.data, ord='fro') / linalg.norm(cov.data, ord='fro')) < 1e-5) + # test with a subset of channels + picks = pick_channels(raw.ch_names, include=raw.ch_names[:5]) + cov = compute_raw_data_covariance(raw, picks=picks) + assert_true(cov_mne.ch_names[:5] == cov.ch_names) + assert_true(linalg.norm(cov.data - cov_mne.data[picks][:, picks], + ord='fro') / linalg.norm(cov.data, ord='fro')) < 1e-6 + def test_cov_estimation_with_triggers(): """Estimate raw with triggers -- 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
