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 11ca75b38a22b7dfa2d48886401b61b2d2af6162 Author: Alexandre Gramfort <[email protected]> Date: Tue Nov 8 22:57:17 2011 -0500 ENH : Sign flip computation for robust label average of signed values --- doc/source/whats_new.rst | 8 +++++ .../plot_compute_mne_inverse_epochs_in_label.py | 13 +++++++- mne/__init__.py | 2 +- mne/label.py | 39 ++++++++++++++++++++++ mne/minimum_norm/tests/test_inverse.py | 21 ++++++++---- 5 files changed, 74 insertions(+), 9 deletions(-) diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst index 35650f3..2ef2950 100644 --- a/doc/source/whats_new.rst +++ b/doc/source/whats_new.rst @@ -3,6 +3,14 @@ What's new .. _changes_0_2: +Current +----------- + +Changelog +~~~~~~~~~ + + - Sign flip computation for robust label average of signed values by `Alex Gramfort`_. + Version 0.2 ----------- diff --git a/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py b/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py index 3975282..473c64e 100644 --- a/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py +++ b/examples/inverse/plot_compute_mne_inverse_epochs_in_label.py @@ -14,6 +14,7 @@ to a brain label. print __doc__ +import numpy as np import pylab as pl import mne from mne.datasets import sample @@ -57,9 +58,19 @@ stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM, label, data = sum(stc.data for stc in stcs) / len(stcs) +# compute sign flip to avoid signal cancelation when averaging signed values +flip = mne.label_sign_flip(label, inverse_operator['src']) + +label_mean = np.mean(data, axis=0) +label_mean_flip = np.mean(flip[:, np.newaxis] * data, axis=0) + ############################################################################### # View activation time-series -pl.plot(1e3 * stcs[0].times, data.T) +h0 = pl.plot(1e3 * stcs[0].times, data.T, 'k') +h1 = pl.plot(1e3 * stcs[0].times, label_mean, 'r', linewidth=3) +h2 = pl.plot(1e3 * stcs[0].times, label_mean_flip, 'g', linewidth=3) +pl.legend((h0[0], h1, h2), ('all dipoles in label', 'mean', + 'mean with sign flip')) pl.xlabel('time (ms)') pl.ylabel('dSPM value') pl.show() diff --git a/mne/__init__.py b/mne/__init__.py index 750bf56..45ee7fb 100644 --- a/mne/__init__.py +++ b/mne/__init__.py @@ -11,7 +11,7 @@ from .source_estimate import read_stc, write_stc, SourceEstimate, morph_data, \ from .surface import read_bem_surfaces, read_surface from .source_space import read_source_spaces from .epochs import Epochs -from .label import label_time_courses, read_label +from .label import label_time_courses, read_label, label_sign_flip from .misc import parse_config, read_reject_parameters from .transforms import transform_coordinates from .proj import read_proj diff --git a/mne/label.py b/mne/label.py index da5fcad..13320fe 100644 --- a/mne/label.py +++ b/mne/label.py @@ -1,4 +1,5 @@ import numpy as np +from scipy import linalg from .source_estimate import read_stc @@ -80,3 +81,41 @@ def label_time_courses(labelfile, stcfile): times = stc['tmin'] + stc['tstep'] * np.arange(stc['data'].shape[1]) return values, times, vertices + + +def label_sign_flip(label, src): + """Compute sign for label averaging + + Parameters + ---------- + label : dict + A label read with the read_label function + src : list of dict + The source space over which is defined the label + + Returns + ------- + flip : array + Sign flip vector (contains 1 or -1) + """ + if len(src) != 2: + raise ValueError('Only source spaces with 2 hemisphers are accepted') + + lh_vertno = src[0]['vertno'] + rh_vertno = src[1]['vertno'] + + # get source orientations + if label['hemi'] == 'lh': + vertno_sel = np.intersect1d(lh_vertno, label['vertices']) + ori = src[0]['nn'][vertno_sel] + elif label['hemi'] == 'rh': + vertno_sel = np.intersect1d(rh_vertno, label['vertices']) + ori = src[1]['nn'][vertno_sel] + else: + raise Exception("Unknown hemisphere type") + + _, _, Vh = linalg.svd(ori, full_matrices=False) + + # Comparing to the direction of the first right singular vector + flip = np.sign(np.dot(ori, Vh[:, 0])) + return flip diff --git a/mne/minimum_norm/tests/test_inverse.py b/mne/minimum_norm/tests/test_inverse.py index d5fe27d..1fdde52 100644 --- a/mne/minimum_norm/tests/test_inverse.py +++ b/mne/minimum_norm/tests/test_inverse.py @@ -4,7 +4,7 @@ from numpy.testing import assert_array_almost_equal, assert_equal from nose.tools import assert_true from ...datasets import sample -from ...label import read_label +from ...label import read_label, label_sign_flip from ...event import read_events from ...epochs import Epochs from ...source_estimate import SourceEstimate @@ -48,7 +48,6 @@ def test_inverse_operator(): With and without precomputed inverse operator. """ evoked = fiff.Evoked(fname_data, setno=0, baseline=(None, 0)) - from copy import deepcopy stc = apply_inverse(evoked, inverse_operator, lambda2, dSPM=False) @@ -113,12 +112,20 @@ def test_apply_mne_inverse_epochs(): reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6) flat = dict(grad=1e-15, mag=1e-15) - events = read_events(fname_event)[:3] + events = read_events(fname_event)[:15] epochs = Epochs(raw, events, event_id, tmin, tmax, picks=picks, baseline=(None, 0), reject=reject, flat=flat) stcs = apply_inverse_epochs(epochs, inverse_operator, lambda2, dSPM, - label=label) + label=label, pick_normal=True) - assert_true(len(stcs) == 1) - assert_true(np.all(stcs[0].data > 0)) - assert_true(np.all(stcs[0].data < 42)) + assert_true(len(stcs) == 4) + assert_true(3 < stcs[0].data.max() < 10) + + + data = sum(stc.data for stc in stcs) / len(stcs) + flip = label_sign_flip(label, inverse_operator['src']) + + label_mean = np.mean(data, axis=0) + label_mean_flip = np.mean(flip[:, np.newaxis] * data, axis=0) + + assert_true(label_mean.max() < label_mean_flip.max()) -- 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
