[med-svn] [python-mne] 92/353: ENH : add support to write BEM surfaces in fif file

Yaroslav Halchenko debian at onerussian.com
Fri Nov 27 17:24:36 UTC 2015


This is an automated email from the git hooks/post-receive script.

yoh pushed a commit to tag 0.4
in repository python-mne.

commit 37d9fb4313d3fd5c237f70feb6023dcd2de50b08
Author: Alexandre Gramfort <alexandre.gramfort at inria.fr>
Date:   Wed Feb 29 17:37:16 2012 +0100

    ENH : add support to write BEM surfaces in fif file
---
 doc/source/whats_new.rst  |  5 ++++
 mne/__init__.py           |  2 +-
 mne/fiff/write.py         | 22 +++++++++++++++++
 mne/surface.py            | 61 ++++++++++++++++++++++++++++++++++++++++-------
 mne/tests/test_surface.py | 12 +++++++---
 5 files changed, 90 insertions(+), 12 deletions(-)

diff --git a/doc/source/whats_new.rst b/doc/source/whats_new.rst
index 0a67d14..277c463 100644
--- a/doc/source/whats_new.rst
+++ b/doc/source/whats_new.rst
@@ -19,6 +19,11 @@ Changelog
 
    - Support for computing sensor space data from a source estimate using an MNE forward solution by `Martin Luessi`_.
 
+   - Support of arithmetic of Covariance by `Alex Gramfort`_.
+
+   - Write BEM surfaces in Python  by `Alex Gramfort`_.
+
+
 Version 0.2
 -----------
 
diff --git a/mne/__init__.py b/mne/__init__.py
index 6dbb5b1..89f389a 100644
--- a/mne/__init__.py
+++ b/mne/__init__.py
@@ -10,7 +10,7 @@ from .source_estimate import read_stc, write_stc, read_w, write_w, \
                              spatio_temporal_src_connectivity, \
                              spatio_temporal_tris_connectivity, \
                              save_stc_as_volume
-from .surface import read_bem_surfaces, read_surface
+from .surface import read_bem_surfaces, read_surface, write_bem_surface
 from .source_space import read_source_spaces
 from .epochs import Epochs
 from .label import label_time_courses, read_label, label_sign_flip, \
diff --git a/mne/fiff/write.py b/mne/fiff/write.py
index 6349ed5..7e7c983 100644
--- a/mne/fiff/write.py
+++ b/mne/fiff/write.py
@@ -87,6 +87,28 @@ def write_float_matrix(fid, kind, mat):
     fid.write(np.array(dims, dtype='>i4').tostring())
 
 
+def write_int_matrix(fid, kind, mat):
+    """Writes a single-precision floating-point matrix tag"""
+    FIFFT_INT = 3
+    FIFFT_MATRIX = 1 << 30
+    FIFFT_MATRIX_INT = FIFFT_INT | FIFFT_MATRIX
+    FIFFV_NEXT_SEQ = 0
+
+    data_size = 4 * mat.size + 4 * 3
+
+    fid.write(np.array(kind, dtype='>i4').tostring())
+    fid.write(np.array(FIFFT_MATRIX_INT, dtype='>i4').tostring())
+    fid.write(np.array(data_size, dtype='>i4').tostring())
+    fid.write(np.array(FIFFV_NEXT_SEQ, dtype='>i4').tostring())
+    fid.write(np.array(mat, dtype='>i4').tostring())
+
+    dims = np.empty(3, dtype=np.int32)
+    dims[0] = mat.shape[1]
+    dims[1] = mat.shape[0]
+    dims[2] = 2
+    fid.write(np.array(dims, dtype='>i4').tostring())
+
+
 def write_id(fid, kind, id_=None):
     """Writes fiff id"""
 
diff --git a/mne/surface.py b/mne/surface.py
index cf948b1..b99b6c7 100644
--- a/mne/surface.py
+++ b/mne/surface.py
@@ -9,6 +9,9 @@ from .fiff.constants import FIFF
 from .fiff.open import fiff_open
 from .fiff.tree import dir_tree_find
 from .fiff.tag import find_tag
+from .fiff.write import write_int, write_float, write_float_matrix, \
+                        write_int_matrix, start_file, end_block, \
+                        start_block, end_file
 
 #
 #   These fiff definitions are not needed elsewhere
@@ -105,27 +108,27 @@ def _read_bem_surface(fid, this, def_coord_frame):
     if tag is None:
         res['id'] = FIFF.FIFFV_BEM_SURF_ID_UNKNOWN
     else:
-        res['id'] = tag.data
+        res['id'] = int(tag.data)
 
     tag = find_tag(fid, this, FIFF_BEM_SIGMA)
     if tag is None:
         res['sigma'] = 1.0
     else:
-        res['sigma'] = tag.data
+        res['sigma'] = float(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['np'] = tag.data
+    res['np'] = int(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
+        res['ntri'] = int(tag.data)
 
     tag = find_tag(fid, this, FIFF.FIFF_MNE_COORD_FRAME)
     if tag is None:
@@ -218,7 +221,8 @@ def _fread3(fobj):
 
 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
+    b1, b2, b3 = np.fromfile(fobj, ">u1",
+                             3 * n).reshape(-1, 3).astype(np.int).T
     return (b1 << 16) + (b2 << 8) + b3
 
 
@@ -244,14 +248,14 @@ def read_surface(filepath):
         if magic == 16777215:  # Quad file
             nvert = _fread3(fobj)
             nquad = _fread3(fobj)
-            coords = np.fromfile(fobj, ">i2", nvert*3).astype(np.float)
+            coords = np.fromfile(fobj, ">i2", nvert * 3).astype(np.float)
             coords = coords.reshape(-1, 3) / 100.0
-            quads = _fread3_many(fobj, nquad*4)
+            quads = _fread3_many(fobj, nquad * 4)
             quads = quads.reshape(nquad, 4)
             #
             #   Face splitting follows
             #
-            faces = np.zeros((2*nquad, 3), dtype=np.int)
+            faces = np.zeros((2 * nquad, 3), dtype=np.int)
             nface = 0
             for quad in quads:
                 if (quad[0] % 2) == 0:
@@ -277,3 +281,44 @@ def read_surface(filepath):
 
     coords = coords.astype(np.float)  # XXX: due to mayavi bug on mac 32bits
     return coords, faces
+
+
+###############################################################################
+# Write
+
+def write_bem_surface(fname, surf):
+    """Read one bem surface
+
+    Parameters
+    ----------
+    fname : string
+        File to write
+
+    surf : dict
+        A surface structued as obtained with read_bem_surfaces
+    """
+
+    # Create the file and save the essentials
+    fid = start_file(fname)
+
+    start_block(fid, FIFFB_BEM)
+    start_block(fid, FIFFB_BEM_SURF)
+
+    write_int(fid, FIFF_BEM_SURF_ID, surf['id'])
+    write_float(fid, FIFF_BEM_SIGMA, surf['sigma'])
+    write_int(fid, FIFF_BEM_SURF_NNODE, surf['np'])
+    write_int(fid, FIFF_BEM_SURF_NTRI, surf['ntri'])
+    write_int(fid, FIFF_BEM_COORD_FRAME, surf['coord_frame'])
+    # write_float_matrix(fid, FIFF_BEM_SURF_NODES, surf['rr'])
+    write_float_matrix(fid, FIFF_BEM_SURF_NODES, surf['rr'])
+
+    if 'nn' in surf and surf['nn'] is not None and len(surf['nn']) > 0:
+        write_float_matrix(fid, FIFF.FIFF_MNE_SOURCE_SPACE_NORMALS, surf['nn'])
+
+    # index start at 0 in Python
+    write_int_matrix(fid, FIFF_BEM_SURF_TRIANGLES, surf['tris'] + 1)
+
+    end_block(fid, FIFFB_BEM_SURF)
+    end_block(fid, FIFFB_BEM)
+
+    end_file(fid)
diff --git a/mne/tests/test_surface.py b/mne/tests/test_surface.py
index 232b7e0..ee1a793 100644
--- a/mne/tests/test_surface.py
+++ b/mne/tests/test_surface.py
@@ -1,9 +1,9 @@
 import os.path as op
 
-# from numpy.testing import assert_array_almost_equal
+from numpy.testing import assert_array_almost_equal
 
 from ..datasets import sample
-from .. import read_bem_surfaces
+from .. import read_bem_surfaces, write_bem_surface
 
 examples_folder = op.join(op.dirname(__file__), '..', '..', 'examples')
 data_path = sample.data_path(examples_folder)
@@ -14,6 +14,12 @@ fname = op.join(data_path, 'subjects', 'sample', 'bem',
 def test_io_bem_surfaces():
     """Testing reading of bem surfaces
     """
-    surf = read_bem_surfaces(fname, add_geom=False)
     surf = read_bem_surfaces(fname, add_geom=True)
+    surf = read_bem_surfaces(fname, add_geom=False)
     print "Number of surfaces : %d" % len(surf)
+
+    write_bem_surface('bem_surf.fif', surf[0])
+    surf_read = read_bem_surfaces('bem_surf.fif', add_geom=False)
+
+    for key in surf[0].keys():
+        assert_array_almost_equal(surf[0][key], surf_read[0][key])

-- 
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