锥束计算机断层扫描图像对应与配准方法

摘要

本发明公布了一种锥束计算机断层扫描图像对应与配准方法，基于混合测度森林，利用非监督聚类森林生成锥束CT图像之间的超体素之间的稠密对应；再利用基于森林的测度估计锥束CT图像之间的超体素之间的稠密对应与锥束CT图像的配准。训练混合测度森林包括：提取超体素特征；训练初始聚类森林，估计超体素之间的相似度；估计弱标注；基于初始聚类森林迭代增强基于森林的测度；估计柔性一致度；在线测试将超体素特征输入到训练好的混合测度森林，估计超体素之间的相似度与对应，并估计柔性一致度；获取锥束CT图像之间的变形参数并实现配准。本发明可以快速建立三维锥束计算机断层扫描图像对应与配准。

说明书

技术领域

本发明涉及口腔临床医学和计算机视觉技术领域，具体涉及锥束计算机断层扫描图像对应与配准方法。

背景技术

锥束计算机断层扫描(锥束CT)图像在口腔正畸临床中广泛使用，用于评价治疗以及度量组织结构的生长发育。高效与可靠的锥束CT图像配准以及图像之间体素的稠密对应是计算机辅助的术中干预与在线的图像属性迁移诸如特征点定位与分割标签迁移的关键。在医学图像处理领域，对三维图像配准已经进行多年的研究。尽管多年的研究对图像配准有了极大的推动，但是高效的锥束CT图像配准以及体素的稠密对应仍具有挑战性。考虑到锥束CT中通常包含上亿个体素，使用一般的图像测度例如互信息与规范相关，并利用大规模的非线性优化通常难以实现在线快速的图像配准。此外当图像配准的初始化较差时，配准易于陷入局部最小。对于诸如计算机辅助的术中干预等在线的交互应用，高效可靠的锥束CT之间的体素稠密对应以及图像配准仍有待解决。

基于缩减采样的技术被用于加速图像配准与体素对应，尽管利用重要性采样加速了相似性测度的梯度计算，但其仍然依赖迭代优化技术并耗时数分钟。而且，在采样过程中，离散样本的采样方式多变并且不能覆盖所有的三维图像。有监督的分类与回归森林具有高效的训练与在线测试能力，并被用于三维医学图像中的体素对应，但是基于森林所获取的对应需要在后处理中利用正则化优化获取一致对应，这也是一个耗时的过程。此外，对于有监督的学习模型，例如分类与回归森林需要大量的训练数据。对三维图像进行标注费时耗力并易于受到个体标注的主观影响。基于超体素分解定义的伪标注可以缓解标注代价，但是，伪标注仅仅来自一个三维图像，难于在有限的标注数据实现具有较好泛化能力强的分类器。

发明内容

为了克服上述现有技术的不足，本发明提供一种锥束计算机断层扫描图像对应与配准方法，基于提出的混合测度森林算法，估计锥束CT图像之间的超体素之间的稠密对应。利用本发明提供的方法，可以快速建立三维锥束计算机断层扫描图像对应与配准，进行计算机辅助的术中干预与在线的特征点定位与分割标签迁移。

本发明的核心是：提出混合测度森林算法，生成锥束CT图像之间超体素之间的稠密对应。其中以迭代森林增强技术获取最优的基于森林的测度，并利用该测度估计锥束CT图像超体素之间的相似度。为了区分并评价非监督聚类森林对锥束CT图像中的超体素划分与相似度估计，利用从聚类随机森林获取的弱标注定义体素以及类别的置信分数，并基于置信分数区分未得到良好聚类的超体素以及相关的超体素类别。在迭代森林增强过程中，基于置信分数定义超体素与对应类别的惩罚权重，定义在节点分裂中的基于弱标注的关于分类器的加权信息增益。同时优先考虑在先前迭代中未获得良好聚类的样本，定义在节点分裂中关于聚类的加权信息增益。结合节点分裂过程中的分类与聚类信息增益，进行混合的信息增益估计，并基于该混合信息增益估计节点分裂的最优参数。在每步迭代中更新的混合测度森林会进一步消除先前迭代中未获得良好聚类的超体素，并更新锥束CT图像之间超体素的相似度矩阵以及超体素之间的稠密对应。为了进一步估计由混合测度森林返回的超体素候选匹配的空间一致性，基于超体素的置信分数定义候选匹配的柔性一致标签。利用回归森林进行一致性评价。在线测试阶段，利用混合测度森林估计由测试图像中分解的超体素之间的稠密对应，进而利用候选匹配的柔性一致标签对稠密对应进行正则化处理。混合测度森林基于非监督学习方法，不需要任何先验的标注信息，通过自学习数据分布以及图像之间一致匹配，在迭代过程中不断增强混合测度森林，并利用基于森林的测度估计锥束CT图像之间超体素的稠密对应与锥束CT图像配准。

本发明提供的技术方案是：

一种锥束计算机断层扫描图像对应与配准方法，基于混合测度森林算法，利用非监督聚类森林生成锥束CT图像之间的超体素之间的稠密对应；利用基于森林的测度，估计锥束CT图像之间的超体素之间的稠密对应与锥束CT图像的配准；包括如下步骤：

1)超体素特征提取：从锥束CT图像中分解超体素，提取超体素特征；

在训练锥束CT图像超体素的混合测度森林中，使用三种类型的超体素特征：灰度特征、空间上下文特征、测地坐标特征。

灰度特征表示为其中，b为超体素中包含体素的灰度直方图；为超体素的近邻超体素的灰度直方图的均值。

基于超体素与包含该超体素的上下文立方体中的超体素之间的灰度特征差异定义上下文特征。定义上下文特征f_c为f_c＝{χ₂(f_a(s)，f_a(s+Δ_k))|s+Δ_k∈P}，即当前超体素灰度特征f_a(s)与上下文超体素灰度特征f_a(s+Δ_k)之间的卡方距离；P为随机采样模式，由一组预先采样的偏移向量Δ_k组成，其中每个向量结合当前超体素的坐标均可确定上下文立方体中的一个超体素；

定义超体素的测地坐标特征f_g为当前超体素s与背景超体素s_g之间的最短测地距离f_g(s)＝min>g|G)；其中，G为在锥束CT图像分解的超体素之间构造的图模型，该图模型的节点为超体素，相邻超体素之间以边连接，并基于相邻超体素的灰度差异定义边上的权重为exp(-ρ||f_a(s_i)-f_a(s_j)||)，其中归一化的系数ρ＝1/max||f_a(s_i)-f_a(s_j)||。

2)训练初始聚类森林，利用非监督的聚类森林估计得到超体素之间的相似度；

给定一组包含M个锥束CT图像的图像集V＝{V_i|i＝0，...，M}与对应锥束CT图像分解的超体素集S＝{s_i|i＝1，...，N}，其中包含N个超体素，利用非监督的聚类森林估计超体素之间的相似度；

21)在非监督聚类森林训练中，利用数据分布不确定性的超椭球体积度量节点分裂的信息增益，表示为式1：

其中，tr为节点分裂后的左子节点l、右子节点r中的超体素数据的协方差矩阵的迹，可有效地避免高维数据协方差矩阵的秩不足问题；系数

22)对于包含n_T棵树的森林，基于聚类森林的相似度定义为式2：

a(s_i，s_j)＝1/n_T∑a_k(s_i，s_j)>

其中，a(s_i，s_j)表示基于聚类森林的测度返回的超体素s_i与s_j的相似度；a_k(s_i，s_j)表示基于聚类森林中的第k棵树返回的超体素s_i与s_j的相似度；对于第k棵树，对应的相似度a_k(s_i，s_j)＝[ι(s_i)＝ι(s_j)]，其中[·]是指示函数；当超体素s_i与s_j被分配到相同的叶子节点ι，即ι(s_i)＝ι(s_j)，则认为超体素s_i与s_j相似；否则不相似；

3)进行估计得到超体素的弱标注；

31)利用基于聚类随机森林的相似度，获取锥束CT图像V_r与V_t之间超体素的稠密对应；

32)在图像V_r与V_t之定义映射函数φ(s_i)＝s_j，其中s_j＝arg>i，s_j)，s_i∈V_r，s_j∈V_t；

33)利用该映射函数，将包含n_s个超体素参照图像V_r的超体素索引集合Y_r＝{y_i|y_i∈{1，...，n_s}}迁移到锥束CT图像V_t，其中，当φ(s_i)＝s_j，y(s_i)＝y(s_j)，即在锥束CT图像V_t中寻找与参照锥束CT图像中超体素s_i对应的超体素s_j，并将参照锥束CT图像中超体素s_i对应的超体索引赋予s_j作为图像V_t的弱标注；

利用从聚类随机森林获取的弱标注定义超体素以及类别的置信分数，并基于置信分数区分未得到良好聚类的超体素以及相关的超体素类别；

根据超体素的弱标注对应的索引迁移，在给定的锥束CT图像集中选择具有最一致的索引迁移的图像作为参照图像，使得锥束CT数据集利用该参照图像获取的对应超体素之间的相似度和最大；即其中，表示聚类森林返回的超体素与的相似度；为避免由于随机选择参照图像造成有偏的索引迁移；

4)训练混合测度森林，所述混合测度森林是基于初始聚类森林，使用迭代的方式增强基于森林的测度；在训练过程中，优先考虑在上一次迭代中未能获得良好聚类的超体素；利用步骤3)得到的超体素索引的弱标注，在节点分裂中通过关于分类的加权信息增益与关于聚类的加权信息增益，寻找最优的节点分裂参数；

41)针对锥束CT图像V^(k)与V^r，定义超体素类别y_i关于第k棵树的置信分数τ(k)为式3：

其中，Q^(k)与Q^r为关于锥束CT图像V^(k)与V^r的n_s×n_s′的距离矩阵，矩阵中元素对应锥束CT图像V^(k)与参照锥束CT图像V^r中标签为Y与Y_s的超体素之间的归一化的欧式距离；Y_s是标签集合Y的一个子集，其元素的个数为n_s′；δ是Heaviside函数；η是预先定义的不一致常数；矩阵Q的第i行是具有标签y_i的超体素与标签为Y_s的超体素之间的空间距离；

42)当图像V^(k)中具有标签y_i的体素与标签为Y_s的体素之间的空间距离关系与参考图像V^r一致时，标签y_i被认为是关于图像V^(k)可靠的标签；

43)利用在所有图像上的置信度定义标签y_i的置信度γ_i为式4：

其中，δ是Heaviside函数；为图像V^(k)上关于类别y_i的置信度；

44)利用标签v的置信度，定义基于弱标注所确定的离散概率分布对应的加权信息增益为式5：

其中，为当前节点数据属于第i类的概率；γ_i为第i类标签的置信分数；a_k为基于集合元素的个数定义系数，

45)利用置信分数区分在先前迭代中聚类较差的样本，并在节点分裂中估计数据分布不确定性时给予惩罚权重；定义超体素的惩罚权重v(s)为式6：

其中，y(s)是超体素s所对应的超体素类别标签，为对于锥束CT图像V^(k)的类别标签y(s)的置信分数；K是惩罚阈值常数；聚类相关的信息增益I_u定义为式7：

其中，是超体素数据集的加权平均；v是超体素的惩罚权重；基于集合元素的个数定义系数

46)估计由弱标注确定的离散概率分布的惩罚权重分类信息增益I_c以及数据分布不确定性的聚类信息增益I_u；混合测度森林在节点分裂时的信息增益I定义为式8：

其中，利用根节点分裂时分类信息增益与聚类的信息增益进行归一化；当更新混合测度森林后，弱标注与惩罚权重也相应地进行更新；使得在迭代过程中混合测度森林节点分裂的信息增益将优先考虑未能得到良好聚类的样本与类别，不断增强基于森林的测度用于估计锥束CT图像之间的对应；

通过上述步骤，得到训练好的混合测度森林；进一步估计柔性一致度。

5)估计柔性一致度，再利用柔性一致度估计在图像上下文约束下的由混合测度森林所确定超体素对应的一致性；

51)基于训练好的混合测度森林，通过式9获取锥束CT图像对超体素之间的对应集C：

其中，φ为图像之间超体素的映射函数；该对应基于混合测度森林测度；其中每一组对应均独立，缺乏每组对应与集合C中其它对应之间的关系；

52)利用回归森林进行柔性一致度估计，再利用柔性一致度估计在图像上下文约束下的由混合测度森林所确定超体素对应的一致性；

将集合C中匹配对标记为z＝(s_i，s_j)，对应的特征为特征f_z的第一项是超体素s_i与s_j特征之间的差异；第二项是对应点对(s_i，s_j)在超体素特征空间中的位置；

基于超体素的置信度引入柔性一致度其中y是超体素s_i与s_j的类别标签；与是超体素类别y关于锥束CT图像与的置信分数；在对应特征f_z与柔性一致度u(z)之间基于回归森林建立回归。

通过上述步骤，得到训练好的混合测度森林。

在线测试阶段，从锥束CT图像中分解超体素，提取超体素特征，并将超体素特征输入到训练好的混合测度森林，估计超体素之间的相似度与超体素之间的对应，并估计对应的超体素的柔性一致度。较高的柔性一致度分数u意味点对中的两个超体素具有图像与图像中的可靠标签，进而z＝(s_i，s_j)与集合C中其余对应一致。本发明具体实施中，柔性一致度u取值为0.1，保留大于0.1的对应。柔性一致度分数较低的超体素对应将被去除，并以该超体素近邻超体素在另一个锥束CT图像中所对应超体素确定新的对应。在获取锥束CT图像之间超体素对应后，即可方便地获取锥束CT图像之间的变形参数例如图像位移场并实现配准。

与现有技术相比，本发明的有益效果是：

本发明提供一种锥束计算机断层扫描图像对应与配准方法，基于混合测度森林，利用非监督聚类森林生成锥束CT图像之间的超体素之间的稠密对应；再利用基于森林的测度估计锥束CT图像之间的超体素之间的稠密对应与锥束CT图像的配准。利用本发明提供的方法，可以快速建立三维锥束计算机断层扫描图像对应与配准，进行计算机辅助的术中干预与在线的特征点定位与分割标签迁移。

附图说明

图1为本发明方法的流程框图。

具体实施方式

下面结合附图，通过实施例进一步描述本发明，但不以任何方式限制本发明的范围。

本发明提供一种锥束计算机断层扫描图像对应与配准方法，基于混合测度森林算法，利用非监督聚类森林生成锥束CT图像之间的超体素之间的稠密对应；利用基于森林的测度，估计锥束CT图像之间的超体素之间的稠密对应与锥束CT图像的配准。本发明基于迭代优化的混合测度森林估计锥束CT图像之间超体素的稠密对应与图像配准，并依据该对应获取自动的图像属性迁移。图1为本发明方法的流程框图，以下实施过程具体针对口腔正畸临床采集的锥束CT图像，利用混合测度森林估计超体素对应以及锥束CT图像配准，具体实施步骤包括：

步骤1：特征提取；

由锥束CT图像分解后得到超体素，再提取超体素特征。在训练锥束CT图像超体素的混合测度森林中，使用三种类型的超体素特征：灰度特征、空间上下文特征、测地坐标特征。基于超体素中包含体素的灰度直方图b定义灰度特征，同时为了抑制局部噪声，引入超体素的一圈近邻超体素的灰度直方图的均值灰度特征表示为基于超体素与包含该超体素的上下文立方体中的超体素之间的灰度特征差异定义上下文特征。在超体素上下文立方体中定义随机采样模式P，该模式由一组预先采样的偏移向量Δ_k组成，其中每个向量结合当前超体素的坐标都可以确定上下文立方体中的一个超体素。上下文特征f_c定义为当前超体素灰度特征f_a(s)与上下文超体素灰度特征f_a(s+Δ_k)之间的卡方距离，f_c＝{χ²(f_a(s)，f_a(s+Δ_k))|s+Δ_k∈P)。超体素的测地坐标f_g定义为当前超体素s与背景超体素s_g之间的最短测地距离。在锥束CT图像分解的超体素之间构造图模型，该图模型的节点为超体素，相邻超体素之间以边连接，并基于相邻超体素的灰度差异定义边上的权重exp(-ρ||f_a(s_i)-f_a(s_j)||)，其中归一化的系数ρ＝1/max||f_a(s_i)-f_a(s_j)||。测地坐标为f_g(s)＝mind(s，s_g|G)。

步骤二：训练初始聚类森林

给定一组包含M个锥束CT图像的图像集V＝{V_i|i＝0，...，M}与对应锥束CT图像分解的超体素集S＝{s_i|i＝1，...，N}，其中包含N个超体素。利用非监督的聚类森林估计超体素之间的相似度。在非监督聚类森林训练中，利用数据分布不确定性的超椭球体积度量节点分裂的信息增益，表示为式1：

I_u由节点分裂后的左子节点、右子节点(分别以下标l与r表示)中的超体素数据的协方差矩阵的迹(tr)定义。基于协方差矩阵的迹的信息增益计算可以有效地避免高维数据协方差矩阵的秩不足问题。基于集合元素的个数定义系数如果超体素s_i与s_j被分配到相同的叶子节点ι，即ι(s_i)＝ι(s_j)，则认为超体素s_i与s_j相似，否则不相似。对于第k棵树，对应的相似度a_k(s_i，s_j)＝[ι(s_i)＝ι(s_i)，其中[]是指示函数。对于包含n_T棵树的森林，基于聚类森林的测度定义为式2：

a(s_i，s_j)＝1/n_T∑a_k(s_i，s_j)>

其中，a(s_i，s_j)表示基于聚类森林的测度返回的超体素s_i与s_j的相似度；a_k(s_i，s_j)表示基于聚类森林中的第k棵树返回的超体素s_i与s_j的相似度。

步骤三：估计弱标注

利用基于聚类随机森林的测度，可获取锥束CT图像V_r与V_t之间超体素的稠密对应，进而在图像V_r与V_t之定义映射函数φ(s_i)＝s_j，其中s_j＝arg>i，s_j)，s_i∈V_r，s_j∈V_t。利用该映射函数，可以将包含n_s个超体素参照图像V_r的超体素索引集合Y_r＝{y_i|y_i∈{1，...，n_s}}迁移到锥束CT图像V_t，其中当φ(s_i)＝s_j，y(s_i)＝y(s_j)，即在锥束CT图像V_t中寻找与参照锥束CT图像中超体素s_i对应的超体素s_j，并将参照锥束CT图像中超体素s_i对应的超体索引赋予s_j作为图像V_t的弱标注。为了避免由于随机选择参照图像造成有偏的索引迁移，在给定的锥束CT图像集中选择具有最一致的索引迁移的图像作为参照图像，使得锥束CT数据集利用该参照图像获取的对应超体素之间的相似度和最大，即其中表示基于聚类森林的测度返回的超体素与的相似度。

步骤四：训练混合测度森林

混合测度森林基于初始聚类森林，使用迭代的方式增强基于森林的测度。在训练过程中，优先考虑在上一次迭代中未能获得良好聚类的超体素。利用超体素索引的弱标注，在节点分裂中结合关于分类的加权信息增益与关于聚类的加权信息增益，寻找最优的节点分裂参数。为了区分上一步迭代中未能得到良好聚类的超体素，定义超体素类别y_i关于第k棵树的置信分数τ(k)为式3：

其中，Q^(k)与Q^r为关于锥束CT图像V^(k)与V^r的n_s×n_s′的距离矩阵，矩阵中元素对应锥束CT图像V^(k)与参照锥束CT图像V^r中标签为Y与Y_s的超体素之间的归一化的欧式距离。Y_s是标签集合Y的一个子集，其元素的个数为n_s′。δ是Heaviside函数。η是预先定义的不一致常数。矩阵Q的第i行是具有标签y_i的超体素与标签为Y_s的超体素之间的空间距离。当图像V^(k)中具有标签y_i的体素与标签为Y_s的体素之间的空间距离关系与参考图像V^r一致时，标签y_i被认为是关于图像V^(k)可靠的标签。利用在所有图像上的置信度定义标签y_i的置信度γ_i：

其中，δ是Heaviside函数。为图像V^(k)上关于类别y_i的置信度。

利用标签y的置信度，基于弱标注所确定的离散概率分布对应的加权信息增益定义为式5：

其中，为当前节点数据属于第i类的概率。γ_i为第i类标签的置信分数。a_k为基于集合元素的个数定义系数，

利用置信分数区分在先前迭代中聚类较差的样本，并在节点分裂中估计数据分布不确定性时给予惩罚权重。超体素的惩罚权重v(s)定义为式6：

其中，y(s)是超体素s所对应的超体素类别标签，为对于锥束CT图像V^(k)的类别标签y(s)的置信分数。K是惩罚阈值常数，聚类相关的信息增益I_u定义为式7：

其中，是超体素数据集的加权平均。v是超体素的惩罚权重。

在训练混合测度森林中，估计由弱标注确定的离散概率分布的惩罚权重分类信息增益I_c以及数据分布不确定性的聚类信息增益I_u。混合测度森林在节点分裂时的信息增益I定义为式8：

其中，利用根节点分裂时分类信息增益与聚类的信息增益进行归一化。当更新混合测度森林后，弱标注与惩罚权重也相应地进行更新。在迭代过程中混合测度森林节点分裂的信息增益将优先考虑未能得到良好聚类的样本与类别，不断增强基于森林的测度用于估计锥束CT图像之间的对应。

步骤五：估计柔性一致度

通过式9，基于混合测度森林可以获取锥束CT图像对超体素之间的对应集C：

其中，φ为图像之间超体素的映射函数。该对应(映射函数φ)基于混合测度森林测度。但是其中的每一组对应都是独立的，即缺乏每组对应与集合C中其它对应之间的关系。系统引入柔性一致度估计在图像上下文约束下的由混合测度森林所确定超体素对应的一致性，并利用回归森林进行柔性一致度估计。将集合C中匹配对标记为z＝(s_i，s_j)。其对应的特征为特征f_z的第一项是超体素s_i与s_j特征之间的差异。第二项是对应点对(s_i，s_j)在超体素特征空间中的位置。基于超体素的置信度引入柔性标注其中y是超体素s_i与s_j的类别标签。与是超体素类别y关于锥束CT图像与的置信分数。在对应特征f_z与柔性一致度u(z)之间基于回归森林建立回归。较高的柔性一致性分数u意味点对中的两个超体素具有图像与图像中的可靠标签，进而对应z＝(s_i，s_j)与集合C中其余对应一致。在线测试阶段，柔性一致性分数较低的超体素对应将被去除，并以该超体素近邻超体素在另一个锥束CT图像中所对应超体素确定新的对应。在获取锥束CT图像之间超体素对应后，即可方便获取锥束CT图像之间的变形参数例如图像位移场并实现配准。

为了验证基于混合测度森林获取锥束CT图像超体素稠密对应，实验测试基于锥束CT图像之间超体素对应的结构分割标签迁移，分别测试了上颌骨、下颌骨分割标签的迁移，其对应DICE相似度系数都在90％以上，具有较好的一致对应。

利用本发明的方法，以完全无监督的方式获取了锥束CT图像之间的超体素稠密对应与图像配准。迭代优化的混合测度森林可以处理先前迭代中未获取良好聚类的超体素，并依靠柔性一致估计获取锥束CT图像之间超体素对应以及图像配准。

需要注意的是，以上实施方式仅用于说明本发明，而并非对本发明的限制，本领域的技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。