CN110473234B - 微分同胚Demons图像配准方法、系统及存储介质 - Google Patents

Abstract

本发明涉及一种微分同胚Demons图像配准方法、系统及存储介质，方法包括以下步骤：S1、选择参考图像s和浮动图像f；S2、设置初始形变位移场t₀；S3、选择Demons力计算更新位移场u；S4、对更新位移场u进行高斯卷积；S5、采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c；S6、对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；S7、采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；S8、判断是否迭代终止，是则输出最新的形变位移场t_n；否则转向步骤S3。

Description

微分同胚Demons图像配准方法、系统及存储介质

技术领域

本发明是关于一种基于自适应步长的微分同胚Demons图像配准方法、系统、计算机程序及存储介质，涉及医学图像处理技术领域。

背景技术

Demons图像配准方法源于光流学模型，是基于图像灰度自动的非刚性图像配准方法，该算法实现简单并且计算效率高，因此吸引了大量研究者对其进行深入研究。针对原始Demons算法存在的一些问题，研究人员分别提出了主动Demons算法、对称Demons算法以及微分同胚Demons算法。其中，主动Demons算法和对称Demons算法主要适用于形变较大图像之间的配准，并且配准精度更高，收敛速度更快。而微分同胚Demons算法则可以生成可逆的变形场，保持图像的拓扑结构。

图像配准实质上是一个优化过程，通过优化找到一个映射关系，使基于浮动图像和参考图像生成的能量函数达到最小，最终使浮动图像经过坐标变换后的体素与参考图像中体素所指向的解剖结构位置在空间一致。研究人员从最小化能量函数(包括相似性测度和正则化项)的角度重新解释了Demons算法，得出微分同胚原始Demons算法、主动Demons算法以及对称Demons算法分别等价于使用Newton-Raphson(牛顿法)、Gauss-Newton以及ESM优化算法对能量函数进行优化。其中，现有的微分同胚Demons图像配准算法的流程如图1所示，实现步骤为：

(1)设置初始的形变位移场t；

(2)选择不同的Demons力(包括原始Demons力、主动Demons力以及对称Demons力，分别对应于原始Demons算法、主动Demons算法以及对称Demons算法)的表达式，计算更新位移场u；

(3)对u进行高斯卷积，从而产生流体映射模型的规则化效果；

(4)通过复合运算更新c←to(α*exp(u))，其中，o代表复合运算，α是步长；

(5)对c进行高斯卷积并更新t，从而产生扩散模型的规则化效果；

(6)判断是否迭代终止，是则输出t；否则转向步骤(2)。

微分同胚Demons算法优化的是每个体素在x、y、z三个方向的位移量，在图1所示微分同胚Demons算法优化过程的关键是确定搜索方向和步长α，其中，搜索方向通过选择的Demons力决定，即更新位移场u；而步长α则往往设置为固定值，这样设置步长有以下缺陷：(1)步长设置为固定值，如果设置过大，则在迭代更新的过程中容易错过极小值；设置过小，则每次更新少，迭代速度慢；(2)相同的步长被应用于所有体素，这样某些体素处位移量优化达到收敛条件时，而其它体素处位移量优化还未收敛，导致迭代收敛慢，配准精度降低；(3)难以设置最佳的初始步长。

发明内容

针对上述问题，本发明的目的是提供一种基于自适应步长的微分同胚Demons图像配准方法，该配准方法中通过对每个体素(三维图像空间的最小单位)设置不同的步长以及优化过程中自适应每个体素的步长，加快收敛速度并提高配准精度，适应于微分同胚原始Demons算法、对称Demons算法以及主动Demons算法。

为实现上述目的，本发明采取以下技术方案：

本发明实施例的第一方面，提供一种微分同胚Demons图像配准方法，包括以下步骤：

S1、选择参考图像s和浮动图像f；

S2、设置初始形变位移场t₀；

S3、选择Demons力计算更新位移场u；

S4、对更新位移场u进行高斯卷积；

S5、采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c；

S6、对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；

S7、采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

S8、判断是否迭代终止，是则输出更新的形变位移场t_n；否则转向步骤S3。

进一步地，初始形变位移场t₀采用与参考图像大小相同的单位矩阵。

进一步地，选择Demons力计算更新位移场u，具体为：

Demons力为原始Demons力、主动Demons力或对称Demons力，分别对应于原始Demons算法、主动Demons算法以及对称Demons算法，选择Demons力计算更新位移场u，其中，

原始Demons力的表达式为：

主动Demons力的表达式为：

对称Demons力的表达式为：

式中，f代表浮动图像，s代表参考图像，(f-s)是两幅图像对应体素灰度值的差异，

和

分别对应浮动图像和参考图像的梯度，u＝(u_x,u_y,u_z)是每个体素在x、y和z方向的位移量，

是归一化因子，u的上限是

进一步地，对更新位移场u进行高斯卷积采用u←K_fluid*u，其中，高斯卷积模板是对连续的三维高斯函数的离散化表示，任意大小的高斯卷积模板都可以通过建立一个(2h+1)×(2h+1)×(2h+1)的M得到，(i,j,k)位置的元素值通过下式确定：

式中，M为三维矩阵，矩阵大小是(2h+1)×(2h+1)×(2h+1)，(i,j,k)是M中体素的位置索引，h是高斯模板的中心，K_fluid即是高斯滤波器的σ值，卷积过程就是在矩阵u上逐点移动卷积模板，使模板中心和u上的点(x,y,z)重合，卷积模板在每一点(x,y,z)的响应是根据模板中的值与u中对应体素的值相乘再求和得到。

进一步地，采用优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c，具体过程为：

S51、基于u更新有偏差的一阶矩和二阶矩m和v，其中：

m_n←β₁·m_n-1+(1-β₁)·u

v_n←β₂·v_n-1+(1-β₂)·u²

式中，β₁是一阶矩估计的指数衰减率，β₂是二阶矩估计的指数衰减率，m和v是中间变换矩阵，其初始设置为零矩阵，在每次迭代中更新，n为迭代次数；

S52、更新无偏差的一阶矩和二阶矩

和

其中：

S53、通过复合运算更新初始位移场t₀得到空间变换c：

式中，ξ为常数。

进一步地，对空间变换c进行高斯卷积得到更新的形变位移场t_n，其中，t_n←K_diff*c，K_fluid即是高斯滤波器的σ值。

进一步地，迭代终止条件设定为3个，满足其中任意一个，迭代终止，3个终止条件分别为：

①当前迭代次数n超过预设的最大迭代次数时，则迭代终止；

②计算参考图像s与变形后的浮动图像f'对应体素灰度差的平方和的平均值，如果该值在三次连续的迭代中持续增大，则迭代终止；

③计算每次迭代后更新的形变位移场t_n与上一次迭代中更新的形变位移场t_n-1之间的差异的平均值，如果该值在连续的三次迭代中持续增大，则迭代终止；

若不符合迭代终止条件，则以当前迭代得到的形变位移场t_n作为初始的形变位移场t₀←t_n，当前迭代得到的变形后的浮动图像作为初始的浮动图像f←f'，迭代次数n←n+1，转向步骤S3继续进行优化，整个优化过程中，参考图像s保持不变。

本发明实施例的第二方面，提供一种微分同胚Demons图像配准系统，该系统包括：

图像选择单元，用于选择参考图像s和浮动图像f；

形变位移场设置单元，用于设置初始形变位移场t₀；

更新位移场计算单元，用于选择Demons力计算更新位移场u；

第一平滑单元，用于对更新位移场u进行高斯卷积；

步长调整单元，用于采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c；

第二平滑单元，用于对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；

空间变换单元，用于采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

判断单元，用于判断是否迭代终止，是则输出更新的形变位移场t_n；否则转向更新位移场计算单元。

本发明实施例的第三方面，提供一种计算机程序，包括计算机程序指令，其中，所述程序指令被处理器执行时用于实现如第一方面所述的微分同胚Demons图像配准方法对应的步骤。

本发明实施例的第四方面，提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序指令，其中，所述程序指令被处理器执行时用于实现如第一方面所述的微分同胚Demons图像配准方法对应的步骤。

本发明由于采取以上技术方案，其具有以下特点：本发明利用adam优化算法自适应调整Demons算法中的步长，在提高配准精度的同时加快收敛速度，减少运行时间，从而拓宽微分同胚Demons算法的应用范围。

附图说明

图1为现有的微分同胚Demons图像配准算法流程图；

图2为本发明实施例1中两幅图像之间的形变位移场；

图3为本发明实施例1中微分同胚Demons图像配准算法流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1：

如图2所示，本实施例针对4DCT不同呼吸时相的图像进行配准，以得到两幅图像之间变形的位移场，该位移场可以用于4D放射治疗计划中器官轮廓的自动推衍以及剂量叠加等。

如图3所示，本实施例提出的基于自适应步长的微分同胚Demons图像配准方法，包括以下步骤：

S1：选择某一呼吸时相的CT图像为参考图像s，另外任一呼吸时相的CT图像为浮动图像f；配准的目的是得到浮动图像上每个体素空间变换到参考图像上对应体素需要的位移量，所有体素处的位移量形成一个位移矢量场。

S2：设置初始的形变位移场t₀；

通常设置每个体素处的初始位移量为1个体素，即位移场t₀为单位矩阵，大小与参考图像大小一致。4DCT不同呼吸时相的图像是通过对患者同一体位扫描得到的数据进行不同时相的重组得到的，因此无须再进行刚体变换等，可以直接设置初始的形变位移场。

S3：选择Demons力(包括原始Demons力、主动Demons力以及对称Demons力，分别对应于原始Demons算法、主动Demons算法以及对称Demons算法)，通过选择的Demons力的表达式计算更新位移场u，其中，

原始Demons力的表达式为：

主动Demons力的表达式为：

对称Demons力的表达式为：

式中，f代表浮动图像，s代表参考图像，(f-s)是两幅图像对应体素灰度值的差异，

和

分别对应浮动图像和参考图像的梯度，u＝(u_x,u_y,u_z)是每个体素在x、y和z方向的位移量，n是迭代次数，

是归一化因子，u的上限是

一般取值0.5～1.0。配准初期变形比较大时，可以使用较小的

这样每次迭代允许的位移比较大，可以加快收敛速度。当算法接近收敛时可以使用较大的

这样，每次迭代允许的位移较小，可以保证配准精度。

S4：对步骤S3中的更新位移场u进行高斯卷积u←K_fluid*u，卷积的目的是为了产生流体映射模型的规则化效果，规则化主要是通过平滑对u进行约束，避免产生不切合实际的变形。

高斯卷积模板是对连续的三维高斯函数的离散化表示。任意大小的高斯卷积模板都可以通过建立一个(2h+1)×(2h+1)×(2h+1)的M得到，其(i,j,k)位置的元素值通过下式确定：

其中，M是一个三维矩阵，矩阵大小是(2h+1)×(2h+1)×(2h+1)，(i,j,k)是M中体素的位置索引，h是高斯模板的中心，K_fluid即是高斯滤波器中的σ值，卷积过程就是在矩阵u上逐点移动卷积模板，使模板中心和u上的点(x,y,z)重合，卷积模板在每一点(x,y,z)的响应是根据模板中的值与u中对应体素的值相乘再求和得到。

S5：使用adam优化算法针对每个体素处更新的位移自适应调整不同的步长，并更新初始的形变位移场，得到未规则化的空间变换c；该步骤得到的c未进行规则化，即未通过其它措施进行约束，使用其对图像进行空间变换可能会导致不符合实际情况的变形。

原始Demons算法中c←to(α*exp(u))，其中t是初始位移场，α是步长，一般设置为固定值1，u是搜索方向。优选地，本实施例使用adam优化算法自适应调整步长，直接得到α*exp(u)对初始位移场t₀进行更新。

Adam优化算法通过对梯度的一阶矩估计(即梯度的均值)和二阶矩估计(即梯度的未中心化的方差)进行综合考虑，计算出更新步长，在adam算法中梯度是搜索方向，而Demons算法中u是搜索方向，于是，使用S4步骤中得到的更新位移场u替代adam算法中的梯度。

1)基于u更新有偏差的一阶矩和二阶矩m和v，其中：

m_n←β₁·m_n-1+(1-β₁)·u

v_n←β₂·v_n-1+(1-β₂)·u²

式中，β₁是一阶矩估计的指数衰减率(一般设置为0.9)，β₂是二阶矩估计的指数衰减率(一般设置为0.999)。m和v是中间变换矩阵，其初始设置为零矩阵，在每次迭代中更新，n为迭代次数；

2)更新无偏差的一阶矩和二阶矩

和

其中：

3)通过复合运算更新初始位移场t₀得到未规则化的空间变换c：

式中，ε是非常小的数，主要是为了防止在实现中除以零(一般设置10E-8)，η是学习率，默认设置为0.001。

即等价于原始Demons算法中的α*exp(u)。

综上，本实施例不再固定步长，而是根据u的一阶矩和二阶矩自适应估计每个体素在每次迭代中的步长。

S6：对未规则化的空间变换c进行高斯卷积得到更新的形变位移场t_n，其中，t_n←K_diff*c，从而产生扩散模型的规则化效果；

高斯卷积模板的实现以及卷积的实现与步骤S4一致，该步骤中K_diff即是高斯滤波器中的σ值。

S7：使用S6中得到的更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

对于浮动图像f上的每个体素，根据每个体素处对应的位移值(t_n中的值)，对每个体素进行位移。空间变换后的体素点经常不再处于整数位置，需要通过灰度插值处理来计算出该输出点的灰度值，此为现有技术，在此不再赘述。

S8：判断是否迭代终止，是则输出t_n；否则转向步骤S3；

具体地，本实施例迭代终止条件设定为3个，满足其中任意一个，迭代终止。迭代终止后，当前迭代中的形变位移场t_n即是最优的形变位移场，输出即可，其中，3个终止条件分别为：

①当前迭代次数n超过预设的最大迭代次数时，则迭代终止；

②计算参考图像s与变形后的浮动图像f'对应体素灰度差的平方和的平均值，如果该值在三次连续的迭代中持续增大，则迭代终止；

③计算每次迭代后更新的形变位移场t_n与上一次迭代中更新的形变位移场t_n-1之间的差异的平均值，如果该值在连续的三次迭代中持续增大，则迭代终止；

若不符合迭代终止条件(即不符合上述迭代终止条件中的任何一个)，则以当前迭代得到的形变位移场t_n作为初始的形变位移场t₀←t_n，当前迭代得到的变形后的浮动图像作为初始的浮动图像f←f'，迭代次数n←n+1，转向步骤S3继续进行优化，整个优化过程中，参考图像s保持不变。

需要说明的是，上述整个实施例过程中，涉及到两个不同的位移场，即形变位移场t和更新位移场u，形变位移场t是整个配准过程优化的目标，在每次迭代中更新，最终用于浮动图像变形到参考图像。而更新位移场u在每次迭代过程中根据浮动图像和参考图像的灰度以及梯度值进行计算，用于对t进行更新，步骤S5中提及的空间变换c，是t更新过程中的一个中间产物，即未进行规则化(平滑)的t。

实施例2：

本实施例还提供一种微分同胚Demons图像配准系统，该系统包括：

图像选择单元，用于选择参考图像s和浮动图像f；

形变位移场设置单元，用于设置初始形变位移场t₀；

更新位移场计算单元，用于选择Demons力计算更新位移场u；

第一平滑单元，用于对更新位移场u进行高斯卷积；

步长调整单元，用于采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c；

第二平滑单元，用于对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；

空间变换单元，用于采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

判断单元，用于判断是否迭代终止，是则输出更新的形变位移场t_n；否则转向更新位移场计算单元。

实施例3：

本实施例还提供一种计算机程序，包括计算机程序指令，其中，程序指令被处理器执行时用于实现实施例1所述的微分同胚Demons图像配准方法对应的步骤。

实施例4：

本实施例还提供一种计算机可读存储介质，计算机可读存储介质上存储有计算机程序指令，其中，所述程序指令被处理器执行时用于实现实施例1所述的微分同胚Demons图像配准方法对应的步骤。

最后应当说明的是：以上实施例仅用于说明本发明的技术方案而非对其保护范围的限制，尽管参照上述实施例对本申请进行了详细的说明,所属领域的普通技术人员应当理解：本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或者等同替换，但这些变更、修改或者等同替换，均在申请待批的权利要求保护范围之内。

Claims (8)

1.一种微分同胚Demons图像配准方法，其特征在于，包括以下步骤：

S1、选择参考图像s和浮动图像f；

S2、设置初始形变位移场t₀；

S3、选择Demons力计算更新位移场u；

S4、对更新位移场u进行高斯卷积；

S5、采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c，具体过程为：

S51、基于u更新有偏差的一阶矩和二阶矩m和v，其中：

m_n←β₁·m_n-1+(1-β₁)·u

v_n←β₂·v_n-1+(1-β₂)·u²

式中，β₁是一阶矩估计的指数衰减率，β₂是二阶矩估计的指数衰减率，m和v是中间变换矩阵，其初始设置为零矩阵，在每次迭代中更新，n为迭代次数；

S52、更新无偏差的一阶矩和二阶矩

和

其中：

S53、通过复合运算更新初始位移场t₀得到空间变换c：

式中，ξ为常数，η是学习率；

S6、对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；

S7、采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

S8、判断是否迭代终止，若是则输出更新的形变位移场t_n；否则转向步骤S3。

2.根据权利要求1所述的微分同胚Demons图像配准方法，其特征在于，初始形变位移场t₀采用与参考图像大小相同的单位矩阵。

3.根据权利要求1所述的微分同胚Demons图像配准方法，其特征在于，选择Demons力计算更新位移场u，具体为：

Demons力为原始Demons力、主动Demons力或对称Demons力，分别对应于原始Demons算法、主动Demons算法以及对称Demons算法，选择Demons力计算更新位移场u，其中，

原始Demons力的表达式为：

主动Demons力的表达式为：

对称Demons力的表达式为：

式中，f代表浮动图像，s代表参考图像，(f-s)是两幅图像对应体素灰度值的差异，▽f和▽s分别对应浮动图像和参考图像的梯度，u＝(u_x,u_y,u_z)是每个体素在x、y和z方向的位移量，

是归一化因子，u的上限是

4.根据权利要求1～3任一项所述的微分同胚Demons图像配准方法，其特征在于，对更新位移场u采用u←K_fluid*u进行高斯卷积，其中，K_fluid即是高斯滤波器的σ值。

5.根据权利要求1～3任一项所述的微分同胚Demons图像配准方法，其特征在于，对空间变换c通过t_n←K_diff*c进行高斯卷积得到更新的形变位移场t_n，其中，K_fluid是高斯滤波器的σ值。

6.根据权利要求1～3任一项所述的微分同胚Demons图像配准方法，其特征在于，迭代终止条件设定为3个，满足其中任意一个，迭代终止，3个终止条件分别为：

①当前迭代次数n超过预设的最大迭代次数时，则迭代终止；

②计算参考图像s与变形后的浮动图像f'对应体素灰度差的平方和的平均值，如果该值在三次连续的迭代中持续增大，则迭代终止；

③计算每次迭代后更新的形变位移场t_n与上一次迭代中更新的形变位移场t_n-1之间的差异的平均值，如果该值在连续的三次迭代中持续增大，则迭代终止；

若不符合迭代终止条件，则以当前迭代得到的形变位移场t_n作为初始的形变位移场t₀←t_n，当前迭代得到的变形后的浮动图像作为初始的浮动图像f←f'，迭代次数n←n+1，转向步骤S3继续进行优化，整个优化过程中，参考图像s保持不变。

7.一种微分同胚Demons图像配准系统，其特征在于，该系统包括：

图像选择单元，用于选择参考图像s和浮动图像f；

形变位移场设置单元，用于设置初始形变位移场t₀；

更新位移场计算单元，用于选择Demons力计算更新位移场u；

第一平滑单元，用于对更新位移场u进行高斯卷积；

步长调整单元，用于采用adam优化算法对浮动图像每个体素处更新的位移自适应调整步长，并更新初始的形变位移场t₀，得到空间变换c，具体过程为：

基于u更新有偏差的一阶矩和二阶矩m和v，其中：

m_n←β₁·m_n-1+(1-β₁)·u

v_n←β₂·v_n-1+(1-β₂)·u²

式中，β₁是一阶矩估计的指数衰减率，β₂是二阶矩估计的指数衰减率，m和v是中间变换矩阵，其初始设置为零矩阵，在每次迭代中更新，n为迭代次数；

更新无偏差的一阶矩和二阶矩

和

其中：

通过复合运算更新初始位移场t₀得到空间变换c：

式中，ξ为常数，η是学习率；

第二平滑单元，用于对空间变换c进行高斯卷积得到更新的形变位移场t_n，n是迭代次数；

空间变换单元，用于采用更新的形变位移场t_n对浮动图像f进行空间变换得到变形后的浮动图像f'；

判断单元，用于判断是否迭代终止，是则输出更新的形变位移场t_n；否则转向更新位移场计算单元。

8.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质上存储有计算机程序指令，其中，所述程序指令被处理器执行时用于实现如权利要求1～6任一项所述的微分同胚Demons图像配准方法对应的步骤。

Images