gifti.py 10.1 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
29
import os.path as op

30
import numpy   as np
31
32
import nibabel as nib

33
34
35
import fsl.utils.path     as fslpath
import fsl.data.constants as constants
import fsl.data.mesh      as fslmesh
36
37


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


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


class GiftiMesh(fslmesh.Mesh):
49
50
51
    """Class which represents a GIFTI surface image. This is essentially
    just a 3D model made of triangles.

52
53
54
    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.
55
56
57
    """


58
    def __init__(self, infile, fixWinding=False, loadAll=False):
59
        """Load the given GIFTI file using ``nibabel``, and extracts surface
60
        data using the  :func:`loadGiftiMesh` function.
61

62
63
        :arg infile:     A GIFTI file (``*..gii``) which contains a surface
                         definition.
64
65
66
67
68
69
70

        :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.
71
72
73

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

76
77
78
        name   = fslpath.removeExt(op.basename(infile), ALLOWED_EXTENSIONS)
        infile = op.abspath(infile)

79
        surfimg, indices, vertices, vdata = loadGiftiMesh(infile)
80

81
82
83
84
        fslmesh.Mesh.__init__(self,
                              indices,
                              name=name,
                              dataSource=infile)
85

86
87
88
89
        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)
90
        self.setMeta(infile, surfimg)
91

92
93
94
        if vdata is not None:
            self.addVertexData(infile, vdata)

95
96
97
98
        # Find and load all other
        # surfaces in the same directory
        # as the specfiied one.
        if loadAll:
99

100
            nvertices = vertices[0].shape[0]
101
            surfFiles = relatedFiles(infile, ALLOWED_EXTENSIONS)
102

103
104
            for sfile in surfFiles:

105
106
107
108
                try:
                    surfimg, _, vertices, _ = loadGiftiMesh(sfile)
                except Exception:
                    continue
109

110
                if vertices[0].shape[0] != nvertices:
111
112
                    continue

113
114
                self.addVertices(vertices[0], sfile, select=False)
                self.setMeta(sfile, surfimg)
115
116


117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
    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

134
        surfimg, _, vertices, _ = loadGiftiMesh(infile)
135

136
137
138
        vertices = self.addVertices(vertices, key, *args, **kwargs)

        self.setMeta(infile, surfimg)
139
140
141
142

        return vertices


143
    def loadVertexData(self, infile, key=None):
144
        """Overrides the :meth:`.Mesh.loadVertexData` method.
145
146

        Attempts to load data associated with each vertex of this
147
148
        ``GiftiMesh`` from the given ``infile``, which may be a GIFTI file or
        a plain text file which contains vertex data.
149
150
        """

151
152
153
        if not infile.endswith('.gii'):
            return fslmesh.Mesh.loadVertexData(self, infile)

154
        infile = op.abspath(infile)
155

156
157
        if key is None:
            key = infile
158

159
160
        vdata = loadGiftiVertexData(infile)[1]
        return self.addVertexData(key, vdata)
161
162


163
def loadGiftiMesh(filename):
164
165
166
167
168
169
    """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).
170
171
      - one or more comprising ``NIFTI_INTENT_POINTSET`` data (the surface
        vertices)
172
173
174

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

175
    :arg filename: Name of a GIFTI file containing surface data.
176
177
178

    :returns:     A tuple containing these values:

179
180
                   - The loaded ``nibabel.gifti.GiftiImage`` instance

181
                   - A ``(M, 3)`` array containing the vertex indices for
182
                     ``M`` triangles.
183

184
185
                   - A list of at least one ``(N, 3)`` arrays containing ``N``
                     vertices.
186

187
188
189
190
                   - A ``(M, N)`` numpy array containing ``N`` data points for
                     ``M`` vertices, or ``None`` if the file does not contain
                     any vertex data.
    """
191

192
    gimg = nib.load(filename)
193

194
195
    pscode  = constants.NIFTI_INTENT_POINTSET
    tricode = constants.NIFTI_INTENT_TRIANGLE
196

197
198
199
    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)]
200

201
    if len(triangles) != 1:
202
203
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one triangle array'.format(filename))
204

205
206
207
208
209
    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]
210
211
    indices  = triangles[0].data

212
213
214
215
    if len(vdata) == 0: vdata = None
    else:               vdata = prepareGiftiVertexData(vdata, filename)

    return gimg, indices, vertices, vdata
216
217
218
219
220


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

221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    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.

237
238
239
    All of the data arrays are concatenated into one ``(M, N)`` array,
    containing ``N`` data points for ``M`` vertices.

240
241
242
    It is assumed that the given file does not contain any
    ``NIFTI_INTENT_POINTSET`` or ``NIFTI_INTENT_TRIANGLE`` data arrays, and
    which contains either:
243

244
245
246
247
248
249
      - 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

250
251
    Returns a ``(M, N)`` numpy array containing ``N`` data points for ``M``
    vertices.
252
253
    """

254
    intents = set([d.intent for d in darrays])
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269

    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.
270
271
272
    if len(darrays) == 1:
        vdata = darrays[0].data
        return vdata.reshape(vdata.shape[0], -1)
273
274
275

    # Otherwise extract and concatenate
    # multiple 1-dimensional arrays
276
    vdata = [d.data for d in darrays]
277
278
279
280

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

282
283
284
    vdata = np.vstack(vdata).T
    vdata = vdata.reshape(vdata.shape[0], -1)

285
    return vdata
286
287


288
def relatedFiles(fname, ftypes=None):
289
290
291
    """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.
292
293
294

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

295
296
    :arg ftype: If provided, only files with suffixes in this list are
                searched for. Defaults to files which contain vertex data.
297
298
    """

299
300
301
    if ftypes is None:
        ftypes = ['.func.gii', '.shape.gii', '.label.gii', '.time.gii']

Paul McCarthy's avatar
Paul McCarthy committed
302
303
304
    # We want to return all files in the same
    # directory which have the following name:

305
    #
306
    # [subj].[hemi].[type].*.[ftype]
Paul McCarthy's avatar
Paul McCarthy committed
307
    #
308

Paul McCarthy's avatar
Paul McCarthy committed
309
    #   where
310
    #     - [subj] is the subject ID, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
311
    #
312
    #     - [hemi] is the hemisphere, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
313
    #
314
    #     - [type] defines the file contents
315
    #
316
    #     - suffix is func, shape, label, time, or `ftype`
Paul McCarthy's avatar
Paul McCarthy committed
317

318
319
    path            = op.abspath(fname)
    dirname, fname  = op.split(path)
Paul McCarthy's avatar
Paul McCarthy committed
320

321
322
323
324
325
    # get the [subj].[hemi] prefix
    try:
        subj, hemi, _ = fname.split('.', 2)
        prefix        = '.'.join((subj, hemi))
    except Exception:
326
        return []
Paul McCarthy's avatar
Paul McCarthy committed
327

328
329
330
    related = []

    for ftype in ftypes:
331
332
        related.extend(
            glob.glob(op.join(dirname, '{}*{}'.format(prefix, ftype))))
333

334
    return [r for r in related if r != path]