Python-MNE-源空间和正模型04：头模型和前向计算

我们知道，在MNE分析中坐标是很重要的，这个前面也提及过了配准的一些方法，总的来说，MNE和freesurfer中使用的配准系统以及他们之间的关系如下图所示：除了传感器坐标之外，所有的坐标系都是笛卡尔坐标系，并且具有“RAS”（x右-y前-z上）的方向。

上图中，在FreeSurfer文件以及MNE中的fif文件中存在的坐标转换信息，以及一些设置为固定值的坐标转换用Tx表示，其中x标识了不同的转换过程。

Head coordinates

这是根据基准坐标（鼻尖和双耳）定义的坐标系。在.fif文件中，EEG电极的位置在这个坐标系中给出，这也就是在MEG、MEG/EEG或EEG采集开始前获取的头部数字化数据。

Device coordinates

这是一个跟MEG设备绑定的坐标系。在MEG测量过程中，通过向3-5个头部位置指示器（HPI）线圈提供电流，根据它们产生的磁场确定它们相对于MEG传感器阵列的位置，这样就能确定MEG设备和头部坐标的关系。

Sensor coordinates

每个MEG传感器都有一个坐标系来定义传感器的方向和位置。利用该坐标系就可以用线圈的几何向量来计算磁场积分。.fif文件中的channel数据就包含各传感器坐标与MEG设备坐标之间转换的信息。

Surface RAS coordinates

FreeSurfer计算出来的脑表面数据用这个坐标系表示。该坐标系统的原点位于FreeSurfer MRI的中心（通常为256 x 256 x 256各向同性的1 mm3体素），坐标轴沿该体素小块的坐标轴定向。BEM表面、信号源在源空间中的位置通常在.fif文件中都用这个坐标系表示。Surface RAS coordinates在MNE中可以被等同于MRI coordinates。

RAS coordinates

该坐标系的轴与Surface RAS coordinates相同，但原点的位置依赖于真正的MRI仪器的定义。在用MNE软件进行分析时，几乎不需要引用这个坐标系，但由于下面讨论的Talairach coordinates是根据它而不是Surface RAS coordinates定义的，因此RAS coordinates参与了Surface RAS coordinates和两个Talairach coordinates之间的转换，在此提及一下。

MNI Talairach coordinates

在https://imaging.mrc-cbu.cam.ac.uk/imaging/MniTalairach中有这个坐标系的定义。

Talairach坐标基于 Talairach and Tournoux (1988) 提出的脑图谱，是一个标准化的大脑坐标系统，Talairach坐标最初是基于一个脑部解剖切片的标准化模板，主要用于小样本（通常是个体脑）的手动配准。而MNI（Montreal Neurological Institute）坐标 是由MNI开发的标准化脑图谱，基于多人的脑部MRI扫描（通常10-152名健康个体平均脑图）。MNI坐标解决了Talairach空间因基于单人脑模板而存在的限制，提供了更适合群体研究的标准化空间。

MNI Talairach坐标用于将EEG/MEG源定位结果与标准大脑模板对齐，方便跨个体、跨实验比较，并与神经解剖学和功能脑区信息对照，二者需要转换工具进行互相映射，这个转换是在FreeSurfer重建过程中确定的。

FreeSurfer Talairach coordinates

MNI Talairach的缺陷是MNI Talairach变换并不完全匹配你自己受试大脑和Talairach的大脑。而因为MNE是基于多个人的，所以MNI的大脑比Talairach的大脑稍大（特别是更高、更深和更长）。从大脑中部向外走得越远，差异就越大。FreeSurfer通过添加一个修正变换来解决这个问题，该变换在MEG/EEG和MRI坐标系中固定用T-和T+表示（下图所示，它是硬编码的，没办法改参数的）。详见https://imaging.mrc-cbu.cam.ac.uk/imaging/MniTalairach（方法2）。

上图的这些转换，放到MNE或者Freesurfer里，其实就是如下图所示：

比如我们之前生成过的，之后也会生成的fwd.fif，inv.fif，-trans.fif，都对应了某个转换过程。

现在，让我们正式开始计算前向算子（forward operator），也就是对应图上的T1转换过程。

计算正算子需要什么

首先导入所需的数据，这里还是使用了自带的sample数据集。

import mne
from mne.datasets import sample

data_path = sample.data_path()

# the raw file containing the channel location + types
sample_dir = data_path / "MEG" / "sample"
raw_fname = sample_dir / "sample_audvis_raw.fif"
# The paths to Freesurfer reconstructions
subjects_dir = data_path / "subjects"
subject = "sample"

计算前向算子，我们需要：

• 一个包含共配准信息的 -trans.fif 文件，上面T2转换过程的一部分
• 一个源空间
• BEM表面

计算和可视化BEM表面

BEM表面是前向计算所需的不同组织之间界面的三角形小面剖分，这些表面有颅骨内表面、颅骨外表面、头皮表面。计算BEM表面需要FreeSurfer，并使用命令行工具mne watershed_bem或mne flash_bem，或相关函数mne.bem.make_watershed_bem()或mne.bem.make_flash_bem()。
一般咱们就用mne watershed_bem，freesurfer生成以后直接跑完了，具体用法如下：

对于EEG，我们使用3层（内颅骨，外颅骨和皮肤），而对于MEG，1层（内颅骨）就足够了，因为磁信号不用考虑电导率，所以不用考虑其他东西

可以用函数mne.viz.plot_bem()来查看所创建的BEM表面。这里我们使用比默认值更小的slices来提高速度。

plot_bem_kwargs = dict(
    subject=subject,
    subjects_dir=subjects_dir,
    brain_surfaces="white",
    orientation="coronal",
    slices=[50, 100, 150, 200],
)

mne.viz.plot_bem(**plot_bem_kwargs)

Using surface:
       /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
       Using surface:
       /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skull.surf
       Using surface:
       /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skin.surf

可视化共配准

共配准是一种操作，它允许在一个共同的坐标系中定位头部和传感器。在MNE软件中，将头部和传感器对齐的转换存储在trans文件中，它是一个以-trans.fif结尾的文件。它可以通过mne.gui.coregistration()或命令行命令mne coreg [options]获得。在这里，我们假设共配已在上一期所说的教程中完成了。

# The transformation file obtained by coregistration
trans = sample_dir / "sample_audvis_raw-trans.fif"

info = mne.io.read_info(raw_fname)
# Here we look at the dense head, which isn't used for BEM computations but
# is useful for coregistration.
mne.viz.plot_alignment(
    info,
    trans,
    subject=subject,
    dig=True,
    meg=["helmet", "sensors"],
    subjects_dir=subjects_dir,
    surfaces="head-dense",
)

Read a total of 3 projection items:
        PCA-v1 (1 x 102)  idle
        PCA-v2 (1 x 102)  idle
        PCA-v3 (1 x 102)  idle
True
Using lh.seghead for head surface.
Getting helmet for system 306m
Channel types:: grad: 203, mag: 102, eeg: 59

计算源空间

源定位是我们的目标结果，而源空间则定义了所要定位的源的位置和方向。有两种类型的源空间：

• 表面源空间：当源在皮层时适用。
• 体积源空间或离散源空间：当源是离散的、或在皮层深部时适用。

表面源空间使用mne.setup_source_space()计算，而体积源空间使用mne.setup_volume_source_space()计算。

现在我们将以'oct4'分辨率计算一个表面源空间。这个分辨率参数事实上就是如何对源空间进行细分，比如ico就是用二十面体，而oct就是用八面体来对源空间进行划分，在正式实验中，我们可以追求更高的分辨率。

src = mne.setup_source_space(
    subject, spacing="oct4", add_dist="patch", subjects_dir=subjects_dir
)
print(src)

etting up the source space with the following parameters:

SUBJECTS_DIR = /home/circleci/mne_data/MNE-sample-data/subjects
Subject      = sample
Surface      = white
Octahedron subdivision grade 4

>>> 1. Creating the source space...

Doing the octahedral vertex picking...
Loading /home/circleci/mne_data/MNE-sample-data/subjects/sample/surf/lh.white...
Mapping lh sample -> oct (4) ...
    Triangle neighbors and vertex normals...
Loading geometry from /home/circleci/mne_data/MNE-sample-data/subjects/sample/surf/lh.sphere...
Setting up the triangulation for the decimated surface...
loaded lh.white 258/155407 selected to source space (oct = 4)

Loading /home/circleci/mne_data/MNE-sample-data/subjects/sample/surf/rh.white...
Mapping rh sample -> oct (4) ...
    Triangle neighbors and vertex normals...
Loading geometry from /home/circleci/mne_data/MNE-sample-data/subjects/sample/surf/rh.sphere...
Setting up the triangulation for the decimated surface...
loaded rh.white 258/156866 selected to source space (oct = 4)

Calculating patch information (limit=0.0 mm)...
    Computing patch statistics...
    Patch information added...
    Computing patch statistics...
    Patch information added...
You are now one step closer to computing the gain matrix
<SourceSpaces: [<surface (lh), n_vertices=155407, n_used=258>, <surface (rh), n_vertices=156866, n_used=258>] MRI (surface RAS) coords, subject 'sample', ~28.7 MiB>

表面源空间src包含两个部分，一个左半球（258个位置）一个右半球（258个位置）。源空间可以在BEM表面上以紫色显示。

mne.viz.plot_bem(src=src, **plot_bem_kwargs)

而体积源空间由半径90mm、中心为（0.0,0.0,40.0）mm的球体内的偶极子网格定义，可以使用以下代码来可视化。从图上看，显然这个球体不完美，因为它有一部分并不局限于大脑，有一部分没有覆盖到脑表面。

sphere = (0.0, 0.0, 0.04, 0.09)
vol_src = mne.setup_volume_source_space(
    subject,
    subjects_dir=subjects_dir,
    sphere=sphere,
    sphere_units="m",
    add_interpolator=False,
)  # just for speed!
print(vol_src)

mne.viz.plot_bem(src=vol_src, **plot_bem_kwargs)

Sphere                : origin at (0.0 0.0 40.0) mm
              radius  : 90.0 mm
grid                  : 5.0 mm
mindist               : 5.0 mm
MRI volume            : /home/circleci/mne_data/MNE-sample-data/subjects/sample/mri/T1.mgz

Reading /home/circleci/mne_data/MNE-sample-data/subjects/sample/mri/T1.mgz...

Setting up the sphere...
Surface CM = (   0.0    0.0   40.0) mm
Surface fits inside a sphere with radius   90.0 mm
Surface extent:
    x =  -90.0 ...   90.0 mm
    y =  -90.0 ...   90.0 mm
    z =  -50.0 ...  130.0 mm
Grid extent:
    x =  -95.0 ...   95.0 mm
    y =  -95.0 ...   95.0 mm
    z =  -50.0 ...  135.0 mm
57798 sources before omitting any.
24365 sources after omitting infeasible sources not within 0.0 - 90.0 mm.
20377 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
    0.005000 0.000000 0.000000     -95.00 mm
    0.000000 0.005000 0.000000     -95.00 mm
    0.000000 0.000000 0.005000     -50.00 mm
    0.000000 0.000000 0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.001000 0.000000 0.000000     128.00 mm
    0.000000 0.000000 0.001000    -128.00 mm
    0.000000 -0.001000 0.000000     128.00 mm
    0.000000 0.000000 0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
    1.000000 -0.000000 -0.000000      -5.27 mm
    -0.000000 1.000000 -0.000000       9.04 mm
    -0.000000 0.000000 1.000000     -27.29 mm
    0.000000 0.000000 0.000000       1.00
<SourceSpaces: [<volume, shape=(np.int64(39), np.int64(39), np.int64(38)), n_used=20377>] MRI (surface RAS) coords, subject 'sample', ~14.3 MiB>
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skull.surf
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skin.surf

所以要解决这个问题，计算基于脑内偶极子网格的体积源空间需要BEM表面，可以使用以下方法：

surface = subjects_dir / subject / "bem" / "inner_skull.surf"
vol_src = mne.setup_volume_source_space(
    subject, subjects_dir=subjects_dir, surface=surface, add_interpolator=False
)  # Just for speed!
print(vol_src)

mne.viz.plot_bem(src=vol_src, **plot_bem_kwargs)

Boundary surface file : /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
grid                  : 5.0 mm
mindist               : 5.0 mm
MRI volume            : /home/circleci/mne_data/MNE-sample-data/subjects/sample/mri/T1.mgz

Reading /home/circleci/mne_data/MNE-sample-data/subjects/sample/mri/T1.mgz...

Loaded bounding surface from /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf (2562 nodes)
Surface CM = (   0.7  -10.0   44.3) mm
Surface fits inside a sphere with radius   91.8 mm
Surface extent:
    x =  -66.7 ...   68.8 mm
    y =  -88.0 ...   79.0 mm
    z =  -44.5 ...  105.8 mm
Grid extent:
    x =  -70.0 ...   70.0 mm
    y =  -90.0 ...   80.0 mm
    z =  -45.0 ...  110.0 mm
32480 sources before omitting any.
22941 sources after omitting infeasible sources not within 0.0 - 91.8 mm.
Source spaces are in MRI coordinates.
Checking that the sources are inside the surface and at least    5.0 mm away (will take a few...)
Checking surface interior status for 22941 points...
    Found  2787/22941 points inside  an interior sphere of radius   43.6 mm
    Found     0/22941 points outside an exterior sphere of radius   91.8 mm
    Found  9357/20154 points outside using surface Qhull
    Found   846/10797 points outside using solid angles
    Total 12738/22941 points inside the surface
Interior check completed in 8836.5 ms
    10203 source space points omitted because they are outside the inner skull surface.
    2362 source space points omitted because of the    5.0-mm distance limit.
10376 sources remaining after excluding the sources outside the surface and less than    5.0 mm inside.
Adjusting the neighborhood info.
Source space : MRI voxel -> MRI (surface RAS)
    0.005000 0.000000 0.000000     -70.00 mm
    0.000000 0.005000 0.000000     -90.00 mm
    0.000000 0.000000 0.005000     -45.00 mm
    0.000000 0.000000 0.000000       1.00
MRI volume : MRI voxel -> MRI (surface RAS)
    -0.001000 0.000000 0.000000     128.00 mm
    0.000000 0.000000 0.001000    -128.00 mm
    0.000000 -0.001000 0.000000     128.00 mm
    0.000000 0.000000 0.000000       1.00
MRI volume : MRI (surface RAS) -> RAS (non-zero origin)
    1.000000 -0.000000 -0.000000      -5.27 mm
    -0.000000 1.000000 -0.000000       9.04 mm
    -0.000000 0.000000 1.000000     -27.29 mm
    0.000000 0.000000 0.000000       1.00
<SourceSpaces: [<volume, shape=(np.int64(29), np.int64(35), np.int64(32)), n_used=10376>] MRI (surface RAS) coords, subject 'sample', ~8.0 MiB>
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/inner_skull.surf
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skull.surf
Using surface: /home/circleci/mne_data/MNE-sample-data/subjects/sample/bem/outer_skin.surf

现在就好了。

现在让我们3D可视化我们生成的表面源空间。

fig = mne.viz.plot_alignment(
    subject=subject,
    subjects_dir=subjects_dir,
    surfaces="white",
    coord_frame="mri",
    src=src,
)
mne.viz.set_3d_view(
    fig,
    azimuth=173.78,
    elevation=101.75,
    distance=0.30,
    focalpoint=(-0.03, -0.01, 0.03),
)

True

计算正算子

现在我们计算正问题的解，也就是正算子。为了减少计算量，我们只计算一个单层BEM，即只计算了可用于MEG的内头骨（EEG要三层）。我们可以使用电导率参数来指定是想要单层的BEM还是三层的BEM。

conductivity = (0.3,)  # for single layer
# conductivity = (0.3, 0.006, 0.3)  # for three layers
model = mne.make_bem_model(
    subject="sample", ico=4, conductivity=conductivity, subjects_dir=subjects_dir
)
bem = mne.make_bem_solution(model)

Creating the BEM geometry...
Going from 4th to 4th subdivision of an icosahedron (n_tri: 5120 -> 5120)
inner skull CM is   0.67 -10.01  44.26 mm
Surfaces passed the basic topology checks.
Complete.

Homogeneous model surface loaded.
Computing the linear collocation solution...
    Matrix coefficients...
        inner skull (2562) -> inner skull (2562) ...
    Inverting the coefficient matrix...
Solution ready.
BEM geometry computations complete.

注意，BEM不涉及任何trans文件的使用。BEM仅取决于头部几何形状和电导率，它独立于MEG的设备坐标和头部坐标。

现在让我们计算正算子，通常被称为增益矩阵或引线场矩阵（gain or leadfield matrix）。我们使用mne.make_forward_solution()的方法，因为它比较重要，所以可以参考这个链接来了解每个参数的详细含义。

fwd = mne.make_forward_solution(
    raw_fname,
    trans=trans,
    src=src,
    bem=bem,
    meg=True,
    eeg=False,
    mindist=5.0,
    n_jobs=None,
    verbose=True,
)
print(fwd)

Source space          : <SourceSpaces: [<surface (lh), n_vertices=155407, n_used=258>, <surface (rh), n_vertices=156866, n_used=258>] MRI (surface RAS) coords, subject 'sample', ~28.7 MiB>
MRI -> head transform : /home/circleci/mne_data/MNE-sample-data/MEG/sample/sample_audvis_raw-trans.fif
Measurement data      : sample_audvis_raw.fif
Conductor model   : instance of ConductorModel
Accurate field computations
Do computations in head coordinates
Free source orientations

Read 2 source spaces a total of 516 active source locations

Coordinate transformation: MRI (surface RAS) -> head
    0.999310 0.009985 -0.035787      -3.17 mm
    0.012759 0.812405 0.582954       6.86 mm
    0.034894 -0.583008 0.811716      28.88 mm
    0.000000 0.000000 0.000000       1.00

Read 306 MEG channels from info
105 coil definitions read
Coordinate transformation: MEG device -> head
    0.991420 -0.039936 -0.124467      -6.13 mm
    0.060661 0.984012 0.167456       0.06 mm
    0.115790 -0.173570 0.977991      64.74 mm
    0.000000 0.000000 0.000000       1.00
MEG coil definitions created in head coordinates.
Source spaces are now in head coordinates.

Employing the head->MRI coordinate transform with the BEM model.
BEM model instance of ConductorModel is now set up

Source spaces are in head coordinates.
Checking that the sources are inside the surface and at least    5.0 mm away (will take a few...)
Checking surface interior status for 258 points...
    Found  54/258 points inside  an interior sphere of radius   43.6 mm
    Found   0/258 points outside an exterior sphere of radius   91.8 mm
    Found   0/204 points outside using surface Qhull
    Found   0/204 points outside using solid angles
    Total 258/258 points inside the surface
Interior check completed in 197.5 ms
    23 source space point omitted because of the    5.0-mm distance limit.
    Computing patch statistics...
    Patch information added...
Checking surface interior status for 258 points...
    Found  51/258 points inside  an interior sphere of radius   43.6 mm
    Found   0/258 points outside an exterior sphere of radius   91.8 mm
    Found   0/207 points outside using surface Qhull
    Found   0/207 points outside using solid angles
    Total 258/258 points inside the surface
Interior check completed in 205.1 ms
    19 source space point omitted because of the    5.0-mm distance limit.
    Computing patch statistics...
    Patch information added...

Checking surface interior status for 306 points...
    Found   0/306 points inside  an interior sphere of radius   43.6 mm
    Found 306/306 points outside an exterior sphere of radius   91.8 mm
    Found   0/  0 points outside using surface Qhull
    Found   0/  0 points outside using solid angles
    Total 0/306 points inside the surface
Interior check completed in 30.5 ms

Composing the field computation matrix...
Computing MEG at 474 source locations (free orientations)...

Finished.
<Forward | MEG channels: 306 | EEG channels: 0 | Source space: Surface with {self['nsource']} vertices | Source orientation: Free>

注意，正算子会去除太靠近颅骨内表面的顶点。例如，这里我们使用的顶点从516个变成了474个。对于许多函数，例如mne.compute_source_morph()，重要的是好好传递fwd['src']或inv['src']参数，以便充分考虑这种自动删除。

print(f"Before: {src}")
print(f'After:  {fwd["src"]}')

Before: <SourceSpaces: [<surface (lh), n_vertices=155407, n_used=258>, <surface (rh), n_vertices=156866, n_used=258>] MRI (surface RAS) coords, subject 'sample', ~28.7 MiB>
After:  <SourceSpaces: [<surface (lh), n_vertices=155407, n_used=235>, <surface (rh), n_vertices=156866, n_used=239>] head coords, subject 'sample', ~28.7 MiB>

我们可以通过访问fwd的内容来了解包含增益矩阵的这个numpy数组：

leadfield = fwd["sol"]["data"]
print(f"Leadfield size : {leadfield.shape[0]} sensors x {leadfield.shape[1]} dipoles")

Leadfield size : 306 sensors x 1422 dipoles

为了在皮质方向约束下提取包含源空间fwd['src']对应的正算子，我们可以使用以下方法：

fwd_fixed = mne.convert_forward_solution(
    fwd, surf_ori=True, force_fixed=True, use_cps=True
)
leadfield = fwd_fixed["sol"]["data"]
print(f"Leadfield size : {leadfield.shape[0]} sensors x {leadfield.shape[1]} dipoles")

    Average patch normals will be employed in the rotation to the local surface coordinates....
    Converting to surface-based source orientations...
    [done]
Leadfield size : 306 sensors x 474 dipoles

这相当于下面的代码：

import numpy as np
n_dipoles = leadfield.shape[1]
vertices = [src_hemi['vertno'] for src_hemi in fwd_fixed['src']]
stc = mne.SourceEstimate(1e-9 * np.eye(n_dipoles), vertices)
leadfield = mne.apply_forward(fwd_fixed, stc, info).data / 1e-9

要将正算子保存到磁盘，可以使用mne.write_forward_solution()并用mne.read_forward_solution()读取。这个文件应该以-fwd.fif结尾。

还有一些关于用模板MRI脑创建EEG正算子的东西在这里不展示了，可以参考链接进行阅读。