Skip to content

Commit ac43da0

Browse files
larsonerlarsoner
committed
ENH: Add helmet deformation for OPMs
1 parent da8c794 commit ac43da0

File tree

4 files changed

+252
-6
lines changed

4 files changed

+252
-6
lines changed

‎mne/data/helmets/Kernel_Flux.fif.gz

-5.89 KB
Binary file not shown.

‎mne/data/helmets/Kernel_Flux_ch_pos.txt

Lines changed: 202 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,202 @@
1+
{
2+
"MA1": [
3+
-0.040249,
4+
0.092195,
5+
0.024061
6+
],
7+
"MA2": [
8+
0.0,
9+
0.10462,
10+
0.016906
11+
],
12+
"MA3": [
13+
0.040249,
14+
0.092195,
15+
0.024061
16+
],
17+
"MA4": [
18+
0.0,
19+
0.083089,
20+
0.053953
21+
],
22+
"MB1": [
23+
-0.040251,
24+
0.063041,
25+
0.05933
26+
],
27+
"MB2": [
28+
0.0,
29+
0.045194,
30+
0.077115
31+
],
32+
"MB3": [
33+
0.040251,
34+
0.063041,
35+
0.05933
36+
],
37+
"MB4": [
38+
-0.036861,
39+
0.026016,
40+
0.076163
41+
],
42+
"MB5": [
43+
0.036861,
44+
0.026016,
45+
0.076163
46+
],
47+
"MC1": [
48+
-0.035716,
49+
-0.017727,
50+
0.077709
51+
],
52+
"MC2": [
53+
0.0,
54+
0.00196,
55+
0.085517
56+
],
57+
"MC3": [
58+
0.035716,
59+
-0.017727,
60+
0.077709
61+
],
62+
"MC4": [
63+
-0.034502,
64+
-0.056871,
65+
0.064752
66+
],
67+
"MC5": [
68+
0.0,
69+
-0.037943,
70+
0.080668
71+
],
72+
"MC6": [
73+
0.034502,
74+
-0.056871,
75+
0.064752
76+
],
77+
"MD1": [
78+
-0.036596,
79+
-0.084035,
80+
0.027244
81+
],
82+
"MD2": [
83+
0,
84+
-0.082423,
85+
0.048749
86+
],
87+
"MD3": [
88+
0.036596,
89+
-0.084035,
90+
0.027244
91+
],
92+
"MD4": [
93+
0,
94+
-0.095211,
95+
0.008834
96+
],
97+
"ME1": [
98+
-0.038329,
99+
-0.084588,
100+
-0.017313
101+
],
102+
"ME2": [
103+
0.0,
104+
-0.086368,
105+
-0.035705
106+
],
107+
"ME3": [
108+
0.038329,
109+
-0.084588,
110+
-0.017313
111+
],
112+
"ME4": [
113+
-0.028545,
114+
-0.071814,
115+
-0.06008
116+
],
117+
"ME5": [
118+
0.028545,
119+
-0.071814,
120+
-0.06008
121+
],
122+
"RA1": [
123+
0.06916,
124+
0.059735,
125+
0.008274
126+
],
127+
"RA2": [
128+
0.070053,
129+
0.032925,
130+
0.039006
131+
],
132+
"RB1": [
133+
0.069394,
134+
-0.009905,
135+
0.045791
136+
],
137+
"RB2": [
138+
0.06596,
139+
-0.048253,
140+
0.029541
141+
],
142+
"RC1": [
143+
0.07757,
144+
0.023536,
145+
-0.014638
146+
],
147+
"RC2": [
148+
0.077612,
149+
-0.013307,
150+
0.002181
151+
],
152+
"RC3": [
153+
0.068927,
154+
-0.049354,
155+
-0.013998
156+
],
157+
"RC4": [
158+
0.06444,
159+
-0.040455,
160+
-0.053382
161+
],
162+
"LA1": [
163+
-0.06916,
164+
0.059735,
165+
0.008274
166+
],
167+
"LA2": [
168+
-0.070053,
169+
0.032925,
170+
0.039006
171+
],
172+
"LB1": [
173+
-0.069394,
174+
-0.009905,
175+
0.045791
176+
],
177+
"LB2": [
178+
-0.06596,
179+
-0.048253,
180+
0.029541
181+
],
182+
"LC1": [
183+
-0.07757,
184+
0.023536,
185+
-0.014638
186+
],
187+
"LC2": [
188+
-0.077612,
189+
-0.013307,
190+
0.002181
191+
],
192+
"LC3": [
193+
-0.068927,
194+
-0.049354,
195+
-0.013998
196+
],
197+
"LC4": [
198+
-0.06444,
199+
-0.040455,
200+
-0.053382
201+
]
202+
}

‎mne/surface.py

Lines changed: 44 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,11 +8,13 @@
88
# Many of the computations in this code were derived from Matti Hämäläinen's
99
# C code.
1010

11+
from collections import OrderedDict
1112
from copy import deepcopy
1213
from functools import partial, lru_cache
13-
from collections import OrderedDict
1414
from glob import glob
15+
import json
1516
from os import path as op
17+
from pathlib import Path
1618
import time
1719
import warnings
1820

@@ -33,6 +35,9 @@
3335
_get_trans,
3436
apply_trans,
3537
Transform,
38+
_fit_matched_points,
39+
_MatchedDisplacementFieldInterpolator,
40+
_angle_between_quats,
3641
)
3742
from .utils import (
3843
logger,
@@ -52,6 +57,8 @@
5257
_import_nibabel,
5358
)
5459

60+
_helmet_path = Path(__file__).parent / "data" / "helmets"
61+
5562

5663
###############################################################################
5764
# AUTOMATED SURFACE FINDING
@@ -205,10 +212,11 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
205212
system, have_helmet = _get_meg_system(info)
206213
if have_helmet:
207214
logger.info("Getting helmet for system %s" % system)
208-
fname = op.join(op.split(__file__)[0], "data", "helmets", system + ".fif.gz")
215+
fname = _helmet_path / f"{system}.fif.gz"
209216
surf = read_bem_surfaces(
210217
fname, False, FIFF.FIFFV_MNE_SURF_MEG_HELMET, verbose=False
211218
)
219+
surf = _scale_helmet_to_sensors(system, surf, info)
212220
else:
213221
rr = np.array(
214222
[
@@ -248,6 +256,40 @@ def get_meg_helmet_surf(info, trans=None, *, verbose=None):
248256
return surf
249257

250258

259+
def _scale_helmet_to_sensors(system, surf, info):
260+
fname = _helmet_path / f"{system}_ch_pos.txt"
261+
if not fname.is_file():
262+
return surf
263+
with open(fname) as fid:
264+
ch_pos_from = json.load(fid)
265+
# find correspondence
266+
fro, to = list(), list()
267+
for key, f_ in ch_pos_from.items():
268+
t_ = [ch["loc"][:3] for ch in info["chs"] if ch["ch_name"].startswith(key)]
269+
if not len(t_):
270+
continue
271+
fro.append(f_)
272+
to.append(np.mean(t_, axis=0))
273+
if len(fro) < 4:
274+
logger.info("Using CAD helmet because fewer than 4 sensors found")
275+
return surf
276+
fro = np.array(fro, float)
277+
to = np.array(to, float)
278+
interp = _MatchedDisplacementFieldInterpolator(fro, to)
279+
new_rr = interp(surf["rr"])
280+
quat, sc = _fit_matched_points(surf["rr"], new_rr)
281+
rot = np.rad2deg(_angle_between_quats(quat[:3]))
282+
tr = 1000 * np.linalg.norm(quat[3:])
283+
logger.info(f" Deforming to match info using {len(fro)} matched points:")
284+
logger.info(f" 1. Affine: {rot:0.1f}°, {tr:0.1f} mm, {sc:0.2f}× scale")
285+
deltas = interp._last_deltas * 1000
286+
mu, mx = np.mean(deltas), np.max(deltas)
287+
logger.info(f" 2. Nonlinear displacement: " f"mean={mu:0.1f}, max={mx:0.1f} mm")
288+
surf["rr"] = new_rr
289+
complete_surface_info(surf, copy=False, verbose=False)
290+
return surf
291+
292+
251293
def _reorder_ccw(rrs, tris):
252294
"""Reorder tris of a convex hull to be wound counter-clockwise."""
253295
# This ensures that rendering with front-/back-face culling works properly

‎mne/transforms.py

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -2095,9 +2095,9 @@ def __init__(self, fro, to):
20952095
assert fro.shape[1] == 3
20962096

20972097
# Prealign using affine + uniform scaling
2098-
trans, scale = _fit_matched_points(fro, to, scale=True)
2099-
trans = _quat_to_affine(trans)
2100-
trans[:3, :3] *= scale
2098+
self._quat, self._scale = _fit_matched_points(fro, to, scale=True)
2099+
trans = _quat_to_affine(self._quat)
2100+
trans[:3, :3] *= self._scale
21012101
self._affine = trans
21022102
fro = apply_trans(trans, fro)
21032103

@@ -2114,6 +2114,8 @@ def __init__(self, fro, to):
21142114
def __call__(self, x):
21152115
assert x.ndim in (1, 2) and x.shape[-1] == 3
21162116
singleton = x.ndim == 1
2117-
out = self._interp(apply_trans(self._affine, x))
2117+
x = apply_trans(self._affine, x)
2118+
out = self._interp(x)
2119+
self._last_deltas = np.linalg.norm(x - out, axis=1)
21182120
out = out[0] if singleton else out
21192121
return out

0 commit comments

Comments
 (0)