gifti.py 13.2 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#!/usr/bin/env python
#
# gifti.py - GIFTI file support.
#
# Author: Paul McCarthy  <pauldmccarthy@gmail.com>
#         Michiel Cottar <michiel.cottaar@ndcn.ox.ac.uk>
#
"""This class provides classes and functions for working with GIFTI files.

The GIFTI file format specification can be found at
http://www.nitrc.org/projects/gifti/

Support is currently very basic - only the following classes/functions
are available:

  .. autosummary::
     :nosignatures:

19
20
21
     GiftiMesh
     loadGiftiMesh
     loadGiftiVertexData
22
     prepareGiftiVertexData
23
     relatedFiles
24
25
26
"""


27
import            glob
28
import            re
29
30
import os.path as op

31
import numpy   as np
32
33
import nibabel as nib

34
import fsl.utils.path     as fslpath
35
import fsl.utils.bids     as bids
36
37
import fsl.data.constants as constants
import fsl.data.mesh      as fslmesh
38
39


40
ALLOWED_EXTENSIONS = ['.surf.gii', '.gii']
41
42
43
44
45
"""List of file extensions that a file containing Gifti surface data
is expected to have.
"""


46
EXTENSION_DESCRIPTIONS = ['GIFTI surface file', 'GIFTI file']
47
48
49
"""A description for each of the :data:`ALLOWED_EXTENSIONS`. """


50
51
52
53
54
55
56
57
58
VERTEX_DATA_EXTENSIONS = ['.func.gii',
                          '.shape.gii',
                          '.label.gii',
                          '.time.gii']
"""File suffixes which are interpreted as GIFTI vertex data files,
containing data values for every vertex in the mesh.
"""


59
class GiftiMesh(fslmesh.Mesh):
60
61
62
    """Class which represents a GIFTI surface image. This is essentially
    just a 3D model made of triangles.

63
64
65
    For each GIFTI surface file that is loaded, the
    ``nibabel.gifti.GiftiImage`` instance is stored in the :class:`.Meta`
    store, with the absolute path to the surface file as the key.
66
67
68
    """


69
    def __init__(self, infile, fixWinding=False, loadAll=False):
70
        """Load the given GIFTI file using ``nibabel``, and extracts surface
71
        data using the  :func:`loadGiftiMesh` function.
72

73
74
75
76
77
78
        If the file contains more than one set of vertices, the additional
        ones are added with keys of the form ``infile_i``, where ``infile`` is
        the absolute path to the file, and ``i`` is an index number, starting
        from 1. See the :meth:`.addVertices` method.

        :arg infile:     A GIFTI file (``*.gii``) which contains a surface
79
                         definition.
80
81
82
83
84
85
86

        :arg fixWinding: Passed through to the :meth:`addVertices` method
                         for the first vertex set.

        :arg loadAll:    If ``True``, the ``infile`` directory is scanned
                         for other surface files which are then loaded
                         as additional vertex sets.
87
88
89

        .. todo:: Allow loading from a ``.topo.gii`` and ``.coord.gii`` file?
                  Maybe.
90
91
        """

92
93
94
        name   = fslpath.removeExt(op.basename(infile), ALLOWED_EXTENSIONS)
        infile = op.abspath(infile)

95
        surfimg, indices, vertices, vdata = loadGiftiMesh(infile)
96

97
98
99
100
        fslmesh.Mesh.__init__(self,
                              indices,
                              name=name,
                              dataSource=infile)
101

102
103
        for i, v in enumerate(vertices):
            if i == 0: key = infile
104
            else:      key = '{}_{}'.format(infile, i)
105
            self.addVertices(v, key, select=(i == 0), fixWinding=fixWinding)
106
        self.setMeta(infile, surfimg)
107

108
109
110
        if vdata is not None:
            self.addVertexData(infile, vdata)

111
112
113
114
        # Find and load all other
        # surfaces in the same directory
        # as the specfiied one.
        if loadAll:
115

116
117
118
119
            # Only attempt to auto-load sensibly
            # named gifti files (i.e. *.surf.gii,
            # rather than *.gii).
            surfFiles = relatedFiles(infile, [ALLOWED_EXTENSIONS[0]])
120
            nvertices = vertices[0].shape[0]
121

122
123
            for sfile in surfFiles:

124
125
126
127
                try:
                    surfimg, _, vertices, _ = loadGiftiMesh(sfile)
                except Exception:
                    continue
128

129
                if vertices[0].shape[0] != nvertices:
130
131
                    continue

132
133
                self.addVertices(vertices[0], sfile, select=False)
                self.setMeta(sfile, surfimg)
134
135


136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
    def loadVertices(self, infile, key=None, *args, **kwargs):
        """Overrides the :meth:`.Mesh.loadVertices` method.

        Attempts to load vertices for this ``GiftiMesh`` from the given
        ``infile``, which may be a GIFTI file or a plain text file containing
        vertices.
        """

        if not infile.endswith('.gii'):
            return fslmesh.Mesh.loadVertices(
                self, infile, key, *args, **kwargs)

        infile = op.abspath(infile)

        if key is None:
            key = infile

153
        surfimg, _, vertices, _ = loadGiftiMesh(infile)
154

155
156
157
158
        for i, v in enumerate(vertices):
            if i == 0: key = infile
            else:      key = '{}_{}'.format(infile, i)
            vertices[i] = self.addVertices(v, key, *args, **kwargs)
159
160

        self.setMeta(infile, surfimg)
161
162
163
164

        return vertices


165
    def loadVertexData(self, infile, key=None):
166
        """Overrides the :meth:`.Mesh.loadVertexData` method.
167
168

        Attempts to load data associated with each vertex of this
169
170
        ``GiftiMesh`` from the given ``infile``, which may be a GIFTI file or
        a plain text file which contains vertex data.
171
172
        """

173
174
175
        if not infile.endswith('.gii'):
            return fslmesh.Mesh.loadVertexData(self, infile)

176
        infile = op.abspath(infile)
177

178
179
        if key is None:
            key = infile
180

181
182
        vdata = loadGiftiVertexData(infile)[1]
        return self.addVertexData(key, vdata)
183
184


185
def loadGiftiMesh(filename):
186
187
188
189
190
191
    """Extracts surface data from the given ``nibabel.gifti.GiftiImage``.

    The image is expected to contain the following``<DataArray>`` elements:

      - one comprising ``NIFTI_INTENT_TRIANGLE`` data (vertex indices
        defining the triangles).
192
193
      - one or more comprising ``NIFTI_INTENT_POINTSET`` data (the surface
        vertices)
194
195
196

    A ``ValueError`` will be raised if this is not the case.

197
    :arg filename: Name of a GIFTI file containing surface data.
198
199
200

    :returns:     A tuple containing these values:

201
202
                   - The loaded ``nibabel.gifti.GiftiImage`` instance

203
                   - A ``(M, 3)`` array containing the vertex indices for
204
                     ``M`` triangles.
205

206
207
                   - A list of at least one ``(N, 3)`` arrays containing ``N``
                     vertices.
208

209
210
211
212
                   - A ``(M, N)`` numpy array containing ``N`` data points for
                     ``M`` vertices, or ``None`` if the file does not contain
                     any vertex data.
    """
213

214
    gimg = nib.load(filename)
215

216
217
    pscode  = constants.NIFTI_INTENT_POINTSET
    tricode = constants.NIFTI_INTENT_TRIANGLE
218

219
220
221
    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)]
222

223
    if len(triangles) != 1:
224
225
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one triangle array'.format(filename))
226

227
228
229
230
    if len(pointsets) == 0:
        raise ValueError('{}: GIFTI surface files must contain '
                         'at least one pointset array'.format(filename))

231
232
    vertices = [ps.data for ps in pointsets]
    indices  = triangles[0].data
233

234
235
236
237
    if len(vdata) == 0: vdata = None
    else:               vdata = prepareGiftiVertexData(vdata, filename)

    return gimg, indices, vertices, vdata
238
239
240
241
242


def loadGiftiVertexData(filename):
    """Loads vertex data from the given GIFTI file.

243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
    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.

259
260
261
    All of the data arrays are concatenated into one ``(M, N)`` array,
    containing ``N`` data points for ``M`` vertices.

262
263
264
    It is assumed that the given file does not contain any
    ``NIFTI_INTENT_POINTSET`` or ``NIFTI_INTENT_TRIANGLE`` data arrays, and
    which contains either:
265

266
267
268
269
270
271
      - One ``(M, N)`` data array containing ``N`` data points for ``M``
        vertices

      - One or more ``(M, 1)`` data arrays each containing a single data point
        for ``M`` vertices, and all with the same intent code

272
273
    Returns a ``(M, N)`` numpy array containing ``N`` data points for ``M``
    vertices.
274
275
    """

276
    intents = {d.intent for d in darrays}
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291

    if len(intents) != 1:
        raise ValueError('{} contains multiple (or no) intents'
                         ': {}'.format(filename, intents))

    intent = intents.pop()

    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.
292
293
294
    if len(darrays) == 1:
        vdata = darrays[0].data
        return vdata.reshape(vdata.shape[0], -1)
295
296
297

    # Otherwise extract and concatenate
    # multiple 1-dimensional arrays
298
    vdata = [d.data for d in darrays]
299
300
301
302

    if any([len(d.shape) != 1 for d in vdata]):
        raise ValueError('{} contains one or more non-vector '
                         'darrays'.format(filename))
303

304
305
306
    vdata = np.vstack(vdata).T
    vdata = vdata.reshape(vdata.shape[0], -1)

307
    return vdata
308
309


310
def relatedFiles(fname, ftypes=None):
311
312
313
    """Given a GIFTI file, returns a list of other GIFTI files in the same
    directory which appear to be related with the given one.  Files which
    share the same prefix are assumed to be related to the given file.
314

315
316
    This function assumes that the GIFTI files are named according to a
    standard convention - the following conventions are supported:
Paul McCarthy's avatar
Paul McCarthy committed
317

318
319
320
321
322
323
324
     - HCP-style, i.e.: ``<subject>.<hemi>.<type>.<space>.<ftype>.gii``
     - BIDS-style, i.e.:
       ``<source_prefix>_hemi-<hemi>[_space-<space>]*_<suffix>.<ftype>.gii``

    If the files are not named according to one of these conventions, this
    function will return an empty list.

325
326
    :arg fname: Name of the file to search for related files for

327
    :arg ftype: If provided, only files with suffixes in this list are
328
                searched for. Defaults to :attr:`VERTEX_DATA_EXTENSIONS`.
329
330
    """

331
    if ftypes is None:
332
        ftypes = VERTEX_DATA_EXTENSIONS
333

334
335
    path           = op.abspath(fname)
    dirname, fname = op.split(path)
Paul McCarthy's avatar
Paul McCarthy committed
336

337
338
339
340
341
    # We want to identify all files in the same
    # directory which are associated with the
    # given file. We assume that the files are
    # named according to one of the following
    # conventions:
342
    #
343
344
    #  - HCP style:
    #      <subject>.<hemi>.<type>.<space>.<ftype>.gii
Paul McCarthy's avatar
Paul McCarthy committed
345
    #
346
347
    #  - BIDS style:
    #      <source_prefix>_hemi-<hemi>[_space-<space>]*.<ftype>.gii
Paul McCarthy's avatar
Paul McCarthy committed
348
    #
349
350
351
352
353
    # We cannot assume consistent ordering of
    # the entities (key-value pairs) within a
    # BIDS style filename, so we cannot simply
    # use a regular expression or glob pattern.
    # Instead, for each style we define:
Paul McCarthy's avatar
Paul McCarthy committed
354
    #
355
356
357
358
    #  - a "matcher" function, which tests
    #    whether the file matches the style,
    #    and returns the important elements
    #    from the file name.
359
    #
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
    #  - a "searcher" function, which takes
    #    the elements of the input file
    #    that were extracted by the matcher,
    #    and searches for other related files

    # HCP style - extract "<subject>.<hemi>"
    # and "<space>".
    def matchhcp(f):
        pat   = r'^(.*\.[LR])\..*\.(.*)\..*\.gii$'
        match = re.match(pat, f)
        if match:
            return match.groups()
        else:
            return None

    def searchhcp(match, ftype):
        prefix, space = match
        template      = '{}.*.{}{}'.format(prefix, space, ftype)
        return glob.glob(op.join(dirname, template))

    # BIDS style - extract all entities (kv
    # pairs), ignoring specific irrelevant
    # ones.
    def matchbids(f):
        try:               match = bids.BIDSFile(f)
        except ValueError: return None
        match.entities.pop('desc', None)
        return match

    def searchbids(match, ftype):
        allfiles = glob.glob(op.join(dirname, '*{}'.format(ftype)))
        for f in allfiles:
            try:               bf = bids.BIDSFile(f)
            except ValueError: continue
            if bf.match(match, False):
                yield f

    # find the first style that matches
    matchers  = [matchhcp,  matchbids]
    searchers = [searchhcp, searchbids]
    for matcher, searcher in zip(matchers, searchers):
        match = matcher(fname)
        if match:
            break

    # Give up if the file does
    # not match any known style.
    else:
408
        return []
Paul McCarthy's avatar
Paul McCarthy committed
409

410
411
    # Build a list of files in the same
    # directory and matching the template
412
413
    related = []
    for ftype in ftypes:
414
415
416
417

        hits = searcher(match, ftype)

        # eliminate dupes
418
        related.extend([h for h in hits if h not in related])
419

420
    # exclude the file itself
421
    return [r for r in related if r != path]