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
|
# Authors: Mark Wronkiewicz <wronk@uw.edu>
# Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD-3-Clause
import os.path as op
import numpy as np
from numpy.testing import assert_allclose
import pytest
from mne import pick_types
from mne.io import read_raw_fif
from mne.datasets import testing
from mne.io.tag import _loc_to_coil_trans
from mne.preprocessing import (read_fine_calibration, write_fine_calibration,
compute_fine_calibration, maxwell_filter)
from mne.preprocessing.tests.test_maxwell import _assert_shielding
from mne.transforms import rot_to_quat, _angle_between_quats
from mne.utils import object_diff
# Define fine calibration filepaths
data_path = testing.data_path(download=False)
fine_cal_fname = op.join(data_path, 'SSS', 'sss_cal_3053.dat')
fine_cal_fname_3d = op.join(data_path, 'SSS', 'sss_cal_3053_3d.dat')
erm_fname = op.join(data_path, 'SSS', '141027_cropped_90Hz_raw.fif')
ctc = op.join(data_path, 'SSS', 'ct_sparse.fif')
cal_mf_fname = op.join(data_path, 'SSS', '141027.dat')
@pytest.mark.parametrize('fname', (cal_mf_fname, fine_cal_fname,
fine_cal_fname_3d))
@testing.requires_testing_data
def test_fine_cal_io(tmp_path, fname):
"""Test round trip reading/writing of fine calibration .dat file."""
temp_fname = op.join(str(tmp_path), 'fine_cal_temp.dat')
# Load fine calibration file
fine_cal_dict = read_fine_calibration(fname)
# Save temp version of fine calibration file
write_fine_calibration(temp_fname, fine_cal_dict)
fine_cal_dict_reload = read_fine_calibration(temp_fname)
# Load temp version of fine calibration file and compare hashes
assert object_diff(fine_cal_dict, fine_cal_dict_reload) == ''
@pytest.mark.slowtest
@testing.requires_testing_data
def test_compute_fine_cal():
"""Test computing fine calibration coefficients."""
raw = read_raw_fif(erm_fname)
want_cal = read_fine_calibration(cal_mf_fname)
got_cal, counts = compute_fine_calibration(
raw, cross_talk=ctc, n_imbalance=1, verbose='debug')
assert counts == 1
assert set(got_cal.keys()) == set(want_cal.keys())
assert got_cal['ch_names'] == want_cal['ch_names']
# in practice these should never be exactly 1.
assert sum([(ic == 1.).any() for ic in want_cal['imb_cals']]) == 0
assert sum([(ic == 1.).any() for ic in got_cal['imb_cals']]) == 0
got_imb = np.array(got_cal['imb_cals'], float)
want_imb = np.array(want_cal['imb_cals'], float)
assert got_imb.shape == want_imb.shape == (306, 1)
got_imb, want_imb = got_imb[:, 0], want_imb[:, 0]
orig_locs = np.array([ch['loc'] for ch in raw.info['chs'][:306]])
want_locs = want_cal['locs']
got_locs = got_cal['locs']
assert want_locs.shape == got_locs.shape
orig_trans = _loc_to_coil_trans(orig_locs)
want_trans = _loc_to_coil_trans(want_locs)
got_trans = _loc_to_coil_trans(got_locs)
dist = np.linalg.norm(got_trans[:, :3, 3] - want_trans[:, :3, 3], axis=1)
assert_allclose(dist, 0., atol=1e-6)
dist = np.linalg.norm(got_trans[:, :3, 3] - orig_trans[:, :3, 3], axis=1)
assert_allclose(dist, 0., atol=1e-6)
orig_quat = rot_to_quat(orig_trans[:, :3, :3])
want_quat = rot_to_quat(want_trans[:, :3, :3])
got_quat = rot_to_quat(got_trans[:, :3, :3])
want_orig_angles = np.rad2deg(_angle_between_quats(want_quat, orig_quat))
got_want_angles = np.rad2deg(_angle_between_quats(got_quat, want_quat))
got_orig_angles = np.rad2deg(_angle_between_quats(got_quat, orig_quat))
for key in ('mag', 'grad'):
# imb_cals value
p = pick_types(raw.info, meg=key, exclude=())
r2 = np.dot(got_imb[p], want_imb[p]) / (
np.linalg.norm(want_imb[p]) * np.linalg.norm(got_imb[p]))
assert 0.99 < r2 <= 1.00001, f'{key}: {r2:0.3f}'
# rotation angles
want_orig_max_angle = want_orig_angles[p].max()
got_orig_max_angle = got_orig_angles[p].max()
got_want_max_angle = got_want_angles[p].max()
if key == 'mag':
assert 8 < want_orig_max_angle < 11, want_orig_max_angle
assert 1 < got_orig_max_angle < 2, got_orig_max_angle
assert 9 < got_want_max_angle < 11, got_want_max_angle
else:
# Some of these angles are large, but mostly this has to do with
# processing a very short (one 10-sec segment), downsampled (90 Hz)
# file
assert 66 < want_orig_max_angle < 68, want_orig_max_angle
assert 56 < got_orig_max_angle < 113, got_orig_max_angle
assert 53 < got_want_max_angle < 60, got_want_max_angle
kwargs = dict(bad_condition='warning', cross_talk=ctc, coord_frame='meg')
raw_sss = maxwell_filter(raw, **kwargs)
raw_sss_mf = maxwell_filter(raw, calibration=cal_mf_fname, **kwargs)
raw_sss_py = maxwell_filter(raw, calibration=got_cal, **kwargs)
_assert_shielding(raw_sss, raw, 26, 27)
_assert_shielding(raw_sss_mf, raw, 61, 63)
_assert_shielding(raw_sss_py, raw, 61, 63)
# redoing with given mag data should yield same result
got_cal_redo, _ = compute_fine_calibration(
raw, cross_talk=ctc, n_imbalance=1, calibration=got_cal,
verbose='debug')
assert got_cal['ch_names'] == got_cal_redo['ch_names']
assert_allclose(got_cal['imb_cals'], got_cal_redo['imb_cals'], atol=5e-5)
assert_allclose(got_cal['locs'], got_cal_redo['locs'], atol=1e-6)
assert sum([(ic == 1.).any() for ic in got_cal['imb_cals']]) == 0
# redoing with 3 imlabance parameters should improve the shielding factor
grad_picks = pick_types(raw.info, meg='grad')
assert len(grad_picks) == 204 and grad_picks[0] == 0
got_grad_imbs = np.array(
[got_cal['imb_cals'][pick] for pick in grad_picks])
assert got_grad_imbs.shape == (204, 1)
got_cal_3, _ = compute_fine_calibration(
raw, cross_talk=ctc, n_imbalance=3, calibration=got_cal,
verbose='debug')
got_grad_3_imbs = np.array([
got_cal_3['imb_cals'][pick] for pick in grad_picks])
assert got_grad_3_imbs.shape == (204, 3)
corr = np.corrcoef(got_grad_3_imbs[:, 0], got_grad_imbs[:, 0])[0, 1]
assert 0.6 < corr < 0.7
raw_sss_py = maxwell_filter(raw, calibration=got_cal_3, **kwargs)
_assert_shielding(raw_sss_py, raw, 68, 70)
|
