Merge pull request #135 from effigies/master

NF: Write capabilities for freesurfer triangle files

Produces files identical to freesurfer up to the point that nibabel
reads. There is additional information after a freesurfer generated
file, but nibabel reads equivalent data and PySurfer displays these
files fine.

Author: matthew-brett

nibabel/freesurfer/__init__.py

@@ -2,5 +2,5 @@
 """
 
 from io import read_geometry, read_morph_data, \
-               read_annot, read_label
+               read_annot, read_label, write_geometry
 from mghformat import load, save, MGHImage

nibabel/freesurfer/io.py

@@ -1,6 +1,8 @@
 from __future__ import with_statement
 
 import numpy as np
+import getpass
+import time
 
 
 def _fread3(fobj):
@@ -93,6 +95,37 @@ def read_geometry(filepath):
     return coords, faces
 
 
+def write_geometry(filepath, coords, faces, create_stamp=None):
+    """Write a triangular format Freesurfer surface mesh.
+
+    Parameters
+    ----------
+    filepath : str
+        Path to surface file
+    coords : numpy array
+        nvtx x 3 array of vertex (x, y, z) coordinates
+    faces : numpy array
+        nfaces x 3 array of defining mesh triangles
+    create_stamp : str
+        User/time stamp (default: "created by <user> on <ctime>")
+    """
+    magic_bytes = np.array([255, 255, 254], dtype=np.uint8)
+
+    if create_stamp is None:
+        create_stamp = "created by %s on %s" % (getpass.getuser(),
+            time.ctime())
+
+    with open(filepath, 'wb') as fobj:
+        magic_bytes.tofile(fobj)
+        fobj.write("%s\n\n" % create_stamp)
+
+        np.array([coords.shape[0], faces.shape[0]], dtype='>i4').tofile(fobj)
+
+        # Coerce types, just to be safe
+        coords.astype('>f4').reshape(-1).tofile(fobj)
+        faces.astype('>i4').reshape(-1).tofile(fobj)
+
+
 def read_morph_data(filepath):
     """Read a Freesurfer morphometry data file.
 

nibabel/freesurfer/tests/test_io.py

@@ -1,11 +1,17 @@
+from __future__ import with_statement
 import os
 from os.path import join as pjoin
+import getpass
+import time
+
+from nibabel.tmpdirs import InTemporaryDirectory
 
 from nose.tools import assert_true
 import numpy as np
 from numpy.testing import assert_equal
 
-from .. import read_geometry, read_morph_data, read_annot, read_label
+from .. import read_geometry, read_morph_data, read_annot, read_label, \
+                write_geometry
 
 
 have_freesurfer = True
@@ -36,6 +42,31 @@ def test_geometry():
     assert_equal(0, faces.min())
     assert_equal(coords.shape[0], faces.max() + 1)
 
+    # Test equivalence of freesurfer- and nibabel-generated triangular files
+    # with respect to read_geometry()
+    with InTemporaryDirectory():
+        surf_path = 'test'
+        create_stamp = "created by %s on %s" % (getpass.getuser(),
+            time.ctime())
+        write_geometry(surf_path, coords, faces, create_stamp)
+
+        coords2, faces2 = read_geometry(surf_path)
+
+        with open(surf_path, 'rb') as fobj:
+            magic = np.fromfile(fobj, ">u1", 3)
+            read_create_stamp = fobj.readline().rstrip('\n')
+
+    assert_equal(create_stamp, read_create_stamp)
+
+    np.testing.assert_array_equal(coords, coords2)
+    np.testing.assert_array_equal(faces, faces2)
+
+    # Validate byte ordering
+    coords_swapped = coords.byteswap().newbyteorder()
+    faces_swapped = faces.byteswap().newbyteorder()
+    np.testing.assert_array_equal(coords_swapped, coords)
+    np.testing.assert_array_equal(faces_swapped, faces)
+
 
 @freesurfer_test
 def test_morph_data():