基于格子波尔兹曼模型的图像去噪方法

专利名称：基于格子波尔兹曼模型的图像去噪方法
技术领域：
本发明涉及一种基于格子波尔兹曼模型(lattice boltzmann model, LBM)来求解非线性扩散方程实现图像去噪的方法，属于图像处理领域。
背景技术：
目前，使用非线性扩散模型进行图像去噪是偏微分方程在图像处理领域中的重要应用。它能在抑制图像噪声的同时，有效保护图像的边缘。使用非线性扩散模型进行图像去噪首先由Perona和Malik于1990年提出，两年后Catte等在模型的理论和实现方法上做出了提高，1998年Weickert通过引入扩散张量改善了纹理图像的平滑效果。然而，由于数字图像固有的不连续性，数学模型所得到的偏微分方程具有非线性，以及图像数据量的庞大等因素，偏微分方程的解析解很难得到或者是根本不存在的，这时有必要借助于数值计算以获取该方程的近似解。在基于偏微分方程理论的图像处理方法中，数值实现问题是限制基于偏微分方程处理图像方法应用的瓶颈。许多非线性扩散方程都是通过显式有限差分(explicit finitedifference, EFD)进行离散实现。这种显式有限差分算法虽然易于实现，但是由于受稳定性的限制，迭代的步长是很小的(如1/4)。这样为了到达预期的扩散时间，就需要很多次迭代，整个计算过程效率低，使得它不适合应用于图像实时处理领域。
格子波尔兹曼模型架起了微观到宏观的桥梁。基于格子波尔兹曼模型的方法求解偏微分方程的出发点是根据系统的微观模型，设计格子波尔兹曼演化方程，在对系统进行模拟的同时实现对偏微分方程的数值求解。作为物理系统的一种数学模型，基于格子波尔兹曼模型的方法具有清晰而明确的物理解释。而作为一种数值解法，它以其清晰的物理思想，简单边界处理和快速的并行计算的优点，广泛的应用于流体力学、化学反应扩散、渗流、交通流等领域。格子波尔兹曼方法为实现处理图像非线性扩散去噪提供了高效的实现途径。

发明内容
本发明的目的在于针对现有技术存在的技术问题，提出一种基于格子波尔兹曼模型的图像去噪方法，该方法不仅能提高图像去噪质量，而且能提高计算效率，尤其适用于图像实时处理。
为了达到上述目的，本发明采用如下技术方案
本发明的基于格子波尔兹曼模型的图像去噪方法包括设立二维格子波尔兹曼模型，它由离散化的计算网格构成，每个网格的节点相当于一个元胞，它的值由粒子的分布函数
4/, (z' = 0,1，…，《)和扩散矢量c' (/ = 0,1,…^决定。
基于格子波尔兹曼模型的方法求解非线性扩散方程实现灰度图像去噪的演化方程为/, (X + C,-, w +1) — A (x，") = (0[/" (x， w) - /, (x， w)]
其中c,为扩散的矢量，/,Oc,w)为迭代次数为n时位于x处具有速度c,的粒子密度分布函
数，式中co为松弛因子，/,"为平衡态分布函数。
基于格子波尔兹曼模型的图像去噪方法的流程如图3所示，本发明提出的方法将图像的边缘截止函数嵌入格子波尔兹曼演化方程的松弛因子中，找到格子波尔兹曼演化方程与宏观方程对应关系来求解非线性扩散方程以实现图像去噪，其步骤如下
(1) 、输入初始图像/(x,0)，节点的值设为对应像素的灰度值；
(2) 、使用二维格子波尔兹曼模型，设置格子波尔兹曼演化方程中各作用方向的初始平衡
态函数/,，X,0);
(3) 、确定格子波尔兹曼演化方程的迭代次数^和迭代步长G
(4) 、遍历图像计算格子波尔兹曼演化方程中的松弛因子CO;
(5) 、计算二维格子波尔兹曼模型的迁移过程/,.(X + C,^) = /,(X,");
(6) 、计算二维格子波尔兹曼模型的作用过程
A (X + C;, W + 1) — A (X,") = C0[/， (X,") - /, (X,")];
(7) 、设"为迭代次数，根据二维格子波尔兹曼模型更新平衡分布函数为/, (X，W);
(8) 、判断是否达到迭代次数yV，若达到A^次时，则输出处理后的图像/(x，AO，若没有达到yV次时，则转步骤(4)，重复步骤(4广(7)，直到达到迭代yV次数后输出处理后的图像/(x，A0。
根据二维格子波尔兹曼模型中离散的扩散矢量Ci的方向数q的不同，可以把二维格子波尔兹曼模型分为D2Q5和D2Q9两类。
上述步骤(2)中的二维格子波尔兹曼模型为D2Q5模型，初始平衡态函数/,w(x,0)为
C=/,(X，0)D2e5=|/(x，0)上述步骤(4)中相关的松弛因子co为2
.+ SC:g其中g为边缘截止函数，C为步长，上述步骤(7)更新平衡分布函数/, (x, w)为
上述步骤(2)中二维格子波尔兹曼模型为D2Q9模型，初始平衡态函数/,"(x,0)为:
/"(x，0)， =/, (x，O)卿+(x,0)上述步骤(4)中相关的松弛因子ro为
2
1 + 3Cg
其中g为边缘截止函数，c为步长。上述步骤(7)更新平衡分布函数/，(x，n)为
r(x,")，4i/,(x,") "=in
基于D2Q5和D2Q9的格子波尔兹曼模型的图像去噪方法在处理图像时使用的边缘截止函数是相同的。上述的步骤(4)中，采用图像梯度的模值来估计边缘位置的方法以保护图像的边缘。在处理灰度图像时，边缘截止函数g为
,t_ t、 f *i3 i 5
其中C^是方差为0"的高斯核，K为平滑阈值。
在处理彩色图像时，综合利用三种颜色(R,G,B)分量图像的信息，在上述步骤(4)中，求取经过与高斯核巻积后图像的梯度系数矩阵的特征值之差作为边缘的估计值，矢量图像的
梯度的平方为即:
其中
<formula>formula see original document page 7</formula>
求系数矩阵A的两个特征值、和、，则矢量图像的边缘截止函数可取为:
<formula>formula see original document page 7</formula>
K为平滑阈值。
本发明与现有技术相比较，具有如下显而易见的突出实质性特点和显著优点上述基于格子波尔兹曼模型的图像去噪方法不仅能实现图像非线性扩散去噪，获得高质量的图像去噪效果，而且在保证算法稳定的情况下能进行大步长的运算，从而有效地提高计算的效率。尤其适用于图像实时处理。


图1是由离散化的计算网格构成的二维格子格子波尔兹曼模型结构示意图2-1是图1中基于D2Q5的格子波尔兹曼模型结构示意图2-2是图1中基于D2Q9的格子波尔兹曼模型结构示意图3是本发明的基于格子波尔兹曼模型的图像去噪方法的流程图4-1是灰度图像的原图4-2灰度图像加噪后的效果图4-3是基于D2Q5模型的格子波尔兹曼方法处理灰度图像去噪效果图；图4-4是基于D2Q9模型的格子波尔兹曼方法处理灰色图像去噪效果图；图5-1是彩色图像的原图5-2是彩色图像加噪后的效果图5-3是基于D2Q5模型的格子波尔兹曼方法处理彩色图像去噪效果图；图5-4是基于D2Q9模型的格子波尔兹曼方法处理彩色图像去噪效果图。
具体实施例方式
下面对本发明的基于格子波尔兹曼模型的图像去噪方法的实施例作详细说明本实施例以本发明的技术方案为前提下进行实施，给出了详细的实施方式，但本发明的保护范围不限于下述的实施例。
本发明的实施例，结合附图详细说明如下
如图l、 2所示，设立二维格子波尔兹曼模型，它由离散化的计算网格构成，每个网格的节点相当于一个元胞，它的值由粒子的分布函数/,力=0，1，-.，^和扩散矢量。(/ = 0,1, —《)决定。
节点值的每次更新(迭代)可分为两个阶段迁移阶段和作用阶段。迁移阶段由邻域节点向中心节点传输粒子，而作用阶段则决定传输的数量。运用格子波尔兹曼模型处理MX N的数字
图像时，每个像素和计算节点可以自然的对应起来，而像素的灰度值则对应节点上的粒子数。图1中节点(圆点)代表图像的每个像素，箭头显示了模型的迁移和作用方向。
根据离散的扩散矢量c,的方向数Q的不同，可以把二维格子波尔兹曼模型DnQq (n维q
方向)分为两类
D2Q5模型，如图2-l所示，它拥有5个离散速度
_ J (0,0), i' = 0
c' _ 1(cos(/ - l);r / 2, sin(/ - / 2)c， / = 1,2,3,4
D2Q9模型，如图2-2所示，它拥有9个离散速度
q = 0 / = 0
;r(/一l);r(/一l)
cos ~^-^, sm ~^-^
2 2
，「 ，7T(f-5)7T、 . ,TT(!'一5) 7T、C;=义，I COS( \ ; + —), sin(~^——^ + —)
4=1 z' = 1,2,3,4义=^/ = 5,6,7,8
基于格子波尔兹曼模型的图像去噪方法的流程如图3所示，本发明提出的方法将图像的边缘截止函数嵌入格子波尔兹曼演化方程的松弛因子中，找到格子波尔兹曼演化方程与宏观方程对应关系来求解非线性扩散方程以实现图像去噪。
8实施例l:灰度图像去噪处理
处理灰度图像的效果如图4所示，4-1、 4-2、 4-3、 4-4依次为原始Lena图像；添加均值为0，方差O.Ol的白噪声的加噪Lena图像；基于D2Q5格子波尔兹曼模型的去噪方法(参数为阈值25，步长5，迭代次数IO)处理后的图像；基于D2Q9格子波尔兹曼模型的去噪方法(参数为阈值25，步长5，迭代次数8)处理后的图像。以基于D2Q9格子波尔兹曼模型的去噪方法为例，其步骤如下
(1) 、输入初始图像/(x，0)，其灰度值设为节点上粒子数，
(2) 、根据模型的邻域方向数，设置各方向初始平衡态函数/r(x,0)，
C(x,0),9=/,(x，0)D2e9=*/(X,0)
(3) 、确定迭代次数为8和迭代步长为5，
(4) 、计算D2Q9模型演化方程中的松弛因子
"(,。*啡D卿=1 + 3"帆"D
其中边缘截止函数为
gOVG *l|>= 1— exp
(5) 、计算格子波尔兹曼模型的迁移过程/,.(x + c,，"^/Xx,")，
(6) 、计算格子波尔兹曼模型的作用过程/,(x + c,，" + l)-/,(x，") = co[/"(X,")-/,(X,")]，
(7) 、设/7为迭代次数，更新平衡分布函数为/r(x,")，
/，(x,")4土/,(x,")
(8) 、判断是否达到迭代次数8，若达到8次时，则输出处理后的图像/(x,8)，若没有达到8次时，则转步骤(4)，重复步骤(4) (7)，直到迭代次数达到8时，输出处理后的图像/(x,8)。
实施例2:彩色图像去噪处理
处理彩色图像的效果如图5所示，5-1、 5-2、 5-3、 5-4依次为p印per原图像；添加均值为0，方差0. 01的白噪声的加噪p印per图像；基于D2Q5格子波尔兹曼模型的去噪方法(参数为阈值IOO，步长5，迭代次数5)处理后的图像；基于D2Q9格子波尔兹曼模型的去噪方法(参数为阈值IOO，步长5，迭代次数4)处理后的图像。以基于D2Q5格子波尔兹曼
模型的去噪方法为例，其步骤如下
(1) 、设彩色图像的各分量图像为/,(x,0)(》1,2,3)
(2) 、分量图像的初始各方向的灰度为/力(X,0)，对分量图像设置各方向的平衡态函数
7(x,0)D2e5=/,,(X，0)D2e5=4/,(X,0)
(3) 、确定迭代次数5和迭代步长5
(4) 、遍历图像计算松弛因子
其中边缘截止函数为
i、 f 3 '3 iSgd — A2|) = 1 — exp|--^-^ 8
\ 、100 ;1
(5) 、计算格子波尔兹曼模型的迁移过程/;,(X + C,,") = /,,(X，")(6)、计算格子波尔兹曼模型的作用过程
/力(x + 、AW + A/) - /力(x， /)=①"-^ (x, /) - /力(x, 0]
(7)、设/7为迭代次数，分别更新各分量图像平衡分布函数为
/7(X，")4h,(X,")
5=0
(8)、判断是否达到迭代次数5，若达到5次时，则输出处理后的图像/7(、5)，若没有达
到5次时，则转步骤(4)，重复步骤(4) (7)，直到迭代次数达到5时，输出处理后的图像/)(x,5)
实验显示，相对以加性算子分裂算法(additive operator splitting, AOS)为代表的有限差分法，基于格子波尔兹曼模型的图像去噪方法可以获得更好的去噪质量，其计算效率也优于A0S算法。采用峰值信噪比(PSNR)客观衡量去噪质量。处理后图像的PSNR越大，说明其越接近理想去噪结果，对应的去噪方法越有效。仿真运行时使用联想旭日410A笔记本(CPUT2080, Memory 1G)，运行环境为MATLAB 2008b。表1为使用本发明的基于D2Q5、 D2Q9格子波尔兹曼模型的方法和AOS方法处理添加了均值为0，方差分别为0. 01和0. 1的高斯白
10噪声的Lena图像的结果。在边缘截止函数中高斯巻积核的方差o统一取为1，阈值为25。选
取PSNR最高的那一次迭代结果，可以看出基于格子波尔兹曼模型的方法比AOS能更加有效的抑制噪声。这里格子波尔兹曼方法的计算时间为串行计算时间，由于格子波尔兹曼模型各个方向的作用和迁移过程是独立的，可以进一步使用并行算法降低计算时间。
去噪方噪声方加噪图像步长迭代次数处理后图像计算时间(s)
法差(db)(db)
D2Q50,0119.98 11028. 050. 5085
D2Q50. 117.03 11019.410. 5224
D2Q90.0119.98 1827.840. 5697
D2Q90. 117.03 1819. 420. 5761
AOS0.0119.98 1326. 860. 5308
AOS0. 117.03 1319.260.5316
表权利要求
1、一种基于格子波尔兹曼模型的图像去噪方法，其特征在于，通过把图像的边缘截止函数嵌入格子波尔兹曼微观演化方程的松弛因子中的方式，在二维格子波尔兹曼模型中找到微观演化方程与宏观非线性扩散方程对应关系来求解非线性扩散方程以实现图像去噪，其步骤如下(1)、输入初始图像I(x，0)，节点的值设为对应像素的灰度值；(2)、使用二维格子波尔兹曼模型，设置格子波尔兹曼微观演化方程中各作用方向的初始平衡态函数Iieq(x，0)；(3)、确定格子波尔兹曼微观演化方程的迭代次数N和迭代步长C；(4)、遍历图像计算格子波尔兹曼演化方程中的松弛因子ω；(5)、计算二维格子波尔兹曼模型的迁移过程Ii(x+ci，n)＝Ii(x，n)；(6)、计算二维格子波尔兹曼模型的作用过程<maths id="math0001" num="0001" ><math><![CDATA[ <mrow><msub> <mi>I</mi> <mi>i</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>+</mo> <msub><mi>c</mi><mi>i</mi> </msub> <mo>,</mo> <mi>n</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo></mrow><mo>-</mo><msub> <mi>I</mi> <mi>i</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo></mrow><mo>=</mo><mi>&omega;</mi><mo>[</mo><msubsup> <mi>I</mi> <mi>i</mi> <mi>eq</mi></msubsup><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo></mrow><mo>-</mo><msub> <mi>I</mi> <mi>i</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo></mrow><mo>]</mo><mo>;</mo> </mrow>]]></math></maths>(7)、设n为迭代次数，根据二维格子波尔兹曼模型更新平衡分布函数为Iieq(x，n)；(8)、判断是否达到迭代次数N，若达到N次时，，则输出处理后的图像I(x，N)，若没有达到N次时，则转步骤(4)，重复步骤(4)~(7)，直到达到迭代N次数后输出处理后的图像I(x，N)。
2、 根据权利要求l所述的基于格子波尔兹曼模型的图像去噪方法，其特征在于上述步 骤(2)中二维格子波尔兹曼模型为D2Q5模型时，初始平衡态函数/"(x，0)为<formula>formula see original document page 2</formula>上述步骤(4)中相关的松弛因子co为<formula>formula see original document page 2</formula>其中g为边缘截止函数，C为步长，上述步骤(7)更新平衡分布函数/"(X，")为<formula>formula see original document page 2</formula>
3、 根据权利要求l所述的基于格子波尔兹曼模型的图像去噪方法，其特征在于上述步骤(2) 二维格子波尔兹曼模型为D2Q9模型时，初始平衡态函数/"(x，0)为: C(x,O)， =/, (x,O), =|/(x，0)上述步骤(4)中相关的松弛因子(D为<formula>formula see original document page 3</formula>其中g为边缘截止函数，C为步长，上述步骤(7)更新平衡分布函数/r(X，")为<formula>formula see original document page 3</formula>
4、根据权利要求1或2或3所述的基于格子波尔兹曼模型的图像去噪方法，其特征在于: 在处理灰度图像时，所述步骤(4)的松弛因子co中边缘截止函数g为<formula>formula see original document page 3</formula>其中C^是方差为a的高斯核，K为平滑阈值。
5、根据权利要求1或2或3所述的基于格子波尔兹曼模型的图像去噪方法，其特征在于: 在处理彩色图像时，所述步骤(4)的松弛因子w中边缘截止函数g为<formula>formula see original document page 3</formula>其中K为平滑阈值，^和^为经过方差o"的高斯核G。与彩色图像巻积后的梯度系数矩阵为<formula>formula see original document page 3</formula>的特征值。
全文摘要
本发明公开了一种基于格子波尔兹曼模型的图像去噪方法，该方法步骤为(1)输入初始图像I(x，0)；(2)设置二维格子波尔兹曼演化方程中各作用方向的初始平衡态函I_i^eq(x，0)；(3)确定格子波尔兹曼演化方程的迭代次数N和迭代步长C；(4)遍历图像计算格子波尔兹曼演化方程中的松弛因子ω；(5)计算格子波尔兹曼模型的迁移过程；(6)计算格子波尔兹曼模型的作用过程；(7)更新平衡分布函数为I_i^eq(x，n)；(8)判断是否达到迭代N次数，若达到N次时，则输出处理后的图像I(x，N)。该方法能在抑制图像噪声的同时，有效的保护图像边缘。不仅能提高图像去噪质量，并且可以实现大步长的迭代计算，从而有效提高了去噪处理的效率。
文档编号G06T5/00GK101673393SQ20091019651
公开日2010年3月17日 申请日期2009年9月25日 优先权日2009年9月25日
发明者严壮志, 王志强, 钱跃竑 申请人:上海大学