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 3cf42669c00367faf3b23358edcbbb52fd864202 Author: Martin Luessi <[email protected]> Date: Tue Jul 31 10:52:12 2012 -0400 ENH: lcmv_epochs() --- doc/source/whats_new.rst | 2 +- mne/beamformer/__init__.py | 2 +- mne/beamformer/_lcmv.py | 152 +++++++++++++++++++++++++++----------- mne/beamformer/tests/test_lcmv.py | 23 +++++- mne/minimum_norm/inverse.py | 27 +------ mne/source_space.py | 43 ++++++++++- 6 files changed, 173 insertions(+), 76 deletions(-) diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 157b60f..7d270c3 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -19,7 +19,7 @@ Changelog - Backporting scipy.signa.firwin2 so filtering works with old scipy by `Alex Gramfort`_. - - LCMV Beamformer for evoked and raw data by `Alex Gramfort`_ and `Martin Luessi`_. + - LCMV Beamformer for evoked data, single trials, and raw data by `Alex Gramfort`_ and `Martin Luessi`_. - Add support for reading named channel selections by `Martin Luessi`_. diff --git a/mne/beamformer/__init__.py b/mne/beamformer/__init__.py index c2b177a..b24e865 100644 --- a/mne/beamformer/__init__.py +++ b/mne/beamformer/__init__.py @@ -1,4 +1,4 @@ """Artifacts finding/correction related functions """ -from ._lcmv import lcmv, lcmv_raw +from ._lcmv import lcmv, lcmv_epochs, lcmv_raw diff --git a/mne/beamformer/_lcmv.py b/mne/beamformer/_lcmv.py index 50e602c..8c8e799 100644 --- a/mne/beamformer/_lcmv.py +++ b/mne/beamformer/_lcmv.py @@ -13,16 +13,19 @@ from ..fiff.proj import make_projector 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 +from ..source_space import label_src_vertno_sel def _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, label=None, picks=None): - """ LCMV beamformer for evoked or raw data + """ LCMV beamformer for evoked data, single epochs, and raw data Parameters ---------- - data : array - Sensor space data + data : array or list / iterable + Sensor space data. If data.ndim == 2 a single observation is assumed + and a single stc is returned. If data.ndim == 3 or if data is + a list / iterable, a list of stc's is returned. info : dict Measurement info tmin : float @@ -42,7 +45,7 @@ def _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, Returns ------- - stc : dict + stc : SourceEstimate (or list of SourceEstimate) Source time courses """ @@ -51,9 +54,6 @@ def _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, if picks is None: picks = pick_types(info, meg=True, eeg=True, exclude=info['bads']) - if len(data) != len(picks): - raise ValueError('data and picks must have the same length') - ch_names = [info['ch_names'][k] for k in picks] # restrict forward solution to selected channels @@ -61,51 +61,33 @@ def _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, # 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") + vertno, src_sel = label_src_vertno_sel(label, forward['src']) 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 + src_sel = 3 * src_sel + src_sel = np.c_[src_sel, src_sel + 1, src_sel + 2] + src_sel = src_sel.ravel() - G = forward['sol']['data'][:, idx_sel] + G = forward['sol']['data'][:, src_sel] else: vertno = _get_vertno(forward['src']) G = forward['sol']['data'] # Handle SSPs proj, ncomp, _ = make_projector(info['projs'], ch_names) - M = np.dot(proj, data) G = np.dot(proj, G) # Handle whitening + data covariance - W, _ = compute_whitener(noise_cov, info, picks) + whitener, _ = compute_whitener(noise_cov, info, picks) - # whiten data and leadfield - M = np.dot(W, M) - G = np.dot(W, G) + # whiten the leadfield + G = np.dot(whitener, 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)) + Cm = np.dot(whitener, np.dot(Cm, whitener.T)) # Cm += reg * np.trace(Cm) / len(Cm) * np.eye(len(Cm)) Cm_inv = linalg.pinv(Cm, reg) @@ -127,19 +109,44 @@ def _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, noise_norm = np.sqrt(noise_norm) - sol = np.dot(W, M) + if isinstance(data, np.ndarray) and data.ndim == 2: + data = [data] + return_single = True + else: + return_single = False + stcs = [] - if is_free_ori: - print 'combining the current components...', - sol = combine_xyz(sol) + for i, M in enumerate(data): + if len(M) != len(picks): + raise ValueError('data and picks must have the same length') - sol /= noise_norm[:, None] + if not return_single: + print "Processing epoch : %d" % (i + 1) - tstep = 1.0 / info['sfreq'] - stc = _make_stc(sol, tmin, tstep, vertno) - print '[done]' + # SSP and whitening + M = np.dot(proj, M) + M = np.dot(whitener, M) - return stc + # project to source space using beamformer weights + 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 / info['sfreq'] + stc = _make_stc(sol, tmin, tstep, vertno) + + if not return_single: + stcs.append(stc) + print '[done]' + + if return_single: + return stc + else: + return stcs def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None): @@ -168,7 +175,7 @@ def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None): Returns ------- - stc : dict + stc : SourceEstimate Source time courses Notes @@ -189,6 +196,61 @@ def lcmv(evoked, forward, noise_cov, data_cov, reg=0.01, label=None): return stc +def lcmv_epochs(epochs, forward, noise_cov, data_cov, reg=0.01, label=None): + """Linearly Constrained Minimum Variance (LCMV) beamformer. + + Compute Linearly Constrained Minimum Variance (LCMV) beamformer + on single trial data. + + NOTE : This implementation is heavilly tested so please + report any issue or suggestions. + + Parameters + ---------- + epochs: Epochs + Single trial epochs + 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 + + Returns + ------- + stc: list of SourceEstimate + The source estimates for all epochs + + 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 + """ + + info = epochs.info + tmin = epochs.times[0] + + # use only the good data channels + def _epochs_data(epochs, picks): + for e in epochs: + yield e[picks, :] + + picks = pick_types(info, meg=True, eeg=True, exclude=info['bads']) + + data = _epochs_data(epochs, picks) + + stcs = _apply_lcmv(data, info, tmin, forward, noise_cov, data_cov, reg, + label) + + return stcs + + 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. @@ -223,7 +285,7 @@ def lcmv_raw(raw, forward, noise_cov, data_cov, reg=0.01, label=None, Returns ------- - stc : dict + stc : SourceEstimate Source time courses Notes diff --git a/mne/beamformer/tests/test_lcmv.py b/mne/beamformer/tests/test_lcmv.py index 69ddd82..c938ecb 100644 --- a/mne/beamformer/tests/test_lcmv.py +++ b/mne/beamformer/tests/test_lcmv.py @@ -6,7 +6,7 @@ from numpy.testing import assert_array_almost_equal import mne from mne.datasets import sample -from mne.beamformer import lcmv, lcmv_raw +from mne.beamformer import lcmv, lcmv_epochs, lcmv_raw examples_folder = op.join(op.dirname(__file__), '..', '..', '..', 'examples') @@ -28,11 +28,13 @@ 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) +forward_fixed = mne.read_forward_solution(fname_fwd, force_fixed=True, + surf_ori=True) events = mne.read_events(fname_event) def test_lcmv(): - """Test LCMV with evoked data + """Test LCMV with evoked data and single trials """ event_id, tmin, tmax = 1, -0.2, 0.2 @@ -67,6 +69,23 @@ def test_lcmv(): assert_true(0.09 < tmax < 0.1) assert_true(2. < np.max(max_stc) < 3.) + # Now test single trial using fixed orientation forward solution + # so we can compare it to the evoked solution + stcs = lcmv_epochs(epochs, forward_fixed, noise_cov, data_cov, reg=0.01) + + epochs.drop_bad_epochs() + assert_true(len(epochs.events) == len(stcs)) + + # average the single trial estimates + stc_avg = np.zeros_like(stc.data) + for stc_single in stcs: + stc_avg += stc_single.data + stc_avg /= len(stcs) + + # compare it to the solution using evoked with fixed orientation + stc_fixed = lcmv(evoked, forward_fixed, noise_cov, data_cov, reg=0.01) + assert_array_almost_equal(stc_avg, stc_fixed.data) + def test_lcmv_raw(): """Test LCMV with raw data diff --git a/mne/minimum_norm/inverse.py b/mne/minimum_norm/inverse.py index c810a5c..db71ace 100644 --- a/mne/minimum_norm/inverse.py +++ b/mne/minimum_norm/inverse.py @@ -28,7 +28,7 @@ from ..forward import compute_depth_prior, compute_depth_prior_fixed, \ is_fixed_orient, compute_orient_prior from ..source_space import read_source_spaces_from_tree, \ find_source_space_hemi, _get_vertno, \ - write_source_spaces + write_source_spaces, label_src_vertno_sel from ..transforms import invert_transform, transform_source_space_to from ..source_estimate import SourceEstimate @@ -559,7 +559,7 @@ def _assemble_kernel(inv, label, method, pick_normal): vertno = _get_vertno(src) if label is not None: - vertno, src_sel = _get_label_sel(label, inv) + vertno, src_sel = label_src_vertno_sel(label, inv['src']) if method != "MNE": noise_norm = noise_norm[src_sel] @@ -623,29 +623,6 @@ def _make_stc(sol, tmin, tstep, vertno): return stc -def _get_label_sel(label, inv): - src = inv['src'] - - if src[0]['type'] != 'surf': - return Exception('Label are only supported with surface source spaces') - - vertno = [src[0]['vertno'], src[1]['vertno']] - - if label['hemi'] == 'lh': - vertno_sel = np.intersect1d(vertno[0], label['vertices']) - src_sel = np.searchsorted(vertno[0], vertno_sel) - vertno[0] = vertno_sel - vertno[1] = np.array([]) - elif label['hemi'] == 'rh': - vertno_sel = np.intersect1d(vertno[1], label['vertices']) - src_sel = np.searchsorted(vertno[1], vertno_sel) + len(vertno[0]) - vertno[0] = np.array([]) - vertno[1] = vertno_sel - else: - raise Exception("Unknown hemisphere type") - - return vertno, src_sel - def _check_method(method, dSPM): if dSPM is not None: diff --git a/mne/source_space.py b/mne/source_space.py index 9b10ea1..a75df95 100644 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -182,7 +182,6 @@ def _read_one_source_space(fid, this): if tag is not None: res['mri_depth'] = int(tag.data) - tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NPOINTS) if tag is None: raise ValueError('Number of vertices not found') @@ -338,6 +337,45 @@ def find_source_space_hemi(src): return hemi +def label_src_vertno_sel(label, src): + """ Find vertex numbers and indices from label + + Parameters + ---------- + label : Label + Source space label + src : dict + Source space + + Returns + ------- + vertno : list of length 2 + Vertex numbers for lh and rh + src_sel : array of int (len(idx) = len(vertno[0]) + len(vertno[1])) + Indices of the selected vertices in sourse space + """ + + if src[0]['type'] != 'surf': + return Exception('Label are only supported with surface source spaces') + + vertno = [src[0]['vertno'], src[1]['vertno']] + + if label['hemi'] == 'lh': + vertno_sel = np.intersect1d(vertno[0], label['vertices']) + src_sel = np.searchsorted(vertno[0], vertno_sel) + vertno[0] = vertno_sel + vertno[1] = np.array([]) + elif label['hemi'] == 'rh': + vertno_sel = np.intersect1d(vertno[1], label['vertices']) + src_sel = np.searchsorted(vertno[1], vertno_sel) + len(vertno[0]) + vertno[0] = np.array([]) + vertno[1] = vertno_sel + else: + raise Exception("Unknown hemisphere type") + + return vertno, src_sel + + def _get_vertno(src): vertno = list() for s in src: @@ -408,7 +446,8 @@ def _write_one_source_space(fid, this): write_float_matrix(fid, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS, this['nn']) if this['ntri'] > 0: - write_int_matrix(fid, FIFF.FIFF_MNE_SOURCE_SPACE_TRIANGLES, this['tris'] + 1) + write_int_matrix(fid, FIFF.FIFF_MNE_SOURCE_SPACE_TRIANGLES, + this['tris'] + 1) # Which vertices are active write_int(fid, FIFF.FIFF_MNE_SOURCE_SPACE_NUSE, this['nuse']) -- 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
