Merge pull request #47 from effigies/master

RF Use nibabel IO functions where available

Author: mwaskom

examples/plot_foci.py

@@ -13,8 +13,8 @@
 import os.path as op
 from numpy import arange
 from numpy.random import permutation
+import nibabel as nib
 from surfer import Brain
-from surfer import io
 
 subject_id = "fsaverage"
 subjects_dir = os.environ["SUBJECTS_DIR"]
@@ -61,7 +61,7 @@
 and find 10 random vertices within the STS.
 """
 annot_path = op.join(subjects_dir, subject_id, "label/lh.aparc.a2009s.annot")
-ids, ctab, names = io.read_annot(annot_path)
+ids, ctab, names = nib.freesurfer.read_annot(annot_path)
 verts = arange(0, len(ids))
 coords = permutation(verts[ids == 74])[:10]
 

examples/plot_parc_values.py

@@ -12,7 +12,8 @@
 import os
 import os.path as op
 import numpy as np
-from surfer import io, Brain
+import nibabel as nib
+from surfer import Brain
 
 subject_id = "fsaverage"
 hemi = "lh"
@@ -31,7 +32,7 @@
 aparc_file = op.join(os.environ["SUBJECTS_DIR"],
                      subject_id, "label",
                      hemi + ".aparc.a2009s.annot")
-labels, ctab, names = io.read_annot(aparc_file)
+labels, ctab, names = nib.freesurfer.read_annot(aparc_file)
 
 """
 Make a random vector of scalar data corresponding to

surfer/io.py

@@ -40,105 +40,6 @@ def _get_subjects_dir(subjects_dir=None):
     return subjects_dir
 
 
-def _fread3(fobj):
-    """Read a 3-byte int from an open binary file object."""
-    b1, b2, b3 = np.fromfile(fobj, ">u1", 3)
-    return (b1 << 16) + (b2 << 8) + b3
-
-
-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
-    return (b1 << 16) + (b2 << 8) + b3
-
-
-def read_geometry(filepath):
-    """Read a triangular format Freesurfer surface mesh.
-
-    Parameters
-    ----------
-    filepath : str
-        Path to surface file
-
-    Returns
-    -------
-    coords : numpy array
-        nvtx x 3 array of vertex (x, y, z) coordinates
-    faces : numpy array
-        nfaces x 3 array of defining mesh triangles
-    """
-    with open(filepath, "rb") as fobj:
-        magic = _fread3(fobj)
-        if magic == 16777215:  # 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
-            quads = _fread3_many(fobj, nquad * 4)
-            quads = quads.reshape(nquad, 4)
-            #
-            #   Face splitting follows
-            #
-            faces = np.zeros((2 * nquad, 3), dtype=np.int)
-            nface = 0
-            for quad in quads:
-                if (quad[0] % 2) == 0:
-                    faces[nface] = quad[0], quad[1], quad[3]
-                    nface += 1
-                    faces[nface] = quad[2], quad[3], quad[1]
-                    nface += 1
-                else:
-                    faces[nface] = quad[0], quad[1], quad[2]
-                    nface += 1
-                    faces[nface] = quad[0], quad[2], quad[3]
-                    nface += 1
-
-        elif magic == 16777214:  # Triangle file
-            create_stamp = fobj.readline()
-            _ = fobj.readline()
-            vnum = np.fromfile(fobj, ">i4", 1)[0]
-            fnum = np.fromfile(fobj, ">i4", 1)[0]
-            coords = np.fromfile(fobj, ">f4", vnum * 3).reshape(vnum, 3)
-            faces = np.fromfile(fobj, ">i4", fnum * 3).reshape(fnum, 3)
-        else:
-            raise ValueError("File does not appear to be a Freesurfer surface")
-
-    coords = coords.astype(np.float)  # XXX: due to mayavi bug on mac 32bits
-    return coords, faces
-
-
-def read_morph_data(filepath):
-    """Read a Freesurfer morphometry data file.
-
-    This function reads in what Freesurfer internally calls "curv" file types,
-    (e.g. ?h. curv, ?h.thickness), but as that has the potential to cause
-    confusion where "curv" also refers to the surface curvature values,
-    we refer to these files as "morphometry" files with PySurfer.
-
-    Parameters
-    ----------
-    filepath : str
-        Path to morphometry file
-
-    Returns
-    -------
-    curv : numpy array
-        Vector representation of surface morpometry values
-
-    """
-    with open(filepath, "rb") as fobj:
-        magic = _fread3(fobj)
-        if magic == 16777215:
-            vnum = np.fromfile(fobj, ">i4", 3)[0]
-            curv = np.fromfile(fobj, ">f4", vnum)
-        else:
-            vnum = magic
-            _ = _fread3(fobj)
-            curv = np.fromfile(fobj, ">i2", vnum) / 100
-    return curv
-
-
 def read_scalar_data(filepath):
     """Load in scalar data from an image.
 
@@ -199,76 +100,6 @@ def read_scalar_data(filepath):
     return scalar_data
 
 
-def read_annot(filepath, orig_ids=False):
-    """Read in a Freesurfer annotation from a .annot file.
-
-    Parameters
-    ----------
-    filepath : str
-        Path to annotation file
-    orig_ids : bool
-        Whether to return the vertex ids as stored in the annotation
-        file or the positional colortable ids
-
-    Returns
-    -------
-    labels : n_vtx numpy array
-        Annotation id at each vertex
-    ctab : numpy array
-        RGBA + label id colortable array
-    names : numpy array
-        Array of region names as stored in the annot file
-
-    """
-    with open(filepath, "rb") as fobj:
-        dt = ">i4"
-        vnum = np.fromfile(fobj, dt, 1)[0]
-        data = np.fromfile(fobj, dt, vnum * 2).reshape(vnum, 2)
-        labels = data[:, 1]
-        ctab_exists = np.fromfile(fobj, dt, 1)[0]
-        if not ctab_exists:
-            raise Exception('Color table not found in annotation file')
-        n_entries = np.fromfile(fobj, dt, 1)[0]
-        if n_entries > 0:
-            length = np.fromfile(fobj, dt, 1)[0]
-            orig_tab = np.fromfile(fobj, '>c', length)
-            orig_tab = orig_tab[:-1]
-
-            names = list()
-            ctab = np.zeros((n_entries, 5), np.int)
-            for i in xrange(n_entries):
-                name_length = np.fromfile(fobj, dt, 1)[0]
-                name = np.fromfile(fobj, "|S%d" % name_length, 1)[0]
-                names.append(name)
-                ctab[i, :4] = np.fromfile(fobj, dt, 4)
-                ctab[i, 4] = (ctab[i, 0] + ctab[i, 1] * (2 ** 8) +
-                              ctab[i, 2] * (2 ** 16) +
-                              ctab[i, 3] * (2 ** 24))
-        else:
-            ctab_version = -n_entries
-            if ctab_version != 2:
-                raise Exception('Color table version not supported')
-            n_entries = np.fromfile(fobj, dt, 1)[0]
-            ctab = np.zeros((n_entries, 5), np.int)
-            length = np.fromfile(fobj, dt, 1)[0]
-            _ = np.fromfile(fobj, "|S%d" % length, 1)[0]  # Orig table path
-            entries_to_read = np.fromfile(fobj, dt, 1)[0]
-            names = list()
-            for i in xrange(entries_to_read):
-                _ = np.fromfile(fobj, dt, 1)[0]  # Structure
-                name_length = np.fromfile(fobj, dt, 1)[0]
-                name = np.fromfile(fobj, "|S%d" % name_length, 1)[0]
-                names.append(name)
-                ctab[i, :4] = np.fromfile(fobj, dt, 4)
-                ctab[i, 4] = (ctab[i, 0] + ctab[i, 1] * (2 ** 8) +
-                                ctab[i, 2] * (2 ** 16))
-        ctab[:, 3] = 255
-    if not orig_ids:
-        ord = np.argsort(ctab[:, -1])
-        labels = ord[np.searchsorted(ctab[ord, -1], labels)]
-    return labels, ctab, names
-
-
 def read_label(filepath, read_scalars=False):
     """Load in a Freesurfer .label file.
 
@@ -492,7 +323,12 @@ def __init__(self, subject_id, hemi, surf, subjects_dir=None):
     def load_geometry(self):
         surf_path = pjoin(self.data_path, "surf",
                           "%s.%s" % (self.hemi, self.surf))
-        self.coords, self.faces = read_geometry(surf_path)
+        self.coords, self.faces = nib.freesurfer.read_geometry(surf_path)
+
+    def save_geometry(self):
+        surf_path = pjoin(self.data_path, "surf",
+                          "%s.%s" % (self.hemi, self.surf))
+        nib.freesurfer.write_geometry(surf_path, self.coords, self.faces)
 
     @property
     def x(self):
@@ -509,7 +345,7 @@ def z(self):
     def load_curvature(self):
         """Load in curvature values from the ?h.curv file."""
         curv_path = pjoin(self.data_path, "surf", "%s.curv" % self.hemi)
-        self.curv = read_morph_data(curv_path)
+        self.curv = nib.freesurfer.read_morph_data(curv_path)
         self.bin_curv = np.array(self.curv > 0, np.int)
 
     def load_label(self, name):
@@ -521,8 +357,8 @@ def load_label(self, name):
         argument.
 
         """
-        label = read_label(pjoin(self.data_path, 'label',
-                                 '%s.%s.label' % (self.hemi, name)))
+        label = nib.freesurfer.read_label(pjoin(self.data_path, 'label',
+                                          '%s.%s.label' % (self.hemi, name)))
         label_array = np.zeros(len(self.x), np.int)
         label_array[label] = 1
         try:

surfer/tests/test_io.py

@@ -15,46 +15,6 @@
 data_path = pjoin(subj_dir, subject_id)
 
 
-def test_geometry():
-    """Test IO of .surf"""
-    surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "inflated"))
-    coords, faces = io.read_geometry(surf_path)
-    assert_equal(0, faces.min())
-    assert_equal(coords.shape[0], faces.max() + 1)
-
-    # Test quad with sphere
-    surf_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "sphere.reg"))
-    coords, faces = io.read_geometry(surf_path)
-    assert_equal(0, faces.min())
-    assert_equal(coords.shape[0], faces.max() + 1)
-
-
-def test_morph_data():
-    """Test IO of morphometry data file (eg. curvature)."""
-    curv_path = pjoin(data_path, "surf", "%s.%s" % ("lh", "curv"))
-    curv = io.read_morph_data(curv_path)
-    assert -1.0 < curv.min() < 0
-    assert 0 < curv.max() < 1.0
-
-
-def test_annot():
-    """Test IO of .annot"""
-    annots = ['aparc', 'aparc.a2005s']
-    for a in annots:
-        annot_path = pjoin(data_path, "label", "%s.%s.annot" % ("lh", a))
-        labels, ctab, names = io.read_annot(annot_path)
-        assert labels.shape == (163842, )
-        assert ctab.shape == (len(names), 5)
-
-
-def test_label():
-    """Test IO of .label"""
-    label_path = pjoin(data_path, "label", "lh.BA1.label")
-    label = io.read_label(label_path)
-    # XXX : test more
-    assert np.all(label > 0)
-
-
 def test_surface():
     """Test of Surface class"""
     for subjects_dir in [None, subj_dir]:

surfer/tests/test_viz.py

@@ -3,6 +3,7 @@
 from os.path import join as pjoin
 from numpy.testing import assert_raises
 from tempfile import mktemp
+import nibabel as nib
 
 from surfer import Brain
 from surfer import io
@@ -99,7 +100,7 @@ def test_foci():
     brain.add_foci(coords, map_surface="white", color="gold")
 
     annot_path = pjoin(subj_dir, subject_id, 'label', 'lh.aparc.a2009s.annot')
-    ids, ctab, names = io.read_annot(annot_path)
+    ids, ctab, names = nib.freesurfer.read_annot(annot_path)
     verts = np.arange(0, len(ids))
     coords = np.random.permutation(verts[ids == 74])[:10]
     scale_factor = 0.7

surfer/viz.py

@@ -7,6 +7,8 @@
 from scipy import ndimage
 from matplotlib.colors import colorConverter
 
+import nibabel as nib
+
 from . import io
 from . import utils
 from .io import Surface, _get_subjects_dir
@@ -598,7 +600,7 @@ def add_annotation(self, annot, borders=True, alpha=1):
                                  % filepath)
 
         # Read in the data
-        labels, cmap, _ = io.read_annot(filepath, orig_ids=True)
+        labels, cmap, _ = nib.freesurfer.read_annot(filepath, orig_ids=True)
 
         # Maybe zero-out the non-border vertices
         if borders:
@@ -828,7 +830,7 @@ def add_morphometry(self, measure, grayscale=False):
         view = mlab.view()
 
         # Read in the morphometric data
-        morph_data = io.read_morph_data(morph_file)
+        morph_data = nib.freesurfer.read_morph_data(morph_file)
 
         # Get a cortex mask for robust range
         self._geo.load_label("cortex")