CN1883386A - 基于机器学习的立体核磁共振脑图像弹性配准方法 - Google Patents

Abstract

本发明涉及一种基于机器学习的立体核磁共振脑图像弹性配准方法，用机器学习的方法得到参考图像中每个点上计算属性向量的最佳尺度，由此得到的最优属性向量和周围各点的属性向量之间保持最大程度的差异性，并和训练样本中对应点上的属性向量保持最大程度的相似性；根据图像中每个点上属性向量的显著性和一致性条件，定义一个评价关键点的标准，通过机器学习的方法自动地层次化地在配准的各个阶段挑选关键点，避免配准过程陷入局部极小值点；最后将基于机器学习的框架和现有的配准算法相结合，完成立体核磁共振脑图像弹性配准。本发明对真实MR图像或模拟MR图像均能提高配准的精确度和鲁棒性，进而为后续的临床应用的可行性和准确性奠定基础。

Description

基于机器学习的立体核磁共振脑图像弹性配准方法

技术领域

本发明涉及一种基于机器学习的立体核磁共振脑图像弹性配准方法，借助机器学习的方法，在立体脑图像的每个点上学习最优属性向量，用以准确地表示该点的特征，并且层次化地挑选图像中的关键点，从而提高弹性配准的精度和鲁棒性。本发明可以为后续的图像融合、精确定位病灶、手术计划的制定及疗效跟踪等临床应用奠定基础，涉及图像弹性配准，机器学习，立体核磁共振脑(MR)图像等领域。

背景技术

医学图像配准具有很重要的临床应用价值。对使用各种不同或相同的成像手段获得的医学图像进行配准不仅可以用于医疗诊断，还可以用于手术计划的制定、放射治疗计划的制定、病理变化的跟踪和治疗效果的评价等各方面。

立体脑图像的非线性配准是近年来医学图像处理领域的研究热点。所谓图像配准实际上是指寻求两幅图像间一对一映射关系，也就是说，要将两幅图像中对应于空间同一位置的点联系起来。这里的映射一般称为变换，在二维空间中表现为二维变换，在三维空间中表现为三维变换。实际配准过程中，根据不同的特点和要求既可以采用简单的刚体变换，也可以采用较复杂的弹性形变。本发明主要针对的是三维图像之间的弹性配准问题。

大多数的配准算法可以分为基于控制点和基于像素的两大类。基于控制点的算法需要人工标记控制点，例如在颅骨中嵌入螺钉，在皮肤上做标记等。这类算法的优点在于简单快速，不需要复杂的优化算法，而且精度较高。而缺点是需要人工干预，很难实现自动化，而且受人的主观影响较大。基于像素的方法可以完全自动化地对两幅图像进行配准，并充分利用图像的统计信息，因此是目前研究的比较多的一类方法。该类方法又可以细分为基于像素和基于属性向量的两个分支。前者仅仅利用到图像的灰度信息，速度快，但是精确度不够。后者在图像的每个点上定义一个属性向量，通过该特征可以准确地匹配在另一幅图像中的对应关系。

HAMMER(Hierarchical Attribute Matching Mechanism for Elastic Registration的缩写)是一种基于属性向量的弹性配准算法(D.Shen and C.Davatzikos.“HAMMER：Hierarchical attribute matching mechanism for elastic registration.”IEEE Trans on Med.Imaging，2002，vol.21，pp.1421-1439)，它的特点是在每个点上定义一个属性向量，着重于解剖结构之间的对应关系，而不像同类型的其他算法那样仅把灰度值作为唯一属性来确定两幅图像之间的点对应关系。这样，HAMMER的配准精度和鲁棒性要明显的优于基于灰度值的算法。HAMMER中的属性向量包括灰度值、边界信息以及各种在该点的邻域内计算出来的能够反映该点解剖结构信息的几何不变矩(Geometric Moment Invariants，GMI)。然而HAMMER的不足之处在于：计算GMI的邻域大小是事先确定的，而且对于图像中的每一个点都在同样大小的邻域内计算GMI。在边界信息丰富的区域，例如大脑皮层(cortical)和脑室(ventricle)的角点附近，在较小的邻域内计算得到的GMI最能反映这些点的解剖结构信息；相反，在白质(White Matter)区域，图像灰度变化缓慢，这时如果还是在小的邻域内计算GMI，就会使得这些点的属性向量非常接近以至难以区分。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于机器学习的立体核磁共振脑图像弹性配准方法，提高配准算法的精确度和鲁棒性，进而为后续的临床应用的可行性和准确性奠定基础。

为实现这一目的，本发明首先提出用机器学习的方法得到参考图像中每个点上计算属性向量的最合适的球形邻域大小，由此得到最优属性向量。最优属性向量不仅要和周围各点的属性向量之间保持最大程度的差异性，而且和训练样本中对应点上的属性向量保持最大程度的相似性。其次，本发明提出根据图像中每个像素点上属性向量的显著性和一致性条件，定义一个评价关键点的标准，通过机器学习的方法自动地层次化地在配准的各个阶段挑选关键点，从而避免了配准过程陷入局部极小值点。最后，在图像配准阶段，本发明把基于机器学习的框架和现有的配准算法(HAMMER)相结合，完成立体核磁共振脑图像弹性配准，提高配准结果的精确度。

本发明的方法包括如下具体步骤：

1、训练数据的准备

这个步骤的主要目的是建立参考图像与各个图像之间的对应关系。首先，采集足够的立体核磁共振脑图像建立训练样本集，训练样本集图像中的每个点上应该包括所有可能在图像配准中使用的属性向量。然后从训练样本集中挑选一个样本作为参考图像。对于训练样本集中的每个个体图像，先将其与参考图像进行线性配准，然后再用现有的弹性配准算法进行非线性配准，得到两幅图像之间的位移场，再通过位移场得到两幅图像点和点之间的对应关系。重复上述过程，得到训练集中所有个体图像与参考图像之间的对应关系。

2、最优属性向量的确定

根据参考图像的局部信息，本发明提出用机器学习的方法确定参考图像中每个点的合适的球形邻域大小，即每个点的最佳尺度，由此计算得到参考图像中每个点上的最优属性向量。其中，最优属性向量必须满足三个条件：1)最优属性向量必须和训练样本中对应点上的属性向量尽可能地保持相似性；2)最优属性向量和周围各点的属性向量之间的差异性为最大；3)最优属性向量所对应的最佳尺度的分布在图像空间中要保证平滑。最后，本发明通过属性向量分布的熵构造了反应上述三个条件的能量函数，并用基于梯度下降的最优化的方法计算得到参考图像中每个点上的最佳尺度。

3、选择关键点

根据参考图像中每个点上的最优属性向量的分布直方图可以确定最优属性向量的显著性和一致性，由此定义一个评价函数M(x)＝Sal(x)/Con(x)。其中，分子Sal(x)为显著性指标，分母Con(x)为一致性指标。在每个点上计算M(x)的值，然后根据M(x)的大小排序，来提取参考图像中的关键点。

4、图像配准

对待配准的个体图像先进行线性配准，使其大致地变形到参考图像，然后在个体图像每个点上根据参考图像中使用过的属性向量，计算得到个体图像中每个点上对应不同球形邻域的所有属性向量，挑选个体图像中的关键点并使个体图像中的关键点在图像配准过程中保持不变；在弹性配准的每次迭代求解过程中，先挑选参考图像中的关键点，并在个体图像中逐一根据属性向量间的欧氏距离进行相似性比较，确定其在个体图像中的对应点；在弹性配准的每次迭代过程的最后，根据参考图像和个体图像中关键点的对应关系生成两幅图像的位移场，根据该位移场对个体图像进行变形，并把插值后的新图像作为下一次迭代过程的输入，直到满足迭代的终止条件，完成立体核磁共振脑图像的弹性配准。

本发明的方法不管是在真实MR图像还是在模拟的MR图像上，配准结果都得到了很大的提高。在对真实MR图像进行的实验中，本发明在大脑皮层附近的配准精度得到了明显的提高；在对模拟数据进行的实验中，本发明的配准误差可以由原来的0.95mm下降到0.66mm，准确性提高了约30.5％。对于特定的脑结构(例如海马趾)，配准后的结果和参考图像的重合度从原来的92.8％提高到93.7％，与此同时体积误差从原来的2.6％下降到1.2％。

附图说明

图1为本发明基于机器学习的弹性配准算法的框架。

图2显示的是从高分辨率到低分辨率的最佳尺度的分布情况。(a)表示是参考图像中的一个切片，(b)-(d)分别表示高分辨率、中分辨率及低分辨率下最佳尺度的分布。(e)中显示的最小的尺度为4(深色)，最大的尺度为24(浅色)。图2最右边的颜色表指示了最佳尺度的大小关系。

图3显示了7个有代表意义的点，分别是脑回的尖部，脑沟的根部、脑室的角点及胼胝体边界上的一个点。(a)-(c)分别表示的是高分辨率、中分辨率及低分辨率。图中圆圈的大小代表相应的最佳尺度，红色实线圆圈表示尺度从4到8的球形邻域；绿色点虚线圆圈表示尺度从8到15的球形邻域；蓝色长虚线圆圈表示尺度为15以上的球形邻域。

图4显示了一个说明最优属性向量的优越性的例子。图4(a)中的红色十字坐标表示参考图像中脑沟根部的一个点，图4(b)中的红色星号和紫色圆点都表示个体图像中脑沟的根部点，但前者是参考图像在个体图像中正确的对应点，而后者则不是。图4(a)和(b)中的实线圆圈和虚线圆圈分别表示小的和大的球形邻域，邻域内的图像内容被抽取到图4(e)中显示。图4(c)和(d)显示的是参考图像中的红色十字坐标点的属性向量和个体图像中所有点的属性向量之间的相似度图，图中矩形框内的深色区域表示相似度大，浅色区域表示相似度小。图4(c)是没有经过本发明改进的结果，而图4(d)是经过本发明优化后的结果。

图5显示了配准过程中关键点随着迭代次数逐渐增加的情况。图5(a)-(d)分别是参考图像在配准初始阶段、中间阶段及最终阶段中的关键点的分布情况，图5(e)显示的是浮动图像中的关键点。

图6显示了本发明的配准结果和没有经过机器学习优化的配准结果。

具体实施方式

以下结合附图和实施例对本发明的技术方案作进一步描述。

本实施例主要包括训练和图像配准两个部分，如图1所示。在训练阶段，主要的任务包括在参考图像的每个点上确定最优属性向量及挑选关键点。在图像配准阶段，把在训练阶段得到的结果和现有医学图像配准算法相结合(例如HAMMER算法)，确定两幅图像之间的对应关系。以下便结合HAMMER算法对本发明的具体实施方法做进一步的详细描述。

实施例采用的MR图像均为T1加权的三维图像。整个发明实现过程如下：

1.训练数据的准备

●采集MR脑图像样本，生成训练集。在本实施例中，样本空间包括18个不同个体的MR图像。

●挑选一个图像作为参考图像，把其余的个体图像线性配准到该参考图像。现有的线性配准算法有FSL，ITK等(均为公开的源代码)。

●用HAMMER算法对线性配准后的个体图像做弹性配准，得到精确的点对应关系。HAMMER算法提供有公开的可执行代码，可在算法的提出者(Dingang Shen)的主页https：//www.rad.upenn.edu/sbia/rsoftware.html下载。

●对参考图像和经过线性配准后的个体图像上的每个点在不同大小的球形邻域内计算 $\mathrm{GMI}\stackrel{\→}{G}\left(S_x\right).$ 其中S＝{i|i＝4*j，j＝1，2，3，...，6}。计算的详细步骤如下：

假设图像已经被分割成白质(White Matter，WM)、灰质(Gray Matter，GM)、脑室(Ventricle)、脑脊髓液(cerebrospinal fluid，CSF)四个部分。先用公式(1)计算图像的三维(p+q+r)阶矩：

$M_{p,q,r}=\underset{{\left(x_1\right)}^2+{\left(x_2\right)}^2+{\left(x_3\right)}^2<S^2}{\&Integral;\&Integral;\&Integral;}{x}_1^p{x}_2^q{x}_3^r{f}_{\mathrm{tissue}}(x_1,x_2,x_3){\mathrm{dx}}_1{\mathrm{dx}}_2{\mathrm{dx}}_3---\left(1\right)$

其中，f_tissue表示组织类别(例如WM，GM，VN，CSF)，s表示计算半径。然后每个点上的属性向量(GMI)就可以用公式(2)计算得到：

GMI₁＝M_0，0，0

GMI₂＝M_2，0，0+M_0，2，0+M_0，0，2

${\mathrm{GMI}}_3=M_{\mathrm{2,0,0}}{M}_{\mathrm{0,2,0}}+M_{\mathrm{2,0,0}}{M}_{\mathrm{0,0,2}}+M_{\mathrm{0,2,0}}{M}_{\mathrm{0,0,2}}-M_{\mathrm{1,0,1}}^2-M_{\mathrm{1,1,0}}^2-M_{\mathrm{0,11}}^2$

${\mathrm{GMI}}_4=M_{\mathrm{2,0,0}}{M}_{\mathrm{0,2,0}}{M}_{\mathrm{0,0,2}}-M_{\mathrm{0,0,2}}{M}_{\mathrm{1,1,0}}^2+{2M}_{\mathrm{1,1,0}}{M}_{\mathrm{1,0,1}}{M}_{\mathrm{0,1,1}}-M_{\mathrm{0,2,0}}{M}_{\mathrm{1,0,1}}^2-M_{\mathrm{2,0,0}}{M}_{\mathrm{0,1,1}}^2---\left(2\right)$

2.最优属性向量的确定

构造能量函数，如公式(3)所示：

$E=\underset{x}{\mathrm{\&Sigma;}}(-E_1(x,S_x)+\mathrm{\&alpha;}{E}_2(x,S_x)+\mathrm{\&beta;}{E}_3(x,S_x))---\left(3\right)$

其中，E₁(x，S_x)表示在大的图像邻域内(半径为5的球形邻域)GMI分布直方图的熵，E₂(x，S_x)表示在小的图像邻域内(半径为2的球形邻域)GMI分布直方图的熵， $E_3(x,S_x)={\mathrm{\&Sigma;}}_{y\&Element;r_x}{(S_x-S_y)}^2$ 表示相邻点上最优尺度的差异。最小化该能量函数，可以得到每个点上的最佳尺度S(x)。本发明采用马尔可夫随机场来进行优化。把公式(3)的前两项(E₁，E₂)当作似然能量，第三项(E₃)作为平滑能量。考虑到三维图像的计算量问题，本发明采用的是局部最优化方法(ICM)来优化马尔可夫随机场。如果有训练样本，则α＝1.0，否则α＝0。β在本实施例中的取值为1.0。

通常为了增加算法的鲁棒性，许多配准算法都使用了多分辨率技术。因此，本发明对每个分辨率都进行了同样的计算，得到了相应的最佳尺度的分布。图2从左到右显示的是从高分辨率到低分辨率的最佳尺度的分布情况。图中最小的尺度为4，最大的尺度为24。图2最右边的颜色表指示了最佳尺度的大小关系。从图中可以看出最佳尺度的分布符合MR脑图像的解剖结构特征。小的尺度都分布在皮层等边界复杂的区域，然后尺度由外向内逐渐增大，在白质等灰度变换平缓的区域达到最大。

图3显示了本发明挑选的7个有代表意义的点，分别是脑回的尖部，脑沟的根部、脑室的角点及胼胝体边界上的一个点。图中圆圈的大小代表相应的最佳尺度，红色实线圆圈表示尺度从4到8的球形邻域；绿色点虚线圆圈表示尺度从8到15的球形邻域；蓝色长虚线圆圈表示尺度为15以上的球形邻域。为了显示方便，低分辨率和中分辨率的图像都上采样到高分辨率图像的大小。

下面给出一个具体的例子说明了本发明在本实施例中的优越性。图4(a)中的红色十字坐标表示参考图像中脑沟根部的一个点，图4(b)中的红色星号和紫色圆点都表示个体图像中脑沟的根部点，但前者是参考图像在个体图像中正确的对应点，而后者则不是。图4(a)和(b)中的实线圆圈和虚线圆圈分别表示小的和大的球形邻域，邻域内的图像内容被抽取到图4(e)中显示。图4(c)和(d)显示的是参考图像中的红色十字坐标点的属性向量和个体图像中所有点的属性向量之间的相似度图，图中矩形框内的深色区域表示相似度大，浅色区域表示相似度小。图4(c)是没有经过本发明改进的结果，而图4(d)是经过本发明优化后的结果。

3.选择关键点

在步骤2中得到最优属性向量 ${\stackrel{\&RightArrow;}{G}}_b=\left\{g_i\right\}$ 后，按照公式(4)计算每个点x上的最优属性向量的显著性值Sal(x)。

$\mathrm{Sal}(s,x)=L(s,x)\&CenterDot;S_x\&CenterDot;\underset{i}{\mathrm{\&Sigma;}}\left|\right|\frac{{\&PartialD;p}_i(s,x)}{\&PartialD;s}{|}_{S_x}\left|\right|---\left(4\right)$

其中 $L(s,x)=-\underset{i}{\mathrm{\&Sigma;}}{p}_i(s,x)\log p_i(s,x),$ p_i(s，x)表示在以点x为中心，s为半径的邻域中，某个属性向量GMI_i出现的概率。需要说明的是公式中出现的s为即为在步骤2中得到的最优属性向量对应的最佳尺度。

最优属性向量的一致性值即为公式(3)中的E₂项，只不过此时的s为最优半径S_b。计算法的方法类似。

最后，在参考图像的每个点上计算M(x)＝Sal(x)/Con(x)，并根据M(x)的大小进行排序。M(x)的值越大，表示该点附近的局部特征越明显，这样的话对应关系也就越容易确定。

图5显示了配准过程中关键点随着迭代次数逐渐增加的情况。图5(a)-(d)分别是参考图像在配准初始阶段、中间阶段及最终阶段中的关键点的分布情况，图5(e)显示的是浮动图像中的关键点。在图中可以看到，初始阶段的关键点分布比较稀疏，集中在脑室、脑沟根部及脑回的顶端，这些恰恰就是局部特征最突出的区域。这些区域中的点的对应关系最容易准确地找到，这在变形的初期确定位移场的大致形状是非常重要的。随着迭代次数的增加，越来越多的边界点变成了关键点。到迭代的最后，所有的边界点都变成了关键点，如图5(d)所示。

4.图像配准

4.1数据准备

对于一幅待配准的图像，首先用线性配准的算法(FSL)把该图像大致变换到参考图像空间。由于本发明的训练结果只针对参考图像，因此对于个体图像来说必须把参考图像空间中所有用到的GMI都计算一遍，即计算 $\mathrm{GMI}\stackrel{\&RightArrow;}{G}\left(S_x\right),$ S＝{i|i＝4*j，j＝1，2，3，...，6}。其次，还需要计算个体图像中每个点上的M(x)值，并按照从大到小的顺序排序，以便在配准的过程中层次化地选择个体图像中的关键点。

4.2弹性配准过程

根据M(x)值的大小挑选个体图像中的关键点{y_i|1≤i≤N_s}，N_s为挑选的关键点个数，约为整个脑图像像素(不包括背景点)个数的15％。个体图像中的关键点在配准过程中保持不变。

选择参考图像中的关键点{x_i|1≤i≤N_T}，N_t为参考图像中的关键点个数。在配准的初始阶段，N_t较小，约为整个脑图像像素(不包括背景点)个数的15％。随着迭代次数的增加，N_t逐渐增大，到迭代收敛时，所有的脑图像像素都可以作为关键点。

4.3确定对应关系并生成位移场

●确定个体图像中的关键点到参考图像的对应关系。即对于任意个体图像中的关键点y_i，在一定的邻域内寻找相似度最高的参考图像中的关键点x_i(在本实施例中，相似度的阈值为0.8)。如果x_i存在，则建立了从h(x_i)到y_i的对应关系。

●对于参考图像中的任意一个关键点x_j，如果在上一步中已经建立了对应关系，则该点发生位移的方向就是沿着h(x_j)到y_i；否则，在x_j的周围搜索个体图像中和x_j相似度最大的点，位移的方向就指向个体图像中最相似的点。注意，如果相似的点不止一个，那么位移的方向为这些力的合力方向。

●用高斯滤波器平滑位移场。

如果迭代终止条件没有满足，重新回到步骤4.3。

本发明的效果

本实施例中，本发明对模拟的MR脑图像及真实图像进行了实验。所有实验都在微机上完成(Pentium IV，3.0GHz)，而且都使用相同的参数。

首先，本发明对18幅老年个体的MR脑图像进行了配准。为了对比实验结果，本发明对同样的18幅图像用原来的HAMMER算法也进行了配准。通过比较发现本发明的配准结果在复杂的大脑皮层附近有明显的提高，如图6所示。在图中的红圈中，本发明还能保持解剖结构的对应关系，但是HAMMER算法的配准结果出现了明显的错误。

其次，本发明还在模拟MR脑图像上进行了对比实验。在本实施例中，先用文献(Z.Xue，D.Shen et al，″Statistical Representation and Simulation ofHigh-Dimensional Deformations：Application to Synthesizing Brain Deformations″.MICCAI，Palm Springs，California，USA，Oct 26～29，2005.)中的算法把参考图像变换成了不同的个体图像，同时得到了参考图像和个体图像之间的位移场。在这个实验中，用本发明把不同个体图像分别变换到参考图像空间。通过验证配准后得到的位移场与事先已知的位移场之间的误差，就可以比较配准结果的精确程度。改进后的算法要明显好于原来的HAMMER算法。在没有训练样本的情况下，配准精度提高了约12.6％；如果所有的模拟数据一半用于训练，另一半用于测试，则在测试集上的精度可以提高30.5％。

Claims (1)

1、一种基于机器学习的立体核磁共振脑图像弹性配准方法，其特征在于包括如下步骤：

1)训练数据的准备：采集足够的立体核磁共振脑图像建立训练样本集，训练样本集图像中的每个点上应该包括所有可能在图像配准中使用的属性向量，然后从训练样本集中挑选一个样本作为参考图像；对于训练样本集中的每个个体图像，先将其与参考图像进行线性配准，然后再用现有的弹性配准算法进行非线性配准，得到两幅图像之间的位移场，再通过位移场得到两幅图像点和点之间的对应关系，由此得到训练集中所有个体图像与参考图像之间的对应关系；

2)最优属性向量的确定：根据参考图像的局部信息，用机器学习的方法确定参考图像中每个点的合适的球形邻域大小，即每个点的最佳尺度，由此计算得到参考图像中每个点上的最优属性向量；其中，最优属性向量必须满足三个条件：1)最优属性向量必须和训练样本中对应点上的属性向量尽可能地保持相似性；2)最优属性向量和周围各点的属性向量之间的差异性为最大；3)最优属性向量所对应的最佳尺度的分布在图像空间中要保证平滑；通过属性向量分布的熵构造反应上述三个条件的能量函数，并用基于梯度下降的最优化方法计算得到参考图像中每个点上的最佳尺度；

3)选择关键点：根据参考图像中每个点上的最优属性向量的分布直方图来确定最优属性向量的显著性和一致性，由此定义一个评价函数M(x)＝Sal(x)/Con(x)，其中，分子Sal(x)为显著性指标，分母Con(x)为一致性指标；在每个点上计算M(x)的值，然后根据M(x)的大小排序，来提取参考图像中的关键点；

4)图像配准：对待配准的个体图像先进行线性配准，使其大致地变形到参考图像，然后在个体图像每个点上根据参考图像中使用过的属性向量，计算得到个体图像中每个点上对应不同球形邻域的所有属性向量，挑选个体图像中的关键点并使个体图像中的关键点在图像配准过程中保持不变；在弹性配准的每次迭代求解过程中，先挑选参考图像中的关键点，并在个体图像中逐一根据属性向量间的欧氏距离进行相似性比较，确定其在个体图像中的对应点；在弹性配准的每次迭代过程的最后，根据参考图像和个体图像中关键点的对应关系生成两幅图像的位移场，根据该位移场对个体图像进行变形，并把插值后的新图像作为下一次迭代过程的输入，直到满足迭代的终止条件，完成立体核磁共振脑图像的弹性配准。

Images