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 <[email protected]> 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 _______________________________________________ debian-med-commit mailing list [email protected] http://lists.alioth.debian.org/cgi-bin/mailman/listinfo/debian-med-commit
