[med-svn] [python-mne] 270/376: ENH : adding support for quad surf files API : no need to pass source space to mne.morph_data
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:23:05 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 b8eefc4482e46ca1fb73f63a6880cdeb9484caef
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Thu May 26 10:59:17 2011 -0400
ENH : adding support for quad surf files
API : no need to pass source space to mne.morph_data
---
examples/inverse/plot_morph_data.py | 6 ++--
mne/source_estimate.py | 25 +++++++++------
mne/surface.py | 62 +++++++++++++++++++++++++++----------
mne/tests/test_source_estimate.py | 11 +------
4 files changed, 64 insertions(+), 40 deletions(-)
diff --git a/examples/inverse/plot_morph_data.py b/examples/inverse/plot_morph_data.py
index b96f7e6..63f9f67 100644
--- a/examples/inverse/plot_morph_data.py
+++ b/examples/inverse/plot_morph_data.py
@@ -25,10 +25,10 @@ subject_to = 'morph'
fname = data_path + '/MEG/sample/sample_audvis-meg'
src_fname = data_path + '/MEG/sample/sample_audvis-meg-oct-6-fwd.fif'
+# Read input stc file
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)
+# Morph
+stc_to = mne.morph_data(subject_from, subject_to, stc_from)
stc_to.save('%s_audvis-meg' % subject_to)
diff --git a/mne/source_estimate.py b/mne/source_estimate.py
index 397c1d8..b263097 100644
--- a/mne/source_estimate.py
+++ b/mne/source_estimate.py
@@ -273,7 +273,8 @@ def mesh_edges(tris):
edges = edges + edges.T
return edges
-def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
+
+def morph_data(subject_from, subject_to, stc_from, grade=5,
subjects_dir=None):
"""Morph a source estimate from one subject to another
@@ -285,8 +286,6 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
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
@@ -307,6 +306,17 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
else:
raise ValueError('SUBJECTS_DIR environment variable not set')
+ tris = list()
+ surf_path_from = os.path.join(subjects_dir, subject_from, 'surf')
+ tris.append(read_surface(os.path.join(surf_path_from, 'lh.sphere.reg'))[1])
+ tris.append(read_surface(os.path.join(surf_path_from, 'rh.sphere.reg'))[1])
+ vertices = [stc_from.lh_vertno, stc_from.rh_vertno]
+
+ 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]
+
maps = read_morph_map(subject_from, subject_to)
lh_data = stc_from.data[:len(stc_from.lh_vertno)]
@@ -315,11 +325,11 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
dmap = [None, None]
for hemi in [0, 1]:
- e = mesh_edges(src_from[hemi]['tris'])
+ e = mesh_edges(tris[hemi])
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
+ idx_use = vertices[hemi]
n_iter = 100 # max nb of smoothing iterations
for k in range(n_iter):
e_use = e[:, idx_use]
@@ -345,11 +355,6 @@ def morph_data(subject_from, subject_to, src_from, stc_from, grade=5,
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]
diff --git a/mne/surface.py b/mne/surface.py
index 085daa5..13522c5 100644
--- a/mne/surface.py
+++ b/mne/surface.py
@@ -210,22 +210,28 @@ def _complete_surface_info(this):
###############################################################################
# Handle freesurfer
-def fread3(fobj):
+def _fread3(fobj):
"""Docstring"""
b1, b2, b3 = np.fromfile(fobj, ">u1", 3)
return (b1 << 16) + (b2 << 8) + b3
+def _fread3_many(fobj, n):
+ """Read 3-byte ints from an open binary file object."""
+ b1, b2, b3 = np.fromfile(fobj, ">u1", 3*n).reshape(-1, 3).astype(np.int).T
+ 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)
+ magic = _fread3(fobj)
if magic == 16777215:
vnum = np.fromfile(fobj, ">i4", 3)[0]
curv = np.fromfile(fobj, ">f4", vnum)
else:
vnum = magic
- fread3(fobj)
+ _fread3(fobj)
curv = np.fromfile(fobj, ">i2", vnum) / 100
bin_curv = 1 - np.array(curv != 0, np.int)
return bin_curv
@@ -234,18 +240,40 @@ def read_curvature(filepath):
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:
+ magic = _fread3(fobj)
+ if magic == 16777215: # Quad file
+ nvert = _fread3(fobj)
+ nquad = _fread3(fobj)
+ coords = np.fromfile(fobj, ">i2", nvert*3).astype(np.float)
+ coords = coords.reshape(-1, 3) / 100.0
+ quads = _fread3_many(fobj, nquad*4)
+ quads = quads.reshape(nquad, 4)
+ #
+ # Face splitting follows
+ #
+ faces = np.zeros((2*nquad, 3), dtype=np.int)
+ nface = 0
+ for quad in quads:
+ if (quad[0] % 2) == 0:
+ faces[nface] = quad[0], quad[1], quad[3]
+ nface += 1
+ faces[nface] = quad[2], quad[3], quad[1]
+ nface += 1
+ else:
+ faces[nface] = quad[0], quad[1], quad[2]
+ nface += 1
+ faces[nface] = quad[0], quad[2], quad[3]
+ nface += 1
+
+ elif magic == 16777214: # Triangle file
+ create_stamp = fobj.readline()
+ _ = fobj.readline()
+ vnum = np.fromfile(fobj, ">i4", 1)[0]
+ fnum = np.fromfile(fobj, ">i4", 1)[0]
+ coords = np.fromfile(fobj, ">f4", vnum * 3).reshape(vnum, 3)
+ faces = np.fromfile(fobj, ">i4", fnum * 3).reshape(fnum, 3)
+ else:
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)
- vertex_coords = vertex_coords.astype(np.float) # XXX : due to mayavi bug
- return vertex_coords, faces
+
+ coords = coords.astype(np.float) # XXX: due to mayavi bug on mac 32bits
+ return coords, faces
diff --git a/mne/tests/test_source_estimate.py b/mne/tests/test_source_estimate.py
index 60cb087..586f04b 100644
--- a/mne/tests/test_source_estimate.py
+++ b/mne/tests/test_source_estimate.py
@@ -2,7 +2,6 @@ import os.path as op
import numpy as np
from numpy.testing import assert_array_almost_equal
-from scipy import linalg
import mne
from mne.datasets import sample
@@ -31,19 +30,11 @@ def test_morph_data():
"""Test morphing of data
"""
import mne
- from mne.datasets import sample
-
subject_from = 'sample'
subject_to = 'morph'
-
fname = op.join(data_path, 'MEG', 'sample', 'sample_audvis-meg')
- src_fname = op.join(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, 3)
+ stc_to = mne.morph_data(subject_from, subject_to, stc_from, 3)
stc_to.save('%s_audvis-meg' % subject_to)
--
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