-
Notifications
You must be signed in to change notification settings - Fork 1.4k
/
Copy pathtest_freesurfer.py
288 lines (264 loc) · 10.2 KB
/
test_freesurfer.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
from pathlib import Path
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_array_equal
import mne
from mne import (
get_volume_labels_from_aseg,
head_to_mni,
read_freesurfer_lut,
read_talxfm,
vertex_to_mni,
)
from mne._freesurfer import (
_check_subject_dir,
_estimate_talxfm_rigid,
_get_mgz_header,
read_lta,
)
from mne.datasets import testing
from mne.transforms import _angle_between_quats, _get_trans, apply_trans, rot_to_quat
data_path = testing.data_path(download=False)
subjects_dir = data_path / "subjects"
fname_mri = data_path / "subjects" / "sample" / "mri" / "T1.mgz"
aseg_fname = data_path / "subjects" / "sample" / "mri" / "aseg.mgz"
trans_fname = data_path / "MEG" / "sample" / "sample_audvis_trunc-trans.fif"
rng = np.random.RandomState(0)
@testing.requires_testing_data
def test_check_subject_dir():
"""Test checking for a Freesurfer recon-all subject directory."""
_check_subject_dir("sample", subjects_dir)
with pytest.raises(ValueError, match="subject folder is incorrect"):
_check_subject_dir("foo", data_path)
@testing.requires_testing_data
def test_mgz_header():
"""Test MGZ header reading."""
nib = pytest.importorskip("nibabel")
header = _get_mgz_header(fname_mri)
mri_hdr = nib.load(fname_mri).header
assert_allclose(mri_hdr.get_data_shape(), header["dims"])
assert_allclose(mri_hdr.get_vox2ras_tkr(), header["vox2ras_tkr"])
assert_allclose(mri_hdr.get_ras2vox(), np.linalg.inv(header["vox2ras"]))
@testing.requires_testing_data
def test_vertex_to_mni():
"""Test conversion of vertices to MNI coordinates."""
pytest.importorskip("nibabel")
# obtained using "tksurfer (sample) (l/r)h white"
vertices = [100960, 7620, 150549, 96761]
coords = np.array(
[
[-60.86, -11.18, -3.19],
[-36.46, -93.18, -2.36],
[-38.00, 50.08, -10.61],
[47.14, 8.01, 46.93],
]
)
hemis = [0, 0, 0, 1]
coords_2 = vertex_to_mni(vertices, hemis, "sample", subjects_dir)
# less than 1mm error
assert_allclose(coords, coords_2, atol=1.0)
@testing.requires_testing_data
def test_head_to_mni():
"""Test conversion of aseg vertices to MNI coordinates."""
# obtained using freeview
coords = (
np.array(
[
[22.52, 11.24, 17.72],
[22.52, 5.46, 21.58],
[16.10, 5.46, 22.23],
[21.24, 8.36, 22.23],
]
)
/ 1000.0
)
xfm = read_talxfm("sample", subjects_dir)
coords_MNI = apply_trans(xfm["trans"], coords) * 1000.0
mri_head_t, _ = _get_trans(trans_fname, "mri", "head", allow_none=False)
# obtained from sample_audvis-meg-oct-6-mixed-fwd.fif
coo_right_amygdala = np.array(
[
[0.01745682, 0.02665809, 0.03281873],
[0.01014125, 0.02496262, 0.04233755],
[0.01713642, 0.02505193, 0.04258181],
[0.01720631, 0.03073877, 0.03850075],
]
)
coords_MNI_2 = head_to_mni(coo_right_amygdala, "sample", mri_head_t, subjects_dir)
# less than 1mm error
assert_allclose(coords_MNI, coords_MNI_2, atol=10.0)
@testing.requires_testing_data
def test_vertex_to_mni_fs_nibabel(monkeypatch):
"""Test equivalence of vert_to_mni for nibabel and freesurfer."""
pytest.importorskip("nibabel")
n_check = 1000
subject = "sample"
vertices = rng.randint(0, 100000, n_check)
hemis = rng.randint(0, 1, n_check)
coords = vertex_to_mni(vertices, hemis, subject, subjects_dir)
read_mri = mne._freesurfer._read_mri_info
monkeypatch.setattr(
mne._freesurfer,
"_read_mri_info",
lambda *args, **kwargs: read_mri(*args, use_nibabel=True, **kwargs),
)
coords_2 = vertex_to_mni(vertices, hemis, subject, subjects_dir)
# less than 0.1 mm error
assert_allclose(coords, coords_2, atol=0.1)
def test_read_lta(tmp_path):
"""Test reading a Freesurfer linear transform array file."""
with open(tmp_path / "test.lta", "w") as fid:
fid.write(
"""type = 0 # LINEAR_VOX_TO_VOX
nxforms = 1
mean = 0.0000 0.0000 0.0000
sigma = 1.0000
1 4 4
0.99221027 -0.05494503 0.11180324 -3.84350586
0.05233596 0.99828744 0.02614108 -9.77523804
-0.11304809 -0.02008611 0.99338663 15.25457001
0 0 0 1
src volume info
valid = 1 # volume info valid
filename = tmp.mgz
volume = 256 256 256
voxelsize = 1 1 1
xras = -1 0 0
yras = 0 0 -1
zras = 0 1 0
cras = -1.19374 -3.31686 3.25835
dst volume info
valid = 1 # volume info valid
filename = tmp.mgz
volume = 256 256 256
voxelsize = 1 1 1
xras = -1 0 0
yras = 0 0 -1
zras = 0 1 0
cras = -1.19374 -3.31686 3.25835"""
)
assert_array_equal(
read_lta(tmp_path / "test.lta"),
np.array(
[
[0.99221027, -0.05494503, 0.11180324, -3.84350586],
[0.05233596, 0.99828744, 0.02614108, -9.77523804],
[-0.11304809, -0.02008611, 0.99338663, 15.25457001],
[0.0, 0.0, 0.0, 1.0],
]
),
)
# test when dst volume != src_volume
with open(tmp_path / "test2.lta", "w") as fid:
fid.write(
"""type = 0 # LINEAR_VOX_TO_VOX
nxforms = 1
mean = 0.0000 0.0000 0.0000
sigma = 1.0000
1 4 4
0.41397345 -0.02919456 -0.00069703 26.37020874
-0.02894894 -0.40985453 -0.06119149 212.38204956
0.00361269 0.0611503 -0.41046342 203.33338928
0 0 0 1
src volume info
valid = 1 # volume info valid
filename = tmp2.mgz
volume = 512 385 512
voxelsize = 0.41499999 0.41541821 0.41499999
xras = -1 0 0
yras = 0 0 1
zras = 0 -1 0
cras = -106.23999786 105.82500458 -79.55259705
dst volume info
valid = 1 # volume info valid
filename = tmp.mgz
volume = 256 256 256
voxelsize = 1 1 1
xras = -1 0 0
yras = 0 0 -1
zras = 0 1 0
cras = -3.68961334 -0.12011719 3.4160614"""
)
assert_allclose(
read_lta(tmp_path / "test2.lta"),
np.array(
[
[0.99752641, -0.07034834, -0.00167959, -236.00043542],
[0.06968626, 0.98660704, 0.14730093, 189.09766694],
[-0.00870528, -0.14735012, 0.98906851, 329.7632126],
[0.0, 0.0, 0.0, 1.0],
]
),
atol=1e-8,
)
@testing.requires_testing_data
@pytest.mark.parametrize(
"fname",
[
None,
Path(mne.__file__).parent / "data" / "FreeSurferColorLUT.txt",
],
)
def test_read_freesurfer_lut(fname, tmp_path):
"""Test reading volume label names."""
pytest.importorskip("nibabel")
atlas_ids, colors = read_freesurfer_lut(fname)
assert list(atlas_ids).count("Brain-Stem") == 1
assert len(colors) == len(atlas_ids) == 1266
label_names, label_colors = get_volume_labels_from_aseg(
aseg_fname, return_colors=True
)
assert isinstance(label_names, list)
assert isinstance(label_colors, list)
assert label_names.count("Brain-Stem") == 1
for c in label_colors:
assert isinstance(c, np.ndarray)
assert c.shape == (4,)
assert len(label_names) == len(label_colors) == 46
with pytest.raises(ValueError, match="must be False"):
get_volume_labels_from_aseg(aseg_fname, return_colors=True, atlas_ids=atlas_ids)
label_names_2 = get_volume_labels_from_aseg(aseg_fname, atlas_ids=atlas_ids)
assert label_names == label_names_2
# long name (only test on one run)
if fname is not None:
return
fname = tmp_path / "long.txt"
names = [
"Anterior_Cingulate_and_Medial_Prefrontal_Cortex-" + hemi
for hemi in ("lh", "rh")
]
ids = np.arange(1, len(names) + 1)
colors = [(id_,) * 4 for id_ in ids]
with open(fname, "w") as fid:
for name, id_, color in zip(names, ids, colors):
out_color = " ".join(f"{x:03}" for x in color)
line = f"{id_} {name} {out_color}\n"
fid.write(line)
lut, got_colors = read_freesurfer_lut(fname)
assert len(lut) == len(got_colors) == len(names) == len(ids)
for name, id_, color in zip(names, ids, colors):
assert name in lut
assert name in got_colors
assert_array_equal(got_colors[name][:3], color[:3])
assert lut[name] == id_
with open(fname, "w") as fid:
for name, id_, color in zip(names, ids, colors):
out_color = " ".join(f"{x:03}" for x in color[:3]) # wrong length!
line = f"{id_} {name} {out_color}\n"
fid.write(line)
with pytest.raises(RuntimeError, match="formatted"):
read_freesurfer_lut(fname)
@testing.requires_testing_data
def test_talxfm_rigid():
"""Test that talxfm_rigid gives reasonable results."""
rigid = _estimate_talxfm_rigid("fsaverage", subjects_dir=subjects_dir)
assert_allclose(rigid, np.eye(4), atol=1e-6)
rigid = _estimate_talxfm_rigid("sample", subjects_dir=subjects_dir)
assert_allclose(np.linalg.norm(rigid[:3, :3], axis=1), 1.0, atol=1e-6)
move = 1000 * np.linalg.norm(rigid[:3, 3])
assert 30 < move < 70
ang = np.rad2deg(_angle_between_quats(rot_to_quat(rigid[:3, :3])))
assert 20 < ang < 25