1
from __future__ import with_statement
7
"""Read a 3-byte int from an open binary file object
19
b1, b2, b3 = np.fromfile(fobj, ">u1", 3)
20
return (b1 << 16) + (b2 << 8) + b3
23
def _fread3_many(fobj, n):
24
"""Read 3-byte ints from an open binary file object.
34
An array of 3 byte int
36
b1, b2, b3 = np.fromfile(fobj, ">u1", 3 * n).reshape(-1,
38
return (b1 << 16) + (b2 << 8) + b3
41
def read_geometry(filepath):
42
"""Read a triangular format Freesurfer surface mesh.
52
nvtx x 3 array of vertex (x, y, z) coordinates
54
nfaces x 3 array of defining mesh triangles
56
with open(filepath, "rb") as fobj:
58
if magic == 16777215: # Quad file
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)
66
# Face splitting follows
68
faces = np.zeros((2 * nquad, 3), dtype=np.int)
71
if (quad[0] % 2) == 0:
72
faces[nface] = quad[0], quad[1], quad[3]
74
faces[nface] = quad[2], quad[3], quad[1]
77
faces[nface] = quad[0], quad[1], quad[2]
79
faces[nface] = quad[0], quad[2], quad[3]
82
elif magic == 16777214: # Triangle file
83
create_stamp = 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)
90
raise ValueError("File does not appear to be a Freesurfer surface")
92
coords = coords.astype(np.float) # XXX: due to mayavi bug on mac 32bits
96
def read_morph_data(filepath):
97
"""Read a Freesurfer morphometry data file.
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.
107
Path to morphometry file
112
Vector representation of surface morpometry values
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)
123
curv = np.fromfile(fobj, ">i2", vnum) / 100
127
def read_annot(filepath, orig_ids=False):
128
"""Read in a Freesurfer annotation from a .annot file.
133
Path to annotation file
135
Whether to return the vertex ids as stored in the annotation
136
file or the positional colortable ids
140
labels : n_vtx numpy array
141
Annotation id at each vertex
143
RGBA + label id colortable array
146
with open(filepath, "rb") as fobj:
148
vnum = np.fromfile(fobj, dt, 1)[0]
149
data = np.fromfile(fobj, dt, vnum * 2).reshape(vnum, 2)
151
ctab_exists = np.fromfile(fobj, dt, 1)[0]
153
raise Exception('Color table not found in annotation file')
154
n_entries = np.fromfile(fobj, dt, 1)[0]
156
length = np.fromfile(fobj, dt, 1)[0]
157
orig_tab = np.fromfile(fobj, '>c', length)
158
orig_tab = orig_tab[:-1]
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]
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))
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]
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]
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))
190
ord = np.argsort(ctab[:, -1])
191
labels = ord[np.searchsorted(ctab[ord, -1], labels)]
192
return labels, ctab, names
195
def read_label(filepath):
196
"""Load in a Freesurfer .label file.
205
label_array : numpy array
206
Array with indices of vertices included in label
209
label_array = np.loadtxt(filepath, dtype=np.int, skiprows=2, usecols=[0])