基于Demons算法的图像去噪算法

摘要

本发明涉及一种基于Demons算法的图像去噪算法，该方法首先在Demons算法的基础上，将扩散过程看作图像配准，建立一种新的基于图像配准的Demons去噪算法，该算法去噪性能较经典的PM算法更优越；其次，考虑到在图像去噪过程中，仅靠梯度信息表征图像的局部特征是不够的，而二阶微分量含有更丰富的信息，将水平集曲率作为控制图像结构的驱动力因素引入到图像配准去噪算法中，建立了一种梯度与曲率双重驱动力的图像去噪算法，双驱动算法；最后，采用加性算子分裂算法(AOS算法)对上述算法进行处理，得到去噪后的图像。本发明不但去噪性能优越，而且图像的整体结构保持完好，去噪后的图像信噪比较其他Demons算法提高了15dB左右，较PM算法提高了25dB左右，清晰度也大幅度提升。

说明书

技术领域

本发明涉及图像处理领域，是基于偏微分方程的图像去噪算法，该算法基于Demons算 法改进。

背景技术

数字图像是很多学科领域获取信息的来源，但图像在采集过程中往往会因为各方面原因 引入噪声。因此，在图像处理和计算机领域，图像去噪是最基本的问题之一。近几十年，偏 微分方程方法开始大量应用于图像处理，在图像的去噪、分割、边缘检测、增强等方面的研 究都取得了显著进展。

目前为止，研究人员已经提出了很多各向异性扩散算法，其中最经典的是令传导系数依 赖于图像梯度的二阶偏微分算法，即PM算法。随着对该技术研究的不断深入,很多实验结果 表明,PM算法存在缺陷:第一,它有时无法正确地区分边缘和噪声,因而处理小尺度空间的噪声 效果不理想；第二,没有理论来验证扩散系数函数的正确性。因此,许多学者对PM算法中扩散 系数进行了改进,试图建立起更有效的保边缘平滑滤波器,并取得了一些进展。然而大部分改 进的算法都是建立在经典算法之上的，没有从根本上解决PM算法的缺陷。

发明内容

本发明旨在处理出更加清晰的图像，以接近原始图像。由于传统算法时效性低，复杂度 较高，且在处理过程中经常会出现过度平滑和平滑不充分的现象，本算法在Demons算法基 础上进行改进，使图像各向异性去噪过程简化成图像配准过程，不需要考虑设计扩散函数、 梯度阈值以及区分边缘等问题，降低了算法复杂度和提高了算法时效性，因此简化了图像去 噪问题，高效的利用了整体结构函数控制图像结构。去噪性能方面，峰值信噪比和清晰度均 大幅提高。

本发明采用以下技术方案：一种基于Demons算法的图像去噪算法,其特征在于，包括以 下步骤：

步骤一、将噪声图像进行Gauss滤波，去除较大噪声；

步骤二、利用Demons算法原理，建立基于Demons算法的图像去噪算法

其中div、▽分别为散度算子和梯度算子，I I₀分别是初始 图像和加噪图像，|▽I|为梯度模值；

步骤三、将水平集曲率作为校正驱动因子引入到Demons图像去 噪模型中，建立结构函数最终建立双驱动去噪算法为 $\frac{\partial I}{\partial t}=\alpha ·\mathrm{div}\left(\frac{\left|{I-I}_0\right|▿I}{{\left|\mathrm{\Delta I}\right|}^2+{\left({I-I}_0\right)}^2}\right)+\beta ·\mathrm{div}\left(\frac{\left|{I-I}_0\right|▿f\left(k\right)}{{|▿f\left(k\right)|}^2+{\left({I-I}_0\right)}^2}\right),$其中|▽I|、|▽f|为梯度模值，α 和β是连贯系数，通过曲线拟合来确定；

步骤四、用半隐式加性算子分裂数值算法对步骤三中双驱动去噪算法进行离散化处理，得 到去噪后图像，过程如下：

经化简，

$\frac{\partial I}{\partial t}=\mathrm{div}(g\left(\right|▿I\left|\right)▿I)=\frac{\partial g}{\partial x}\frac{\partial I}{\partial x}+g\frac{{\partial }^{2}I}{{\partial }^{2}x}+\frac{\partial g}{\partial y}\frac{\partial I}{\partial y}+g\frac{{\partial }^{2}I}{{\partial }^{2}y}$

$\frac{\partial f}{\partial t}=\mathrm{div}(G\left(\right|▿f\left|\right)▿f)=\frac{\partial G}{\partial x}\frac{\partial f}{\partial x}+G\frac{{\partial }^{2}f}{{\partial }^{2}x}+\frac{\partial G}{\partial y}\frac{\partial f}{\partial y}+G\frac{{\partial }^{2}f}{{\partial }^{2}y}$

用表示I(x_i,y_j；t_n),f(x_i,y_j；t_n)，其隐式格式为

$I^{n+1}{=I}^n+\mathrm{\Delta t}\underset{l\in \{x,y\}}{\Sigma }{A}_l\left(I^n\right)I^{n+1}$

$f^{n+1}{=f}^n+\mathrm{\Delta t}\underset{l\in \{x,y\}}{\Sigma }{A}_l^{\prime }\left(f^n\right)f^{n+1}$

其中Δt为时间步长，A_x、A_y分别为对x和y方向进行差分对应的系数矩阵，

求解Iⁿ⁺¹，fⁿ⁺¹：

$\left(\begin{array}{c}{I}^{n+1}=I^n+2\mathrm{\Delta t}{A}_x\left(I^n\right)I^{n+1}\\ I^{n+1}=I^n+{2\mathrm{\Delta tA}}_y\left(I^n\right)I^{n+1}\end{array}\right)$

$\left(\begin{array}{c}{f}^{n+1}=f^n+2\mathrm{\Delta t}{A}_x\left(f^n\right)f^{n+1}\\ f^{n+1}=f^n+{2\mathrm{\Delta tA}}_y\left(f^n\right)f^{n+1}\end{array}\right)$

变形可得

$\left(\begin{array}{c}{I}^{n+1}={(I-2\Delta {\mathrm{tA}}_x\left(I^n\right))}^{-1}{I}^n\\ I^{n+1}={(I-{2\mathrm{\Delta tA}}_y\left(I^n\right))}^{-1}{I}^n\end{array}\right)$

$\left(\begin{array}{c}{f}^{n+1}={(f-2\Delta {\mathrm{tA}}_x\left(f^n\right))}^{-1}{f}^n\\ I^{n+1}={(f-{2\mathrm{\Delta tA}}_y\left(f^n\right))}^{-1}{f}^n\end{array}\right)$

1)当i＝1，…，N时，计算和三个对角线上的元素：$({\alpha }_k^{\left(i\right)},k=1,···,N),$$({\beta }_k^{\left(i\right)},k=1,···,N-1),({\gamma }_k^{\left(i\right)},k=2,···,N)$和$({\alpha }_{1k}^{\left(i\right)},k=1,···,N),({\beta }_{1k}^{\left(i\right)},k=1,···,N-1),({\gamma }_{1k}^{\left(i\right)},k=2,···,N)$并采用追赶法求解$\left({I-2\mathrm{\tau A}}_{x,i}^n\right)I_{1i}^{n+1}=I_{1i}^n,\left({f-2\mathrm{\tau A}}_{x,i}^n\right)f_{1i}^{n+1}=f_{1i}^n,$得到

2)当j＝1,…,M时，同样计算和的三个对角线上的元素，求解 $\left({I-2\mathrm{\tau A}}_{y,j}^n\right)I_{2j}^{n+1}=I_{2j}^n,\left({f-2\mathrm{\tau A}}_{y,j}^n\right)f_{2j}^{n+1}=f_{2j}^n,$得到

3)计算$I^{n+1}+f^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1})+\frac{1}{2}(f_1^{n+1}+f_2^{n+1});$

4)重复1)～3)，经过多次迭代，得到清晰图像。

上述Demons算法即其中F是参考图像，M是浮动图像， 为(x,y)处从浮动图像向参考图像形变所需的位移量，div、▽分别为散度算子和梯 度算子。

本发明达到的有益效果：①在算法的复杂度方面，需要的信息量少，方法简单，不需要着 重分析图像的局部细节信息，只需从图像的整体结构考虑，建立整体结构函数来控制图像结 构，实现图像去噪，使受污染图像更接近原始图像；②在算法的时效性方面，因为本发明需 要的信息量少，实施的复杂度低，从而降低了算法的处理时间；③在算法的去噪性能方面， 图像的峰值信噪比和清晰度均大幅提高，受噪声污染的图像经本算法处理后更加接近原始图 像。

附图说明

图1为建立基于Demons算法的图像去噪算法的流程图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发 明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用 于限定本发明。

如图1所示，将噪声图像进行预处理，即进行Gauss滤波，去除较大噪声，接下来判断 图像结构是否变形，即考虑含噪图像的梯度|▽I|做边缘检测时，由于受噪声干扰，使得对图 像边缘、纹理的检测不够准确，在去噪过程中破坏图像本身的边缘纹理等细节信息，而使图 像的结构发生变形。

(1)在图像结构没有变形的情况下，接着检测图像的所有像素是否直通，如果所有像素 直通，则得到并输出清晰的图像；如果不是所有像素直通，则回到噪声图像预处理阶段进行 处理。

(2)在图像结构有变形的情况下，由于图像的水平集曲率在有噪声的情况下，能够有效 检测图像的边缘纹理等细节信息，故建立水平集曲率的结构控制函数，用来维持图像的结构， 使之不发生变形。将水平集曲率作为校正驱动因子引入到Demons 图像去噪模型$\frac{\partial I}{\partial t}=\mathrm{div}\left(\frac{\left|{I-I}_0\right|▿I}{{|▿I|}^2+{\left({I-I}_0\right)}^2}\right)$中，建立结构函数$f\left(\kappa \right)=1-\exp (-\frac{{\kappa }^2}{100}),$最终建立双 驱动去噪算法为$\frac{\partial I}{\partial t}=\alpha ·\mathrm{div}\left(\frac{\left|{I-I}_0\right|▿I}{{|▿I|}^2+{\left({I-I}_0\right)}^2}\right)+\beta ·\mathrm{div}\left(\frac{\left|{I-I}_0\right|▿f\left(k\right)}{{|▿f\left(k\right)|}^2+{\left({I-I}_0\right)}^2}\right).$

接下来用半隐式加性算子分裂数值算法(AOS算法)对双驱动去噪算法进行离散化处理。 经化简，得到$\frac{\partial I}{\partial t}=\mathrm{div}(g\left(\right|▿I\left|\right)▿I)=\frac{\partial g}{\partial x}\frac{\partial I}{\partial x}+g\frac{{\partial }^{2}I}{{\partial }^{2}x}+\frac{\partial g}{\partial y}\frac{\partial I}{\partial y}+g\frac{{\partial }^{2}I}{{\partial }^{2}y}$和

$\frac{\partial f}{\partial t}=\mathrm{div}(G\left(\right|▿f\left|\right)▿f)=\frac{\partial G}{\partial x}\frac{\partial \mathrm{fI}}{\partial x}+G\frac{{\partial }^{2}f}{{\partial }^{2}x}+\frac{\partial G}{\partial y}\frac{\partial f}{\partial y}+G\frac{{\partial }^{2}f}{{\partial }^{2}y},$用表示I(x_i,y_j；t_n),f(x_i,y_j；t_n)，其 隐式格式分别为$I^{n+1}=I^n+\mathrm{\Delta t}\underset{l\in \{x,y\}}{\Sigma }{A}_l\left(I^n\right)I^{n+1},f^{n+1}=f^n+\mathrm{\Delta t}\underset{l\in \{x,y\}}{\Sigma }{A}_l^{\prime }\left(f^n\right)f^{n+1},$其中Δt为时间步长，A_x、 A_y分别为对x和y方向进行差分对应的系数矩阵。

求解Iⁿ⁺¹，fⁿ⁺¹：

$\left(\begin{array}{c}{I}^{n+1}=I^n+2\mathrm{\Delta t}{A}_x\left(I^n\right)I^{n+1}\\ I^{n+1}=I^n+{2\mathrm{\Delta tA}}_y\left(I^n\right)I^{n+1}\end{array}\right)$

$\left(\begin{array}{c}{f}^{n+1}=f^n+2\mathrm{\Delta t}{A}_x\left(f^n\right)f^{n+1}\\ f^{n+1}=f^n+{2\mathrm{\Delta tA}}_y\left(f^n\right)f^{n+1}\end{array}\right)$

变形可得

$\left(\begin{array}{c}{I}^{n+1}={(I-2\Delta {\mathrm{tA}}_x\left(I^n\right))}^{-1}{I}^n\\ I^{n+1}={(I-{2\mathrm{\Delta tA}}_y\left(I^n\right))}^{-1}{I}^n\end{array}\right)$

$\left(\begin{array}{c}{f}^{n+1}={(f-2\Delta {\mathrm{tA}}_x\left(f^n\right))}^{-1}{f}^n\\ I^{n+1}={(f-{2\mathrm{\Delta tA}}_y\left(f^n\right))}^{-1}{f}^n\end{array}\right)$

1)当i＝1，…，N时，计算和三个对角线上的元素：$({\alpha }_k^{\left(i\right)},k=1,···,N),$$({\beta }_k^{\left(i\right)},k=1,···,N-1),({\gamma }_k^{\left(i\right)},k=2,···,N)$和$({\alpha }_{1k}^{\left(i\right)},k=1,···,N),({\beta }_{1k}^{\left(i\right)},k=1,···,N-1),({\gamma }_{1k}^{\left(i\right)},k=2,···,N)$求解$\left({I-2\mathrm{\tau A}}_{x,i}^n\right)I_{1i}^{n+1}=I_{1i}^n,\left({f-2\mathrm{\tau A}}_{x,i}^n\right)f_{1i}^{n+1}=f_{1i}^n,$得到

2)当j＝1,…,M时，同样计算和的三个对角线上的元素，求解 $\left({I-2\mathrm{\tau A}}_{y,j}^n\right)I_{2j}^{n+1}=I_{2j}^n,\left({f-2\mathrm{\tau A}}_{y,j}^n\right)f_{2j}^{n+1}=f_{2j}^n,$得到

3)计算$I^{n+1}+f^{n+1}=\frac{1}{2}(I_1^{n+1}+I_2^{n+1})+\frac{1}{2}(f_1^{n+1}+f_2^{n+1});$

4)重复1)～3)，经过多次迭代，直至得到清晰图像。

最后检测图像中所有像素是否直通，是则输出清晰图像，否则重回到噪声图像预处理步 骤。

本发明不但去噪性能优越，而且图像的整体结构保持完好，去噪后的图像信噪比较 Demons去噪其他算法提高了15dB左右，较PM算法提高了25dB左右，清晰度也大幅度提升。

以上是本发明的较佳实施方式，但本发明的保护范围不限于此。任何熟悉本领域的技术 人员在本发明所揭露的技术范围内，未经创造性劳动想到的变换或替换，都应涵盖在本发明 的保护范围之内。因此本发明的保护范围应以权利要求所限定的保护范围为准。