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 8502ef8b5cd4c842b3f5ebbf46cbc3bfbe00db16 Author: Alexandre Gramfort <[email protected]> Date: Fri Apr 22 11:49:34 2011 -0400 ENH : adding the possibility to morph source estimates FIX : some fix to read sparse data in fif FIX : fix failing stats test --- examples/plot_morph_data.py | 44 +++++++ mne/__init__.py | 2 +- mne/bem_surfaces.py | 2 +- mne/fiff/tag.py | 5 +- mne/source_space.py | 1 - mne/stats/cluster_level.py | 2 +- mne/stats/tests/test_cluster_level.py | 8 +- mne/stc.py | 218 ++++++++++++++++++++++++++++++++++ mne/surfer.py | 50 ++++++++ 9 files changed, 322 insertions(+), 10 deletions(-) diff --git a/examples/plot_morph_data.py b/examples/plot_morph_data.py new file mode 100644 index 0000000..bb1471f --- /dev/null +++ b/examples/plot_morph_data.py @@ -0,0 +1,44 @@ +""" +========================================================== +Morph source estimates from one subject to another subject +========================================================== + +A source estimate from a given subject 'sample' is morphed +to the anatomy of another subject 'morph'. The output +is a source estimate defined on the anatomy of 'morph' + +""" +# Author: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +print __doc__ + +import mne +from mne.datasets import sample + +data_path = sample.data_path('.') + +subject_from = 'sample' +subject_to = 'morph' + +fname = data_path + '/MEG/sample/sample_audvis-meg' +fname = data_path + '/MEG/sample/sample_audvis-meg' +src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif' + +stc_from = mne.SourceEstimate(fname) +src_from = mne.read_source_spaces(src_fname) + +stc_to = mne.morph_data(subject_from, subject_to, src_from, stc_from) + +stc_to.save('%s_audvis-meg' % subject_to) + +# View source activations +import pylab as pl +pl.plot(stc_from.times, stc_from.data.mean(axis=0), 'r', label='from') +pl.plot(stc_to.times, stc_to.data.mean(axis=0), 'b', label='to') +pl.xlabel('time (ms)') +pl.ylabel('Mean Source amplitude') +pl.legend() +pl.show() + diff --git a/mne/__init__.py b/mne/__init__.py index 289ef1d..58ec38a 100755 --- a/mne/__init__.py +++ b/mne/__init__.py @@ -4,7 +4,7 @@ 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 .stc import read_stc, write_stc, SourceEstimate +from .stc import read_stc, write_stc, SourceEstimate, morph_data from .bem_surfaces import read_bem_surfaces from .source_space import read_source_spaces from .epochs import Epochs diff --git a/mne/bem_surfaces.py b/mne/bem_surfaces.py index 2c62811..eb2c18a 100755 --- a/mne/bem_surfaces.py +++ b/mne/bem_surfaces.py @@ -101,6 +101,7 @@ def _read_bem_surface(fid, this, def_coord_frame): # Read all the interesting stuff # tag = find_tag(fid, this, FIFF_BEM_SURF_ID) + if tag is None: res['id'] = FIFF.FIFFV_BEM_SURF_ID_UNKNOWN else: @@ -117,7 +118,6 @@ def _read_bem_surface(fid, this, def_coord_frame): fid.close() raise ValueError('Number of vertices not found') - res = dict() res['np'] = tag.data tag = find_tag(fid, this, FIFF_BEM_SURF_NTRI) diff --git a/mne/fiff/tag.py b/mne/fiff/tag.py index e05d1ef..a2c156b 100755 --- a/mne/fiff/tag.py +++ b/mne/fiff/tag.py @@ -172,19 +172,20 @@ def read_tag(fid, pos=None): nrow = dims[1] ncol = dims[2] sparse_data = np.fromfile(fid, dtype='>f4', count=nnz) + shape = (dims[1], dims[2]) if matrix_coding == matrix_coding_CCS: # CCS sparse.csc_matrix() sparse_indices = np.fromfile(fid, dtype='>i4', count=nnz) sparse_ptrs = np.fromfile(fid, dtype='>i4', count=ncol + 1) tag.data = sparse.csc_matrix((sparse_data, sparse_indices, - sparse_ptrs), shape=dims) + sparse_ptrs), shape=shape) else: # RCS sparse_indices = np.fromfile(fid, dtype='>i4', count=nnz) sparse_ptrs = np.fromfile(fid, dtype='>i4', count=nrow + 1) tag.data = sparse.csr_matrix((sparse_data, sparse_indices, - sparse_ptrs), shape=dims) + sparse_ptrs), shape=shape) else: raise ValueError('Cannot handle other than dense or sparse ' 'matrices yet') diff --git a/mne/source_space.py b/mne/source_space.py index 48322a1..5beaa28 100755 --- a/mne/source_space.py +++ b/mne/source_space.py @@ -3,7 +3,6 @@ # # License: BSD (3-clause) -from math import sqrt import numpy as np from .fiff.constants import FIFF diff --git a/mne/stats/cluster_level.py b/mne/stats/cluster_level.py index c00ea1e..aa82aec 100755 --- a/mne/stats/cluster_level.py +++ b/mne/stats/cluster_level.py @@ -291,7 +291,7 @@ permutation_cluster_1samp_test.__test__ = False # # X, threshold, tail = condition1, 1.67, 1 # # X, threshold, tail = -condition1, -1.67, -1 # # # X, threshold, tail = condition1, 1.67, 0 -# # fs, clusters, cluster_p_values, histogram = permutation_cluster_t_test( +# # fs, clusters, cluster_p_values, histogram = permutation_cluster_1samp_test( # # condition1, n_permutations=500, tail=tail, # # threshold=threshold) # diff --git a/mne/stats/tests/test_cluster_level.py b/mne/stats/tests/test_cluster_level.py index 5c7d78a..a92130f 100755 --- a/mne/stats/tests/test_cluster_level.py +++ b/mne/stats/tests/test_cluster_level.py @@ -2,7 +2,7 @@ import numpy as np from numpy.testing import assert_equal, assert_array_equal from ..cluster_level import permutation_cluster_test, \ - permutation_cluster_t_test + permutation_cluster_1samp_test noiselevel = 20 @@ -41,15 +41,15 @@ def test_cluster_permutation_t_test(): """ my_condition1 = condition1[:,:,None] # to test 2D also - T_obs, clusters, cluster_p_values, hist = permutation_cluster_t_test( + T_obs, clusters, cluster_p_values, hist = permutation_cluster_1samp_test( my_condition1, n_permutations=500, tail=0) assert_equal(np.sum(cluster_p_values < 0.05), 1) - T_obs_pos, _, cluster_p_values_pos, _ = permutation_cluster_t_test( + T_obs_pos, _, cluster_p_values_pos, _ = permutation_cluster_1samp_test( my_condition1, n_permutations=500, tail=1, threshold=1.67) - T_obs_neg, _, cluster_p_values_neg, _ = permutation_cluster_t_test( + T_obs_neg, _, cluster_p_values_neg, _ = permutation_cluster_1samp_test( -my_condition1, n_permutations=500, tail=-1, threshold=-1.67) assert_array_equal(T_obs_pos, -T_obs_neg) diff --git a/mne/stc.py b/mne/stc.py index 13919d9..71a620a 100755 --- a/mne/stc.py +++ b/mne/stc.py @@ -3,6 +3,8 @@ # # License: BSD (3-clause) +import os +import copy import numpy as np @@ -144,3 +146,219 @@ class SourceEstimate(object): write_stc(fname + '-rh.stc', tmin=self.tmin, tstep=self.tstep, vertices=self.rh_vertno, data=rh_data) print '[done]' + + def __repr__(self): + s = "%d vertices" % (len(self.lh_vertno) + len(self.rh_vertno)) + s += ", tmin : %s (ms)" % (1e3 * self.tmin) + s += ", tmax : %s (ms)" % (1e3 * self.times[-1]) + s += ", tstep : %s (ms)" % (1e3 * self.tstep) + s += ", data size : %s x %s" % self.data.shape + return "SourceEstimate (%s)" % s + + +############################################################################### +# Morphing + +from .fiff.constants import FIFF +from .fiff.tag import find_tag +from .fiff.open import fiff_open +from .fiff.tree import dir_tree_find +from .bem_surfaces import read_bem_surfaces +from .surfer import read_surface + + +def read_morph_map(subject_from, subject_to, subjects_dir=None): + """Read morph map generated with mne_make_morph_maps + + Parameters + ---------- + subject_from : string + Name of the original subject as named in the SUBJECTS_DIR + subject_to : string + Name of the subject on which to morph as named in the SUBJECTS_DIR + subjects_dir : string + Path to SUBJECTS_DIR is not set in the environment + + Returns + ------- + maps : dict + The morph maps for the 2 hemisphere + """ + + 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') + + # Does the file exist + name = '%s/morph-maps/%s-%s-morph.fif' % (subjects_dir, subject_from, + subject_to) + if not os.path.exists(name): + name = '%s/morph-maps/%s-%s-morph.fif' % (subjects_dir, subject_to, + subject_from) + if not os.path.exists(name): + raise ValueError('The requested morph map does not exist') + + fid, tree, _ = fiff_open(name) + + # Locate all maps + maps = dir_tree_find(tree, FIFF.FIFFB_MNE_MORPH_MAP) + if len(maps) == 0: + fid.close() + raise ValueError('Morphing map data not found') + + # Find the correct ones + leftmap = None + rightmap = None + for m in maps: + tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP_FROM) + if tag.data == subject_from: + tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP_TO) + if tag.data == subject_to: + # Names match: which hemishere is this? + tag = find_tag(fid, m, FIFF.FIFF_MNE_HEMI) + if tag.data == FIFF.FIFFV_MNE_SURF_LEFT_HEMI: + tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP) + leftmap = tag.data + print '\tLeft-hemisphere map read.' + elif tag.data == FIFF.FIFFV_MNE_SURF_RIGHT_HEMI: + tag = find_tag(fid, m, FIFF.FIFF_MNE_MORPH_MAP) + rightmap = tag.data + print '\tRight-hemisphere map read.' + + fid.close() + if leftmap is None: + raise ValueError('Left hemisphere map not found in %s' % name) + + if rightmap is None: + raise ValueError('Left hemisphere map not found in %s' % name) + + return leftmap, rightmap + + +def mesh_edges(tris): + """Returns sparse matrix with edges as an adjacency matrix + + Parameters + ---------- + tris : array of shape [n_triangles x 3] + The triangles + + Returns + ------- + edges : sparse matrix + The adjacency matrix + """ + from scipy import sparse + npoints = np.max(tris) + 1 + ntris = len(tris) + a, b, c = tris.T + edges = sparse.coo_matrix((np.ones(ntris), (a, b)), + shape=(npoints, npoints)) + edges = edges + sparse.coo_matrix((np.ones(ntris), (b, c)), + shape=(npoints, npoints)) + edges = edges + sparse.coo_matrix((np.ones(ntris), (c, a)), + shape=(npoints, npoints)) + edges = edges.tocsr() + edges = edges + edges.T + return edges + + +def morph_data(subject_from, subject_to, src_from, stc_from, grade=5, + subjects_dir=None): + """Morph a source estimate from one subject to another + + The functions requires to set MNE_ROOT and SUBJECTS_DIR variables. + + Parameters + ---------- + subject_from : string + Name of the original subject as named in the SUBJECTS_DIR + subject_to : string + Name of the subject on which to morph as named in the SUBJECTS_DIR + src_from : dict + Source space for original "from" subject + stc_from : SourceEstimate + Source estimates for subject "from" to morph + grade : int + Resolution of the icosahedral mesh (typically 5) + subjects_dir : string + Path to SUBJECTS_DIR is not set in the environment + + Returns + ------- + stc_to : SourceEstimate + Source estimate for the destination subject. + """ + from scipy import sparse + + 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') + + maps = read_morph_map(subject_from, subject_to) + + lh_data = stc_from.data[:len(stc_from.lh_vertno)] + rh_data = stc_from.data[-len(stc_from.rh_vertno):] + data = [lh_data, rh_data] + dmap = [None, None] + + for hemi in [0, 1]: + e = mesh_edges(src_from[hemi]['tris']) + e.data[e.data == 2] = 1 + n_vertices = e.shape[0] + e = e + sparse.eye(n_vertices, n_vertices) + idx_use = np.where(src_from[hemi]['inuse'])[0] # XXX + n_iter = 100 # max nb of smoothing iterations + for k in range(n_iter): + data1 = e[:, idx_use] * np.ones(len(idx_use)) + data[hemi] = e[:, idx_use] * data[hemi] + idx_use = np.where(data1)[0] + if (k == (n_iter - 1)) or (len(idx_use) >= n_vertices): + data[hemi][idx_use, :] /= data1[idx_use][:, None] + break + else: + data[hemi] = data[hemi][idx_use, :] / data1[idx_use][:, None] + + print '\t%d smooth iterations done.' % k + dmap[hemi] = maps[hemi] * data[hemi] + + ico_file_name = os.path.join(os.environ['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 + + sphere = os.path.join(subjects_dir, subject_to, 'surf', 'lh.sphere.reg') + lhs = read_surface(sphere)[0] + sphere = os.path.join(subjects_dir, subject_to, 'surf', 'rh.sphere.reg') + rhs = read_surface(sphere)[0] + + nearest = np.zeros((2, ico['np']), dtype=np.int) + + lhs /= np.sqrt(np.sum(lhs ** 2, axis=1))[:, None] + rhs /= np.sqrt(np.sum(rhs ** 2, axis=1))[:, None] + + rr = ico['rr'] + dr = 16 + for k in range(0, len(rr), dr): + dots = np.dot(rr[k:k + dr], lhs.T) + nearest[0, k:k + dr] = np.argmax(dots, axis=1) + dots = np.dot(rr[k:k + dr], rhs.T) + nearest[1, k:k + dr] = np.argmax(dots, axis=1) + + stc_to = copy.copy(stc_from) + stc_to.lh_vertno = nearest[0] + stc_to.rh_vertno = nearest[1] + stc_to.data = np.r_[dmap[0][nearest[0], :], dmap[1][nearest[1], :]] + + print '[done]' + + return stc_to diff --git a/mne/surfer.py b/mne/surfer.py new file mode 100644 index 0000000..80727cf --- /dev/null +++ b/mne/surfer.py @@ -0,0 +1,50 @@ +"""Set of tools to interact with Freesurfer data +""" + +# Authors: Alexandre Gramfort <[email protected]> +# +# License: BSD (3-clause) + +import numpy as np + + +def fread3(fobj): + """Docstring""" + b1 = np.fromfile(fobj, ">u1", 1)[0] + b2 = np.fromfile(fobj, ">u1", 1)[0] + b3 = np.fromfile(fobj, ">u1", 1)[0] + return (b1 << 16) + (b2 << 8) + b3 + + +def read_curvature(filepath): + """Load in curavature values from the ?h.curv file.""" + with open(filepath, "rb") as fobj: + magic = fread3(fobj) + if magic == 16777215: + vnum = np.fromfile(fobj, ">i4", 3)[0] + curv = np.fromfile(fobj, ">f4", vnum) + else: + vnum = magic + fread3(fobj) + curv = np.fromfile(fobj, ">i2", vnum) / 100 + bin_curv = 1 - np.array(curv != 0, np.int) + return bin_curv + + +def read_surface(filepath): + """Load in a Freesurfer surface mesh in triangular format.""" + with open(filepath, "rb") as fobj: + magic = fread3(fobj) + if magic == 16777215: + raise NotImplementedError("Quadrangle surface format reading not " + "implemented") + elif magic != 16777214: + raise ValueError("File does not appear to be a Freesurfer surface") + create_stamp = fobj.readline() + blankline = fobj.readline() + del blankline + vnum = np.fromfile(fobj, ">i4", 1)[0] + fnum = np.fromfile(fobj, ">i4", 1)[0] + vertex_coords = np.fromfile(fobj, ">f4", vnum*3).reshape(vnum, 3) + faces = np.fromfile(fobj, ">i4", fnum*3).reshape(fnum, 3) + return vertex_coords, faces -- 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
