diff --git a/CHANGELOG.rst b/CHANGELOG.rst index b792deadc8444062fbfc6f81e12555ba3901982f..ae032ca0801613df43503e37426aaf164d1a272e 100644 --- a/CHANGELOG.rst +++ b/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 diff --git a/fsl/data/gifti.py b/fsl/data/gifti.py index b0b018c20e5c1955597fb75f6c740fc796ba861e..c4d8ce9c79b8062b413fceaf14a5581b68503437 100644 --- a/fsl/data/gifti.py +++ b/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): diff --git a/tests/test_gifti.py b/tests/test_gifti.py index b01635d1f85fc768dab8f6808c429a150d430f63..e39e05b6e54532ed36e4d93086335c0726468b31 100644 --- a/tests/test_gifti.py +++ b/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)