This is an automated email from the git hooks/post-receive script. yoh pushed a commit to annotated tag v0.1 in repository python-mne.
commit 33b704536cf9e61c0131d5b6d7171d47741d9ebd Author: Alexandre Gramfort <[email protected]> Date: Tue Jun 7 12:21:35 2011 -0400 ENH: adding tools to generate source space connectivity --- mne/__init__.py | 5 ++- mne/source_estimate.py | 73 ++++++++++++++++++++++++++++++++++++++- mne/stats/cluster_level.py | 34 ++++++++++-------- mne/tests/test_source_estimate.py | 28 ++++++++++++++- 4 files changed, 123 insertions(+), 17 deletions(-) diff --git a/mne/__init__.py b/mne/__init__.py index 03b5cb5..79fa389 100644 --- a/mne/__init__.py +++ b/mne/__init__.py @@ -4,7 +4,9 @@ from .cov import read_cov, write_cov, write_cov_file, Covariance, \ compute_raw_data_covariance, compute_covariance from .event import read_events, write_events, find_events from .forward import read_forward_solution -from .source_estimate import read_stc, write_stc, SourceEstimate, morph_data +from .source_estimate import read_stc, write_stc, SourceEstimate, morph_data, \ + spatio_temporal_src_connectivity, \ + spatio_temporal_tris_connectivity from .surface import read_bem_surfaces, read_surface from .source_space import read_source_spaces from .epochs import Epochs @@ -12,3 +14,4 @@ from .label import label_time_courses, read_label from .misc import parse_config, read_reject_parameters import fiff import artifacts +import stats diff --git a/mne/source_estimate.py b/mne/source_estimate.py index 02e61ad..b96f353 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -6,6 +6,7 @@ import os import copy import numpy as np +from scipy import sparse def read_stc(filename): @@ -259,7 +260,6 @@ def mesh_edges(tris): edges : sparse matrix The adjacency matrix """ - from scipy import sparse npoints = np.max(tris) + 1 ntris = len(tris) a, b, c = tris.T @@ -384,3 +384,74 @@ def morph_data(subject_from, subject_to, stc_from, grade=5, smooth=None, print '[done]' return stc_to + + +def spatio_temporal_src_connectivity(src, n_times): + """Compute connectivity for a source space activation over time + + Parameters + ---------- + src : source space + The source space. + + n_times : int + Number of time instants + + Returns + ------- + connectivity : sparse COO matrix + The connectivity matrix describing the spatio-temporal + graph structure. If N is the number of vertices in the + source space, the N first nodes in the graph are the + vertices are time 1, the nodes from 2 to 2N are the vertices + during time 2, etc. + + """ + if src[0]['use_tris'] is None: + raise Exception("The source space does not appear to be an ico " + "surface. Connectivity cannot be extracted from " + "non-ico source spaces.") + lh_tris = np.searchsorted(np.unique(src[0]['use_tris']), + src[0]['use_tris']) + rh_tris = np.searchsorted(np.unique(src[1]['use_tris']), + src[1]['use_tris']) + tris = np.concatenate((lh_tris, rh_tris + np.max(lh_tris) + 1)) + return spatio_temporal_tris_connectivity(tris, n_times) + + +def spatio_temporal_tris_connectivity(tris, n_times): + """Compute connectivity from triangles and time instants""" + edges = mesh_edges(tris).tocoo() + n_vertices = edges.shape[0] + print "-- number of connected vertices : %d" % n_vertices + nnz = edges.col.size + aux = n_vertices * np.arange(n_times)[:, None] * np.ones((1, nnz), np.int) + col = (edges.col[None, :] + aux).ravel() + row = (edges.row[None, :] + aux).ravel() + if n_times > 1: # add temporal edges + o = (n_vertices * np.arange(n_times - 1)[:, None] + + np.arange(n_vertices)[None, :]).ravel() + d = (n_vertices * np.arange(1, n_times)[:, None] + + np.arange(n_vertices)[None, :]).ravel() + row = np.concatenate((row, o, d)) + col = np.concatenate((col, d, o)) + data = np.ones(edges.data.size * n_times + 2 * n_vertices * (n_times - 1), + dtype=np.int) + connectivity = sparse.coo_matrix((data, (row, col)), + shape=(n_times * n_vertices, ) * 2) + return connectivity + + +def _get_ico_tris(grade): + """Get triangles for ico surface.""" + mne_root = os.environ.get('MNE_ROOT') + if mne_root is None: + raise Exception('Please set MNE_ROOT environment variable.') + ico_file_name = os.path.join(mne_root, 'share', 'mne', 'icos.fif') + surfaces = read_bem_surfaces(ico_file_name) + for s in surfaces: + if s['id'] == (9000 + grade): + ico = s + break + return ico['tris'] + \ No newline at end of file diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index 9ae9f4c..ae08ae7 100644 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -12,6 +12,24 @@ from scipy import stats, sparse, ndimage from .parametric import f_oneway +def _get_components(x_in, connectivity): + """get connected components from a mask and a connectivity matrix""" + from scikits.learn.utils._csgraph import cs_graph_components + mask = np.logical_and(x_in[connectivity.row], x_in[connectivity.col]) + data = connectivity.data[mask] + row = connectivity.row[mask] + col = connectivity.col[mask] + shape = connectivity.shape + idx = np.where(x_in)[0] + row = np.concatenate((row, idx)) + col = np.concatenate((col, idx)) + data = np.concatenate((data, np.ones(len(idx), dtype=data.dtype))) + connectivity = sparse.coo_matrix((data, (row, col)), shape=shape) + _, components = cs_graph_components(connectivity) + # print "-- number of components : %d" % np.unique(components).size + return components + + def _find_clusters(x, threshold, tail=0, connectivity=None): """For a given 1d-array (test statistic), find all clusters which are above/below a certain threshold. Returns a list of 2-tuples. @@ -68,21 +86,9 @@ def _find_clusters(x, threshold, tail=0, connectivity=None): if x.ndim > 1: raise Exception("Data should be 1D when using a connectivity " "to define clusters.") - from scikits.learn.utils._csgraph import cs_graph_components - mask = np.logical_and(x_in[connectivity.row], x_in[connectivity.col]) - if np.sum(mask) == 0: + if np.sum(x_in) == 0: return [], np.empty(0) - mask = np.logical_and(x_in[connectivity.row], x_in[connectivity.col]) - data = connectivity.data[mask] - row = connectivity.row[mask] - col = connectivity.col[mask] - shape = connectivity.shape - idx = np.where(x_in)[0] - row = np.concatenate((row, idx)) - col = np.concatenate((col, idx)) - data = np.concatenate((data, np.ones(len(idx), dtype=data.dtype))) - connectivity = sparse.coo_matrix((data, (row, col)), shape=shape) - _, components = cs_graph_components(connectivity) + components = _get_components(x_in, connectivity) labels = np.unique(components) clusters = list() sums = list() diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py index a84eb30..ef75230 100644 --- a/mne/tests/test_source_estimate.py +++ b/mne/tests/test_source_estimate.py @@ -1,10 +1,14 @@ import os.path as op import numpy as np -from numpy.testing import assert_array_almost_equal +from numpy.testing import assert_array_almost_equal, assert_array_equal import mne from mne.datasets import sample +from mne import stats +from mne.source_estimate import spatio_temporal_tris_connectivity, \ + spatio_temporal_src_connectivity + examples_folder = op.join(op.dirname(__file__), '..', '..', 'examples') data_path = sample.data_path(examples_folder) @@ -42,3 +46,25 @@ def test_morph_data(): mean_from = stc_from.data.mean(axis=0) mean_to = stc_to.data.mean(axis=0) assert np.corrcoef(mean_to, mean_from).min() > 0.99 + + +def test_spatio_temporal_tris_connectivity(): + """Test spatio-temporal connectivity""" + tris = np.array([[0, 1, 2], [3, 4, 5]]) + connectivity = spatio_temporal_tris_connectivity(tris, 2) + x = [1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1] + components = stats.cluster_level._get_components(np.array(x), connectivity) + assert_array_equal(components, + [0, 0, -2, -2, -2, -2, 0, -2, -2, -2, -2, 1]) + + +def test_spatio_temporal_src_connectivity(): + """Test spatio-temporal connectivity""" + tris = np.array([[0, 1, 2], [3, 4, 5]]) + src = [dict(), dict()] + connectivity = spatio_temporal_tris_connectivity(tris, 2) + src[0]['use_tris'] = np.array([[0, 1, 2]]) + src[1]['use_tris'] = np.array([[0, 1, 2]]) + connectivity2 = spatio_temporal_src_connectivity(src, 2) + assert_array_equal(connectivity.todense(), connectivity2.todense()) + -- 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
