gifti.py 9.16 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
     relatedFiles
23
24
25
"""


26
import            glob
27
28
import os.path as op

29
import numpy   as np
30
31
import nibabel as nib

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


37
38
39
# We include '.gii' here because not all surface
# GIFTIs follow the file suffix convention.
ALLOWED_EXTENSIONS = ['.surf.gii', '.gii']
40
41
42
43
44
"""List of file extensions that a file containing Gifti surface data
is expected to have.
"""


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


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

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


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

63
64
65
66
67
68
69
70
        :arg infile:     A GIFTI surface file (``*.surf.gii``).

        :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, vertices, indices = loadGiftiMesh(infile)
80

81

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

87
88
        self.addVertices(vertices, infile, fixWinding=fixWinding)
        self.setMeta(infile, surfimg)
89

90
91
92
93
        # Find and load all other
        # surfaces in the same directory
        # as the specfiied one.
        if loadAll:
94

95
            nvertices = vertices.shape[0]
96
            surfFiles = relatedFiles(infile, ALLOWED_EXTENSIONS)
97

98
99
            for sfile in surfFiles:

100
                surfimg, vertices, _ = loadGiftiMesh(sfile)
101
102
103
104
105
106
107
108

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

                self.addVertices(vertices, sfile, select=False)
                self.setMeta(    sfile, surfimg)


109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
    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

        surfimg, vertices, _ = loadGiftiMesh(infile)

128
129
130
        vertices = self.addVertices(vertices, key, *args, **kwargs)

        self.setMeta(infile, surfimg)
131
132
133
134

        return vertices


135
    def loadVertexData(self, infile, key=None):
136
        """Overrides the :meth:`.Mesh.loadVertexData` method.
137
138

        Attempts to load data associated with each vertex of this
139
140
        ``GiftiMesh`` from the given ``infile``, which may be a GIFTI file or
        a plain text file which contains vertex data.
141
142
        """

143
144
145
        if not infile.endswith('.gii'):
            return fslmesh.Mesh.loadVertexData(self, infile)

146
        infile = op.abspath(infile)
147

148
149
        if key is None:
            key = infile
150

151
152
        vdata = loadGiftiVertexData(infile)[1]
        return self.addVertexData(key, vdata)
153
154


155
def loadGiftiMesh(filename):
156
157
158
159
160
161
162
163
164
165
    """Extracts surface data from the given ``nibabel.gifti.GiftiImage``.

    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).

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

166
    :arg filename: Name of a GIFTI file containing surface data.
167
168
169

    :returns:     A tuple containing these values:

170
171
                   - The loaded ``nibabel.gifti.GiftiImage`` instance

172
                   - A ``(N, 3)`` array containing ``N`` vertices.
173

174
175
                   - A ``(M, 3))`` array containing the vertex indices for
                     ``M`` triangles.
176
177
    """

178
    gimg = nib.load(filename)
179

180
181
182
183
184
185
186
    pointsetCode = constants.NIFTI_INTENT_POINTSET
    triangleCode = constants.NIFTI_INTENT_TRIANGLE

    pointsets = [d for d in gimg.darrays if d.intent == pointsetCode]
    triangles = [d for d in gimg.darrays if d.intent == triangleCode]

    if len(gimg.darrays) != 2:
187
188
189
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one pointset array and one '
                         'triangle array'.format(filename))
190

191
    if len(pointsets) != 1:
192
193
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one pointset array'.format(filename))
194

195
    if len(triangles) != 1:
196
197
        raise ValueError('{}: GIFTI surface files must contain '
                         'exactly one triangle array'.format(filename))
198
199
200
201
202
203
204
205
206
207
208
209
210

    vertices = pointsets[0].data
    indices  = triangles[0].data

    return gimg, vertices, indices


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

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

212
213
214
215
216
217
218
      - 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

    Returns a tuple containing:
219

220
      - The loaded ``nibabel.gifti.GiftiImage`` object
221

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

    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):
238

239
240
241
242
243
244
245
        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:
246
247
        vdata = gimg.darrays[0].data
        return gimg, vdata.reshape(vdata.shape[0], -1)
248
249
250
251
252
253
254
255

    # Otherwise extract and concatenate
    # multiple 1-dimensional arrays
    vdata = [d.data for d in gimg.darrays]

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

257
258
259
260
    vdata = np.vstack(vdata).T
    vdata = vdata.reshape(vdata.shape[0], -1)

    return gimg, vdata
261
262


263
def relatedFiles(fname, ftypes=None):
264
265
266
    """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.
267
268
269

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

270
271
    :arg ftype: If provided, only files with suffixes in this list are
                searched for. Defaults to files which contain vertex data.
272
273
    """

274
275
276
    if ftypes is None:
        ftypes = ['.func.gii', '.shape.gii', '.label.gii', '.time.gii']

Paul McCarthy's avatar
Paul McCarthy committed
277
278
279
    # We want to return all files in the same
    # directory which have the following name:

280
    #
281
    # [subj].[hemi].[type].*.[ftype]
Paul McCarthy's avatar
Paul McCarthy committed
282
    #
283

Paul McCarthy's avatar
Paul McCarthy committed
284
    #   where
285
    #     - [subj] is the subject ID, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
286
    #
287
    #     - [hemi] is the hemisphere, and matches fname
Paul McCarthy's avatar
Paul McCarthy committed
288
    #
289
    #     - [type] defines the file contents
290
    #
291
    #     - suffix is func, shape, label, time, or `ftype`
Paul McCarthy's avatar
Paul McCarthy committed
292

293
294
    path            = op.abspath(fname)
    dirname, fname  = op.split(path)
Paul McCarthy's avatar
Paul McCarthy committed
295

296
297
298
299
300
    # get the [subj].[hemi] prefix
    try:
        subj, hemi, _ = fname.split('.', 2)
        prefix        = '.'.join((subj, hemi))
    except Exception:
301
        return []
Paul McCarthy's avatar
Paul McCarthy committed
302

303
304
305
    related = []

    for ftype in ftypes:
306
307
        related.extend(
            glob.glob(op.join(dirname, '{}*{}'.format(prefix, ftype))))
308

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