Medical Image Non-rigid Registration Based on Adaptive Fractional Order
-
摘要: 现有的医学图像配准算法对于灰度均匀、弱边缘以及弱纹理图像易陷入局部最优从而导致配准精度低下、收敛速度缓慢. 分数阶主动Demons (Fractional active Demons, FAD)算法是解决该问题的有效方法, 并且适用于图像的非刚性配准. 但FAD中的最佳分数阶阶次是人工交互选取, 并且对整幅图像都是固定不变的. 为了解决该问题, 提出一种阶次自适应的主动Demons算法并将其应用到医学图像的非刚性配准中. 算法首先根据图像的局部特征建立分数阶阶次自适应的数学模型, 并逐像素计算最优阶次, 基于该阶次构造Riemann-Liouvill (R-L)分数阶微分动态模板; 然后将自适应R-L分数阶微分引入到Active Demons算法, 在一定程度上缓解了图像配准在弱边缘和弱纹理区域易陷入局部最优问题, 从而提高了配准精度. 通过在两个医学图像库上进行实验验证, 实验结果表明该方法可以处理灰度均匀、弱纹理和弱边缘的医学图像非刚性配准, 配准精度得到较大提升.
-
关键词:
- 自适应分数阶 /
- 主动Demons算法 /
- 自适应模型 /
- 非刚性配准 /
- 医学图像
Abstract: The existing medical image registration methods have limitation in registration images with intensity uniformity, weak edges and weak texture, troubled by inclining to local minimum, which always result in low registration accuracy and efficiency. Active Demons algorithm with fractional differential has been proved to be effective for non-rigid image registration. However, it searches the optimal order of fractional differential manually, lack of order adaptive in images registration. To address the problem, this paper applies multi-resolution and adaptive fractional differential to active Demons, and proposes a novel images registration method for 3D medical images registration. Firstly, a mathematical model of adaptive order for fractional differential is constructed, which adopts local image features such as gradient magnitude and information entropy, therefore the optimal order and differential dynamic template are adjusted adaptively; Secondly, multi-resolution strategy is introduced to adaptive fractional differential active Demons algorithm, optimization falling into local minimum is avoided, therefore the registration accuracy is improved once more. Lastly, extensive experiments show that the proposed algorithm is capable of registration images with intensity uniformity and weak texture. And the optimal order of fractional differential can be calculated adaptively. Furthermore, the presented methods is capable of avoiding falling into local optimum, thus the registration accuracy can be improved greatly. -
医学图像配准[1]是将不同传感器或不同视点或不同时间段获得的同一场景的两幅或多幅图像进行匹配, 目的是寻求一种最优变换(或最优形变场), 使一幅图像与另一幅图像在空间上对齐, 用以纠正图像的形变. 随着计算机和现代医学成像技术的高速发展, 出现了包含各种不同类型信息的图像, 如CT (Computer tomography), MRI (Magnetic resonance imaging)等, 单模态图像为医生提供单一片面的信息, 要得到更加完整且互补的图像信息, 需要将包括不同类型信息的多模态图像进行配准; 此外, 在不同的视点或不同时间段对同一个人的某个器官进行拍摄, 其图像也存在差异, 为了做出更加准确可靠的判断, 也需要将不同视点或不同时段的图像进行配准. 图像配准是医学图像理解与分析的基础, 医学疾病的诊断和治疗、模拟手术导航等应用都需要图像的精确配准作为前提条件. 目前医学图像配准还存在许多急需要解决的难题, 如图像局部存在对应缺失和较大形变, 器官自身的不规则生理运动, 图像边缘模糊、纹理结构不清晰等特点, 传统的方法对配准该类图像的效果不够理想.
基于光流场模型的方法由于检测精度较高和稳定性较好得到学者们的重视. 其中, Demons光流场的配准方法[2]具有快速和高效的特点, 所以基于Demons光流场的配准方法在医学图像上得到了广泛的应用. 经典Demons算法的基本思想是将配准看作浮动图像像素在参考图像像素灰度梯度信息驱动下向参考图像逐步扩散的过程. 但是利用参考图像的梯度作为驱动力驱动变形图像时, 当梯度信息不足, 容易出现匹配错误, 并且只适合配准较小形变的图像. Wang等[3]在Demons算法的基础上提出了主动Demons算法, 允许参考图像和变形图像的梯度共同驱动像素点朝着对方对应的像素点移动, 也即同时使用参考图像和变形图像的梯度作为驱动力. 文献[3-6]的实验结果表明, 将经典Demons中参考图像的单向驱动力转换为双向力, 能处理较大形变的图像配准, 即使参考图像的梯度很小时, 也能得到较高的配准精度. Vercauteren等[7]提出将Demons算法与微分同胚相结合, 保证了变形场的可逆性、可微性和空间点的一一对应, 阻止了变形空间的折叠. Hao等[8]针对较大变形和多模态问题, 将传统Demons 算法和局部结构张量相结合, 获得了较好的配准效果. 针对肺部器官自身的不规则生理运动造成的肺部图像大形变问题, Lu 等[9]在传统Demons算法的基础上进行改进, 提出一种新的精度更高的Demons算法, 实验结果表明, 改进的配准算法能有效配准较大变形图像. 闫德勤等[10]将流形学习的思想引入到微分同胚的Demons 中, 提出了适合较大形变的图像配准算法, 图像的拓扑结构得到较好的保持, 配准精度也得到提高. 薛鹏等[11]将平衡系数引入到主动 Demons 算法中, 用弹性系数与平衡系数共同来调节驱动力的大小, 从而提高配准精度和配准效率, 该算法可以较好配准一般的医学图像, 但是对于灰度均匀和弱纹理图像的配准不够理想. Lu等[12]针对螺旋CT图像存在较大噪声问题, 使用了微分同胚的Demons算法, 实验结果表明该算法可较完整地得到清晰的肺通气功能图. 但是, 上述基于扩散理论的Demons 算法的驱动力均来自于图像的灰度梯度, 当图像的局部区域存在低对比度造成的灰度均匀、弱边缘和弱纹理, 优化易陷入局部极值导致配准精度低下. 针对该问题, 张桂梅等[13]提出将分数阶微积分引入到Active Demons 中, 提出了分数阶梯度驱动的Active Demons算法, 增强了图像的梯度驱动力, 提高了图像配准的精度和效率, 适合于灰度均匀、弱纹理和弱边缘的图像配准效果, 但是该方法的分数阶阶次需要通过人工交互性选取. 目前, 随着深度学习的兴起, 有学者将深度学习的方法运用到变形医学图像的配准. Hu等[14]利用图像的局部分割块来训练ConvNets, 得到相应的网络模型来完成全局和局部的图像配准任务. Balakrishnan等[15]提出了基于学习的算法完成3D变形医学图像对的配准任务, 该方法将配准任务转换为参数化函数, 针对感兴趣的区域, 对参数函数进行最佳化处理. 但是基于块的无监督配准方法, 需要后加工进行处理, 并且后处理不能在卷积神经网络中直接进行, 同时它的配准质量评价仍依赖于其他基于特征的方法. 基于深度学习的图像配准方法虽然在不同的数据集上获得了较好配准性能. 但配准模型依赖于对较完善的配准样本的训练, 大部分用于配准的训练样本都是合成转换参数得到的, 并且还需要手工标注出相应的信息. 对于实际的医学图像, 由于结构复杂, 合成的训练集很难达到与真实图像接近, 从而使得训练出的配准模型精度不够理想.
分数阶微积分是在传统整数阶微积分的基础上延伸和发展出来的, 因为它具有长记忆性, 非局部性和弱奇异性, 能在保持信号中、高频成分的同时, 非线性增强信号的低频成分, 并且阶次更灵活、连续可调, 因此在图像处理如图像增强、图像去噪、图像边缘提取和图像分割等领域均得到较广泛的应用. Pu等[16]和Chen等[17]将分数阶微积分应用到图像去噪中, 不仅提高了图像的峰值信噪比, 而且纹理保持效果也得到有效提升. Mathieu等[18]将分数阶微分用于边缘检测, 实验结果表明, 分数阶微分的边缘掩模算子的结果优于整数阶微分掩膜算子. Ren等[19]提出了一种自适应分数阶能量驱动的活动轮廓模型, 并将其应用在图像分割中, 但该方法仅将分数阶微分拟合项增加到现有的拟合方程中, 分数阶最佳阶次仍是通过实验凭经验选取. 虽然分数阶微积分的算法在处理灰度均匀, 弱边缘、弱纹理图像具有一定的优势, 但是, 其最佳分数阶微积分的阶次需要通过多次实验人工选取, 费时费力. 张桂梅等[20]提出了基于自适应分数阶的可变区域拟合(Region-scalable fitting, RSF)模型, 解决了RSF模型在分割弱纹理、弱边缘图像时, 演化曲线易陷入局部极值的问题. Yu等[21]针对图像在去噪的同时易使纹理细节丢失并使边缘模糊的问题, 基于图像的局部特征信息构造了分数阶阶次自适应的分段函数, 实验结果表明, 该算法在去除噪声的同时较好地保存了图像的纹理和边缘.
针对医学图像非刚性配准中的灰度均匀、弱纹理和弱边缘以及阶次的自适应问题, 本文在文献[13]的基础上, 提出了一种新的医学图像非刚性配准方法. 本文主要工作有:
1)根据图像的局部信息熵和梯度模值建立了阶次自适应的分数阶模型, 自动计算最佳的分数阶阶次, 同时基于最佳阶次构造了自适应分数阶阶次的动态模板, 解决了人工寻找最佳阶次缺乏阶次自适应的问题;
2)将自适应分数阶引入到主动Demons算法, 在一定程度上缓解了在图像弱纹理和弱边缘区域的配准陷入局部极小值, 提高了图像配准的精度;
3)将提出的算法应用于医学图像的非刚性配准, 实验结果表明本文的算法可以处理灰度均匀、弱纹理和弱边缘的医学图像配准, 配准精度得到有效提升.
1. 背景知识
1.1 主动Demons算法
图像配准可以近似为从变形图像在驱动力的作用下流到参考图像的过程. 因此, 由光流场模型得到的速度场可以作为图像配准的位移场. Demons算法及各种改进的Demons算法仅利用图像的梯度信息驱动变形图像与参考图像接近, 但当图像中的某些区域灰度接近时, 则该区域的梯度值将很小, 甚至为零, 此时驱动不足易导致配准错误, 此外, 传统的Demons 算法只适合于存在小变形的图像, 当变形较大时图像配准的精度则下降较多. 针对该问题, Wang等[3]提出了Active-Demons算法, 主要是将两个绝对梯度图像中的两个单向驱动力转换为双向驱动力, 以此来提高配准的精度和收敛的速度, 表达式为
$$ \begin{split} u =\;& (M-F)\left(\frac{\nabla F}{|\nabla F|^{2}+(M-F)^{2}}+\right.\\ &\left.\frac{\nabla M}{|\nabla M|^{2}+(M-F)^{2}}\right) \end{split} $$ (1) 其中,
$ F $ 和$ M $ 分别为参考图像和变形图像的灰度,$ \nabla F $ 和$ \nabla M $ 则分别为它们的灰度梯度. Active-Demons算法的驱动力是双向的, 分别来自参考图像和变形图像的梯度以及变形图像和参考图像的梯度, 故可较好地配准大形变的图像. 为便于调整驱动力大小, 在式(1) 中加入了均化系数$ \beta , $ 则有$$ \begin{split} u =\;& (M-F)\left(\frac{\nabla F}{|\nabla F|^{2}+\beta^{2}(M-F)^{2}}+\right.\\ &\left.\frac{\nabla M}{|\nabla M|^{2}+\beta^{2}(M-F)^{2}}\right) \end{split} $$ (2) 通过调整
$ \beta $ , 不仅可提高配准精度, 且可加快收敛速度.1.2 分数阶微分定义
分数阶微积分理论是分形学的基础之一, 它在图像处理领域受到了广泛的关注. 对任意函数
$ f(x) $ 在区间$ (a, x) $ 上的$ \alpha $ 阶微积分可表示为$$ _aD_x^\alpha f(x) = \left\{ \begin{array}{l} \dfrac{{{{\rm{d}}^\alpha }}}{{{\rm{d}}{{{\tau}}^\alpha }}}f(x),\;\;\;\;\;\;\;\;\;\;\;\alpha > 0\\ f(x),\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\alpha = 0\\ \displaystyle\int_a^x f (\tau ){\rm{d}}\tau ,\;\;\;\;\;\;\;\;\alpha < {\rm{0}} \end{array} \right. $$ (3) 其中,
$ \alpha $ 表示阶次, 当$ \alpha>0 $ 时,$ {}_aD_x^\alpha $ 是微分算子, 当$ \alpha<0 $ 时, 则$ {}_aD_x^\alpha $ 是积分算子.通常式(3)中的
$ \alpha $ 取整数, 若当$ \alpha $ 为分数或小数时, 则由整数阶微积分演变为分数阶微积分. 分数阶微积分的三种经典定义有Grünwald-Letnikov (G-L)定义、Riemann-Liouville (R-L)定义和Caputo定义[16]. R-L定义与G-L定义可通过卷积计算得到, 故广泛用于信号处理领域, 但R-L定义的计算简单, 表达式更清晰, 故我们选用R-L定义, 构造分数阶微分模板, 并应用大形变的医学图像非刚性配准中.2. 本文算法
2.1 结合分数阶的Active Demons算法
2.1.1 R-L微分模板的构造
函数
$ f(x) $ 在区间上$ (a, x) $ 的整数一阶、二阶积分分别如式(4)和式(5)所示$$ \frac{{\rm{d}}^{-1}}{{\rm{d}}x^{-1}}f(x) = \int_a^x f(\tau){\rm{d}} \tau \qquad\qquad\qquad\qquad\qquad\;$$ (4) $$ \frac{{\rm{d}}^{-2}}{{\rm{d}}x^{-2}}f(x) = \iint_D{f(\tau)}{\rm{d}}\tau {\rm{d}}x = \int_a^x f(\tau)(x-\tau){\rm{d}} \tau $$ (5) 同理,
$ f(x) $ 的$ n $ 阶积分为$$ \frac{{\rm{d}}^{-n}}{{\rm{d}}x^{-n}}f(x) = \frac{1}{(n-1)!}\int_a^x f(\tau)(x-\tau)^{n-1}{\rm{d}} \tau $$ (6) 将整数
$ n $ 扩展为分数, 并引入Gamma函数, 则可推导出$ \alpha $ 阶分数阶积分的R-L定义为$$ {}_aD_x^{-\alpha}f(x) = \frac{1}{\Gamma(\alpha)}\int_a^x f(\tau)(x-\tau)^{\alpha-1}{\rm{d}} \tau $$ (7) 其中
$$ \Gamma(\alpha) = \int_0^\infty {\rm e}^{-t} t^{\alpha-1}{{\rm{d}}{t}} $$ (8) 若正实数
$ \alpha $ 满足条件:$ 0\leq(n-1)<\alpha\leq n, $ 其中$ n $ 为不小于$ \alpha $ 的最大整数, 则$ f(x) $ 的$ \alpha $ 阶微分的R-L 定义为$$ {}_aD_x^{\alpha}f(x) = \frac{1}{\Gamma(n-\alpha)}\frac{{\rm{d}}^{{n}}}{{\rm{d}}x^n}\int_a^x (x-\tau)^{n-\alpha-1}f(\tau){\rm{d}} \tau $$ (9) 特别地, 当微分阶次
$ \alpha $ 在区间$ (0, 1) $ 取值时,$ n = 1, $ 则式(9)变为$$ {}_aD_x^{\alpha}f(x) = \frac{1}{\Gamma(1-\alpha)}\frac{\rm{d}}{{\rm{d}}x}\int_a^x (x-\tau)^{-\alpha}f(\tau){\rm{d}} \tau $$ (10) 将式(10)改写为卷积分形式
$$ {}_aD_x^{\alpha}f(x) = \frac{x^{-\alpha}*f(x)}{\Gamma(1-\alpha)} = h(x, \alpha)*f(x) $$ (11) 其中
$$ h(x, \alpha) = \frac{x^{-\alpha}}{\Gamma(1-\alpha)} \nonumber $$ 式(11)中,
$ h(x, \alpha) $ 为通项式系数.将通项式系数扩展到二维, 则有
$$ h(x, y, \alpha) = \frac{({x^2+y^2})^{-\frac{\alpha}{2}}}{\Gamma(1-\alpha)}, \qquad (x, y)\neq(0, 0) $$ (12) 式(12)分别对
$ x, y $ 求一阶导数, 那么可获得图像在$ x, y $ 的分数阶微分模板, 如下:$$ \begin{split} &H_x(x, y, \alpha) = -\frac{\alpha x}{\Gamma(1-\alpha)}({x^2+y^2})^{-\frac{\alpha}{2}-1}\\ &H_y(x, y, \alpha) = -\frac{\alpha y}{\Gamma(1-\alpha)}({x^2+y^2})^{-\frac{\alpha}{2}-1} \end{split} $$ (13) 分别将
$ H_x(x, y, \alpha) $ 和$ H_y(x, y, \alpha) $ 离散化, 则:$$ \begin{split} &H_x(x, y, \alpha) = -\frac{\alpha x_M}{\Gamma(1-\alpha)}({x_M^2+y_M^2})^{-\frac{\alpha}{2}-1}\\& H_y(x, y, \alpha) = -\frac{\alpha y_M}{\Gamma(1-\alpha)}({x_M^2+y_M^2})^{-\frac{\alpha}{2}-1} \end{split} $$ (14) 其中,
$x_M = -K, -K+1, \cdots, K-1, K,\ y_M = -L,$ $-L +1, \cdots,L-1, L.$ 若模板大小为$ M\times N,K, L $ 分别为不大于$ M/2 $ 和$ L/2 $ 的最大整数. 即:$ H_x(x, y, \alpha) $ 和$ H_y(x, y, \alpha) $ 是大小为$ (2K+1)\times(2L+1) $ 的微分模板. 如模板大小为$ 5\times 5, $ 则$ K, L $ 分别为2. 分别将$ H_x(x, y, \alpha) $ 和$ H_y(x, y, \alpha) $ 与图像进行卷积运算则可获得x和y方向的分数阶梯度, 从而得到x和y轴方向的形变场.2.1.2 分数阶主动Demons算法
现有的各种Demons方法均以图像的灰度梯度作为驱动力, 当图像中某区域灰度均匀、存在弱纹理和弱边缘时, 配准优化易陷入局部最优. 但是, 分数阶微分不仅可增强图像的细节信息, 且可保留图像的平滑区域信息[22]. 因此我们用分数阶梯度替代图像中的整数阶梯度, 则有
$$ \begin{split} u_x =\;& (M-F)\Biggr(\frac{H_x(x, y, \alpha)* F}{{|H_x(x, y, \alpha)* F}|^2+\beta ^2(M-F)^2}+\\ &\frac{H_x(x, y, \alpha)* M}{{|H_x(x, y, \alpha)* M}|^2+\beta ^2 (M-F)^2}\Biggr)\\ u_y =\; &(M-F)\Biggr(\frac{H_y(x, y, \alpha)* F}{{|H_y(x, y, \alpha)* F}|^2+\beta ^2(M-F)^2}+\\ &\frac{H_y(x, y, \alpha)* M}{{|H_y(x, y, \alpha)* M}|^2+\beta ^2 (M-F)^2}\Biggr) \\[-20pt] \end{split} $$ (15) 文献[13]的研究证明了分数阶Active Demons算法能同时提升配准的精度和效率. 但是对不同的图像, 其最佳的阶次均不相同, 需经过多次实验寻找最佳阶次, 费时费力, 且缺乏自适应性.
2.2 分数阶自适应数学模型的构造
2.2.1 图像中局部特征的选择
根据文献[13]的研究结果可知, 若图像中存在灰度均匀、弱纹理和弱边缘时, 主动Demons算法在配准过程中易陷入局部最优, 但将分数阶引入到Active Demons算法中, 分数阶阶次的大小将影响配准的结果, 该文献通过多次实验手动调整分数阶的最佳阶次, 缺乏阶次调整的自适应性. 基于此, 本文尝试构造自适应阶次的分数阶模型. 为了构造该模型, 首先需要寻求能反映图像局部信息的特征量, 并得到图像的局部特征与分数阶阶次的关系. 由于图像梯度反映了图像灰度在空间上的变化率, 表现在图像上为梯度大的地方, 其像素灰度变化明显; 梯度变化小的地方, 像素灰度则变化平缓; 灰度相同区域, 则梯度接近于零. 图像梯度可以认为是纹理的量化反映. 图像的信息熵[23]反映了图像纹理信息的丰富程度, 在图像的边缘具有更大的图像信息熵, 而在图像平滑和纹理区域信息熵则更小. 故我们选用图像局部信息熵和梯度模值作为局部信息特征, 它们的概念和计算公式如下.
定义1. 图像的梯度反映了图像中灰度的变化情况, 在图像的纹理或边缘区域, 灰度变化较其他地方更大, 也即梯度值较大, 在图像的平滑区域, 灰度变化小, 即图像梯度值较小. 图像梯度可表示为
$$ G[f(x, y)] = \begin{bmatrix} G_x \\ G_y \end{bmatrix} = \begin{bmatrix} {\dfrac{\partial I}{\partial x}} \\ {\dfrac{\partial I}{\partial y}}\end{bmatrix} $$ (16) 梯度模值为
$$ mag(G[f(x, y)]) = |G| = \sqrt{G_x^2+G_y^2} $$ (17) 定义2. 信息熵也能反映图像的局部特征信息, 在图像的纹理或边缘区域, 信息熵通常比较大, 在图像的平滑区域, 信息熵通常比较小, 其定义为
$$ S = -\sum\limits_{{I_{ij}}\in\omega}P_{ij}\log_n(P_{ij}) $$ (18) 其中,
$ S $ 是信息熵,$ I_{ij} $ 为在像素点$ (i, j) $ 处的灰度值,$ \omega $ 为模板,$ P_{ij} $ 为在模板内具有相同灰度值的概率.2.2.2 建立分数阶阶次自适应的数学模型
文献[13]的研究结果表明, 对于不同的图像, 其最佳阶次不同, 也即是图像的最佳阶次与图像的特征有关, 图像的局部信息熵和梯度模值能较好地反映图像的局部特征, 故它们对分数阶最佳阶次的确定有较大影响. 图像某区域的局部信息熵、梯度模值越大, 则该区域为边缘区域可能性越大; 反之该区域为弱纹理或平滑区域的可能性越大; 对于数字图像, 弱纹理和平滑区域对应于信号的中低频; 边缘和噪声对应于信号的高频, 也即图像的某局部区域上述两个特征信息值越大, 该区域越可能为高频区域, 反之越可能为中低频区域. 根据分数阶微积分的特性可知, 分数阶微分算子在阶次选择适当时可以在增强信号的高频成分的同时非线性地保留信号的中低频; 且提升高频信号的幅度随着分数阶阶次的增加急速增长, 即阶次越大, 提升的幅度就越大; 在保留信号中的中低频成分时, 随着分数阶阶次的增大保留的幅度越小. 总之, 要增强高频信号, 需要选择较大的阶次; 要保留中低频信号, 则需要选择较小的阶次. 反映到弱纹理和弱边缘区域的图像配准上, 需要增强弱纹理、弱边缘和平滑区域的局部梯度驱动力, 需要选择较小的阶次, 相对而言, 图像的边缘区域则需要选择较大的阶次; 由此表明图像某区域处的梯度值和信息熵越大, 该区域应选择较大的阶次, 反之则应选择较小的阶次进行处理. 由此分析可知, 分数阶阶次与图像的局部信息相关联, 故应根据图像的局部信息选取不同的阶次, 而不是整幅图像均选择一个固定的阶次. 由于图像的信息熵和梯度模值均可反映出图像中纹理的丰富程度, 故我们根据图像这两个局部信息特征构造阶次与图像局部特征自适应变化的模型, 从而逐像素地计算最优阶次.
将图像的信息熵和梯度模值进行归一化, 使其满足
$|G|\in[0, 1] , |S|\in[0, 1] ,$ 然后再将它们进行融合, 如下式:$$ f(G,S) = m\left| G \right| + n\left| S \right|,\qquad 0 \le f \le 1 $$ (19) 其中,
$ m $ 和$ n $ 为权值, 根据多次实验, 选择$m = n =$ $0.5$ .从分数阶的幅频特性曲线[24]可知: 当阶次为0时, 对信号无增强效果; 当阶次大于1时, 对信号的中低频成分提升的幅度偏低, 也即对图像的弱纹理和平滑区域增强不明显. 由于本文的目的是需要增强图像中弱纹理和平滑区域的梯度, 所以分数阶的阶次应位于[0, 1]区间. 所需构造的自适应函数必须单调有界, 也即是通过自适应函数计算得到的阶次要满足两个条件, 首先阶次位于[0, 1]区间, 其次随着图像的局部特征的增加而单调递增. 而如图1所示的反正切函数正好能满足这个要求, 选择反正切函数的右侧作为原型函数, 因为它是单调递增并有界. 所以我们以反正切函数为原型函数, 建立自适应函数:
$$ \alpha = k\times\arctan{(f)}+b $$ (20) 其中,
$ \alpha $ 为阶次,$ k ,b $ 为待定的系数. 在信息熵及图像的梯度模值较大处,$ f $ 应取较大值, 此时对图像高频区域的增强程度较大, 所选取的分数阶阶次应该较大; 同理, 在信息熵和图像梯度模值较小处,$ f $ 应该取较小值, 此时对低频区域信息的保留程度应该较大, 所选取的分数阶阶次应该较小. 那么, 令$ f = 1 $ 时,$ \alpha = 1; $ 令$ f = 0 $ 时,$ \alpha = 0. $ 于是有$$ \left\{ \begin{array}{l} k \times \arctan \left( 1 \right) + b = 1\\ b = 0 \end{array} \right. $$ (21) 从式(21)可以得出
$ k = 4/\pi , b = 0 ,$ 并将其代入式(20), 有$$ \alpha = \dfrac{4}{\pi}\times\arctan{(f)} $$ (22) 式(22)即为自适应分数阶阶次的数学模型.
2.3 多分辨率策略的图像配准
构造自适应分数阶阶次数学模型后, 图像的每个像素点都对应一个最佳分数阶阶次, 由于分数阶微分模板是动态计算得到, 并且计算分数阶微分模板比整数阶微分模板更为复杂, 所以提出的模型相较于固定阶次模型的配准时间会有所增加. 但是, 将多分辨率策略引入到自适应分数阶的主动Demons中, 可以适当减少配准时间.
本文利用下采样的多分辨率策略对图像进行分层配准, 这是因为下采样的多分辨率策略相比较于小波变换的多分辨率策略具有如下优点: 1) 能提高信号的峰值信噪比; 2) 下采样是采用隔点采样的方式, 对图像信息进行了压缩存储, 从而可以提高采用效率, 降低时间成本. 我们使用下采样方法分解输入的图像, 得到两层不同分辨率的变形图像和参考图像, 首先在低分辨率层实现图像的初配准, 并将该配准结果作为下一层的输入, 然后在高分辨率层实现图像的精配准.
2.4 本文的算法步骤
步骤1. 导入变形图像和参考图像;
步骤2. 使用下采样方法分解导入的图像, 得到两层不同分辨率的变形图像和参考图像;
步骤3. 对于相同分辨率层的变形图像和参考图像, 基于图像的信息熵和梯度模值建立阶次自适应的分数阶模型, 并逐像素计算最优阶次, 基于该阶次构造其R-L分数阶微分动态模板;
步骤4. 由式(15)计算变形图像与参考图像之间的形变场
$ u $ ;步骤5. 用高斯函数平滑处理配准后的变形图像与参考图像间的位移形变场;
步骤6. 如果收敛, 转步骤7, 否则转步骤3;
步骤7. 如果是最后一层, 结束配准, 否则转步骤3, 继续配准下一层图像.
3. 实验
3.1 实验图像库
为了更公正地评价各算法的性能, 通常都在标准医学图像库中进行测试对比. 本文使用了两种不同的医学图像库进行实验验证, 分别是Whole Brain Atlas图像库和BrainWeb图像库. 本文主要是针对脑部3D图像不同层的切片图像进行配准研究. 我们选取了3D脑部图像中的3个有代表性、对比度较高的切片层, 以及脑部图像的冠状面、矢状面和横截面图像进行配准实验. 对于第1个图像库中图像的配准测试我们选用了不同模态和不同尺度的切片图像, 第2个图像库采用了其不同模态下的冠状面、矢状面和横截面图像.
3.2 计算机环境及性能评价标准
本文计算机环境如下: Win7 32位, MATLAB R2010a, 计算机内存2 GB, CPU 主频3.2 GB.
相似性度量是判断图像配准质量的常见方法之一. 相似性度量的评价指标较多, 但每个度量均有其自身的特点, 故相似性度量的选择要以图像的特征和性质为依据. 本文采用均方误差(Mean square error, MSE)和重叠比(Dice ratio)作为图像配准结果的评价指标, 另外我们还展示了主观视觉对比结果. 其中定量评价指标中的MSE表示参考图像和配准后的图像中的体素差异, 用于评估图像对之间的灰度分布的相似性和局部的对比性, Dice ratio表示配准后的图像与参考图像在同一个结构区域中的像素位置的重叠比, 其表达式分别如式(23)和式(24)所示:
$$ MSE(f, m(\phi)) = \frac{1}{|\Omega|}\sum\limits_{p\in\Omega}[f(p)-[m(\phi)](p)]^2 $$ (23) 其中,
$ f $ 表示参考图像,$ m $ 表示变形图像,$ m(\phi) $ 表示配准后的图像;$ \phi $ 表示变形图像与参考图像之间对应的形变场,$ p $ 表示单个体素,$ \Omega $ 表示图像对应的空间像素域.$$ Dice(A, B) = \frac{|A|\bigcap|B|}{|A|\bigcup|B|} $$ (24) 其中,
$ A $ 表示参考图像中的像素标签位置,$ B $ 表示配准后的图像中与参考图像位置对应的像素标签位置.本文迭代收敛的条件为MSE随着迭代次数增加不再变小, Dice ratio 随着迭代次数增加不再变大, 即达到一个稳定值.
3.3 实验结果及分析
3.3.1 BrainWeb图像库
为测试提出算法的有效性, 我们选用了BrainWeb的图像分别在不同模态及不同尺度进行实验验证, 并与文献[7]和文献[13]的方法进行对比.
本次实验中, 我们选择了具有代表性和比较典型的3个切片层图像进行实验比较, 如以MR-T1模态, 1 mm厚不同层的切片图像I、II、III为参考图像, 其大小为
$181\times217\times181,$ 如图2(a)所示; 以MR-PD模态、切片厚度为3 mm的不同层的切片图像I、II、III 为变形图像、其大小为$181\times 217\times$ $181 ,$ 如图2(b)所示. 用文献[7]、文献[13]和本文方法分别进行实验, 配准结果如图2所示. 其中图2(c)、图2(d)、图2(e)分别为文献[7], 文献[13]和本文方法的配准结果图. 从图2(c)可以看出, 经过配准纠正后的图像与变形图像相比有所变化, 但是与参考图像的形状、灰度均不接近, 也即是经过配准后, 变形图像没有得到有效矫正. 观察图2(d)可以发现, 配准后的图像底部变尖了, 与参考图像还是不够接近, 并且图像的灰度与参考图像相差较大, 形状与参考图像也不接近, 这表明用文献[7]和文献[13]的方法对输入的参考图像和变形图像配准失败. 再观察图2(e)可知, 配准后图像底部的那条线变清晰了, 图像的亮度、形状与参考图像均较接近, 从而表明本文方法配准效果较好. 对比实验结果的定量分析如表1和表2所示, 从表1和表2可知, 采用文献[7]和文献[13]的算法配准后的图像MSE有所降低, Dice ratio有所增加, 但是性能提升幅度不大, 而本文方法得到的配准后图像与参考图像之间的$ MSE $ 有较大减少, Dice ratio有较大增加, 从而表明本文算法的配准性能要优于文献[7]和文献[13]的算法. 这是由于本文将自适应分数阶引入到主动Demons算法, 并自适应计算了图像每个像素点的最佳分数阶阶次, 从而在一定程度上缓解了图像配准优化陷入局部极值, 提高了配准精度. 所以, 从定性分析和定量结果都可以看出, 本文提出的算法对不同模态医学图像的非刚性配准的效果较好.表 1 均方误差比较Table 1 Mean square error comparison
图 2 不同切片层的BrainWeb图像配准结果Fig. 2 Registration result of BrainWeb image of different slice layers表 2 Dice ratio比较Table 2 Dice ratio comparison3.3.2 Whole Brain Atlas图像库
为了验证本文方法的有效性, 我们进一步用Whole Brain Atlas 图像库中的图像进行实验, 本次实验中, 我们选择了2个不同模态的不同端面的图像进行实验, 如分别选择模态MR-T1和模态MR-T2的冠状面、矢状面和横截面图像进行对比实验, 同时与文献[7]和文献[13]的算法进行对比.
3.3.2.1. 冠状面图像配准
我们选不同模态的脑部冠状面图像进行实验, 结果如图3所示. 其中图3(c)、图3(d)和图3(e)分别是文献[7]、文献[13]和本文方法的配准结果. 观察图3可知, 图3(c)和图3(d)都存在变形, 脑部图像的部分细节未得到正确配准, 而本文方法的配准效果较好, 除了图像的亮度稍暗外, 脑部图像的细节也得到成功匹配, 并且与参考图像的形状、轮廓均接近. 我们再进行了相应的定量分析实验, 结果如表3, 本文算法得到的两个评价指标相比文献[7]、文献[13]均有一定程度的改善, 如MSE得到有效减少, Dice ratio得到较大幅度提高, 故本文算法配准效果最好. 这是因为我们在主动Demons算法中融入了分数阶微积分, 并根据图像的局部特征, 计算出了图像各个像素点的最佳分数阶阶次.
表 3 冠状面配准精度对比Table 3 Comparison of registration accuracy of coronal plane3.3.2.2. 矢状面图像配准
我们选择了不同模态的脑部矢状面图像进行实验, 结果如图4, 其中图4(a)是参考图像, 图4(b)是变形图像, 图4(c)、图4(d)和图4(e)分别是文献[7]、文献[13]和本文方法的配准结果. 从图4可以看出, 图4(c)的结果与参考图像相差较远, 图4(d)次之, 图4(e)与参考图像很接近. 这说明文献[13]的算法要优于文献[7], 这是因为文献[13]在调节驱动力的时候引入了分数阶梯度驱动, 但是由于文献[13]中的分数阶阶次对整幅图像都是一个固定值, 而本文方法在文献[13]的基础加入了分数阶阶次的自适应, 即能根据图像的局部信息特征自动计算各个像素点的最佳阶次, 所以效果最好. 该实验的定量分析结果展示在表4中, 观察表4可知, 本文方法的两个评价指标均是最优的. 这是因为本文算法引入了自适应分数阶, 增强了图像的梯度驱动力, 图像配准的精度得到提高, 从而证明了提出的方法对于矢状面图像配准是有效的.
表 4 矢状面配准精度对比Table 4 Comparison of registration accuracy of sagittal plane3.3.2.3. 横切面图像实验
我们选择了不同模态的脑部横切面图像进行配准实验, 实验结果如图5, 其中图5(c)、图5(d)和图5(e)分别是文献[7]、文献[13]和本文方法的配准结果. 从图5可以看出, 文献[7]和文献[13]的算法都没有配准成功, 而本文算法的配准结果除了图像亮度稍有差别外, 其他的都与参考图像比较接近. 再从表5中可以观察到, 本文方法的两个评价指标均是最优的. 这是因为我们建立了图像局部信息与分数阶阶次之间的关系, 逐像素计算了最优阶次, 此外因为分数阶具有弱导数的特点, 对于图像灰度均匀区域的信息有较好的保留, 而对图像弱边缘和弱纹理区域的信息有较大的增强.
表 5 横切面配准精度对比Table 5 Comparison of registration accuracy of transverse plane3.3.3 配准时间实验
本文用配准时间对配准的效率进行比较, 使用BrainWeb图像库中的脑部图像进行测试. 由于文献[13]中采用固定的分数阶阶次对图像进行处理, 并且寻找最佳的分数阶阶次是通过人工多次实验选取, 本文在该文献的基础上建立了分数阶阶次自适应的数学模型, 本节将比较文献[13]和本文算法的配准时间, 实验结果如表6所示. 从表6可以看出, 本文提出算法的配准时间较多, 这是因为本文算法需要逐像素计算最佳阶次, 并基于每个最佳阶次构造分数阶微分模板; 而文献[13]的算法只是通过交互式选择了1个固定的阶次, 整幅图像也只对应1 个分数阶微分模板, 故本文算法的时间开销更大. 但是使用文献[13]的算法寻找最佳阶次时需要通过多次实验, 也即是每次都需要更换1个分数阶阶次重复进行实验. 例如阶次从0.1~0.9需要做9次实验, 这样9次实验的总时间则比本文算法配准时间更多. 总之, 本文算法相对于文献[13]的单次实验时间开销更大, 但对于多次实验的时间总和则可能具有一定程度的优势, 并且不需要通过多次实验寻找最佳阶次, 实现了分数阶阶次的自适应计算.
表 6 不同算法的时间对比(s)Table 6 Time comparison of two methods (s)不同切片层的图像 文献 [13] 的方法 (不同的阶次) 本文方法 $ \alpha $= 0.1 $ \alpha $= 0.2 $ \alpha $= 0.3 $ \alpha $= 0.4 $ \alpha $= 0.5 $ \alpha $= 0.6 $ \alpha $= 0.7 $ \alpha $= 0.8 $ \alpha $= 0.9 总计时间 I 3.21 3.14 2.98 2.76 3.03 2.89 2.92 2.67 2.58 26.18 17.69 II 3.72 3.65 3.37 3.54 3.68 3.03 3.29 3.21 3.16 30.65 19.08 III 4.64 4.61 4.53 4.57 4.65 4.06 4.52 4.18 4.49 40.25 26.83 由于本文引入了自适应分数阶微分, 需要逐像素计算分数阶最佳阶次, 从而增加了算法的计算量, 基于此, 我们加入多分辨率策略. 同时, 我们也对比了用多分辨率策略和不采用多分辨率策略进行配准所花费的时间, 仍使用BrainWeb图像库中的脑部图像进行测试. 实验结果如表7所示, 从表7可知, 采用多分率策略进行配准后, 时间开销有了较大的减少, 这也是大部分配准算法采用多分辨率策略配准的原因.
表 7 两种策略的时间对比(s)Table 7 Time comparison of two strategies (s)图像 不采用多分辨率 采用多分辨率 I 30.25 17.69 II 28.04 19.08 III 40.63 26.83 4. 结论
将自适应R-L分数阶微分引入到主动 Demons算法中, 解决了医学图像中存在灰度均匀、弱纹理和弱边缘区域的配准精度低下问题. 基于图像的局部信息建立了分数阶阶次自适应的数学模型, 逐像素计算了最优阶次, 并根据该阶次建立了动态的微分模板, 实现了分数阶阶次的自适应性; 再则, 将自适应分数阶微分引入到Active Demons算法中, 并采用多分辨率策略对不同分辨率的图像分别进行了配准, 在一定程度上缓解了图像配准陷入局部最优, 从而提高了配准精度. 用医学库的图像进行了大量实验, 结果表明本文提出的算法较好地实现了医学图像的非刚性配准, 建立了图像的局部特征与最佳分数阶阶次的对应关系, 实现了阶次的自适应, 配准精度得到有效提高. 本文主要是针对脑部3D图像不同层的切片图像进行配准研究, 后续我们将考虑研究对整体的3D图像进行配准. 此外还将研究更多类型的医学图像非刚性配准.
-
图 2 不同切片层的BrainWeb图像配准结果
Fig. 2 Registration result of BrainWeb image of different slice layers
表 3 冠状面配准精度对比
Table 3 Comparison of registration accuracy of coronal plane
下载: 导出CSV
表 4 矢状面配准精度对比
Table 4 Comparison of registration accuracy of sagittal plane
下载: 导出CSV
表 5 横切面配准精度对比
Table 5 Comparison of registration accuracy of transverse plane
下载: 导出CSV
表 6 不同算法的时间对比(s)
Table 6 Time comparison of two methods (s)
不同切片层的图像 文献 [13] 的方法 (不同的阶次) 本文方法 $ \alpha $ = 0.1$ \alpha $ = 0.2$ \alpha $ = 0.3$ \alpha $ = 0.4$ \alpha $ = 0.5$ \alpha $ = 0.6$ \alpha $ = 0.7$ \alpha $ = 0.8$ \alpha $ = 0.9总计时间 I 3.21 3.14 2.98 2.76 3.03 2.89 2.92 2.67 2.58 26.18 17.69 II 3.72 3.65 3.37 3.54 3.68 3.03 3.29 3.21 3.16 30.65 19.08 III 4.64 4.61 4.53 4.57 4.65 4.06 4.52 4.18 4.49 40.25 26.83
下载: 导出CSV
表 7 两种策略的时间对比(s)
Table 7 Time comparison of two strategies (s)
图像 不采用多分辨率 采用多分辨率 I 30.25 17.69 II 28.04 19.08 III 40.63 26.83
下载: 导出CSV
-
[1] 申艳平. 医学图像配准技术. 中国医学物理杂志, 2013, 30(1): 3385−3389Shen Yan-Ping. Review of image registration methods for medical images. Chinese Journal of Medical Physics, 2013, 30(1): 3385−3389 [2] Thirion J P. Image matching as a diffusion process: An analogy with Maxwell's demons. Medical Image Analysis, 1998, 2(3): 243−260 doi: 10.1016/S1361-8415(98)80022-4 [3] Wang H, Dong L, O′Daniel J, Mohan R, Garden A S, Ang K K, et al. Validation of an accelerated “demons” algorithm fordeformable image registration in radiation therapy. Physics in Medicine and Biology, 2005, 50(12): 2887−2905 doi: 10.1088/0031-9155/50/12/011 [4] Ying S H, Li D, Xiao B, Peng Y X, Do S Y, Xu M F. Nonlinear image registration with bidirectional metric and reciprocal regularization. PloS One, 2017, 12(2): e0172432 doi: 10.1371/journal.pone.0172432 [5] Du S Y, Guo Y R, Sanroma G, Ni D, Wu G R, Shen D G. Building dynamic population graph for accurate correspondence detection. Medical Image Analysis, 2015, 26(1): 256−267 doi: 10.1016/j.media.2015.10.001 [6] Du S Y, Zhang C J, Wu Z Z, Liu J, Xue J R. Robust isotropic sffigcaling ICP algorithm with bidirectional distance and bounded rotation angle. Neurocomputing, 2016, 215: 160−168 [7] Vercauteren T, Pennec X, Perchant A, Ayache N. Symmetric log-domain diffeomorphic registration: A demons-based approach. Med Image Comput Comput Assist Interv, 2008, 11(Pt1): 754−761 [8] 郝培博, 陈震, 江少锋, 汪洋. 基于Demons算法的多模态医学图像非刚性配准研究. 生物医学工程学杂志, 2014, 1: 161−165Hao P B, Chen Z, Jiang S F, Wang Y. Research on non-rigid registration of multi-modal medical image based on demons algorithm. Journal of Biomedical Engineering, 2014, 1: 161−165 [9] Lu X Q, Yu H F, Zhao Y, Hou H, Li Y H. Three-dimensional lung medical image registration based on improved demons algorithm. Optik–International Journal for Light and Electron Optics, 2015, 127(4): 1893−1899 [10] 闫德勤, 刘彩凤, 刘胜蓝, 刘德山. 大形变微分同胚图像配准快速算法. 自动化学报, 2015, 41(8): 1461−1470Yan De-Qin, Liu Cai-Feng, Liu Sheng-Lan, Liu De-Shan. A fast image registration algorithm for diffeomorphic image with large deformation. Acta Automatica Sinica, 2015, 41(8): 1461−1470 [11] 薛鹏, 杨佩, 曹祝楼, 贾大宇, 董恩清. 基于平衡系数的 Active demons非刚性配准算法. 自动化学报, 2016, 42(9): 1389−1400Xue Peng, Yang Pei, Cao Zhu-Lou, Jia Da-Yu, Dong En-Qing. Active demons non-rigid registration algorithm based on balance coefficient. Acta Automatica Sinica, 2016, 42(9): 1389−1400 [12] 路玉昆, 巩贯忠, 虞刚, 李登旺, 尹勇. 基于微分同胚Demons形变配准算法获取肺通气图. 中国医学物理学杂志, 2017, 34(7): 666−670Lu Y K, Gong G Z, Yu G, Li D W, Yin Y. Diffeomorphic demons registration algorithm for obtaining lung ventilation image. Chinese Journal of Medical Physics, 2017, 34(7): 666−670 [13] 张桂梅, 曹红洋, 陈阳泉, 刘建新. 基于分数阶梯度驱动的主动Demons算法研究. 电子学报, 2016, 44(12): 2834−2841 doi: 10.3969/j.issn.0372-2112.2016.12.004Zhang Gui-Mei, Cao Hong-Yang, Chen Yang-Quan, Liu Jian-Xin. Research on active demons based on fractional differentiation gradient driving. Acta Eletronica Sinica, 2016, 44(12): 2834−2841 doi: 10.3969/j.issn.0372-2112.2016.12.004 [14] Hu Y P, Modat M, Gibson E, Ghavami N, Bonmati E, Moore C M, Emberton M, Noble J A, Barratt D C, Vercauteren T. Label-driven weakly-supervised learning for multimodal deformable image registration. In: Proceedings of the 15th IEEE International Symposium on Biomedical Imaging, 2018. 1070–1074 [15] Balakrishnan G, Zhao A, Sabuncu M R, Guttag J, Dalca A V. An unsupervised learning model for deformable medical image registration. In: Proceedings of the 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition. Salt Lake City, UT, USA: IEEE, 2018. 235–242 [16] Pu Y F, Siarry P, Zhou J L, Zhang N. A fractional partial differential equation based multiscale denoising model for texture image. Mathematical Methods in the Applied Sciences, 2014, 37(12): 1784−1806 doi: 10.1002/mma.2935 [17] Chen D L, Chen Y Q, Xue D Y. Fractional-order total variation image denoising based on proximity algorithm. Applied Mathematics and Computation, 2015, 257: 537−545 [18] Mathieu B, Melchior P, Oustaloup A, Ceyral C. Fractional differentiation for edge detection. Signal Processing, 2003, 83(11): 2421−2432 doi: 10.1016/S0165-1684(03)00194-4 [19] Ren Z M. Adaptive active contour model driven by fractional order fitting energy. Signal Processing, 2015, 117: 138−150 [20] 张桂梅, 徐继元, 刘建新. 一种新的基于自适应分数阶的活动轮廓模型. 计算机研究与发展, 2017, 54(5): 1045−1056 doi: 10.7544/issn1000-1239.2017.20160301Zhang Gui-Mei, Xu Ji-Yuan, Liu Jian-Xin. A new active contour model based on adaptive fractional order. Journal of Computer Research and Development, 2017, 54(5): 1045−1056 doi: 10.7544/issn1000-1239.2017.20160301 [21] Yu J M, Tan L J, Zhou S B, Wang L P, Siddique M A. Image denoising algorithm based on entropy and adaptive fractional order calculus operator. IEEE Access, 2017, 5: 12275−12285 doi: 10.1109/ACCESS.2017.2718558 [22] 张桂梅, 郭黎娟, 熊邦书, 储珺. 基于多分辨率和自适应分数阶的Active Demons算法. 计算机研究与发展, 2018, 55(12): 2753−2763 doi: 10.7544/issn1000-1239.2018.20170523Zhang Gui-Mei, Guo Li-Juan, Xiong Bang-Shu, Chu Jun. Active demons algorithm based on multi-resolution and adaptive fractional differential. Journal of Computer Research and Development, 2018, 55(12): 2753−2763 doi: 10.7544/issn1000-1239.2018.20170523 [23] Du S Y, Xu G L, Zhang S R, Zhang X T, Gao Y, Chen B D. Robust rigid registration algorithm based on pointwise correspondence and correntropy. Pattern Recognition Letters, 2020, 132: 91−98 期刊类型引用(4)
1. 任会峰，伏建雄，鄢锋，李富. 基于旋转不变非线性局部模糊编码的膝骨关节炎辅助诊断. 软件导刊. 2022(04): 25-31 . 
百度学术2. 何宇，邓小龙，罗琦. 基于分数阶广义积分器的电网同步技术. 南京信息工程大学学报(自然科学版). 2022(03): 348-360 . 百度学术3. 杜晓刚，王玉琪，王福海，雷涛，张学军. 形变引导正则化的医学图像Demons快速配准算法. 小型微型计算机系统. 2022(12): 2580-2590 . 百度学术4. 郑伟，周杨，李文华，刘帅奇，张晓丹，马泽鹏. 基于变惯性系数active demons算法的DTI多通道配准研究. 河北大学学报(自然科学版). 2021(04): 436-442 . 百度学术其他类型引用(10)
-
下载:
计量
- 文章访问数: 1552
- HTML全文浏览量: 148
- PDF下载量: 198
- 被引次数: 14
