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 48d2bfa82a437c806ddd01634a6993e54bb8857f Merge: a075c40 8a3accd Author: Martin Luessi <[email protected]> Date: Mon Jul 16 14:12:25 2012 -0400 merged branch mluessi/simulation mne/__init__.py | 1 + mne/simulation/__init__.py | 3 +- mne/simulation/source.py | 177 +++++++++++++++++++++++++++++++++++- mne/simulation/tests/__init__.py | 0 mne/simulation/tests/test_source.py | 18 ++++ 5 files changed, 196 insertions(+), 3 deletions(-) diff --cc mne/simulation/__init__.py index 6bb7dde,0648a66..e0d887b --- a/mne/simulation/__init__.py +++ b/mne/simulation/__init__.py @@@ -1,5 -1,2 +1,6 @@@ +"""Data simulation code +""" + +from .source import select_source_in_label, generate_sparse_stc, \ - generate_evoked ++ generate_evoked, circular_source_labels + -from source import circular_source_labels diff --cc mne/simulation/source.py index b314d7b,799aae0..de416b4 --- a/mne/simulation/source.py +++ b/mne/simulation/source.py @@@ -1,196 -1,180 +1,369 @@@ -# Authors: Martin Luessi <[email protected]> +# Authors: Alexandre Gramfort <[email protected]> ++# Martin Luessi <[email protected]> +# Matti Hamalainen <[email protected]> # # License: BSD (3-clause) - + import os ++import copy import numpy as np +from scipy import signal - - import copy + from scipy.sparse import csr_matrix +from ..fiff.pick import pick_channels_cov +from ..minimum_norm.inverse import _make_stc +from ..utils import check_random_state +from ..forward import apply_forward + from ..surface import read_surface + from ..source_estimate import mesh_edges +def generate_evoked(fwd, stc, evoked, cov, snr=3, tmin=None, tmax=None, + fir_filter=None, random_state=None): + """Generate noisy evoked data + + Parameters + ---------- + fwd : dict + a forward solution + stc : SourceEstimate object + The source time courses + evoked : Evoked object + An instance of evoked used as template + cov : Covariance object + The noise covariance + snr : float + signal to noise ratio in dB. It corresponds to + 10 * log10( var(signal) / var(noise) ) + tmin : float | None + start of time interval to estimate SNR. If None first time point + is used. + tmax : float + start of time interval to estimate SNR. If None last time point + is used. + fir_filter : None | array + FIR filter coefficients e.g. [1, -1, 0.2] + random_state : None | int | np.random.RandomState + To specify the random generator state. + + Returns + ------- + evoked : Evoked object + The simulated evoked data + """ + evoked = apply_forward(fwd, stc, evoked) + noise = generate_noise_evoked(evoked, cov, fir_filter) + evoked_noise = add_noise_evoked(evoked, noise, snr, tmin=tmin, tmax=tmax) + return evoked_noise + + +def generate_noise_evoked(evoked, noise_cov, fir_filter=None, + random_state=None): + """Creates noise as a multivariate Gaussian + + The spatial covariance of the noise is given from the cov matrix. + + Parameters + ---------- + evoked : evoked object + an instance of evoked used as template + cov : Covariance object + The noise covariance + fir_filter : None | array + FIR filter coefficients + random_state : None | int | np.random.RandomState + To specify the random generator state. + + Returns + ------- + noise : evoked object + an instance of evoked + """ + noise = copy.deepcopy(evoked) + noise_cov = pick_channels_cov(noise_cov, include=noise.info['ch_names']) + rng = check_random_state(random_state) + n_channels = np.zeros(noise.info['nchan']) + n_samples = evoked.data.shape[1] + noise.data = rng.multivariate_normal(n_channels, noise_cov.data, + n_samples).T + if fir_filter is not None: + noise.data = signal.lfilter([1], fir_filter, noise.data, axis=-1) + return noise + + +def add_noise_evoked(evoked, noise, snr, tmin=None, tmax=None): + """Adds noise to evoked object with specified SNR. + + SNR is computed in the interval from tmin to tmax. + + Parameters + ---------- + evoked : Evoked object + An instance of evoked with signal + noise : Evoked object + An instance of evoked with noise + snr : float + signal to noise ratio in dB. It corresponds to + 10 * log10( var(signal) / var(noise) ) + tmin : float + start time before event + tmax : float + end time after event + + Returns + ------- + evoked_noise : Evoked object + An instance of evoked corrupted by noise + """ + evoked = copy.deepcopy(evoked) + times = evoked.times + if tmin is None: + tmin = np.min(times) + if tmax is None: + tmax = np.max(times) + tmask = (times >= tmin) & (times <= tmax) + tmp = np.mean((evoked.data[:, tmask] ** 2).ravel()) / \ + np.mean((noise.data ** 2).ravel()) + tmp = 10 * np.log10(tmp) + noise.data = 10 ** ((tmp - float(snr)) / 20) * noise.data + evoked.data += noise.data + return evoked + + +def select_source_in_label(fwd, label, random_state=None): + """Select source positions using a label + + Parameters + ---------- + fwd : dict + a forward solution + label : dict + the label (read with mne.read_label) + random_state : None | int | np.random.RandomState + To specify the random generator state. + + Returns + ------- + lh_vertno : list + selected source coefficients on the left hemisphere + rh_vertno : list + selected source coefficients on the right hemisphere + """ + lh_vertno = list() + rh_vertno = list() + + rng = check_random_state(random_state) + + if label['hemi'] == 'lh': + src_sel_lh = np.intersect1d(fwd['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']) + 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): + """Generate sources time courses from waveforms and labels + + Parameters + ---------- + fwd : dict + a forward solution + labels : list of dict + The labels + stc_data : array + The waveforms + tmin : float + The beginning of the timeseries + tstep : float + The time step (1 / sampling frequency) + random_state : None | int | np.random.RandomState + To specify the random generator state. + + Returns + ------- + stc : SourceEstimate + The generated source time courses. + """ + rng = check_random_state(random_state) + vertno = [[], []] + for label in labels: + lh_vertno, rh_vertno = select_source_in_label(fwd, label, rng) + vertno[0] += lh_vertno + vertno[1] += rh_vertno + vertno = map(np.array, vertno) + 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 + -- 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
