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 4c12dd9b60ec08690cb7b98a1bb15bc9a97b27c0 Author: Alexandre Gramfort <[email protected]> Date: Tue Jul 17 10:58:15 2012 +0200 ENH : simulation code refactoring --- examples/plot_simulate_evoked_data.py | 3 +- mne/__init__.py | 2 +- mne/label.py | 148 ++++++++++++++++++++++- mne/simulation/__init__.py | 4 +- mne/simulation/evoked.py | 4 +- mne/simulation/source.py | 217 ++++------------------------------ mne/source_estimate.py | 30 +++++ 7 files changed, 207 insertions(+), 201 deletions(-) diff --git a/examples/plot_simulate_evoked_data.py b/examples/plot_simulate_evoked_data.py index 7d4c9b5..e3d2506 100644 --- a/examples/plot_simulate_evoked_data.py +++ b/examples/plot_simulate_evoked_data.py @@ -63,7 +63,8 @@ stc_data *= 100 * 1e-9 # use nAm as unit # time translation stc_data[1] = np.roll(stc_data[1], 80) -stc = generate_sparse_stc(fwd, labels, stc_data, tmin, tstep, random_state=0) +stc = generate_sparse_stc(fwd['src'], labels, stc_data, tmin, tstep, + random_state=0) ############################################################################### # Generate noisy evoked data diff --git a/mne/__init__.py b/mne/__init__.py index 0deffd2..352c2ca 100644 --- a/mne/__init__.py +++ b/mne/__init__.py @@ -17,7 +17,7 @@ from .surface import read_bem_surfaces, read_surface, write_bem_surface from .source_space import read_source_spaces from .epochs import Epochs from .label import label_time_courses, read_label, label_sign_flip, \ - write_label, stc_to_label + write_label, stc_to_label, grow_labels from .misc import parse_config, read_reject_parameters from .transforms import transform_coordinates from .proj import read_proj, write_proj, compute_proj_epochs, \ diff --git a/mne/label.py b/mne/label.py index 40eae0d..40d470c 100644 --- a/mne/label.py +++ b/mne/label.py @@ -1,8 +1,13 @@ +# Authors: Alexandre Gramfort <[email protected]> +# Martin Luessi <[email protected]> +# +# License: BSD (3-clause) + import os import numpy as np from scipy import linalg -from .source_estimate import read_stc, mesh_edges +from .source_estimate import read_stc, mesh_edges, mesh_dist from .surface import read_surface @@ -224,3 +229,144 @@ def stc_to_label(stc, src, smooth=5): labels.append(label) return labels + + +def _verts_within_dist(graph, source, max_dist): + """Find all vertices wihin a maximum geodesic distance from source + + Parameters: + ----------- + graph : scipy.sparse.csr_matrix + Sparse matrix with distances between adjacent vertices + source : int + Source vertex + max_dist: float + Maximum geodesic distance + + Returns: + -------- + verts : array + Vertices within max_dist + dist : array + Distances from source vertex + """ + dist_map = {} + dist_map[source] = 0 + verts_added_last = [source] + # add neighbors until no more neighbors within max_dist can be found + while len(verts_added_last) > 0: + verts_added = [] + for i in verts_added_last: + v_dist = dist_map[i] + row = graph[i, :] + neighbor_vert = row.indices + neighbor_dist = row.data + for j, d in zip(neighbor_vert, neighbor_dist): + n_dist = v_dist + d + if j in dist_map: + if n_dist < dist_map[j]: + dist_map[j] = n_dist + else: + if n_dist <= max_dist: + dist_map[j] = n_dist + # we found a new vertex within max_dist + verts_added.append(j) + verts_added_last = verts_added + + verts = np.sort(np.array(dist_map.keys(), dtype=np.int)) + dist = np.array([dist_map[v] for v in verts]) + + return verts, dist + + +def grow_labels(subject, seeds, extents, hemis, subjects_dir=None): + """Generate circular labels in source space with region growing + + This function generates a number of labels in source space by growing + regions starting from the vertices defined in "seeds". For each seed, a + label is generated containing all vertices within a maximum geodesic + distance on the white matter surface from the seed. + + Note: "extents" and "hemis" can either be arrays with the same length as + seeds, which allows using a different extent and hemisphere for each + label, or integers, in which case the same extent and hemisphere is + used for each label. + + Parameters: + ---------- + subject : string + Name of the subject as in SUBJECTS_DIR + seeds : array or int + Seed vertex numbers + extents : array or float + Extents (radius in mm) of the labels + hemis : array or int + Hemispheres to use for the labels (0: left, 1: right) + subjects_dir : string + Path to SUBJECTS_DIR if not set in the environment + + Returns: + -------- + labels : list + List with lables. Each label is a dictionary with keys: + comment Comment with seed vertex and extent + hemi Hemisphere of the label ('lh', or 'rh') + vertices Vertex indices (0 based) + pos Locations in millimeters + values Distances in millimeters from seed + """ + if subjects_dir is None: + if 'SUBJECTS_DIR' in os.environ: + subjects_dir = os.environ['SUBJECTS_DIR'] + else: + raise ValueError('SUBJECTS_DIR environment variable not set') + + # make sure the inputs are arrays + seeds = np.atleast_1d(seeds) + extents = np.atleast_1d(extents) + hemis = np.atleast_1d(hemis) + + n_seeds = len(seeds) + + if len(extents) != 1 and len(extents) != n_seeds: + raise ValueError('The extents parameter has to be of length 1 or ' + 'len(seeds)') + + if len(hemis) != 1 and len(hemis) != n_seeds: + raise ValueError('The hemis parameter has to be of length 1 or ' + 'len(seeds)') + + # make the arrays the same length as seeds + if len(extents) == 1: + extents = np.tile(extents, n_seeds) + + if len(hemis) == 1: + hemis = np.tile(hemis, n_seeds) + + hemis = ['lh' if h == 0 else 'rh' for h in hemis] + + # load the surfaces and create the distance graphs + tris, vert, dist = {}, {}, {} + for hemi in set(hemis): + surf_fname = os.path.join(subjects_dir, subject, 'surf', + hemi + '.white') + vert[hemi], tris[hemi] = read_surface(surf_fname) + dist[hemi] = mesh_dist(tris[hemi], vert[hemi]) + + # create the patches + labels = [] + for seed, extent, hemi in zip(seeds, extents, hemis): + label_verts, label_dist = _verts_within_dist(dist[hemi], seed, extent) + + # create a label + label = dict() + label['comment'] = 'Circular label: seed=%d, extent=%0.1fmm' %\ + (seed, extent) + label['vertices'] = label_verts + label['pos'] = vert[hemi][label_verts] + label['values'] = label_dist + label['hemi'] = hemi + + labels.append(label) + + return labels diff --git a/mne/simulation/__init__.py b/mne/simulation/__init__.py index 70c5eba..7076a90 100644 --- a/mne/simulation/__init__.py +++ b/mne/simulation/__init__.py @@ -3,6 +3,4 @@ from .evoked import generate_evoked -from .source import select_source_in_label, generate_sparse_stc, \ - circular_source_labels - +from .source import select_source_in_label, generate_sparse_stc diff --git a/mne/simulation/evoked.py b/mne/simulation/evoked.py index cdecb32..09813fa 100644 --- a/mne/simulation/evoked.py +++ b/mne/simulation/evoked.py @@ -1,6 +1,6 @@ # Authors: Alexandre Gramfort <[email protected]> +# Daniel Strohmeier <[email protected]> # Martin Luessi <[email protected]> -# Matti Hamalainen <[email protected]> # # License: BSD (3-clause) import copy @@ -123,5 +123,3 @@ def add_noise_evoked(evoked, noise, snr, tmin=None, tmax=None): noise.data = 10 ** ((tmp - float(snr)) / 20) * noise.data evoked.data += noise.data return evoked - - diff --git a/mne/simulation/source.py b/mne/simulation/source.py index ef54171..bddb5ed 100644 --- a/mne/simulation/source.py +++ b/mne/simulation/source.py @@ -1,26 +1,21 @@ # Authors: Alexandre Gramfort <[email protected]> # Martin Luessi <[email protected]> -# Matti Hamalainen <[email protected]> +# Daniel Strohmeier <[email protected]> # # License: BSD (3-clause) -import os import numpy as np -from scipy.sparse import csr_matrix - from ..minimum_norm.inverse import _make_stc from ..utils import check_random_state -from ..surface import read_surface -from ..source_estimate import mesh_edges -def select_source_in_label(fwd, label, random_state=None): +def select_source_in_label(src, label, random_state=None): """Select source positions using a label Parameters ---------- - fwd : dict - A forward solution + src : list of dict + The source space label : dict the label (read with mne.read_label) random_state : None | int | np.random.RandomState @@ -39,18 +34,18 @@ def select_source_in_label(fwd, label, random_state=None): rng = check_random_state(random_state) if label['hemi'] == 'lh': - src_sel_lh = np.intersect1d(fwd['src'][0]['vertno'], label['vertices']) + src_sel_lh = np.intersect1d(src[0]['vertno'], label['vertices']) idx_select = rng.randint(0, len(src_sel_lh), 1) lh_vertno.append(src_sel_lh[idx_select][0]) else: - src_sel_rh = np.intersect1d(fwd['src'][1]['vertno'], label['vertices']) + src_sel_rh = np.intersect1d(src[1]['vertno'], label['vertices']) idx_select = rng.randint(0, len(src_sel_rh), 1) rh_vertno.append(src_sel_rh[idx_select][0]) return lh_vertno, rh_vertno -def generate_sparse_stc(fwd, labels, stc_data, tmin, tstep, random_state=0): +def generate_sparse_stc(src, labels, stc_data, tmin, tstep, random_state=0): """Generate sparse sources time courses from waveforms and labels This function randomly selects a single vertex in each label and assigns @@ -58,8 +53,8 @@ def generate_sparse_stc(fwd, labels, stc_data, tmin, tstep, random_state=0): Parameters ---------- - fwd : dict - A forward solution + src : list of dict + The source space labels : list of dict The labels stc_data : array (shape: len(labels) x n_times) @@ -81,18 +76,25 @@ def generate_sparse_stc(fwd, labels, stc_data, tmin, tstep, random_state=0): rng = check_random_state(random_state) vertno = [[], []] - for label in labels: - lh_vertno, rh_vertno = select_source_in_label(fwd, label, rng) + lh_data = list() + rh_data = list() + for label_data, label in zip(stc_data, labels): + lh_vertno, rh_vertno = select_source_in_label(src, label, rng) vertno[0] += lh_vertno vertno[1] += rh_vertno + if len(lh_vertno) != 0: + lh_data.append(label_data) + elif len(rh_vertno) != 0: + rh_data.append(label_data) + else: + raise ValueError('No vertno found.') vertno = map(np.array, vertno) - # XXX there is a bug here, the ordering of the stc_data is assumed - # to be the same as in vertno (not the same as labels, left comes first) - stc = _make_stc(stc_data, tmin, tstep, vertno) + data = np.r_[lh_data, rh_data] + stc = _make_stc(data, tmin, tstep, vertno) return stc -def generate_stc(fwd, labels, stc_data, tmin, tstep, value_fun=None): +def generate_stc(src, labels, stc_data, tmin, tstep, value_fun=None): """Generate sources time courses from waveforms and labels This function generates a source estimate with extended sources by @@ -111,8 +113,8 @@ def generate_stc(fwd, labels, stc_data, tmin, tstep, value_fun=None): Parameters ---------- - fwd : dict - A forward solution + src : list of dict + The source space labels : list of dict The labels stc_data : array (shape: len(labels) x n_times) @@ -139,7 +141,7 @@ def generate_stc(fwd, labels, stc_data, tmin, tstep, value_fun=None): hemi_to_ind['lh'], hemi_to_ind['rh'] = 0, 1 for i, label in enumerate(labels): hemi_ind = hemi_to_ind[label['hemi']] - src_sel = np.intersect1d(fwd['src'][hemi_ind]['vertno'], + src_sel = np.intersect1d(src[hemi_ind]['vertno'], label['vertices']) if value_fun is not None: idx_sel = np.searchsorted(label['vertices'], src_sel) @@ -179,172 +181,3 @@ def generate_stc(fwd, labels, stc_data, tmin, tstep, value_fun=None): stc = _make_stc(stc_data, tmin, tstep, vertno) return stc - - -def mesh_dist(tris, vert): - """ Generate an adjacency matrix where the entries are the distances - between neighboring vertices - - Parameters: - ----------- - tris : array (n_tris x 3) - Mesh triangulation - vert : array (n_vert x 3) - Vertex locations - - Returns: - -------- - dist_matrix : scipy.sparse.csr_matrix - Sparse matrix with distances between adjacent vertices - """ - edges = mesh_edges(tris).tocoo() - - # Euclidean distances between neighboring vertices - dist = np.sqrt(np.sum((vert[edges.row, :] - vert[edges.col, :]) ** 2, - axis=1)) - - dist_matrix = csr_matrix((dist, (edges.row, edges.col)), shape=edges.shape) - - return dist_matrix - - -def _verts_within_dist(graph, source, max_dist): - """ Find all vertices wihin a maximum geodesic distance from source - - Parameters: - ----------- - graph : scipy.sparse.csr_matrix - Sparse matrix with distances between adjacent vertices - source : int - Source vertex - max_dist: float - Maximum geodesic distance - - Returns: - -------- - verts : array - Vertices within max_dist - dist : array - Distances from source vertex - """ - dist_map = {} - dist_map[source] = 0 - verts_added_last = [source] - # add neighbors until no more neighbors within max_dist can be found - while len(verts_added_last) > 0: - verts_added = [] - for i in verts_added_last: - v_dist = dist_map[i] - row = graph[i, :] - neighbor_vert = row.indices - neighbor_dist = row.data - for j, d in zip(neighbor_vert, neighbor_dist): - n_dist = v_dist + d - if j in dist_map: - if n_dist < dist_map[j]: - dist_map[j] = n_dist - else: - if n_dist <= max_dist: - dist_map[j] = n_dist - # we found a new vertex within max_dist - verts_added.append(j) - verts_added_last = verts_added - - verts = np.sort(np.array(dist_map.keys(), dtype=np.int)) - dist = np.array([dist_map[v] for v in verts]) - - return verts, dist - - -def circular_source_labels(subject, seeds, extents, hemis, subjects_dir=None): - """ Generate circular labels in source space - - This function generates a number of labels in source space by growing - regions starting from the vertices defined in "seeds". For each seed, a - label is generated containing all vertices within a maximum geodesic - distance on the white matter surface from the seed. - - Note: "extents" and "hemis" can either be arrays with the same length as - seeds, which allows using a different extent and hemisphere for each - label, or integers, in which case the same extent and hemisphere is - used for each label. - - Parameters: - ---------- - subject : string - Name of the subject as in SUBJECTS_DIR - seeds : array or int - Seed vertex numbers - extents : array or float - Extents (radius in mm) of the labels - hemis : array or int - Hemispheres to use for the labels (0: left, 1: right) - subjects_dir : string - Path to SUBJECTS_DIR if not set in the environment - - Returns: - -------- - labels : list - List with lables. Each label is a dictionary with keys: - comment Comment with seed vertex and extent - hemi Hemisphere of the label ('lh', or 'rh') - vertices Vertex indices (0 based) - pos Locations in millimeters - values Distances in millimeters from seed - """ - if subjects_dir is None: - if 'SUBJECTS_DIR' in os.environ: - subjects_dir = os.environ['SUBJECTS_DIR'] - else: - raise ValueError('SUBJECTS_DIR environment variable not set') - - # make sure the inputs are arrays - seeds = np.atleast_1d(seeds) - extents = np.atleast_1d(extents) - hemis = np.atleast_1d(hemis) - - n_seeds = len(seeds) - - if len(extents) != 1 and len(extents) != n_seeds: - raise ValueError('The extents parameter has to be of length 1 or ' - 'len(seeds)') - - if len(hemis) != 1 and len(hemis) != n_seeds: - raise ValueError('The hemis parameter has to be of length 1 or ' - 'len(seeds)') - - # make the arrays the same length as seeds - if len(extents) == 1: - extents = np.tile(extents, n_seeds) - - if len(hemis) == 1: - hemis = np.tile(hemis, n_seeds) - - hemis = ['lh' if h == 0 else 'rh' for h in hemis] - - # load the surfaces and create the distance graphs - tris, vert, dist = {}, {}, {} - for hemi in set(hemis): - surf_fname = os.path.join(subjects_dir, subject, 'surf', - hemi + '.white') - vert[hemi], tris[hemi] = read_surface(surf_fname) - dist[hemi] = mesh_dist(tris[hemi], vert[hemi]) - - # create the patches - labels = [] - for seed, extent, hemi in zip(seeds, extents, hemis): - label_verts, label_dist = _verts_within_dist(dist[hemi], seed, extent) - - # create a label - label = dict() - label['comment'] = 'Circular label: seed=%d, extent=%0.1fmm' %\ - (seed, extent) - label['vertices'] = label_verts - label['pos'] = vert[hemi][label_verts] - label['values'] = label_dist - label['hemi'] = hemi - - labels.append(label) - - return labels - diff --git a/mne/source_estimate.py b/mne/source_estimate.py index cccb049..4bf741f 100644 --- a/mne/source_estimate.py +++ b/mne/source_estimate.py @@ -9,6 +9,7 @@ import copy from math import ceil import numpy as np from scipy import sparse +from scipy.sparse import csr_matrix from .parallel import parallel_func @@ -552,6 +553,35 @@ def mesh_edges(tris): return edges +def mesh_dist(tris, vert): + """Compute adjacency matrix weighted by distances + + It generates an adjacency matrix where the entries are the distances + between neighboring vertices. + + Parameters: + ----------- + tris : array (n_tris x 3) + Mesh triangulation + vert : array (n_vert x 3) + Vertex locations + + Returns: + -------- + dist_matrix : scipy.sparse.csr_matrix + Sparse matrix with distances between adjacent vertices + """ + edges = mesh_edges(tris).tocoo() + + # Euclidean distances between neighboring vertices + dist = np.sqrt(np.sum((vert[edges.row, :] - vert[edges.col, :]) ** 2, + axis=1)) + + dist_matrix = csr_matrix((dist, (edges.row, edges.col)), shape=edges.shape) + + return dist_matrix + + def _morph_buffer(data, idx_use, e, smooth, n_vertices, nearest, maps): n_iter = 100 # max nb of smoothing iterations for k in range(n_iter): -- 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
