具有全数据或部分数据的多模态医学图像配准

李正伟

（成都理工大学工程技术学院，四川乐山 614007）

图像分析涉及到图像分割、重构、分类和检测等多种任务[1]。通过比较或融合两个或多个图像可以寻求重要信息，图像之间不可避免的错位使得图像配准成为研究重点。图像配准的目的是在同一场景的两个或多个图像之间找到最佳的空间对齐，涉及图像配准的图像可以使用相同的(单模态)或不同的(多模态)成像技术捕获[2]。在计算机断层扫描(CT)和磁共振成像(MRI)的情况下，在轴平面上拍摄器官的完整横截面图像。因此，将新扫描与原始扫描进行比较将需要具有部分重叠的图像配准。由于部分重叠的两幅图像之间缺乏一一对应关系，因此使用最常用的配准方法很难获得精确的配准[3]。

单模态医学图像配准方法的目的是对同一成像技术(如CT)得到的图像进行配准。文献[4]提出了基于快速傅里叶变换(FFT)技术，该技术可寻求参考图像和感测图像之间的最佳平移、旋转和缩放参数。这种相位相关技术的扩展后来发展为基于Fourier-Mellin 变换(FMT)的图像配准方法[5]。文献[6]提出了改进的FMT 方法，消除了直角坐标到对数极坐标的转换。评估并验证了该方法对噪声的准确性和鲁棒性。

多模态图像配准算法侧重于寻找使用各种成像方式(如CT 和MRI)生成图像之间的对应关系，并从不同医学成像方式的融合中提供密集的视觉信息。文献[7]实现了一种基于图像强度全局和局部变化的通用图像配准方法。该算法在全数据的多模态图像和部分数据的单模态图像上进行了实验，但由于采用弹性模型估计变形场而存在特殊的缺陷。基于信息论的度量，如互信息(MI)[8]及其变化，已广泛用于基于强度的模态间配准，但其存在收敛速度慢、精度低和对实现决策敏感等缺陷。文献[9]提出了标准化互信息(NMI)度量，当边缘熵之和的增长速度快于联合熵时，可以处理由于失准而导致的MI 值的增长。

考虑到CT和MRI是两种广泛使用的成像方式，可以产生具有可比空间分辨率的横截面图像，因此MRI具有更高的对比度分辨率，而CT扫描则相对较低。常用的多模态配准方法往往由于缺失部分数据导致配准无效，降维算法（即流形学习）常用于图像和数据处理领域，其目的是研究嵌入高维空间的低维流形[10]。

该文改进了文献[11]中利用Laplacian 特征映射来表示多模态图像配准的方法，实现了多模态到单模态的转换。在多模态医学图像上直接应用单模态配准技术以及恢复强放缩、旋转和平移，并且对完全重叠和部分重叠两种类型的图像都适用。此外，还提出了一种快速且易于使用的对齐技术，以此补偿由主向量反射引起的随机符号歧义。对齐技术有助于更复杂的优化算法以及基于强度度量的应用。提出的强度变换解决了部分重叠的多模态医学图像配准问题。从实验的角度出发，提出了一套定性和定量的分析方法来检验系统的性能和准确性，其中临床输入图像同时受到多个自由度（平移、旋转、放缩）的约束，在不同的多模态配准中获得了较低的平均绝对误差。

1 数据降维

将高维数据映射到低维空间的变换降维技术可分为线性和非线性方法。该文引入主成分分析(PCA)[12]和Laplacian 特征映射[11]分别用于线性和非线性技术。光谱嵌入可作为寻找与特定矩阵的上（或下）特征值相关联的一组特征向量。

考虑集合x1,…,xn∈M的流形M表面上有n个采样数据点，其中xi∈RD。寻找一组对应点y1,…,yn，其中yi∈RD(d＜＜D)，用于保留X=[x1,...,xn]的特定相关信息。

1.1 主成分分析(PCA)

PCA 的目的是通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量[12]。在PCA 中，假设数据集的方差是对数据集中存储信息量的度量，则方差越大，存储的信息越多。

给定矩阵X∈RD×n，其列为D维数据点，则矩阵Xc∈RD×n定义为中心数据点的矩阵（每行的样本均值为零）。投影由Y=AXc给出，其中A为d×n矩阵。PCA 优化目标函数为：

其中，A为d×D正交矩阵。Y的协方差矩阵为：

利用Lagrange 乘数技术，将最大化问题简化为求解协方差矩阵V的特征值问题：

根据线性代数，任何MMT形式的矩阵都是对称矩阵。因此，V为对称矩阵，由此得出结论：根据谱定理，V为正交可对角化矩阵，并且只有实的非负特征值。正交可对角化矩阵V具有正交的非零实特征向量。

矩阵V具有最大数量D个特征值和大小为D×1 的D个特征向量。设λ1≥λ2≥···≥λD≥0（按降序排列），相应的正交特征向量为。矩阵V的特征向量是数据点的主成分(或主方向)，数据点沿主分量(或主方向)方差最大。图1 给出了样本合成数据集的主要组成部分。这些主成分构成一组标准基向量，将高维数据点投影到低维空间，同时保留大部分信息。所需的线性变换矩阵A由与前d个最大特征值相关的前d个主成分组成[13]。

图1 样本合成数据集的主成分

1.2 Laplacian特征映射

Laplacian 特征映射的目的是找到高维流形的低维嵌入作为结构表示。在保留局部信息的同时，提取结构信息并将其嵌入到低维空间中。结构表示得益于保留局部性和结构等价性。保留局部性是指同一流形的两个相似块映射到新的坐标后，其结构表示是相似的。结构等价性是指不同流形的两个相似块的结构表示相似。

Laplacian 特征映射基于Laplacian 图，其算法包括3 个主要步骤[11]：1）构造邻接图；2）选择权值；3）寻找特征映射。每个步骤都涉及选择适当的实现方式并为多个参数确定适当的值。由Laplacian 特征映射所保证的保留局部性和结构等价性是该文方法的关键决定因素。

2 多模态图像配准

该文提出的多模态到单模态变换的目的是将多模态输入图像变换成相似强度的坐标系。通过提取比单个像素强度更多的信息来研究小图像块的结构。由于来自不同模态的图像的内部结构相似，因此，转换将导致多模态图像的相似结构表示。多模态到单模态变换的图像配准如图2 所示。

图2 多模态到单模态变换的图像配准

2.1 构建高维空间

在图像配准问题中，只存在参考和感测两种图像。在高维空间中，通过从每个输入图像中提取小块来构造点云。高维空间的每个维度专用于每个补丁的一个图像像素。考虑大小为s×s的每个补丁都比前一个补丁移动了一个像素，对于大小nr×nc的图像，高维流形在s2=D维空间中包括nr×nc=N个点。使用点的坐标以矩阵格式表示这些点集，对于大小为200×200、块大小为3×3的图像，矩阵表示为40 000×9的矩阵，使用图像的小块和相应的点云矩阵表示构建高维点云的示例，如图3 所示。对于大小为N像素、块为D个像素的灰度图像，矩阵大小为N×D。

图3 高维空间中点云的矩阵表示

D的选择取决于图像的应用和大小。当选择相对较小的D时，嵌入的流形类似于原始图像。当D选择较大的值会导致图像模糊。不同D值时T2-MRI图像切片及其嵌入流形的图像示例，如图4所示。

图4 不同D值时T2-MRI图像切片及其嵌入流形的图像

图4(a)给出了原始T2-MRI 图像的切片。利用Laplacian 特征映射对图像进行了研究。图4(b)～(d)说明了嵌入流形后参数D的选定值对渲染图像清晰度的影响。随着D的增加，嵌入流形的清晰度降低。从图4(d)可以看出，使用15×15 的补丁会导致大脑结构内部细节的丢失。

2.2 流形学习

Laplacian 特征映射降维需要构造一个无向邻接图G=(V,E)，图中包含高维空间RD中的一组点xi∈V和一组指定邻域连通度的边E。对于给定的点集，k近邻通过在高维空间中找到每个点的k近邻并指定该点与其每个近邻之间的边来构造图[14]。通过构造有向图，当节点i属于节点j的k近邻集合，而节点j不在节点i的k近邻集合时，在Laplacian 特征映射方法中，通过在有单向边的地方添加更多边（称为对称kNN 图）或删除单向边（称为相互kNN 图）来生成无向图[15]。

为了构造邻接图，该文将k近邻与对称kNN 图相结合。在降维过程中，参数j的选择只是内存速度与确保具有单个连接之间的权衡。由于医学扫描使用平滑的强度级别范围，从而获得平滑的高维点云，因此，对于200×200 像素左右的医学扫描，能够获得k=20 的单个连接分量。如果输入图像有噪点，可能导致多个连接组件的图，因此使用低通滤波器或增加k。

选择每条边的权重决定了连接的能力，该文借鉴文献[16]中图像的Laplacian 算子与流形上的Laplacian Beltrami(L-B)算子的可比性。为了得到L-B算子的最佳近似，该文使用以下热核加权方式：

其中，σ∈R为唯一的参数。如果两个节点i和j通过一条边连接，则wij=1，否则wij=0。

相邻点之间的距离是局部信息的度量方法。因此，通过选择合适的参数σ值来保持各向异性指数加权方案。在式（4）的指数加权方案中，邻近点的权重高于远处点的权重。当σ与||xi-xj||相比相对较大时，所有连接的权值都接近于1，从而产生无权邻接图。相反，相对较小的σ将导致不显著的边缘，进而使得特征值的收敛失败。该文提出构造与σ相关的权重矩阵S(σ2)为：

对于多个σ值，绘制S曲线与对数刻度(σ2)的关系图。S曲线在σ=0 和σ=∞处有两个渐近线，并用于选择σ，其中图形在总范围的上半部分呈线性。为每幅图像计算S曲线增加了额外的处理时间。该文选择热核的方差等于所有边的最大平方距离：

所有边的权重都限制在e-1/2≈0.606 和1 之间。在减少计算时间的同时，嵌入流形和配准结果均有效且可靠。

最优嵌入在保留局部信息的同时将点xi映射到点yi，即连接点xi和xj在分别映射到yi和yj之后将保持尽可能靠近。最小化以下目标函数：

对于给定的具有wij项的对称权值矩阵W，该文定义了对角度量矩阵D，使得。度量矩阵D的项是W的列（或行）之和。该文还将Laplacian 矩阵定义为L=D-W。如果集合V包含N个点，则W以及D和L均为稀疏N×N对称矩阵。

通过添加正交约束yTD1=0 来消除平凡解，并约束yTDy=1 来消除嵌入中的任意比例因子，将式(7)的最小化问题简化为：

由广义特征值问题的最小特征值解得到式（8）的解向量y：

式（9）为L-B 算子的特征分解或谱分解，将矩阵用其特征值和特征向量表示，特征值的集合为L的谱。L-B 算子特征值和特征向量具有等距不变量性质，即如果流形没有拉伸（例如弯曲成额外的维度），光谱值不会改变。因此，两个方向不同的流形将具有相似的谱表示。式（4）和（7）表明，L为实对称且半正定矩阵。因此，本征函数1 和本征值λ=0是式（9）的平凡解。特征值为0 的多重性与图的连通分量的数目相关，对于具有多个连通分量的图，L为块对角矩阵，其中每个块表示连通分量的Laplacian矩阵。

考虑到Laplacian 矩阵的性质，需要去掉与特征值为0 的所有特征向量，并使用下一个具有最小非零特征值的d个特征向量来将流形嵌入到低d维欧几里得空间中。嵌入流形的形式为N×d矩阵Y=[y1,…,yd]，其中第i行显示第i点的嵌入坐标。

该文使用ARPACK 软件求解广义特征值问题的d个最小特征值，ARPACK 代表ARnoldi 包可用于逼近几个特征值和对应的特征向量。由于矩阵L为对称矩阵，因此，将Arnoldi 过程简化为Lanczos 过程的变体，即隐式重新启动Lanczos 方法(IRLM)[17]。由特征解算器进行的计算近似使得最小特征值不完全等于零。与大于0 的最小特征值相对应的第一个d特征向量是在低维空间嵌入流形的正交基。此时，尺寸为N×1 的d个特征向量在d维空间中定位嵌入流形的N个点。

输入图像的结构表示是通过重新排序特征向量到原始输入图像的大小。每个特征向量表示具有特定派生度量的输入图像。特征向量的指数越高，结构表示的导数较高。使用特征图像1、6、11 和18 的人脑T1-MRI 的结构表示，如图5 所示。Laplacian 特征映射降维的必要步骤，如算法1 所示。

图5 用Laplacian特征图表示T1-MRI的结构

算法1：Laplacian 特征映射降维（结构表示）

输入：高维空间中N个点的集合X=(x1,…,xN)T∈RN×D

输出：低维空间中N个点的集合Y=(y1,…,yN)T∈RN×d

步骤1：计算每两个数据点之间的距离；

步骤2：选择k和对称kNN 图，构造k近邻构造邻接图；

步骤3：对于任意的xi,xj∈V，定义σ2=max(||xi-xj||2)；

步骤4：采用热核加权方案为每两个邻居点之间的每个边分配权重；

步骤5：构造稀疏实对称N×N矩阵W,D,L；

步骤6：从L中找出连接组件的数量nc；

步骤7：求解(nc+d) 个特征值和相应特征向量的广义特征值问题Ly=λDy；

步骤8：对特征向量进行排序，以递增顺序表示特征值；

步骤9：去掉第一个nc个特征向量；

步骤10：返回剩余的d个特征向量。

2.3 流形对齐

Lanczos 算法以随机初始向量开始初始化，虽然Laplacian 特征映射保留了局部信息，但不能确保嵌入流形在低维空间的对齐。因此，需要进行流形对齐[18]。

PCA 提供由平移和旋转组成刚体变换的参数。由于任何特征向量都是有效的特征向量，因此两个嵌入流形之间的映射可能需要旋转。采用PCA 的刚体变换无法处理这种随机符号模糊问题。如果两个流形的主成分遵循由向量的叉积在垂直主向量之间定义的右手或左手规则，PCA 将使流形具有足够的精度确保对齐。两个嵌入流形的右手规则检查，如图6 所示。其中，参考流形和感测流形的前3 个主方向分别遵循左手和右手规则。

图6 两个嵌入流形的右手规则检查

该文提出了一种快速而直接的方法来检查参考流形和感测流形是否符合右手规则，每两幅对应的特征图像所对应的协方差的符号是它们关联方向的度量。如果它们有负关联，则将-1 与感知嵌入流形的特征图像相乘。当使用该文方法补偿反射误差时，PCA计算刚体变换矩阵R，其对齐中心嵌入流形为：

其中，矩阵A和矩阵B分别为感测流形和参考流形的主方向。由于矩阵A为正交矩阵，根据PCA的性质，旋转矩阵简化为：

通过将两个流形的中心移动到原点，然后旋转感测流形的主方向以对准参考流形的主方向来实现对齐。利用PCA 对感测嵌入流形进行刚体变换，使其与参考嵌入流形对齐。在检查右侧规则符合性后，可表示为：

其中，cs和cr分别为的感测点和参考点。

图7 给出了流形对齐完成预期强度转换的必要性。不同模态使用不同的相对强度水平来显示大脑的同一部分。T1-MRI 和T2-MRI 图像中的脑脊液、白质和灰质表现不同。通过流形学习，还没有达到必要的强度调整。然而，流形学习和流形对齐在两幅图像上执行强度变换，使得它们使用可比较的强度级别。

图7 T1-MRI特征图像1和T2-MRI强度变换特征图像1之间的颜色相似性

图7 将多模态扫描转换为相同强度的坐标系，虽然原始的T1-MRI 和T2-MRI 存在强度变化，但T2-MRI 强度转换后的特征图像与T1-MRI 特征图像使用可比较的相对强度水平，在图像中产生必要的对比度。

2.4 图像配准

该文使用流形学习和流形对齐的多模态到单模态转换提供了多模态扫描的结构表示。这种新的图像表示代替了原始图像，适用于使用单模态配准技术寻找配准参数。然后，使用估计的变换矩阵对原始的多模态扫描进行对齐。该文分别使用基于强度的图像配准和基于FMT[5]的图像配准来求解具有全部和部分数据的图像配准问题。基于强度的配准迭代利用规则步长梯度下降优化器来调整变换参数，使优化沿着均方梯度的极值方向进行。

3 实验分析

为了验证该文提出的强度变换方法的性能，使用Laplacian 特征映射和提出的流形对齐技术，在完全或部分数据可用的情况下，来改善和促进多模态医学扫描配准。

3.1 实验装置

该文实验在Intel(R) Xeon(R) E3-1245 CPU，3.50 GHz，32 Gb 内存的个人计算机上运行Matlab R2013a 环境进行。为了验证所提出方法的有效性，利用Elastix toolbox 进行基于NMI[9]的配准。在实验中，NMI 和MI 度量使用50 个具有直方图的所有像素进行计算。

对BrainWeb 数据库中的T1、T2 和PD 加权MR图像进行了实验。实验中的图像含有3%的噪声和20%的强度不均匀性。同时，还使用了CT 和Vanderbilt 数据库中的不同MR 模式，该数据库从RIRE 数据库进行收集。此外，RIRE 数据库提供了利用对静场不均匀性进行校正而生成的校正图像。所有图像的大小都调整为200×200，以便节省处理时间和内存使用。

3.2 多模态到单模态转换

为了直观地验证所提出的强度变换性能，使用与参考特征图像中相似的相对强度级别来重新绘制感测特征图像。使用来自不同数据集的图像来检验该方法的有效性，该文分别研究了每幅图像的流形，使用大小为3×3 的块得到了9 维空间中包含40 000个点的流形，构造了k=20 近邻的Laplacian 图，并将流形嵌入到三维空间中。以T1 加权MRI 为参考，将其余部分与之对齐。多模态到单模态转换的结果如图8 所示。其中，第一行为人脑的原始MRI 和CT 扫描，第二行流形学习后每种模式的第一幅特征图像；第三行流形对齐后的相同特征图像。T2、PD-MRI和CT 的特征图像的强度映射与T1-MRI 的特征图像配准。在经过多种流形对齐之后，所有模态都使用了与T1-MRI 特征图像1 类似的强度映射。

图8 多模态医学扫描的强度变换

3.3 全数据的多模态配准

为了验证该文方法在全数据多模态图像配准中的有效性，利用RIRE 数据集进行了两组实验。由于目前的方法只适用于二维图像，因此获得了切片对齐。首先，使用相同的参数D=9，k=10 和d=3 分别对每幅图像的高维流形进行研究。然后，采用流形对齐技术将多模态扫描转换为单模态强度坐标系。其次，利用各模态的第一幅特征图像估计刚体配准参数。为此，该文采用了正则梯度下降优化器和均方误差度量作为单模态配准的方法。最后，在原始图像集上引入了配准参数。

使用所提出的方法检查了用多种MR 模式配准CT 扫描的性能。由于没有RIRE 测试数据集的真实信息，还需测量配准前后图像对之间的相似性度量MI。MI量化了参考图像和感测图像之间的统计相关性程度。当前数据集包含10 名患者的CT、T1、T2 和PD-MRI扫描。其中，只有6名患者有经过校正的MRI。在统计误差p＜0.5%内，每对CT-MR 扫描配准前后测量的MI 的均值和标准差（单位：mm），如表1 所示。

表1 每对CT-MR扫描配准前后测量的MI均值和标准差

RIRE 数据库中CT 和PD-MRI 扫描样本的对比如图9 所示。结果表明，即使在配准过程不进行MI值的计算和优化，该方法仍能提高MI 值。

图9 RIRE数据库中CT和PD-MRI扫描样本的对比

为了进一步定量地评价该文方法的性能，利用RIRE 训练数据集，在[-π/4,π/4]范围内随机生成旋转角度，并沿每个轴进行平移以此约束大脑区域保持在坐标系内，进而创建了一个CT（感测图像）的非对齐切片。采用3 种不同的方法对失真图像进行配准：1）文献[11]中提出的基于Laplacian 图的配准方法；2）文献[8]中提出的MI 度量的多模态配准方法；3）该文方法。利用5 个随机重定位点的平均绝对误差(MAE)计算配准误差。由于RIRE 训练数据集只有患者的多模态图像，该文对每对模态重复实验30 次。距离误差的平均值和标准差结果（单位：mm），如表2所示。

表2 每对CT-MR模态扫描的距离误差的平均值和标准差

在表2 中，粗体数字表示每个模态对获得的最佳结果。除带星号数据外，其余结果均在统计误差p＜0.1%内。CT-T1 校正扫描的配准误差均在统计误差p＜0.3%内。CT-T2 校正扫描的配准误差高于MI，但无统计学意义。

3.4 部分数据的多模态配准

在多模态医学扫描与部分数据的配准问题中，从不同模态下获取图像时，感测图像（部分数据）只覆盖部分参考图像（全部数据），因此，需要在较大数据空间内定位小数据。

利用随机参数对遥感图像进行旋转、平移和裁剪，生成局部数据的合成实例。在实验中，在感测图像中加入了额外的放缩。使用所提出的方法对两幅图像进行研究，然后将每幅图像的特征图像传递给FMT 算法。FMT 算法将相应地计算转换参数(平移、旋转和放缩)，利用估计的变换参数对原始的全图像和局部图像进行联合配准。原始图像和配准图像对比结果如图10 所示。

图10 部分数据的多模态图像配准对比

将患者曾经扫描的T1-MRI 与覆盖患者大脑部分的新捕获的T2-MRI（感测）进行配准。两幅图像均在轴向采集，在进行新的扫描（T2）时，患者的成像位置与T1 不同。参考图像和感测图像以及多种方法的配准结果，如图11 所示。

图11 感测图像中的部分数据对合成实例进行配准

4 结束语

该文将所提出的多模态到单模态的变换作为预处理步骤，而不必考虑图像的重叠，这有助于恢复平移、旋转和放缩。在医学图像配准策略的背景下，该方法可作为参数化方法，但不需要对参数进行精确的匹配和赋值。为进一步调整高维流形的光滑性提供了自由度，最终将对特征值问题的收敛速度产生积极的影响。为了实现精确的结构表示，利用所提出的多模态到单模态变换，生成多个特征图像作为结构表示的输出。虽然该文只考虑了人脑图像的刚性变换（除了考虑仿射变换的部分数据的配准外）来评估所提出的变换性能，但对其他图像的配准以及非刚性配准都可以借助此方法。