[med-svn] [python-mne] 20/376: Reading BEM surfaces
Yaroslav Halchenko
debian at onerussian.com
Fri Nov 27 17:21:57 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 84354567a78870dccfc61c709c4fcadde5adcef8
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date: Thu Dec 30 14:52:37 2010 -0500
Reading BEM surfaces
---
examples/read_bem_surfaces.py | 20 ++++
fiff/__init__.py | 1 +
fiff/bem_surfaces.py | 206 ++++++++++++++++++++++++++++++++++++++++++
3 files changed, 227 insertions(+)
diff --git a/examples/read_bem_surfaces.py b/examples/read_bem_surfaces.py
new file mode 100644
index 0000000..d6567b0
--- /dev/null
+++ b/examples/read_bem_surfaces.py
@@ -0,0 +1,20 @@
+"""Reading BEM surfaces
+"""
+print __doc__
+
+import fiff
+
+fname = 'MNE-sample-data/subjects/sample/bem/sample-5120-bem-sol.fif'
+
+surf = fiff.read_bem_surfaces(fname)
+
+print "Number of surfaces : %d" % len(surf)
+
+###############################################################################
+# Show result
+
+# 3D source space
+points = surf[0]['rr']
+faces = surf[0]['tris']
+from enthought.mayavi import mlab
+mlab.triangular_mesh(points[:,0], points[:,1], points[:,2], faces)
diff --git a/fiff/__init__.py b/fiff/__init__.py
index 10be9cc..bfc53a9 100644
--- a/fiff/__init__.py
+++ b/fiff/__init__.py
@@ -8,4 +8,5 @@ from .raw import setup_read_raw, read_raw_segment, read_raw_segment_times
from .event import read_events
from .forward import read_forward_solution
from .stc import read_stc
+from .bem_surfaces import read_bem_surfaces
diff --git a/fiff/bem_surfaces.py b/fiff/bem_surfaces.py
new file mode 100644
index 0000000..2e0fe4c
--- /dev/null
+++ b/fiff/bem_surfaces.py
@@ -0,0 +1,206 @@
+import numpy as np
+
+from .constants import FIFF
+from .open import fiff_open
+from .tree import dir_tree_find
+from .tag import find_tag
+from scipy import linalg
+
+#
+# These fiff definitions are not needed elsewhere
+#
+FIFFB_BEM = 310 # BEM data
+FIFFB_BEM_SURF = 311 # One of the surfaces
+#
+FIFF_BEM_SURF_ID = 3101 # int surface number
+FIFF_BEM_SURF_NAME = 3102 # string surface name
+FIFF_BEM_SURF_NNODE = 3103 # int # of nodes on a surface
+FIFF_BEM_SURF_NTRI = 3104 # int # number of triangles on a surface
+FIFF_BEM_SURF_NODES = 3105 # float surface nodes (nnode,3)
+FIFF_BEM_SURF_TRIANGLES = 3106 # int surface triangles (ntri,3)
+FIFF_BEM_SURF_NORMALS = 3107 # float surface node normal unit vectors (nnode,3)
+FIFF_BEM_COORD_FRAME = 3112 # The coordinate frame of the mode
+FIFF_BEM_SIGMA = 3113 # Conductivity of a compartment
+
+
+def read_bem_surfaces(fname, add_geom=False):
+ """
+ #
+ # [surf] = mne_read_bem_surfaces(fname, add_geom)
+ #
+ # Reads source spaces from a fif file
+ #
+ # fname - The name of the file or an open file id
+ # add_geom - Add geometry information to the surfaces
+ #
+ """
+ #
+ # Default coordinate frame
+ #
+ coord_frame = FIFF.FIFFV_COORD_MRI
+ #
+ # Open the file, create directory
+ #
+ fid, tree, _ = fiff_open(fname)
+ #
+ # Find BEM
+ #
+ bem = dir_tree_find(tree, FIFFB_BEM)
+ if bem is None:
+ fid.close()
+ raise ValueError, 'BEM data not found'
+
+ bem = bem[0]
+ #
+ # Locate all surfaces
+ #
+ bemsurf = dir_tree_find(bem, FIFFB_BEM_SURF)
+ if bemsurf is None:
+ fid.close()
+ raise ValueError, 'BEM surface data not found'
+
+ print '\t%d BEM surfaces found' % len(bemsurf)
+ #
+ # Coordinate frame possibly at the top level
+ #
+ tag = find_tag(fid, bem, FIFF_BEM_COORD_FRAME)
+ if tag is not None:
+ coord_frame = tag.data
+ #
+ # Read all surfaces
+ #
+ surf = []
+ for bsurf in bemsurf:
+ print '\tReading a surface...'
+ this = read_bem_surface(fid, bsurf, coord_frame)
+ print '[done]'
+ if add_geom:
+ complete_surface_info()
+ surf.append(this)
+
+ print '\t%d BEM surfaces read' % len(surf)
+
+ fid.close()
+
+ return surf
+
+
+def read_bem_surface(fid, this, def_coord_frame):
+ """
+ """
+ res = dict()
+ #
+ # 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:
+ res['id'] = tag.data
+
+ tag = find_tag(fid, this, FIFF_BEM_SIGMA)
+ if tag is None:
+ res['sigma'] = 1.0
+ else:
+ res['sigma'] = tag.data
+
+ tag = find_tag(fid, this, FIFF_BEM_SURF_NNODE)
+ if tag is None:
+ fid.close()
+ raise ValueError, 'Number of vertices not found'
+
+ res = dict()
+ res['np'] = tag.data
+
+ tag = find_tag(fid, this,FIFF_BEM_SURF_NTRI)
+ if tag is None:
+ fid.close()
+ raise ValueError, 'Number of triangles not found'
+ else:
+ res['ntri'] = tag.data
+
+ tag = find_tag(fid, this, FIFF.FIFF_MNE_COORD_FRAME)
+ if tag is None:
+ tag = find_tag(fid, this, FIFF_BEM_COORD_FRAME)
+ if tag is None:
+ res['coord_frame'] = def_coord_frame
+ else:
+ res['coord_frame'] = tag.data
+ else:
+ res['coord_frame'] = tag.data
+ #
+ # Vertices, normals, and triangles
+ #
+ tag = find_tag(fid, this, FIFF_BEM_SURF_NODES)
+ if tag is None:
+ fid.close()
+ raise ValueError, 'Vertex data not found'
+
+ res['rr'] = tag.data.astype(np.float) # XXX : double because of mayavi bug
+ if res['rr'].shape[0] != res['np']:
+ fid.close()
+ raise ValueError, 'Vertex information is incorrect'
+
+ tag = find_tag(fid, this, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS)
+ if tag is None:
+ res['nn'] = []
+ else:
+ res['nn'] = tag.data
+ if res['nn'].shape[0] != res['np']:
+ fid.close()
+ raise ValueError, 'Vertex normal information is incorrect'
+
+ tag = find_tag(fid, this, FIFF_BEM_SURF_TRIANGLES)
+ if tag is None:
+ fid.close()
+ raise ValueError, 'Triangulation not found'
+
+ res['tris'] = tag.data - 1 # index start at 0 in Python
+ if res['tris'].shape[0] != res['ntri']:
+ fid.close()
+ raise ValueError, 'Triangulation information is incorrect'
+
+ return res
+
+
+def complete_surface_info(this):
+ """ XXX : should be factorize with complete_source_space_info
+ """
+ #
+ # Main triangulation
+ #
+ print '\tCompleting triangulation info...'
+ print 'triangle normals...'
+ this.tri_area = np.zeros(this['ntri'])
+ r1 = this['rr'][this['tris'][:,0],:]
+ r2 = this['rr'][this['tris'][:,1],:]
+ r3 = this['rr'][this['tris'][:,2],:]
+ this.tri_cent = (r1+r2+r3) /3.0
+ this['tri_nn'] = np.cross((r2-r1), (r3-r1))
+ #
+ # Triangle normals and areas
+ #
+ for p in range(this['ntri']):
+ size = linalg.norm(this['tri_nn'][p,:])
+ this.tri_area[p] = size / 2.0
+ if size > 0.0:
+ this['tri_nn'][p,:] = this['tri_nn'][p,:] / size
+ #
+ # Accumulate the vertex normals
+ #
+ print 'vertex normals...'
+ this['nn'] = np.zeros((this['np'], 3))
+ for p in range(this['ntri']):
+ this['nn'][this['tris'][p,:],:] = this['nn'][this['tris'][p,:],:] \
+ + np.kron(np.ones((3, 1)), this['tri_nn'][p,:])
+ #
+ # Compute the lengths of the vertex normals and scale
+ #
+ print 'normalize...'
+ for p in range(this['np']):
+ size = linalg.norm(this['nn'][p,:])
+ if size > 0:
+ this['nn'][p,:] = this['nn'][p,:] / size
+
+ print '[done]\n'
+ return this
--
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