Merge pull request #461 from jaeilepp/fix-surf

MRG: Fix to surface reading with new quad.

Follow up to PR #460 which implemented read of surface data.

The data format is actually different when reading new quad files.

Use a file from freesurfer that has this format (added in the nibabel-data submodules) to test.

Author: matthew-brett

nibabel/freesurfer/io.py

@@ -46,6 +46,7 @@ def _fread3_many(fobj, n):
 
 
 def _read_volume_info(fobj):
+    """Helper for reading the footer from a surface file."""
     volume_info = OrderedDict()
     head = np.fromfile(fobj, '>i4', 1)
     if not np.array_equal(head, [20]):  # Read two bytes more
@@ -106,13 +107,17 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False):
     """
     volume_info = OrderedDict()
 
+    TRIANGLE_MAGIC = 16777214
+    QUAD_MAGIC = 16777215
+    NEW_QUAD_MAGIC = 16777213
     with open(filepath, "rb") as fobj:
         magic = _fread3(fobj)
-        if magic in (16777215, 16777213):  # Quad file
+        if magic in (QUAD_MAGIC, NEW_QUAD_MAGIC):  # 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
+            (fmt, div) = (">i2", 100.) if magic == QUAD_MAGIC else (">f4", 1.)
+            coords = np.fromfile(fobj, fmt, nvert * 3).astype(np.float) / div
+            coords = coords.reshape(-1, 3)
             quads = _fread3_many(fobj, nquad * 4)
             quads = quads.reshape(nquad, 4)
             #
@@ -132,7 +137,7 @@ def read_geometry(filepath, read_metadata=False, read_stamp=False):
                     faces[nface] = quad[0], quad[2], quad[3]
                     nface += 1
 
-        elif magic == 16777214:  # Triangle file
+        elif magic == TRIANGLE_MAGIC:  # Triangle file
             create_stamp = fobj.readline().rstrip(b'\n').decode('utf-8')
             fobj.readline()
             vnum = np.fromfile(fobj, ">i4", 1)[0]

nibabel/freesurfer/tests/test_io.py

@@ -16,7 +16,7 @@
 from .. import (read_geometry, read_morph_data, read_annot, read_label,
                 write_geometry, write_morph_data, write_annot)
 
-from ...tests.nibabel_data import get_nibabel_data
+from ...tests.nibabel_data import get_nibabel_data, needs_nibabel_data
 from ...fileslice import strided_scalar
 from ...testing import clear_and_catch_warnings
 
@@ -114,6 +114,23 @@ def test_geometry():
     np.testing.assert_array_equal(faces_swapped, faces)
 
 
+@freesurfer_test
+@needs_nibabel_data('nitest-freesurfer')
+def test_quad_geometry():
+    """Test IO of freesurfer quad files."""
+    new_quad = pjoin(get_nibabel_data(), 'nitest-freesurfer', 'subjects',
+                     'bert', 'surf', 'lh.inflated.nofix')
+    coords, faces = read_geometry(new_quad)
+    assert_equal(0, faces.min())
+    assert_equal(coords.shape[0], faces.max() + 1)
+    with InTemporaryDirectory():
+        new_path = 'test'
+        write_geometry(new_path, coords, faces)
+        coords2, faces2 = read_geometry(new_path)
+        assert_equal(coords, coords2)
+        assert_equal(faces, faces2)
+
+
 @freesurfer_test
 def test_morph_data():
     """Test IO of morphometry data file (eg. curvature)."""