[7f9fb8]: / mne / io / fil / sensors.py

Download this file

146 lines (120 with data), 4.5 kB

  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
# Authors: The MNE-Python contributors.
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.
from copy import deepcopy
import numpy as np
from ...utils import logger
def _refine_sensor_orientation(chanin):
"""Improve orientation matrices based on multiaxis measures.
The ex and ey elements from _convert_channel_info were oriented not
based on the physical orientation of the sensor.
It doesn't have to be this way, we can use (if available) the orientation
information from mulit-axis recordings to refine these elements.
"""
logger.info("Refining sensor orientations...")
chanout = deepcopy(chanin)
tmpname = list()
for ii in range(len(chanin)):
tmpname.append(chanin[ii]["ch_name"])
for ii in range(len(chanin)):
tmploc = deepcopy(chanin[ii]["loc"])
tmploc = tmploc.reshape(3, 4, order="F")
if np.isnan(tmploc.sum()) is False:
target, flipFlag = _guess_other_chan_axis(tmpname, ii)
if np.isnan(target) is False:
targetloc = deepcopy(chanin[target]["loc"])
if np.isnan(targetloc.sum()) is False:
targetloc = targetloc.reshape(3, 4, order="F")
tmploc[:, 2] = targetloc[:, 3]
tmploc[:, 1] = flipFlag * np.cross(tmploc[:, 2], tmploc[:, 3])
chanout[ii]["loc"] = tmploc.reshape(12, order="F")
logger.info("[done]")
return chanout
def _guess_other_chan_axis(tmpname, seedID):
"""Try to guess the name of another axis of a multiaxis sensor."""
# see if its using the old RAD/TAN convention first, otherwise use XYZ
if tmpname[seedID][-3:] == "RAD":
prefix1 = "RAD"
prefix2 = "TAN"
flipflag = 1.0
elif tmpname[seedID][-3:] == "TAN":
prefix1 = "TAN"
prefix2 = "RAD"
flipflag = -1.0
elif tmpname[seedID][-1:] == "Z" or tmpname[seedID][-3:] == "[Z]":
prefix1 = "Z"
prefix2 = "Y"
flipflag = -1.0
elif tmpname[seedID][-1:] == "Y" or tmpname[seedID][-3:] == "[Y]":
prefix1 = "Y"
prefix2 = "Z"
flipflag = 1.0
elif tmpname[seedID][-1:] == "X" or tmpname[seedID][-3:] == "[X]":
prefix1 = "X"
prefix2 = "Y"
flipflag = 1.0
else:
prefix1 = "?"
prefix2 = "?"
flipflag = 1.0
target_name = tmpname[seedID][: -len(prefix1)] + prefix2
target_id = np.where([t == target_name for t in tmpname])[0]
target_id = target_id[0] if len(target_id) else np.nan
return target_id, flipflag
def _get_pos_units(pos):
"""Get the units of a point cloud.
Determines the units a point cloud of sensor positions, provides the
scale factor required to ensure the units can be converted to meters.
"""
# get rid of None elements
nppos = np.empty((0, 3))
for ii in range(0, len(pos)):
if pos[ii] is not None and sum(np.isnan(pos[ii])) == 0:
nppos = np.vstack((nppos, pos[ii]))
idrange = np.empty(shape=(0, 3))
for ii in range(0, 3):
q90, q10 = np.percentile(nppos[:, ii], [90, 10])
idrange = np.append(idrange, q90 - q10)
size = np.linalg.norm(idrange)
unit, sf = _size2units(size)
return unit, sf
def _size2units(size):
"""Convert the size returned from _get_pos_units into a physical unit."""
if size >= 0.050 and size < 0.500:
unit = "m"
sf = 1
elif size >= 0.50 and size < 5:
unit = "dm"
sf = 10
elif size >= 5 and size < 50:
unit = "cm"
sf = 100
elif size >= 50 and size < 500:
unit = "mm"
sf = 1000
else:
unit = "unknown"
sf = 1
return unit, sf
def _get_plane_vectors(ez):
"""Get two orthogonal vectors orthogonal to ez (ez will be modified).
Note: the ex and ey positions will not be realistic, this can be fixed
using _refine_sensor_orientation.
"""
assert ez.shape == (3,)
ez_len = np.sqrt(np.sum(ez * ez))
if ez_len == 0:
raise RuntimeError("Zero length normal. Cannot proceed.")
if np.abs(ez_len - np.abs(ez[2])) < 1e-5: # ez already in z-direction
ex = np.array([1.0, 0.0, 0.0])
else:
ex = np.zeros(3)
if ez[1] < ez[2]:
ex[0 if ez[0] < ez[1] else 1] = 1.0
else:
ex[0 if ez[0] < ez[2] else 2] = 1.0
ez /= ez_len
ex -= np.dot(ez, ex) * ez
ex /= np.sqrt(np.sum(ex * ex))
ey = np.cross(ez, ex)
return ex, ey