Merge branch 'enh/flexi_gifti' into 'master'

Enh/flexi gifti

See merge request fsl/fslpy!96

CHANGELOG.rst

@@ -23,6 +23,8 @@ Changed
 * Minimum required version of ``nibabel`` is now 2.3.
 * The :class:`.Image` class now fully delegates to ``nibabel`` for managing
   file handles.
+* The :class:`.GiftiMesh` class can now load surface files which contain
+  vertex data.
 
 
 Removed

fsl/data/gifti.py

@@ -19,6 +19,7 @@ are available:
      GiftiMesh
      loadGiftiMesh
      loadGiftiVertexData
+     prepareGiftiVertexData
      relatedFiles
 """
 
@@ -34,15 +35,13 @@ import fsl.data.constants as constants
 import fsl.data.mesh      as fslmesh
 
 
-# We include '.gii' here because not all surface
-# GIFTIs follow the file suffix convention.
 ALLOWED_EXTENSIONS = ['.surf.gii', '.gii']
 """List of file extensions that a file containing Gifti surface data
 is expected to have.
 """
 
 
-EXTENSION_DESCRIPTIONS = ['GIFTI surface file', 'GIFTI file']
+EXTENSION_DESCRIPTIONS = ['GIFTII surface file', 'GIFTI file']
 """A description for each of the :data:`ALLOWED_EXTENSIONS`. """
 
 
@@ -60,7 +59,8 @@ class GiftiMesh(fslmesh.Mesh):
         """Load the given GIFTI file using ``nibabel``, and extracts surface
         data using the  :func:`loadGiftiMesh` function.
 
-        :arg infile:     A GIFTI surface file (``*.surf.gii``).
+        :arg infile:     A GIFTI file (``*..gii``) which contains a surface
+                         definition.
 
         :arg fixWinding: Passed through to the :meth:`addVertices` method
                          for the first vertex set.
@@ -76,34 +76,42 @@ class GiftiMesh(fslmesh.Mesh):
         name   = fslpath.removeExt(op.basename(infile), ALLOWED_EXTENSIONS)
         infile = op.abspath(infile)
 
-        surfimg, vertices, indices = loadGiftiMesh(infile)
-
+        surfimg, indices, vertices, vdata = loadGiftiMesh(infile)
 
         fslmesh.Mesh.__init__(self,
                               indices,
                               name=name,
                               dataSource=infile)
 
-        self.addVertices(vertices, infile, fixWinding=fixWinding)
+        for i, v in enumerate(vertices):
+            if i == 0: key = infile
+            else:      key = '{} [{}]'.format(infile, i)
+            self.addVertices(v, key, select=(i == 0), fixWinding=fixWinding)
         self.setMeta(infile, surfimg)
 
+        if vdata is not None:
+            self.addVertexData(infile, vdata)
+
         # Find and load all other
         # surfaces in the same directory
         # as the specfiied one.
         if loadAll:
 
-            nvertices = vertices.shape[0]
+            nvertices = vertices[0].shape[0]
             surfFiles = relatedFiles(infile, ALLOWED_EXTENSIONS)
 
             for sfile in surfFiles:
 
-                surfimg, vertices, _ = loadGiftiMesh(sfile)
+                try:
+                    surfimg, _, vertices, _ = loadGiftiMesh(sfile)
+                except Exception:
+                    continue
 
-                if vertices.shape[0] != nvertices:
+                if vertices[0].shape[0] != nvertices:
                     continue
 
-                self.addVertices(vertices, sfile, select=False)
-                self.setMeta(    sfile, surfimg)
+                self.addVertices(vertices[0], sfile, select=False)
+                self.setMeta(sfile, surfimg)
 
 
     def loadVertices(self, infile, key=None, *args, **kwargs):
@@ -123,7 +131,7 @@ class GiftiMesh(fslmesh.Mesh):
         if key is None:
             key = infile
 
-        surfimg, vertices, _ = loadGiftiMesh(infile)
+        surfimg, _, vertices, _ = loadGiftiMesh(infile)
 
         vertices = self.addVertices(vertices, key, *args, **kwargs)
 
@@ -157,9 +165,10 @@ def loadGiftiMesh(filename):
 
     The image is expected to contain the following``<DataArray>`` elements:
 
-      - one comprising ``NIFTI_INTENT_POINTSET`` data (the surface vertices)
       - one comprising ``NIFTI_INTENT_TRIANGLE`` data (vertex indices
         defining the triangles).
+      - one or more comprising ``NIFTI_INTENT_POINTSET`` data (the surface
+        vertices)
 
     A ``ValueError`` will be raised if this is not the case.
 
@@ -169,42 +178,65 @@ def loadGiftiMesh(filename):
 
                    - The loaded ``nibabel.gifti.GiftiImage`` instance
 
-                   - A ``(N, 3)`` array containing ``N`` vertices.
-
-                   - A ``(M, 3))`` array containing the vertex indices for
+                   - A ``(M, 3)`` array containing the vertex indices for
                      ``M`` triangles.
-    """
 
-    gimg = nib.load(filename)
+                   - A list of at least one ``(N, 3)`` arrays containing ``N``
+                     vertices.
 
-    pointsetCode = constants.NIFTI_INTENT_POINTSET
-    triangleCode = constants.NIFTI_INTENT_TRIANGLE
+                   - A ``(M, N)`` numpy array containing ``N`` data points for
+                     ``M`` vertices, or ``None`` if the file does not contain
+                     any vertex data.
+    """
 
-    pointsets = [d for d in gimg.darrays if d.intent == pointsetCode]
-    triangles = [d for d in gimg.darrays if d.intent == triangleCode]
+    gimg = nib.load(filename)
 
-    if len(gimg.darrays) != 2:
-        raise ValueError('{}: GIFTI surface files must contain '
-                         'exactly one pointset array and one '
-                         'triangle array'.format(filename))
+    pscode  = constants.NIFTI_INTENT_POINTSET
+    tricode = constants.NIFTI_INTENT_TRIANGLE
 
-    if len(pointsets) != 1:
-        raise ValueError('{}: GIFTI surface files must contain '
-                         'exactly one pointset array'.format(filename))
+    pointsets = [d for d in gimg.darrays if d.intent == pscode]
+    triangles = [d for d in gimg.darrays if d.intent == tricode]
+    vdata     = [d for d in gimg.darrays if d.intent not in (pscode, tricode)]
 
     if len(triangles) != 1:
         raise ValueError('{}: GIFTI surface files must contain '
                          'exactly one triangle array'.format(filename))
 
-    vertices = pointsets[0].data
+    if len(pointsets) == 0:
+        raise ValueError('{}: GIFTI surface files must contain '
+                         'at least one pointset array'.format(filename))
+
+    vertices = [ps.data for ps in pointsets]
     indices  = triangles[0].data
 
-    return gimg, vertices, indices
+    if len(vdata) == 0: vdata = None
+    else:               vdata = prepareGiftiVertexData(vdata, filename)
+
+    return gimg, indices, vertices, vdata
 
 
 def loadGiftiVertexData(filename):
     """Loads vertex data from the given GIFTI file.
 
+    See :func:`prepareGiftiVertexData`.
+
+    Returns a tuple containing:
+
+      - The loaded ``nibabel.gifti.GiftiImage`` object
+
+      - A ``(M, N)`` numpy array containing ``N`` data points for ``M``
+        vertices
+    """
+    gimg = nib.load(filename)
+    return gimg, prepareGiftiVertexData(gimg.darrays, filename)
+
+
+def prepareGiftiVertexData(darrays, filename=None):
+    """Prepares vertex data from the given list of GIFTI data arrays.
+
+    All of the data arrays are concatenated into one ``(M, N)`` array,
+    containing ``N`` data points for ``M`` vertices.
+
     It is assumed that the given file does not contain any
     ``NIFTI_INTENT_POINTSET`` or ``NIFTI_INTENT_TRIANGLE`` data arrays, and
     which contains either:
@@ -215,17 +247,11 @@ def loadGiftiVertexData(filename):
       - One or more ``(M, 1)`` data arrays each containing a single data point
         for ``M`` vertices, and all with the same intent code
 
-    Returns a tuple containing:
-
-      - The loaded ``nibabel.gifti.GiftiImage`` object
-
-      - A ``(M, N)`` numpy array containing ``N`` data points for ``M``
-        vertices
+    Returns a ``(M, N)`` numpy array containing ``N`` data points for ``M``
+    vertices.
     """
 
-    gimg = nib.load(filename)
-
-    intents = set([d.intent for d in gimg.darrays])
+    intents = set([d.intent for d in darrays])
 
     if len(intents) != 1:
         raise ValueError('{} contains multiple (or no) intents'
@@ -235,20 +261,19 @@ def loadGiftiVertexData(filename):
 
     if intent in (constants.NIFTI_INTENT_POINTSET,
                   constants.NIFTI_INTENT_TRIANGLE):
-
         raise ValueError('{} contains surface data'.format(filename))
 
     # Just a single array - return it as-is.
     # n.b. Storing (M, N) data in a single
     # DataArray goes against the GIFTI spec,
     # but hey, it happens.
-    if len(gimg.darrays) == 1:
-        vdata = gimg.darrays[0].data
-        return gimg, vdata.reshape(vdata.shape[0], -1)
+    if len(darrays) == 1:
+        vdata = darrays[0].data
+        return vdata.reshape(vdata.shape[0], -1)
 
     # Otherwise extract and concatenate
     # multiple 1-dimensional arrays
-    vdata = [d.data for d in gimg.darrays]
+    vdata = [d.data for d in darrays]
 
     if any([len(d.shape) != 1 for d in vdata]):
         raise ValueError('{} contains one or more non-vector '
@@ -257,7 +282,7 @@ def loadGiftiVertexData(filename):
     vdata = np.vstack(vdata).T
     vdata = vdata.reshape(vdata.shape[0], -1)
 
-    return gimg, vdata
+    return vdata
 
 
 def relatedFiles(fname, ftypes=None):

tests/test_gifti.py

@@ -70,11 +70,12 @@ def test_loadGiftiMesh():
     testdir  = op.join(op.dirname(__file__), 'testdata')
     testfile = op.join(testdir, 'example.surf.gii')
 
-    gimg, verts, idxs = gifti.loadGiftiMesh(testfile)
+    gimg, idxs, verts, _ = gifti.loadGiftiMesh(testfile)
 
     assert isinstance(gimg, nib.gifti.GiftiImage)
-    assert tuple(verts.shape) == (642,  3)
-    assert tuple(idxs.shape)  == (1280, 3)
+    assert len(verts)            == 1
+    assert tuple(verts[0].shape) == (642,  3)
+    assert tuple(idxs.shape)     == (1280, 3)
 
     badfiles = glob.glob(op.join(testdir, 'example_bad*surf.gii'))
 
@@ -234,3 +235,67 @@ def test_relatedFiles():
         for s in rsurfaces:
             result = gifti.relatedFiles(s)
             assert sorted(rrelated) == sorted(result)
+
+TEST_VERTS = np.array([
+    [0, 0, 0],
+    [1, 0, 0],
+    [0, 1, 0],
+    [0, 0, 1]])
+TEST_IDXS = np.array([
+    [0, 1, 2],
+    [0, 3, 1],
+    [0, 2, 3],
+    [1, 3, 2]])
+
+TEST_VERT_ARRAY = nib.gifti.GiftiDataArray(
+    TEST_VERTS, intent='NIFTI_INTENT_POINTSET')
+TEST_IDX_ARRAY  = nib.gifti.GiftiDataArray(
+    TEST_IDXS, intent='NIFTI_INTENT_TRIANGLE')
+
+def test_GiftiMesh_surface_and_data():
+
+    data1   = np.random.randint(0, 10, len(TEST_VERTS))
+    data2   = np.random.randint(0, 10, len(TEST_VERTS))
+    expdata = np.vstack([data1, data2]).T
+    verts   = TEST_VERT_ARRAY
+    tris    = TEST_IDX_ARRAY
+    data1   = nib.gifti.GiftiDataArray(data1, intent='NIFTI_INTENT_SHAPE')
+    data2   = nib.gifti.GiftiDataArray(data2, intent='NIFTI_INTENT_SHAPE')
+    gimg    = nib.gifti.GiftiImage(darrays=[verts, tris, data1, data2])
+
+    with tempdir():
+        fname = op.abspath('test.gii')
+        gimg.to_filename(fname)
+        surf = gifti.GiftiMesh(fname)
+
+        assert np.all(surf.vertices  == TEST_VERTS)
+        assert np.all(surf.indices   == TEST_IDXS)
+        assert surf.vertexDataSets() == [fname]
+        assert np.all(surf.getVertexData(fname) == expdata)
+
+
+
+def test_GiftiMesh_multiple_vertices():
+
+    tris   = TEST_IDX_ARRAY
+    verts1 = TEST_VERT_ARRAY
+    verts2 = nib.gifti.GiftiDataArray(
+        TEST_VERTS * 5, intent='NIFTI_INTENT_POINTSET')
+
+    gimg  = nib.gifti.GiftiImage(darrays=[verts1, verts2, tris])
+
+    with tempdir():
+        fname = op.abspath('test.gii')
+        gimg.to_filename(fname)
+        surf  = gifti.GiftiMesh(fname)
+
+        expvsets = [fname,
+                    '{} [{}]'.format(fname, 1)]
+
+        assert np.all(surf.vertices == TEST_VERTS)
+        assert np.all(surf.indices  == TEST_IDXS)
+        assert  surf.vertexSets()   == expvsets
+
+        surf.vertices = expvsets[1]
+
+        assert np.all(surf.vertices == TEST_VERTS * 5)