[med-svn] [python-mne] 07/353: ENH : Sign flip computation for robust label average of signed values

[Yaroslav Halchenko]

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 <alexandre.gramfort at inria.fr>
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