gifti.py 10.7 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 = ['GIFTI 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
64
65
66
67
        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
68
                         definition.
69
70
71
72
73
74
75

        :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.
76
77
78

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

81
82
83
        name   = fslpath.removeExt(op.basename(infile), ALLOWED_EXTENSIONS)
        infile = op.abspath(infile)

84
        surfimg, indices, vertices, vdata = loadGiftiMesh(infile)
85

86
87
88
89
        fslmesh.Mesh.__init__(self,
                              indices,
                              name=name,
                              dataSource=infile)
90

91
92
        for i, v in enumerate(vertices):
            if i == 0: key = infile
93
            else:      key = '{}_{}'.format(infile, i)
94
            self.addVertices(v, key, select=(i == 0), fixWinding=fixWinding)
95
        self.setMeta(infile, surfimg)
96

97
98
99
        if vdata is not None:
            self.addVertexData(infile, vdata)

100
101
102
103
        # Find and load all other
        # surfaces in the same directory
        # as the specfiied one.
        if loadAll:
104

105
106
107
108
            # Only attempt to auto-load sensibly
            # named gifti files (i.e. *.surf.gii,
            # rather than *.gii).
            surfFiles = relatedFiles(infile, [ALLOWED_EXTENSIONS[0]])
109
            nvertices = vertices[0].shape[0]
110

111
112
            for sfile in surfFiles:

113
114
115
116
                try:
                    surfimg, _, vertices, _ = loadGiftiMesh(sfile)
                except Exception:
                    continue
117

118
                if vertices[0].shape[0] != nvertices:
119
120
                    continue

121
122
                self.addVertices(vertices[0], sfile, select=False)
                self.setMeta(sfile, surfimg)
123
124


125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    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

142
        surfimg, _, vertices, _ = loadGiftiMesh(infile)
143

144
145
146
147
        for i, v in enumerate(vertices):
            if i == 0: key = infile
            else:      key = '{}_{}'.format(infile, i)
            vertices[i] = self.addVertices(v, key, *args, **kwargs)
148
149

        self.setMeta(infile, surfimg)
150
151
152
153

        return vertices


154
    def loadVertexData(self, infile, key=None):
155
        """Overrides the :meth:`.Mesh.loadVertexData` method.
156
157

        Attempts to load data associated with each vertex of this
158
159
        ``GiftiMesh`` from the given ``infile``, which may be a GIFTI file or
        a plain text file which contains vertex data.
160
161
        """

162
163
164
        if not infile.endswith('.gii'):
            return fslmesh.Mesh.loadVertexData(self, infile)

165
        infile = op.abspath(infile)
166

167
168
        if key is None:
            key = infile
169

170
171
        vdata = loadGiftiVertexData(infile)[1]
        return self.addVertexData(key, vdata)
172
173


174
def loadGiftiMesh(filename):
175
176
177
178
179
180
    """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).
181
182
      - one or more comprising ``NIFTI_INTENT_POINTSET`` data (the surface
        vertices)
183
184
185

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

186
    :arg filename: Name of a GIFTI file containing surface data.
187
188
189

    :returns:     A tuple containing these values:

190
191
                   - The loaded ``nibabel.gifti.GiftiImage`` instance

192
                   - A ``(M, 3)`` array containing the vertex indices for
193
                     ``M`` triangles.
194

195
196
                   - A list of at least one ``(N, 3)`` arrays containing ``N``
                     vertices.
197

198
199
200
201
                   - A ``(M, N)`` numpy array containing ``N`` data points for
                     ``M`` vertices, or ``None`` if the file does not contain
                     any vertex data.
    """
202

203
    gimg = nib.load(filename)
204

205
206
    pscode  = constants.NIFTI_INTENT_POINTSET
    tricode = constants.NIFTI_INTENT_TRIANGLE
207

208
209
210
    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)]
211

212
    if len(triangles) != 1:
213
214
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one triangle array'.format(filename))
215

216
217
218
219
    if len(pointsets) == 0:
        raise ValueError('{}: GIFTI surface files must contain '
                         'at least one pointset array'.format(filename))

220
221
    vertices = [ps.data for ps in pointsets]
    indices  = triangles[0].data
222

223
224
225
226
    if len(vdata) == 0: vdata = None
    else:               vdata = prepareGiftiVertexData(vdata, filename)

    return gimg, indices, vertices, vdata
227
228
229
230
231


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

232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
    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.

248
249
250
    All of the data arrays are concatenated into one ``(M, N)`` array,
    containing ``N`` data points for ``M`` vertices.

251
252
253
    It is assumed that the given file does not contain any
    ``NIFTI_INTENT_POINTSET`` or ``NIFTI_INTENT_TRIANGLE`` data arrays, and
    which contains either:
254

255
256
257
258
259
260
      - 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

261
262
    Returns a ``(M, N)`` numpy array containing ``N`` data points for ``M``
    vertices.
263
264
    """

265
    intents = set([d.intent for d in darrays])
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280

    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.
281
282
283
    if len(darrays) == 1:
        vdata = darrays[0].data
        return vdata.reshape(vdata.shape[0], -1)
284
285
286

    # Otherwise extract and concatenate
    # multiple 1-dimensional arrays
287
    vdata = [d.data for d in darrays]
288
289
290
291

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

293
294
295
    vdata = np.vstack(vdata).T
    vdata = vdata.reshape(vdata.shape[0], -1)

296
    return vdata
297
298


299
def relatedFiles(fname, ftypes=None):
300
301
302
    """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.
303
304
305

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

306
307
    :arg ftype: If provided, only files with suffixes in this list are
                searched for. Defaults to files which contain vertex data.
308
309
    """

310
311
312
    if ftypes is None:
        ftypes = ['.func.gii', '.shape.gii', '.label.gii', '.time.gii']

Paul McCarthy's avatar
Paul McCarthy committed
313
314
315
    # We want to return all files in the same
    # directory which have the following name:

316
    #
317
    # [subj].[hemi].[type].*.[ftype]
Paul McCarthy's avatar
Paul McCarthy committed
318
    #
319

Paul McCarthy's avatar
Paul McCarthy committed
320
    #   where
321
    #     - [subj] is the subject ID, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
322
    #
323
    #     - [hemi] is the hemisphere, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
324
    #
325
    #     - [type] defines the file contents
326
    #
327
    #     - suffix is func, shape, label, time, or `ftype`
Paul McCarthy's avatar
Paul McCarthy committed
328

329
330
    path            = op.abspath(fname)
    dirname, fname  = op.split(path)
Paul McCarthy's avatar
Paul McCarthy committed
331

332
333
334
335
336
    # get the [subj].[hemi] prefix
    try:
        subj, hemi, _ = fname.split('.', 2)
        prefix        = '.'.join((subj, hemi))
    except Exception:
337
        return []
Paul McCarthy's avatar
Paul McCarthy committed
338

339
340
341
    related = []

    for ftype in ftypes:
342
343
        hits = glob.glob(op.join(dirname, '{}*{}'.format(prefix, ftype)))
        related.extend([h for h in hits if h not in related])
344

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