[med-svn] [python-mne] 326/353: ENH: LCMV for raw data
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:25:25 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 db09465a5d83070d4aa1ef8c0dc8788fbeb8f3d9
Author: Martin Luessi <mluessi at nmr.mgh.harvard.edu>
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
More information about the debian-med-commit
mailing list