Merge pull request #414 from effigies/write_morph

ENH: Add write_morph_data to freesurfer.io

The first commit was in the midst of #249 and causing PEP8 warnings. Not otherwise relevant to Cifti, so I'm submitting as its own PR. Fixed the style and brought a little more in line with freesurfer.io functions.

Author: matthew-brett

nibabel/freesurfer/__init__.py

@@ -1,6 +1,6 @@
 """Reading functions for freesurfer files
 """
 
-from .io import read_geometry, read_morph_data, \
+from .io import read_geometry, read_morph_data, write_morph_data, \
     read_annot, read_label, write_geometry, write_annot
 from .mghformat import load, save, MGHImage

nibabel/freesurfer/io.py

@@ -6,6 +6,7 @@
 
 
 from .. externals.six.moves import xrange
+from ..openers import Opener
 
 
 def _fread3(fobj):
@@ -160,6 +161,52 @@ def read_morph_data(filepath):
     return curv
 
 
+def write_morph_data(file_like, values, fnum=0):
+    """Write Freesurfer morphometry data `values` to file-like `file_like`
+
+    Equivalent to FreeSurfer's `write_curv.m`_
+
+    See also:
+    http://www.grahamwideman.com/gw/brain/fs/surfacefileformats.htm#CurvNew
+
+    .. _write_curv.m: \
+    https://github.com/neurodebian/freesurfer/blob/debian-sloppy/matlab/write_curv.m
+
+    Parameters
+    ----------
+    file_like : file-like
+        String containing path of file to be written, or file-like object, open
+        in binary write (`'wb'` mode, implementing the `write` method)
+    values : array-like
+        Surface morphometry values
+
+        Shape must be (N,), (N, 1), (1, N) or (N, 1, 1)
+    fnum : int, optional
+        Number of faces in the associated surface
+    """
+    magic_bytes = np.array([255, 255, 255], dtype=np.uint8)
+
+    vector = np.asarray(values)
+    vnum = np.prod(vector.shape)
+    if vector.shape not in ((vnum,), (vnum, 1), (1, vnum), (vnum, 1, 1)):
+        raise ValueError("Invalid shape: argument values must be a vector")
+
+    i4info = np.iinfo('i4')
+    if vnum > i4info.max:
+        raise ValueError("Too many values for morphometry file")
+    if not i4info.min <= fnum <= i4info.max:
+        raise ValueError("Argument fnum must be between {0} and {1}".format(
+                         i4info.min, i4info.max))
+
+    with Opener(file_like, 'wb') as fobj:
+        fobj.write(magic_bytes)
+
+        # vertex count, face count (unused), vals per vertex (only 1 supported)
+        fobj.write(np.array([vnum, fnum, 1], dtype='>i4'))
+
+        fobj.write(vector.astype('>f4'))
+
+
 def read_annot(filepath, orig_ids=False):
     """Read in a Freesurfer annotation from a .annot file.
 

nibabel/freesurfer/tests/test_io.py

@@ -10,12 +10,13 @@
 
 from nose.tools import assert_true
 import numpy as np
-from numpy.testing import assert_equal, dec
+from numpy.testing import assert_equal, assert_raises, dec
 
 from .. import (read_geometry, read_morph_data, read_annot, read_label,
-                write_geometry, write_annot)
+                write_geometry, write_morph_data, write_annot)
 
 from ...tests.nibabel_data import get_nibabel_data
+from ...fileslice import strided_scalar
 
 
 DATA_SDIR = 'fsaverage'
@@ -92,6 +93,33 @@ def test_morph_data():
     curv = read_morph_data(curv_path)
     assert_true(-1.0 < curv.min() < 0)
     assert_true(0 < curv.max() < 1.0)
+    with InTemporaryDirectory():
+        new_path = 'test'
+        write_morph_data(new_path, curv)
+        curv2 = read_morph_data(new_path)
+        assert_equal(curv2, curv)
+
+
+def test_write_morph_data():
+    """Test write_morph_data edge cases"""
+    values = np.arange(20, dtype='>f4')
+    okay_shapes = [(20,), (20, 1), (20, 1, 1), (1, 20)]
+    bad_shapes = [(10, 2), (1, 1, 20, 1, 1)]
+    big_num = np.iinfo('i4').max + 1
+    with InTemporaryDirectory():
+        for shape in okay_shapes:
+            write_morph_data('test.curv', values.reshape(shape))
+            # Check ordering is preserved, regardless of shape
+            assert_equal(values, read_morph_data('test.curv'))
+        assert_raises(ValueError, write_morph_data, 'test.curv',
+                      np.zeros(shape), big_num)
+        # Windows 32-bit overflows Python int
+        if np.dtype(np.int) != np.dtype(np.int32):
+            assert_raises(ValueError, write_morph_data, 'test.curv',
+                          strided_scalar((big_num,)))
+        for shape in bad_shapes:
+            assert_raises(ValueError, write_morph_data, 'test.curv',
+                          values.reshape(shape))
 
 
 @freesurfer_test