~ubuntu-branches/ubuntu/wily/nibabel/wily-proposed

« back to all changes in this revision

Viewing changes to nibabel/freesurfer/io.py

  • Committer: Package Import Robot
  • Author(s): Yaroslav Halchenko
  • Date: 2012-05-06 12:58:22 UTC
  • mfrom: (1.1.3)
  • Revision ID: package-import@ubuntu.com-20120506125822-3jiwjkmdqcxkrior
Tags: 1.2.0-1
* New upstream release: bugfixes, support for new formats
  (Freesurfer, ECAT, etc), nib-ls tool.
* debian/copyright
  - corrected reference to the format specs
  - points to github's repository as the origin of sources
  - removed lightunit license/copyright -- removed upstream
  - added netcdf module license/copyright terms
* debian/control
  - added python-fuse to Recommends and Build-Depends for dicomfs
  - boosted policy compliance to 3.9.3 (no further changes)
* debian/watch
  - adjusted to download numbered tag

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
from __future__ import with_statement
 
2
 
 
3
import numpy as np
 
4
 
 
5
 
 
6
def _fread3(fobj):
 
7
    """Read a 3-byte int from an open binary file object
 
8
 
 
9
    Parameters
 
10
    ----------
 
11
    fobj : file
 
12
        File descriptor
 
13
 
 
14
    Returns
 
15
    -------
 
16
    n : int
 
17
        A 3 byte int
 
18
    """
 
19
    b1, b2, b3 = np.fromfile(fobj, ">u1", 3)
 
20
    return (b1 << 16) + (b2 << 8) + b3
 
21
 
 
22
 
 
23
def _fread3_many(fobj, n):
 
24
    """Read 3-byte ints from an open binary file object.
 
25
 
 
26
    Parameters
 
27
    ----------
 
28
    fobj : file
 
29
        File descriptor
 
30
 
 
31
    Returns
 
32
    -------
 
33
    out : 1D array
 
34
        An array of 3 byte int
 
35
    """
 
36
    b1, b2, b3 = np.fromfile(fobj, ">u1", 3 * n).reshape(-1,
 
37
                                                    3).astype(np.int).T
 
38
    return (b1 << 16) + (b2 << 8) + b3
 
39
 
 
40
 
 
41
def read_geometry(filepath):
 
42
    """Read a triangular format Freesurfer surface mesh.
 
43
 
 
44
    Parameters
 
45
    ----------
 
46
    filepath : str
 
47
        Path to surface file
 
48
 
 
49
    Returns
 
50
    -------
 
51
    coords : numpy array
 
52
        nvtx x 3 array of vertex (x, y, z) coordinates
 
53
    faces : numpy array
 
54
        nfaces x 3 array of defining mesh triangles
 
55
    """
 
56
    with open(filepath, "rb") as fobj:
 
57
        magic = _fread3(fobj)
 
58
        if magic == 16777215:  # Quad file
 
59
            nvert = _fread3(fobj)
 
60
            nquad = _fread3(fobj)
 
61
            coords = np.fromfile(fobj, ">i2", nvert * 3).astype(np.float)
 
62
            coords = coords.reshape(-1, 3) / 100.0
 
63
            quads = _fread3_many(fobj, nquad * 4)
 
64
            quads = quads.reshape(nquad, 4)
 
65
            #
 
66
            #   Face splitting follows
 
67
            #
 
68
            faces = np.zeros((2 * nquad, 3), dtype=np.int)
 
69
            nface = 0
 
70
            for quad in quads:
 
71
                if (quad[0] % 2) == 0:
 
72
                    faces[nface] = quad[0], quad[1], quad[3]
 
73
                    nface += 1
 
74
                    faces[nface] = quad[2], quad[3], quad[1]
 
75
                    nface += 1
 
76
                else:
 
77
                    faces[nface] = quad[0], quad[1], quad[2]
 
78
                    nface += 1
 
79
                    faces[nface] = quad[0], quad[2], quad[3]
 
80
                    nface += 1
 
81
 
 
82
        elif magic == 16777214:  # Triangle file
 
83
            create_stamp = fobj.readline()
 
84
            _ = fobj.readline()
 
85
            vnum = np.fromfile(fobj, ">i4", 1)[0]
 
86
            fnum = np.fromfile(fobj, ">i4", 1)[0]
 
87
            coords = np.fromfile(fobj, ">f4", vnum * 3).reshape(vnum, 3)
 
88
            faces = np.fromfile(fobj, ">i4", fnum * 3).reshape(fnum, 3)
 
89
        else:
 
90
            raise ValueError("File does not appear to be a Freesurfer surface")
 
91
 
 
92
    coords = coords.astype(np.float)  # XXX: due to mayavi bug on mac 32bits
 
93
    return coords, faces
 
94
 
 
95
 
 
96
def read_morph_data(filepath):
 
97
    """Read a Freesurfer morphometry data file.
 
98
 
 
99
    This function reads in what Freesurfer internally calls "curv" file types,
 
100
    (e.g. ?h. curv, ?h.thickness), but as that has the potential to cause
 
101
    confusion where "curv" also refers to the surface curvature values,
 
102
    we refer to these files as "morphometry" files with PySurfer.
 
103
 
 
104
    Parameters
 
105
    ----------
 
106
    filepath : str
 
107
        Path to morphometry file
 
108
 
 
109
    Returns
 
110
    -------
 
111
    curv : numpy array
 
112
        Vector representation of surface morpometry values
 
113
 
 
114
    """
 
115
    with open(filepath, "rb") as fobj:
 
116
        magic = _fread3(fobj)
 
117
        if magic == 16777215:
 
118
            vnum = np.fromfile(fobj, ">i4", 3)[0]
 
119
            curv = np.fromfile(fobj, ">f4", vnum)
 
120
        else:
 
121
            vnum = magic
 
122
            _ = _fread3(fobj)
 
123
            curv = np.fromfile(fobj, ">i2", vnum) / 100
 
124
    return curv
 
125
 
 
126
 
 
127
def read_annot(filepath, orig_ids=False):
 
128
    """Read in a Freesurfer annotation from a .annot file.
 
129
 
 
130
    Parameters
 
131
    ----------
 
132
    filepath : str
 
133
        Path to annotation file
 
134
    orig_ids : bool
 
135
        Whether to return the vertex ids as stored in the annotation
 
136
        file or the positional colortable ids
 
137
 
 
138
    Returns
 
139
    -------
 
140
    labels : n_vtx numpy array
 
141
        Annotation id at each vertex
 
142
    ctab : numpy array
 
143
        RGBA + label id colortable array
 
144
 
 
145
    """
 
146
    with open(filepath, "rb") as fobj:
 
147
        dt = ">i4"
 
148
        vnum = np.fromfile(fobj, dt, 1)[0]
 
149
        data = np.fromfile(fobj, dt, vnum * 2).reshape(vnum, 2)
 
150
        labels = data[:, 1]
 
151
        ctab_exists = np.fromfile(fobj, dt, 1)[0]
 
152
        if not ctab_exists:
 
153
            raise Exception('Color table not found in annotation file')
 
154
        n_entries = np.fromfile(fobj, dt, 1)[0]
 
155
        if n_entries > 0:
 
156
            length = np.fromfile(fobj, dt, 1)[0]
 
157
            orig_tab = np.fromfile(fobj, '>c', length)
 
158
            orig_tab = orig_tab[:-1]
 
159
 
 
160
            names = list()
 
161
            ctab = np.zeros((n_entries, 5), np.int)
 
162
            for i in xrange(n_entries):
 
163
                name_length = np.fromfile(fobj, dt, 1)[0]
 
164
                name = np.fromfile(fobj, "|S%d" % name_length, 1)[0]
 
165
                names.append(name)
 
166
                ctab[i, :4] = np.fromfile(fobj, dt, 4)
 
167
                ctab[i, 4] = (ctab[i, 0] + ctab[i, 1] * (2 ** 8) +
 
168
                              ctab[i, 2] * (2 ** 16) +
 
169
                              ctab[i, 3] * (2 ** 24))
 
170
        else:
 
171
            ctab_version = -n_entries
 
172
            if ctab_version != 2:
 
173
                raise Exception('Color table version not supported')
 
174
            n_entries = np.fromfile(fobj, dt, 1)[0]
 
175
            ctab = np.zeros((n_entries, 5), np.int)
 
176
            length = np.fromfile(fobj, dt, 1)[0]
 
177
            _ = np.fromfile(fobj, "|S%d" % length, 1)[0]  # Orig table path
 
178
            entries_to_read = np.fromfile(fobj, dt, 1)[0]
 
179
            names = list()
 
180
            for i in xrange(entries_to_read):
 
181
                _ = np.fromfile(fobj, dt, 1)[0]  # Structure
 
182
                name_length = np.fromfile(fobj, dt, 1)[0]
 
183
                name = np.fromfile(fobj, "|S%d" % name_length, 1)[0]
 
184
                names.append(name)
 
185
                ctab[i, :4] = np.fromfile(fobj, dt, 4)
 
186
                ctab[i, 4] = (ctab[i, 0] + ctab[i, 1] * (2 ** 8) +
 
187
                                ctab[i, 2] * (2 ** 16))
 
188
        ctab[:, 3] = 255
 
189
    if not orig_ids:
 
190
        ord = np.argsort(ctab[:, -1])
 
191
        labels = ord[np.searchsorted(ctab[ord, -1], labels)]
 
192
    return labels, ctab, names
 
193
 
 
194
 
 
195
def read_label(filepath):
 
196
    """Load in a Freesurfer .label file.
 
197
 
 
198
    Parameters
 
199
    ----------
 
200
    filepath : str
 
201
        Path to label file
 
202
 
 
203
    Returns
 
204
    -------
 
205
    label_array : numpy array
 
206
        Array with indices of vertices included in label
 
207
 
 
208
    """
 
209
    label_array = np.loadtxt(filepath, dtype=np.int, skiprows=2, usecols=[0])
 
210
    return label_array