[med-svn] [python-mne] 210/376: ENH : adding the possibility to morph source estimates FIX : some fix to read sparse data in fif FIX : fix failing stats test
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:22:39 UTC 2015
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 <alexandre.gramfort at inria.fr>
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 <gramfort at nmr.mgh.harvard.edu>
+#
+# 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 <gramfort at nmr.mgh.harvard.edu>
+#
+# 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
More information about the debian-med-commit
mailing list