[med-svn] [python-mne] 333/353: ENH: lcmv_epochs()
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:26 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 3cf42669c00367faf3b23358edcbbb52fd864202
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
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
More information about the debian-med-commit
mailing list