CN112734814B - 三维颅面锥形束ct图像配准方法 - Google Patents

Abstract

本发明公布了一种三维颅面锥形束CT图像配准方法，建立基于非监督学习的配准全卷积神经网络，通过成对CBCT图像之间稠密一致的体素级别对应求解，有效实现三维颅面CBCT的图像配准，并满足可逆一致的要求，实现三维颅面锥形束CT(CBCT)图像的高效的在线配准，大幅减少配准用时。

Description

三维颅面锥形束CT图像配准方法

技术领域

本发明涉及口腔临床医学和计算机视觉技术领域，具体涉及一种三维颅面锥形束计算机断层扫描(Cone Beam Computed Tomography，简称CBCT)图像配准方法。

背景技术

三维图像稠密体素对应的求解是医学图像分析中的关键性技术，被广泛应用于统计形状分析，特征点与分割迁移，以及术前术后疗效评价等任务中。

现有传统的三维图像非刚性配准方法通过在线的非线性迭代优化减小三维图像稠密体素级别的表观差异，常用的几种测度包括均方差、归一化交叉相关、互信息等。非刚性配准中的形变模型包括弹性体模型、demons、微分同胚映射、B样条自由形变模型等。考虑到三维图像配准中大规模的待求解参数，传统配准方法的优化求解通常有很大的计算时间代价并且易于陷入局部极小。一些方法使用降采样和统计形变模型来减少问题求解的参数，但仍依赖于迭代的参数优化求解。基于数值求解的非刚性配准方法则需要手工的关键点和分割图标注。保持拓扑、变形场没有折叠、可逆平滑的一一映射是解剖图像配准应该具有的重要性质，为此包括LDDMM、SyN、微分同胚demons等利用微分同胚变形模型的配准方法被用于解剖结构配准中。

近年来，回归模型被引入到三维图像非刚性配准中，基于学习的三维图像非刚性配准方法通过寻找输入图像和配准变形场的回归模型来避免在线的迭代优化，比如支持向量机、随机森林和深度神经网络被引入到非刚性配准中，相比于传统的优化方法，在配准耗时上有数量级的减少。有监督的基于深度学习的三维图像配准方法在训练过程中需要关键点和变形场的真值，而三维图像的手工标注难度会限制训练数据的获取。非监督的基于深度学习的三维图像配准方法使用空间变换网络来减小体素级别的表观差异。已有的基于深度学习的三维图像配准方法通过三维图像分块配准减少所需的内存占用，或者只求解一张图像到参考图像的单向映射。为了求解可逆的一一映射，已有的使用微分同胚映射的基于深度学习的三维图像配准方法需要重参数化层将神经网络预测出的速度场转化成双向的变形场。一些方法利用了可逆约束直接推理双向的变形场，避免了中间层的估计，但依赖于手工提取特征。因此，利用神经网络直接预测可逆一致的变形场目前尚未出现有效的解决方案。

发明内容

本发明提供一种三维颅面锥形束CT(CBCT)图像配准方法，目的是克服上述缺陷，对三维图像进行高效的在线配准，并满足可逆一致的要求。

为了达到上述目的，本发明建立一种基于非监督学习的配准框架利用全卷积神经网络，实现了成对CBCT图像之间稠密一致的体素级别对应求解，相比传统迭代优化配准方法大幅减少配准用时。本发明使用卷积神经网络回归输入图像V和模板图像T之间的体素级别的双向映射，通过引入组合可逆约束保证正向和反向变形场满足可逆一致的要求，并在网络训练过程中使用数据集中的任意二元组和三元组通过对变形场的约束加强变形场的闭环一致。在配准网络训练中，***允许用户指定感兴趣结构(SOI)，增强配准难度较高的解剖结构的配准效果，比如骨密度较低、与周围软组织边界模糊的颅底。对于任意一对给定图像的在线配准，训练好的卷积神经网络分别返回每张图像与模板图像T的正向和方向变形场，通过变形场的组合得到一对图像可逆一致的正向和反向变形场。

本发明中，计算机断层扫描：CT(Computed Tomography)；锥形束计算机断层扫描：Cone Beam Computed Tomography(CBCT)。

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

一种三维颅面锥形束CT(CBCT)图像配准方法，建立基于非监督学习的配准全卷积神经网络，通过成对CBCT图像之间稠密一致的体素级别对应求解，有效地实现三维颅面锥形束CT(CBCT)的图像配准，可大幅减少配准用时；包括如下步骤：

1)定义模板图像，经过平均变形场变形后得到模板图像T；

将一组三维图像数据集v＝{V_i|i＝1，...，N}配准到一个随机选取的参考图像V_r，得到变形场φ_i，N是数据集的图像数量。模板图像T是经过平均变形场变形后得到的图像。

公式1中，表示利用变形场/>对图像V_r进行变形。

2)配准网络的搭建与预训练，包括如下过程：

21)使用基于固定模板图像的双向配准网络求解三维图像V和模板图像T之间体素级别的平滑一一映射，V和T是三维单通道的灰度图像；

22)使用卷积神经网络回归输入的图像和输出的双向变形场：

h_Θ，T：V→[φ_f，φ_b]

其中，φ_f为输入图像V变形到模板图像T的变形场，φ_b为模板图像T变形到输入图像V的变形场，Θ是配准网络h的可学习参数(例如卷积神经网络的卷积核权重)。设定图像对 以及/>具有相似的解剖表观(可通过公式4来衡量相似程度)。

23)卷积神经网络使用对称的带有长残差连接的编解码器结构，根据输入的单通道三维图像V推测六通道的正向和反向变形场[φ_f，φ_b]。编解码器由编码器和解码器两部分构成，解码阶段的特征图像由前一层解码层输出的特征图像和编码阶段的同尺寸特征图像组成，以促进特征的传播和训练收敛。

本发明使用对称的带有长残差连接的编解码器结构根据输入的一张单通道三维图像V推测六通道的正向和反向变形场[φ_f，φ_b]。编码器包含六个卷积核大小3×3×3步长为1的卷积层，每个卷积层后为实例正则化(instance normalization)、Leaky ReLU激活函数和2×2×2最大池化层。解码器包含六个卷积核大小3×3×3步长为1的卷积层。每个卷积层之前为×2上采样层，之后为实例正则化和ReLU激活函数层。

24)使用合成数据对卷积神经网络进行预训练，使用随机采样的B样条变形参数生成变形场，并确保生成的变形场没有折叠；利用生成的变形场对模板图像T进行变形，生成合成图像V_syn，并计算出反向变形场。使用合成图像和对应的正向和反向变形场对神经网络进行预训练。

3)采用无监督方式训练步骤2)构建的配准网络，得到训练好的配准网络；包括如下过程：

对于具有分布的训练数据集v＝{V_i|i＝1，...，N}，N是训练数据集的图像数量，通过寻找最优的网络可优化参数Θ使损失函数最小。

其中，h为配准网络；Θ是配准网络h的可学习参数，V是训练图像，T是模板图像，是训练数据集的分布。

配准网络训练的损失函数为：

其中，L_sim是相似性度量损失，L_reg是变形场正则化项，L_inv是变形场组合可逆约束项，L_str是感兴趣结构(SO1)配准损失，Θ是卷积神经网络可优化参数，α_sim、α_reg、α_inv、α_str分别是各损失项的权重。

下面依次对每项损失函数进行说明：

31)配准网络预测出正向和反向变形场[φ_f，φ_b]，相似性测度L_sim度量利用正向变形场φ_f，i变形后的输入图像与模板图像T在体素级别的差异，以及利用反向变形场φ_b，i变形后的模板图像/>与输入图像V_i的差异，m是一个批数据的图像数量。相似性度量损失L_sim表示为公式4：

32)变形场正则化项L_reg用来保证正向和反向变形场的平滑性，使用变形场的空间梯度作为正则化项，表示为公式5。

其中，φ_f，i，φ_b，i为正向和反向变形场，||·||_F是Frobenius范数，m是一个批数据的图像数量。

33)变形场组合可逆约束项L_inv：

本发明使用组合可逆约束使正向和反向变形场满足正向逆向一致与闭环一致，变形场组合可逆约束项L_inv表示为公式6：

L_inv的第一项要求双向变换[φ_f，φ_b]是互逆的，m是训练中一批图像的数量。||·||_F是Frobenius范数。L_inv的第二项，要求k个图像{V_i|i＝0，...，k-1}的变形场沿着闭合环c的级联是一个恒等变换。对于有N张三维图像的训练集，长度大于等于2的闭合环的个数为/>为了简化训练，可以只考虑两个或三个体数据形成的闭环，K设置为3。L_inv的第二项要求任意两张图像的正向和反向变形场互逆，任意三张图像的变形场满足闭环一致。

公式6中，φ_f，i，φ_b，i为正向和反向变形场，||·||_F是Frobenius范数，m是一个批数据的图像数量，为k个变形场沿着闭合环c的组合。对于变形场的组合，利用标准坐标场Q经过多个变形场变形后变为Q′，Q′-Q即为变形场的组合结果。

34)感兴趣结构(SOI)配准损失L_str：

在三维颅面图像的配准中，一些解剖结构比如颅底由于较低的骨密度以及设备相关的噪声，其与周围的软组织的边界相对模糊，配准难度较高，本发明在网络训练中使用感兴趣结构(SOI)配准损失函数L_str来改善精细结构的配准，表示为公式7。

其中，M为感兴趣区域(SOI)的掩模，其中体素值为1表示感兴趣区域，其余体素值为0。m为一个批数据的图像数量，V_i为网络的输入图像，T为模板图像，φ_f，i为网络输出的正向变形场，φ_b，i为网络输出的反向变形场。运算符★表示逐个元素的矩阵乘法。

4)对待配准的三维图像(成对CBCT图像)进行在线配准：

在线配准阶段，对于任意一对浮动图像V_M和固定图像V_F，记作(V_M，V_F)，通过步骤3)训练好的配准网络，分别返回V_M和V_F的正向和反向变形场并通过变形场组合得到最终的配准变形场φ_M→F和φ_F→M：

其中，利用φ_M→F表示可以将浮动图像变形到固定图像，利用φ_F→M可以将固定图像变形到浮动图像，应与V_M具有相似的表观特征，/>应与V_F具有相似的表观特征。得到变形场即配准好图像。浮动图像利用配准生成的变形场做变形后，得到的图像与固定图像表观越相似，则配准效果越好。

通过上述步骤，实现三维颅面锥形束CT(CBCT)图像的配准。

本发明的有益效果：

利用本发明提供的三维颅面锥形束CT(CBCT)图像配准方法，可以求解一对三维图像体素级别的稠密对应，并满足可逆一致的要求，实现三维颅面锥形束CT(CBCT)图像的高效的在线配准。

附图说明

图1是本发明方法对应的配准网络训练框架图；

图2是本发明的在线配准流程框图。

图3是本发明实施例中一对图像的配准变形场和配准结果示意；

其中，(a)为输入图像对(V_M，V_F)；(b)和(c)分别是V_M和V_F的正向和反向变形场；(d)是正向变形场φ_M→F和反向变形场φ_F→M(e)变形后的图像；

图4是本发明实施例中一对三维图像配准前以及经过正向配准和经过反向配准的的结果在矢状面、冠状面、轴状面的重叠示意；

其中，(a)为矢状面；(b)为冠状面；(c)轴状面重叠。(d)变形场φ_f和φ_b的切片。

具体实施方式

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

本发明提供的三维颅面锥形束CT(CBCT)图像配准方法，使用卷积神经网络回归输入图像和双向变形场，并利用空间变换网络减小目标图像和变形后图像的体素级别的表观差异，通过组合可逆约束保证正向和反向变形场的可逆一致性，在线配准时避免了现有传统配准方法的迭代优化，可大幅减少配准用时。

下面，参照附图和实施例对本发明具体实施方式做进一步说明。

步骤1：定义模板图像

将一组采集到的口腔颅面三维锥形束CT图像数据集v＝{V_i|i＝1，...，N}被配准到一个随机选取的参考图像V_r，得到变形场φ_i，N是数据集的图像数量。模板图像T是经过平均变形场变形后得到的图像。

步骤2：配准网络搭建与预训练

本发明使用基于固定模板图像的双向配准网络求解三维图像V和模板图像T之间体素级别的平滑一一映射，V和T是三维单通道的灰度图像，使用卷积神经网络回归输入的图像和输出的双向变形场h_Θ，T：V→[φ_f，φ_b]，φ_f：V→T为输入图像V变形到模板图像T的变形场，φ_b：T→V为模板图像T变形到输入图像V的变形场，Θ是配准网络h的可学习参数，例如卷积神经网络的卷积核权重。图像对以及/>被期望具有相似的解剖表观。在训练中通过度量前向和反向两个方面的配准相似性得到最优的参数Θ。

本发明使用对称的带有长残差连接的编解码器结构根据输入的一张单通道三维图像V推测六通道的正向和反向变形场[φ_f，φ_b]。编码器包含六个卷积核大小3×3×3步长为1的卷积层，每个卷积层后为实例正则化(instance normalization)、Leaky ReLU激活函数和2×2×2最大池化层。解码器包含六个卷积核大小3×3×3步长为1的卷积层。每个卷积层之前为×2上采样层，之后为实例正则化和ReLU激活函数层。解码阶段的特征图像由前一层解码层输出的特征图像和编码阶段的同尺寸特征图像组成，以促进特征的传播和训练收敛。

本发明使用合成数据对卷积神经网络进行预训练，使用随机采样的B样条变形参数生成变形场，并确保生成的变形场没有折叠，利用生成的变形场对模板图像T进行变形生成合成图像V，并计算出反向变形场。使用合成图像和对应的正向和反向变形场对神经网络进行预训练。

步骤3：无监督方式训练配准网络

本发明使用组合的可逆约束来保证变形场的可逆一致和闭环一致性。对于具有分布D_v的训练数据集v＝{V_i|i＝1，...，N}，N是训练数据集的图像数量，通过寻找最优的网络可优化参数Θ使损失函数最小。

其中，Θ是配准网络h的可学习参数，V是训练图像，T是模板图像，是训练数据集的分布。

配准网络训练的损失函数为：

其中，L_sim是相似性度量损失，L_reg是变形场正则化项，L_inv是变形场组合可逆约束项，L_str是感兴趣结构(SOI)配准损失，Θ是卷积神经网络可优化参数，α_sim、α_reg、α_inv、α_str分别是各损失项的权重。

下面依次对每项损失函数进行说明：

1)配准网络预测出正向和反向变形场[φ_f，φ_b]，相似性测度L_sim度量利用正向变形场φ_f，i变形后的输入图像与模板图像T在体素级别的差异，以及利用反向变形场φ_b，i变形后的模板图像/>与输入图像V_i的差异，m是一个批数据的图像数量。

2)正则化项L_reg用来保证正向和反向变形场的平滑性，使用变形场的空间梯度作为正则化项。

其中，φ_f，i，φ_b，i为正向和反向变形场，||·||_F是Frobenius范数，m是一个批数据的图像数量。

3)本发明使用组合可逆约束使正向和反向变形场满足正向逆向一致与闭环一致，L_inv的第一项要求双向变换[φ_f，φ_b]是互逆的，m是训练中一批图像的数量。||·||_F是Frobenius范数。L_inv的第二项，要求k个图像{V_i|i＝0，...，k-1}的变形场沿着闭合环c的级联是一个恒等变换。对于有N张三维图像的训练集，长度大于等于2的闭合环的个数为/>为了简化训练，可以只考虑两个或三个体数据形成的闭环，K设置为3。L_inv的第二项要求任意两张图像的正向和反向变形场互逆，任意三张图像的变形场满足闭环一致。

其中，φ_f，i，φ_b，i为正向和反向变形场，||·||_F是Frobenius范数，m是一个批数据的图像数量，为k个变形场沿着闭合环c的组合。对于变形场的组合，利用标准坐标场Q经过多个变形场变形后变为Q′，Q′-Q即为变形场的组合结果。

4)在三维颅面图像的配准中，一些解剖结构比如颅底由于较低的骨密度以及设备相关的噪声，其与周围的软组织的边界相对模糊，配准难度较高，本发明在网络训练中使用感兴趣结构(SOI)配准损失函数L_str来改善精细结构的配准。

其中，M为感兴趣区域(SOI)的掩模，其中体素值为1表示感兴趣区域，其余体素值为0。m为一个批数据的图像数量，V_i为网络的输入图像，T为模板图像，φ_f，i为网络输出的正向变形场，φ_b，i为网络输出的反向变形场。运算符★表示逐个元素的矩阵乘法。

步骤4：成对CBCT图像在线配准

在线配准阶段，给定任意一对浮动图像和固定图像(V_M，V_F)，配准网络分别返回V_M和V_F的正向和反向变形场并通过变形场组合得到最终的变形场

为了验证配准网络对成对CBCT图像的配准精度，对于合成数据，计算网络预测的正向和反向变形场与变形场真实值的均方距离，平均误差小于0.25mm，同时变形场具有较好的可逆性，雅可比行列式小于等于0的体素个数平均值为7.44，表明本发明的配准方法预测出的变形场具有较好的可逆性。

利用本发明的方法，可以快速有效地对成对的图像(如锥形束CT图像)进行非刚性配准。***利用卷积神经网络回归输入图像与双向变形场，避免了传统配准方法的迭代求解，大幅减少配准用时。***利用组合可逆约束，进一步要求对于训练数据集中任意二元组和三元组，变形场满足闭合环一致的要求，使配准网络预测出的正向和反向变形场具有较好的可逆性。***利用感兴趣区域(SOI)配准损失增强精细解剖结构的配准。该***有效克服了传统非刚性配准方法计算耗时、易陷入局部最小的问题，并满足变形场可逆一致的要求。

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

Claims (10)

1.一种三维颅面锥形束CT图像配准方法，所述锥形束CT简称为CBCT；所述配准方法建立基于非监督学习的配准全卷积神经网络，通过成对CBCT图像之间稠密一致的体素级别对应求解，有效实现三维颅面CBCT的图像配准；包括如下步骤：

1)定义模板图像，经过平均变形场变形后得到模板图像T；

将一组三维图像数据集v＝{V_i|i＝1，...，N}配准到随机选取的参考图像V_r，得到变形场φ_i；其中，N是数据集的图像数量；

模板图像T是经过平均变形场变形后得到的图像，用公式1表示：

公式1中，表示利用变形场/>对图像V_r进行变形；

2)搭建配准网络并进行预训练，包括如下过程：

21)使用基于固定模板图像的双向配准网络求解三维图像和模板图像之间体素级别的平滑一一映射；三维图像和模板图像均为三维单通道的灰度图像；

22)使用卷积神经网络回归输入的图像和输出的双向变形场：

h_Θ，T：V→[φ_f，φ_b]

其中，V为输入图像；φ_f为输入图像V变形到模板图像T的变形场；φ_b为模板图像T变形到输入图像V的变形场；Θ是配准网络h的可学习参数；

23)使用对称的带有长残差连接的编解码器结构，根据单通道的输入图像V推测得到六通道的正向和反向变形场[φ_f，φ_b]；解码阶段的特征图像由前一层解码层输出的特征图像和编码阶段的同尺寸特征图像组成；

24)使用合成图像和对应的正向和反向变形场数据对卷积神经网络进行预训练，使用随机采样的B样条变形参数生成变形场，并使得生成的变形场没有折叠；利用生成的变形场对模板图像T进行变形，生成合成图像V_syn，并计算出反向变形场；

3)采用无监督方式训练步骤2)的配准网络，得到训练好的配准网络；包括如下过程：

对于具有分布的训练数据集v＝{V_i|i＝1，...，N}，N是训练数据集的图像数量，通过寻找最优的网络可优化参数Θ使损失函数最小，表示为公式2：

其中，h为配准网络；Θ是配准网络h的可学习参数；V是训练图像；T是模板图像；是训练数据集的分布；

配准网络训练的损失函数表示为公式3：

其中，L_sim为相似性度量损失；L_reg是变形场正则化项；L_inv是变形场组合可逆约束项；L_str为感兴趣结构SOI配准损失；Θ是卷积神经网络可优化参数；α_sim、α_reg、α_inv、α_str分别是相应损失项的权重；

4)对待配准的三维CBCT图像进行在线配准：

对于任意一对浮动图像V_M和固定图像V_F，记作(V_M，V_F)，通过步骤3)训练好的配准网络，分别返回V_M和V_F的正向和反向变形场并通过变形场组合得到最终的配准变形场φ_M→F和φ_F→M，表示为：

其中，利用φ_M→F表示可以将浮动图像变形到固定图像，利用φ_F→M可以将固定图像变形到浮动图像，与V_M具有相似的表观特征，/>与V_F具有相似的表观特征；

得到变形场即实现图像配准；浮动图像利用配准生成的变形场做变形后得到的图像即为配准后的图像；

通过上述步骤，即实现三维颅面CBCT图像的配准。

2.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤22)中，配准网络h的可学习参数Θ为卷积神经网络的卷积核权重。

3.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤22)中，具体设定图像对与/>具有相似的解剖表观，在训练中通过度量前向和反向的配准相似性得到最优的参数Θ。

4.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤23)中，编码器包含六个卷积核大小为3×3×3、步长为1的卷积层，每个卷积层后为实例正则化、Leaky ReLU激活函数和2×2×2最大池化层。

5.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤23)中，解码器包含六个卷积核大小为3×3×3、步长为1的卷积层；每个卷积层之前为×2上采样层，之后为实例正则化和ReLU激活函数层。

6.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤3)配准网络训练的损失函数中，相似性度量损失L_sim具体表示为公式4：

其中，m是图像数量；φ_f，i为正向变形场；φ_b，i为反向变形场；表示利用正向变形场φ_f，i变形后的输入图像；/>表示利用反向变形场φ_b，i变形后的模板图像。

7.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤3)配准网络训练的损失函数中，变形场正则化项L_reg用于保证正向和反向变形场的平滑性，使用变形场的空间梯度作为正则化项，表示为公式5：

其中，φ_f，i、φ_b，i为正向和反向变形场，||·||_F是Frobenius范数，m是图像数量。

8.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤3)配准网络训练的损失函数中，变形场组合可逆约束项L_inv表示为公式6：

L_inv的第一项要求双向变换[φ_f，φ_b]是互逆的，m是训练中一批图像的数量；||·||_F是Frobenius范数；

L_inv的第二项中，为k个变形场沿着闭合环c的组合；k个图像{V_i|i＝0，...，k-1}的变形场沿着闭合环c的级联/>是一个恒等变换，对于有N张三维图像的训练集，长度大于等于2的闭合环的个数为/>

9.如权利要求8所述三维颅面锥形束CT图像配准方法，其特征是，具体只考虑两个或三个体数据形成的闭环，k取值为3。

10.如权利要求1所述三维颅面锥形束CT图像配准方法，其特征是，步骤3)配准网络训练的损失函数中，感兴趣结构SOI配准损失L_str用于改善精细结构的配准，表示为公式7：

其中，M为感兴趣区域SOI的掩模，其中体素值为1表示感兴趣区域，其余体素值为0；m为图像数量；V_i为网络的输入图像；T为模板图像；φ_f，i为网络输出的正向变形场；φ_b，i为网络输出的反向变形场；运算符★表示逐个元素的矩阵乘法。