脑磁图的探头位置坐标说明

本文对脑磁图的探头位置坐标及其提取和干预方式进行简单说明。

脑磁图的探头位置坐标说明

定位点及标准空间

一般来说，脑电或脑磁实验前，都需要对被试头部位置和设备探头的位置进行空间匹配。这个匹配不是对齐，而是想个办法让计算机知道被试大脑与设备之间的位置关系。从实用角度来讲，这个工作可以依靠在被试的眉间和双耳位置贴上三个位置传感器来实现。

因为这三个点天然地不共线，而不共线的三个点可以唯一地确定一个三维平面。

更好的是，从大脑的构造来讲，它可以比较准确地通过这三个点来定位。（眼睛保健操要揉这三个点还是有一定科学依据的，划掉）三个控制点的位置如下图，三点与大脑的位置关系图实在是比较大，因为放在最后的附录中。

三点位置

位置提取及手动干预

在实验过程中，设备是能够“知道”这三个控制点的位置的，再加上设备“知道”自己的探头在哪里。因此，这些数据都会存储在实验数据中。在实验结束后，我们当然可以通过数据分析的方法将这些数据提取出来。提取数据的工具有很多，我目前用着比较顺手的是 MNE 工具。

MNE - MNE 1.1.1 documentation^[1]

下面，我们顺着工作流，逐步进行到位置提取这一步，首先进行数据读取

# Read raw data, prevent preload to save time
raw = mne.io.read_raw_ctf(file_path, preload=False)

全部数据都存储在这个变量当中，它的非常有用的属性是 info 属性，包含实验的全部信息，比如实验设备有 273 个磁导联、 64 个脑电导联、28个参考定位导联，采样频率为 1200 Hz，初始滤波参数为 600 Hz 的带通滤波等。

我们当然能够从中获取转移矩阵，这个矩阵将用来进行空间平移、旋转和缩放等仿射变换，使用它的方法见之后的函数，这是我从 MNE 包中抄下来的，因为怕它改版之后又双叒叕变了

# Fetch dev_head_t from info, the info contains necessary informations besides data
dev_head_t = raw.info['dev_head_t']

def apply_trans(trans, pts, move=True):
    """Apply a transform matrix to an array of points.

    Parameters
    ----------
    trans : array, shape = (4, 4) | instance of Transform
        Transform matrix.
    pts : array, shape = (3,) | (n, 3)
        Array with coordinates for one or n points.
    move : bool
        If True (default), apply translation.

    Returns
    -------
    transformed_pts : shape = (3,) | (n, 3)
        Transformed point(s).
    """

    if isinstance(trans, dict):
        trans = trans['trans']
    pts = np.array(pts)
    if pts.size == 0:
        return pts.copy()

    # apply rotation & scale
    out_pts = np.dot(pts, trans[:3, :3].T)
    # apply translation
    if move:
        out_pts += trans[:3, 3]

    return out_pts

可以想见的，我们需要的全部定位点位置信息也在这个数据包里，获取方法如下。另外，有了正变换就应该对应反变换的方法，因此提供反变换的思路如下。当然，通过变换还可以在变换中途进行一些人为操作，比如将左耳的控制点降低 1 厘米

# Read info for position of the points
dps = raw.info['dig'][:3]
dev_head_t = raw.info['dev_head_t']
display(np.array([e['r'] for e in dps]), dev_head_t, [e for e in dps])

# Apply transform and inverse-transform
trans = dev_head_t['trans'].copy()
pt = np.array(dps[1]['r'])
# Transform
pt1 = np.dot(pt, trans[:3, :3].T) + trans[:3, 3]
# Downward the LPA by 1 cm
pt1[2] -= 0.01
# Inverse transform
pt2 = np.dot(pt1 - trans[:3, 3], np.linalg.inv(trans[:3, :3].T))
# Compare different
display(pt, pt1, pt2)

# Write back to the info
dps[1]['r'] = pt2

不人为操作的结果，用于检验正确性

进行人为操作的结果，用于验证功能

最后，我们需要整理一下代码，让它可以正确地找到三样东西，一是三个定位点的位置，另一个定位点的位置，这些都是系统内置的定位信息，另外还需要找到的探头位置，也就是真正采集到神经信号的位置。相应代码如下

# All digital sensors position, used for localization
dig = raw.info['dig']
df = pd.DataFrame(
    [
        (e['ident'], e['kind']._name, e['r'])
        for e in dig
    ],
    columns=['ident', 'kind', 'loc']
)

# All signal channels position
chs = raw.info['chs']
df = pd.DataFrame(
    [
        (e['scanno'], e['logno'], e['ch_name'], e['kind']._name, e['loc'])
        for e in chs
        if not np.isnan(e['loc'][0]) # Drop not-channel channels
    ],
    columns=['scanno', 'logno', 'chname', 'kind', 'loc']
)