CN106296649A - 一种基于水平集模型的纹理图像分割方法 - Google Patents

Abstract

一种基于水平集模型的纹理图像分割方法，包括以下步骤：一、对输入图像进行Log‑Gabor滤波，输出滤波后的图像；二、对经过滤波后的图像提取LSS描述符，得到由描述符分量构成的图像；三、在由描述符分量构成的图像上构造纹理能量项；四、在输入图像上构造LGDF能量项；五、综合三和四中的两种能量项，组成本方法能量函数；六、通过最小化五中的能量函数，求解能量函数的分割曲线，并使用格子波尔曼法加速求解过程。本发明基于水平集模型，利用LSS描述符设计了一个新的纹理能量项，与传统LGDF能量项结合后对各种纹理图像都有良好的分割效果，对噪声和灰度不均具有鲁棒性，并通过格子玻尔兹曼法优化提升了图像分割效率。

Description

一种基于水平集模型的纹理图像分割方法

技术领域

本发明涉及一种基于水平集模型的纹理图像分割方法，属于图像分割以及图像处理技术领域。

背景技术

尽管图像分割是图像处理领域最基本的问题之一，仅处于D.Marr所提出的视觉计算理论中图像处理与图像分析层次之间，但是图像分割至今仍为一项极为困难的工作。主要原因其一是因为图像分割的“病态”特性，其分割结果的正确性并不确定，具有二义性，与人的主观感受和心理作用有密切关系；其二是因为图像中的信息十分丰富，常常同时拥有灰度、色彩、纹理等多种成分，难以通过一个统一的模型进行概括和衡量。所以科研人员对于图像分割的研究方式主要是在一个大的框架下，对具体的图像类型进行分析，提出有针对性的分割方法。在这些图像分割方法中，基于水平集模型(Level set model)的方法能够将图像的相关属性和先验知识融入统一的能量函数构造过程，具有良好的可扩展性，是图像分割领域的研究热点之一。

水平集模型是一种引入了水平集函数，将分割曲线表示为高一维的标量函数的模型。水平集模型能够使得分割曲线自适应地改变拓补结构，并且它是一种无参数(Parameterization free)模型，分割曲线在演化过程中不受其自身参数表示的影响，具有数值稳定性。纹理特征是基于水平集模型的纹理图像分割方法中最重要的方面。对于纹理图像，需要分析其纹理特征，将纹理特征融入能量函数构造过程中，从而得到有针对性的分割方法。局部自相似(Local self-similarity，LSS)描述符是一种常用于数值化纹理特征的描述符。它具有独立于灰度和色彩的理想性质，因此可在能量函数中引入LSS纹理能量项，提高纹理图像的分割效果。

与本发明相关的文章和专利有以下4篇，下文将分别对其进行分析：

文章(1)：《IEEE Transactions on Image Processing》2001年10卷第2期，题目为：“Active contours without edges”，提出了一种最经典的基于水平集的图像分割方法——CV(Chan-Vese)算法，但是CV算法对灰度不均和局部有噪声图像无法做出准确的分割。

文章(2)：《Signal Processing》2009年89卷第12期，题目为：“Active contoursdriven by local Gaussian distribution fitting energy”，通过最大后验概率(Maximum a posteriori probability，MAP)思想提出了一种基于水平集的图像分割方法——LGDF(Local Gaussian Distribution Fitting)算法，LGDF对灰度不均和局部有噪声图像能够有较好的处理，但是LGDF不能有效地分割纹理图像。

文章(3)：《Pattern Recognition》2015年48卷第4期，题目为：An Intensity-Texture model based level set method for image segmentation，通过ASLVD(Adaptive Scale Local Variation Degree)滤波提取纹理特征，并将其融入能量函数构造中，提出了一种基于水平集模型的纹理图像分割方法，但是该方法分割结果有一定瑕疵，未实现曲线演化加速功能。

专利(1)：申请号2013103713364，标题为“基于两通道纹理分割主动轮廓模型的图像分割方法”，利用图像的灰度和梯度构造两个纹理特征通道，分别在两个通道上应用水平集模型进行纹理图像分割，但是该方法更容易受到图像灰度和梯度不均匀的影响，并且增加了通道数，计算复杂度较高。

上述已有的基于水平集模型的图像分割方法虽然在纹理图像分割上取得了一定成果，然而，由于受限于纹理特征表示方法，分割结果并不理想，并且在复杂图像上缺乏曲线演化加速过程，演化所需的迭代数和时间较长。本发明旨在解决现有方法对纹理图像分割效果不理想的问题，同时对曲线演化进行加速，提出了一种对灰度不均和噪声具有鲁棒性的基于水平集模型的纹理图像分割方法。

发明内容

本发明的目的是为了解决现有方法对纹理图像分割困难以及对灰度不均和噪声敏感的缺陷，提出了一种基于水平集模型的纹理图像分割方法。

本发明的核心技术思想是通过计算图像的LSS描述符，构建图像的纹理能量项，将其融入水平集模型能量函数，然后再利用曲线演化求解能量函数，并且在该过程中应用格子波尔曼方法(Lattice Boltzmann Method，LBM)加速演化，减少演化所需的迭代次数，提高效率；

一种基于水平集模型的纹理图像分割方法，简称本方法；

具体包括以下步骤：

步骤1，对输入的一幅图像进行Log-Gabor滤波，以增强图像的纹理特征，输出Log-Gabor滤波后的图像；

其中，记输入图像为I，Log-Gabor滤波后输出的图像为I_G。Log-Gabor滤波需在图像的频率域上进行卷积操作，则输出I_G与输入I的关系为：

表示傅里叶变换，*为卷积符号，G为Log-Gabor的传递函数，即：

G(w；σ_f,σ_θ,θ,S)＝exp(-log(w/w₀)/2(log(k/w₀))²) (2)

其中，exp表示以常数e为底的指数函数，log表示自然对数；w为图像I的像素点在频率域上对应的值，w₀为在频率域上滤波中心的像素点的值，k为控制传递函数形状的参数，通常取k∈{k|0.55w₀≤k≤0.85w₀}，以便使得滤波频率参数σ_f＝k/w₀∈[0.55,0.85]，σ_θ为角度参数，通常取0.41～0.74，θ为方向参数，取0～360，S为小波尺度，通常为1～3；

步骤2，对经过步骤1Log-Gabor滤波后的图像提取LSS描述符，得到由描述符分量构成的图像；

步骤2.1，对于I_G的每一个像素p，以p为中心选取N×N像素的区域(N通常取21～61的奇数)，若图像边界处区域不足N×N大小，则复制图像边界处的值来进行扩展，使得区域大小为N×N；例如，假设图像左上角边界处如式(3)所示，扩展2像素，则扩展后的图像为式(4)：

在每个N×N区域中，以p为中心选取n×n像素的小块(n通常取3～9的奇数)，记为t_p；对于N×N区域的每一个像素x(可以为p也可以非p)，以x为中心选取n×n像素的小块，记为t_x，若在区域边界处小块不足n×n大小，同样复制区域边界处的值来进行扩展；

计算t_p和t_x的平方差(Sum of Square Differences，SSD)，即：

$SSD(t_p,t_x)=\left|\right|t_p-t_x|{|}_F^2---\left(5\right)$

||·||_F为F-范数；

步骤2.2，根据步骤2.1输出的SSD计算N×N区域上的相关矩阵C_p(C_p的大小为N×N)：

$C_p\left(x\right)=\exp (-\frac{SSD(t_p,t_x)}{\mathrm{\α}})---\left(6\right)$

α为自适应参数，可取x所在的N×N区域中最大的SSD值；

步骤2.3，将C_p从直角坐标转换到对数极(Log-polar)坐标上；

步骤2.4，按照Log-polar坐标(ρ,θ)将C_p划分为M块，再对每个块取最大值，将这些最大值拼接起来即LSS描述符L_p：

$L_p(\mathrm{\ρ},\mathrm{\θ})=\underset{x\∈BIN(\mathrm{\ρ},\mathrm{\θ})}{max}\left\{C_p\left(x\right)\right\}---\left(7\right)$

其中，每个块记为BIN(ρ,θ)，是L_p的第i个描述符分量，将按照像素p的空间域顺序构成则就是步骤2输出的由LSS描述符分量构成的图像，i的取值范围是i∈{1,2,...,M}；

步骤3，在步骤2输出的由LSS描述符分量构成的图像上构造纹理能量项，具体为：

记图像的每个像素点为x，x所具有的纹理能量为：

$\begin{array}{c}{E}_x^{LSS}={\mathrm{\λ}}_1\mathrm{\ω}(x-y)\∫\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-in}^i\left(x\right)|H\left(\mathrm{\φ}\right(y\left)\right)dy\\ +{\mathrm{\λ}}_2\mathrm{\ω}(x-y)\∫\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-out}^i\left(x\right)|(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy\end{array}---\left(8\right)$

y为x邻域内的像素，φ为水平集函数，λ₁和λ₂分别为曲线内部平衡参数和曲线外部平衡参数，通常取0.5～2.0，H为Heaviside函数，即：

$H\left(\mathrm{\φ}\right)=\frac{1}{2}\left(1+\frac{2}{\mathrm{\π}}arctan\left(\frac{\mathrm{\φ}}{\mathrm{\ϵ}}\right)\right)---\left(9\right)$

ε为Heavisde函数参数，通常取0.1～10；ω(d)为权重函数，其中d＝x-y：

$\mathrm{\ω}\left(d\right)=\left\{\begin{array}{cc}\frac{1}{a}\exp \left(-\frac{|d{|}^2}{2{\mathrm{\σ}}^2}\right)& if\left|d\right|\≤\mathrm{\ρ}\\ 0& if\left|d\right|>\mathrm{\ρ}\end{array}\right.---\left(10\right)$

a和σ为权重函数参数，ρ为距离阈值；和分别为演化曲线内部和外部的均值：

$L_{p-in}^i\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)L_p^i\left(y\right)H\left(\mathrm{\φ}\left(y\right)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\left(y\right)\right)dy}---\left(11\right)$

$L_{p-out}^i\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)L_p^i\left(y\right)\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)dy}{\∫\mathrm{\ω}(y-x)\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)dy}---\left(12\right)$

则纹理能量项E^LSS如式(13)：

$E^{LSS}={\∫}_{\mathrm{\Ω}}{E}_x^{LSS}dx---\left(13\right)$

步骤4，在步骤1中的输入图像上构造LGDF能量项；

其中，LGDF能量项记为E^LGDF，通过式(14)计算：

$\begin{array}{c}{E}^{LGDF}=-{\mathrm{\λ}}_1\∫\∫\mathrm{\ω}(x-y)\ln p_{1,x}\left(I\right(y\left)H\right(\mathrm{\φ}\left(y\right)\left)\right)dydx\\ -{\mathrm{\λ}}_2\∫\∫\mathrm{\ω}(x-y)\ln p_{2,x}\left(I\left(y\right)\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)\right)dydx\end{array}---\left(14\right)$

p_j,x(I(y))(其中j∈{1,2})为概率分布：

$p_{j,x}\left(I\right(y\left)\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j\left(x\right)}\exp \left(-\frac{{\left(u_j\right(x)-I(y\left)\right)}^2}{2{\mathrm{\σ}}_j{\left(x\right)}^2}\right)---\left(15\right)$

其中u₁、u₂、σ₁(x)²和σ₂(x)²分别为：

$u_1\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)I\left(y\right)H\left(\mathrm{\φ}\right(y\left)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\right(y\left)\right)dy}---\left(16\right)$

$u_2\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)I\left(y\right)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}{\∫\mathrm{\ω}(y-x)\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)dy}---\left(17\right)$

${\mathrm{\σ}}_1{\left(x\right)}^2=\frac{\∫\mathrm{\ω}(y-x){\left(u_1\right(x)-I(y\left)\right)}^{2}H\left(\mathrm{\φ}\right(y\left)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\left(y\right)\right)dy}---\left(18\right)$

${\mathrm{\σ}}_2{\left(x\right)}^2=\frac{\∫\mathrm{\ω}(y-x){\left(u_2\left(x\right)-I\left(y\right)\right)}^2\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)dy}{\∫\mathrm{\ω}(y-x)(1-H\left(\mathrm{\φ}\right(y\left)\right))dy}---\left(19\right)$

步骤5，综合步骤4输出的LGDF能量项和步骤3输出的LSS纹理能量项，组成本方法的能量函数；

其中，记本方法的能量函数为E^LL，则：

E^LL＝ξ₁E^LSS+ξ₂E^LGDF+ν·Length(C)+μ·P(x) (20)

ξ₁和ξ₂分别为LSS权重参数和LGDF权重参数，通常取1～40，Length(C)＝∫|▽H(φ(x))|dx为曲线长度，为防止重新初始化项，避免了曲线演化过程中的重新初始化过程，ν和μ分别为曲线长度和防止重新初始化项的系数，ν通常取0.0001*255²～0.01*255²，μ通常取0.1～1.5；

步骤6，通过最小化步骤5所得的能量函数求解能量函数的分割曲线；同时，在此过程中使用LBM加速曲线演化，减少迭代次数；

步骤6.1，通过步骤5所得能量函数获取梯度下降流公式。其中，式(20)能量函数对应的梯度下降流公式为：

$\begin{array}{c}\frac{\∂\mathrm{\φ}}{\∂t}=-{\mathrm{\ξ}}_1\mathrm{\δ}\left(\mathrm{\φ}\right)(e_{11}-e_{12})-{\mathrm{\ξ}}_2\mathrm{\δ}\left(\mathrm{\φ}\right)(e_{21}-e_{22})\\ +v\mathrm{\δ}\left(\mathrm{\φ}\right)div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\mathrm{\μ}\left({\▿}^2\mathrm{\φ}-div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)\right)\end{array}---\left(21\right)$

δ为Dirac函数，即式(9)中Heaviside函数的导数，e₁₁、e₁₂、e₂₁和e₂₂分别为：

$e_{11}\left(x\right)=\∫\mathrm{\ω}(x-y)\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-in}^i\left(x\right)|dy---\left(22\right)$

$e_{11}\left(x\right)=\∫\mathrm{\ω}(x-y)\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-out}^i\left(x\right)|dy---\left(23\right)$

$e_{21}\left(x\right)=\∫\mathrm{\ω}(x-y)\[ln\left({\mathrm{\σ}}_1\right(y\left)\right)+\frac{{\left(u_1\right(y)-I(x\left)\right)}^2}{{\mathrm{\σ}}_1{\left(y\right)}^2}\]---\left(24\right)$

$e_{22}\left(x\right)=\∫\mathrm{\ω}(x-y)\[ln\left({\mathrm{\σ}}_2\right(y\left)\right)+\frac{{\left(u_2\right(y)-I(x\left)\right)}^2}{{\mathrm{\σ}}_2{\left(y\right)}^2}\]---\left(25\right)$

步骤6.2，将梯度下降流公式转化为LBM形式，以便更新水平集函数；

其中，根据LBM，梯度下降流公式的流体扩散方程形式为：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\γ}div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+F---\left(26\right)$

$\mathrm{\γ}=\frac{2}{9}(2\mathrm{\τ}-1)---\left(27\right)$

联立梯度下降流公式(15)与LBM公式(20)、式(21)，解得LBM的松弛时间τ和内在引导力F为：

$\mathrm{\τ}=\frac{9}{4}\left(\mathrm{\δ}\right(\mathrm{\φ})-\mathrm{\μ})+\frac{1}{2}---\left(28\right)$

F＝-ξ₁δ(φ)(e₁₁-e₁₂)-ξ₂δ(φ)(e₂₁-e₂₂)+μ▽²φ (29)

对于每个LBM中的晶格(在本方法中即像素)，它们的碰撞、对流和反弹结果分别为：

$\begin{array}{c}{f}_c(x^{*},t+\mathrm{\Δ}t)=f_c(x,t)+\frac{1}{\mathrm{\τ}}(A_c\mathrm{\φ}-f_c(x,t\left)\right)\\ +3\×\frac{2\mathrm{\τ}-1}{2\mathrm{\τ}}{A}_c(F\·v_c)\end{array}---\left(30\right)$

f_c(x+v_cΔt,t+Δt)＝f_c(x^*,t+Δt) (31)

f_c(x+v_cΔt,t+Δt)＝f_c(x+v_-cΔt,t+Δt) (32)

其中A_c为LBM中的方向参数，v_c为晶格的速度向量，v_-c为与v_c方向相反的速度向量。叠加以上三种结果，利用式(33)更新水平集函数：

$\mathrm{\φ}=\underset{c}{\Σ}{f}_c---\left(33\right)$

步骤6.3，重复步骤6.2的水平集函数更新过程，直到水平集函数收敛或误差达到预设的阈值内，即可得到曲线演化结果，即纹理图像分割曲线为：

C(t+Δt):＝{x|φ(x)＝0} (34)

至此，从步骤1到步骤6，完成了一种基于水平集模型的纹理图像分割方法。

有益效果

一种基于水平集模型的纹理图像分割方法，对比现有技术，具有如下有益效果：

1.本发明所提出的在水平集模型中融入纹理特征的方法，能够有效地对纹理图像进行分割，结果优于现有的纹理分割方法；

由于LSS描述符独立于图像灰度和色彩，只表示图像纹理特征的优良性质，利用LSS构造纹理能量项能够充分体现图像不同区域的纹理特征，因而本方法能够有效区分和分割各种纹理图像，在各种纹理图像上的分割结果也验证了本方法的有效性；

2.本发明所提出的能量函数包含了LGDF能量项，对灰度不均和噪声纹理图像具有鲁棒性；

3.本发明引入了LBM加速曲线演化过程，有效减少了演化所需迭代次数，提高了图像分割效率。

附图说明

图1为本发明“一种基于水平集模型的纹理图像分割方法”及实施例1中的流程示意图；

图2为本发明“一种基于水平集模型的纹理图像分割方法”实施例2中，对于纹理图像的分割结果及其与现有方法的对比；

图3为本发明“一种基于水平集模型的纹理图像分割方法”实施例2中，对于灰度不均或有噪声纹理图像的分割结果及其与现有方法的对比。

具体实施方式

为使发明目的、技术方案及优点更加清楚明白，下面将结合附图对本发明的实施例进行详细描述。以下实施例以本发明的技术方案为前提进行实施，给出详细实施方式和具体操作过程，但本发明的保护范围不限于以下实施例。

实施例1

本实施例阐述了将本发明“一种基于水平集模型的纹理图像分割方法”应用于图像“Leopard1.bmp”的流程：

图1为本方法及本实施例的算法流程，从图1可以看出，本方法包括以下步骤：

步骤A：Log-Gabor滤波；

具体到本实施例中，输入的图像I为“Leopard1.bmp”，并对I进行Log-Gabor滤波得到I_G；

本实施例中Log-Gabor参数取：滤波频率参数σ_f＝k/w₀＝0.65，角度参数σ_θ＝0.52，方向参数θ为0，小波尺度S为3；

步骤B：提取LSS描述符；

具体过程可参考步骤2，区域大小N取41，小块大小n取5，块数M取40(ρ均分为4个区间，θ均分为10个区间)。对每个像素p计算其对应的LSS描述符L_p，获取每个描述符L_p的40个分量再按照像素p的空间域顺序，将描述符分量构成一系列新图像

步骤C：构造纹理能量项；

具体到实施例中，详细步骤与步骤3相同。其中，参数λ₁＝1.0，λ₂＝1.0，Heaviside函数的参数ε＝1.0，权重函数的参数σ＝5.0，距离阈值ρ＝11；

步骤D：构造LGDF能量项；

具体到本实施例中，从输入图像I中构造LGDF能量项的参数分别为：λ₁＝1.0，λ₂＝1.0，Heaviside函数参数ε＝1.0，权重函数参数σ＝5.0，距离阈值ρ＝11；

步骤E：合并步骤C和步骤D的能量项，构造能量函数；

本实施例取ξ₁＝5.0，ξ₂＝10.0，ν＝0.0002*255²，μ＝1.0；

步骤F：使用LBM加速求解能量函数；

具体方法见步骤6，取LBM方向参数A_c(c∈{0,1,2,3,4,5,6,7,8})分别为：A₀＝4/9，A_1,2,3,4＝1/9，A_5,6,7,8＝1/36，速度向量v_c取单位向量，从而可对能量函数的梯度下降流公式进行LBM快速演化，经过20次迭代获得最终的分割曲线；

至此，从步骤A到F，完成了本实施例一种基于水平集模型的纹理图像分割方法。

实施例2

本实施例具体阐述了在6个纹理图像上执行本发明步骤1到步骤6所得的分割结果(分割结果用白色线条表示)，同时与现有的图像分割方法进行比较。

图2和图3分别为3张纹理图像，以及3张有噪声/灰度不均纹理图像的分割结果。此6张纹理图像即步骤1中的输入图像；

图2分为5行3列，每行代表使用一种图像分割方法所得的结果(前两行分别为原始图像和人工分割的标准结果)，每列代表不同的原始图像。其中，对比的方法有：LGDF算法、ASLVD算法和本方法；共有3张原始图像，从左至右分别为：“Leopard1.bmp”、“Zebra1.bmp”和“Zebra2.png”；

图3分为4行3列，每行代表使用一种图像分割方法所得的结果(第一行为原始图像)，每列代表不同的原始图像。其中，对比的方法有：CV算法、LGDF算法和本方法；共有3张原始图像，从左至右分别为：“MAP.jpg”、“LSS.jpg”和“star.jpg”。

本实施例中进一步利用JS(Jaccard Similarity)相似性来度量CV、LGDF、ASLVD算法和本方法图像分割的准确性，JS相似性定义为：

$JS(S_R,S_{GT})=\frac{|S_R\∩S_{GT}|}{|S_R\∪S_{GT}|}---\left(35\right)$

其中S_R为图像分割方法所得的结果，S_GT为同一图像人工分割的标准结果，|·|表示分割曲线所围成的区域。

从图2能够明显看出LGDF算法和ASLVD算法所得结果和本方法相比具有一定瑕疵(例如在“Leopard1.bmp”上，LGDF算法未能够分割“耳朵”部分，ASLVD算法对物体边缘明显变化的区域处理不够精确；其他两幅图像相比本方法也各有不足)。它们与标准结果的JS相似性分别为：

表1不同方法相分割结果相对于人工分割标准结果的JS相似性

本方法在3张图像上的分割结果都具有最高的JS相似性，说明本方法最接近于人工分割的标准结果，在纹理图像分割上优于其他两种现有方法。

从图3中能够看出CV算法明显受限于“MAP.jpg”中的噪声和“LSS.jpg”与“star.jpg”中的灰度不均；而LGDF算法虽然在“MAP.jpg”上取得了一定效果，但是相比本方法所得结果具有较强的边缘抖动效应，效果仍不理想，且在其他两幅图像上分割结果较差。这证明了在灰度不均和噪声图像上，本方法相比现有方法具有一定优势。

实施例3

本实施例为本方法分割图像所需的曲线演化迭代次数，与去除LBM加速单元情况下的迭代次数之对比。两种情况下，处理不同图像时所需的迭代次数分别如表2所示：

表2本方法与去除LBM单元情况下的迭代次数对比

能够看出去除LBM加速单元的情况下，迭代次数明显增加；在6幅图像上，加速后的效率比加速前提升了2～6倍，这表明了LBM加速的有效性。引入LBM明显减少了本方法的曲线演化迭代次数，减少了演化所需时间，提高了图像分割效率。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡是在本发明的精神和原则内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于水平集模型的纹理图像分割方法，其特征在于：

一种基于水平集模型的纹理图像分割方法，简称本方法；

具体包括以下步骤：

步骤1，对输入的一幅图像进行Log-Gabor滤波，以增强图像的纹理特征，输出Log-Gabor滤波后的图像；

步骤2，对经过步骤1Log-Gabor滤波后的图像提取LSS描述符，得到由描述符分量构成的图像；

步骤3，在步骤2输出的由LSS描述符分量构成的图像上构造纹理能量项；

步骤4，在步骤1中的输入图像上构造LGDF能量项；

步骤5，综合步骤4输出的LGDF能量项和步骤3输出的LSS纹理能量项，组成本方法的能量函数；

步骤6，通过最小化步骤5所得的能量函数求解能量函数的分割曲线；同时，在此过程中使用LBM加速曲线演化，减少迭代次数；

至此，从步骤1到步骤6，完成了一种基于水平集模型的纹理图像分割方法。

2.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤1中，记输入图像为I，Log-Gabor滤波后输出的图像为I_G。Log-Gabor滤波需在图像的频率域上进行卷积操作，则输出I_G与输入I的关系为：

表示傅里叶变换，*为卷积符号，G为Log-Gabor的传递函数，即：

G(w；σ_f，σ_θ，θ,S)＝exp(-log(w/w₀)/2(log(k/w₀))²) (2)

其中，exp表示以常数e为底的指数函数，log表示自然对数；w为图像I的像素点在频率域上对应的值，w₀为在频率域上滤波中心的像素点的值，k为控制传递函数形状的参数，通常取k∈{k|0.55w₀≤k≤0.85w₀}，以便使得滤波频率参数σ_f＝k/w₀∈[0.55,0.85]，σ_θ为角度参数，通常取0.41～0.74，θ为方向参数，取0～360，S为小波尺度，通常为1～3。

3.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤2，具体为：

步骤2.1，对于I_G的每一个像素p，以p为中心选取N×N像素的区域(N通常取21～61的奇数)，若图像边界处区域不足N×N大小，则复制图像边界处的值来进行扩展，使得区域大小为N×N；

例如，假设图像左上角边界处如式(3)所示，扩展2像素，则扩展后的图像为式(4)：

在每个N×N区域中，以p为中心选取n×n像素的小块(n通常取3～9的奇数)，记为t_p；对于N×N区域的每一个像素x(可以为p也可以非p)，以x为中心选取n×n像素的小块，记为t_x，若在区域边界处小块不足n×n大小，同样复制区域边界处的值来进行扩展；

计算t_p和t_x的平方差(Sum of Square Differences，SSD)，即：

$SSD(t_p,t_x)=\left|\right|t_p-t_x|{|}_F^2---\left(5\right)$

||·||_F为F-范数；

步骤2.2，根据步骤2.1输出的SSD计算N×N区域上的相关矩阵C_p(C_p的大小为N×N)：

$C_p\left(x\right)=\exp (-\frac{SSD(t_p,t_x)}{\mathrm{\α}})---\left(6\right)$

α为自适应参数，可取x所在的N×N区域中最大的SSD值；

步骤2.3，将C_p从直角坐标转换到对数极(Log-polar)坐标上；

步骤2.4，按照Log-polar坐标(ρ,θ)将C_p划分为M块，再对每个块取最大值，将这些最大值拼接起来即LSS描述符L_p：

$L_p(\mathrm{\ρ},\mathrm{\θ})=\underset{x\∈BIN(\mathrm{\ρ},\mathrm{\θ})}{max}\left\{C_p\left(x\right)\right\}---\left(7\right)$

其中，每个块记为BIN(ρ,θ)；是L_p的第i个描述符分量，将按照像素p的空间域顺序构成则就是步骤2输出的由LSS描述符分量构成的图像，i的取值范围是i∈{1,2,...,M}。

4.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤3，具体为：

记图像的每个像素点为x，x所具有的纹理能量为：

$\begin{array}{c}{E}_x^{LSS}={\mathrm{\λ}}_1\mathrm{\ω}\left(x-y\right)\∫\underset{i=1}{\overset{M}{\mathrm{\Σ}}}|L_p^i\left(y\right)-L_{p-in}^i\left(x\right)|H\left(\mathrm{\φ}\left(y\right)\right)dy\\ +{\mathrm{\λ}}_2\mathrm{\ω}\left(x-y\right)\∫\underset{i=1}{\overset{M}{\mathrm{\Σ}}}|L_p^i\left(y\right)-L_{p-out}^i\left(x\right)|\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)dy\end{array}---\left(8\right)$

y为x邻域内的像素，φ为水平集函数，λ₁和λ₂分别为曲线内部平衡参数和曲线外部平衡参数，通常取0.5～2.0，H为Heaviside函数，即：

$H\left(\mathrm{\φ}\right)=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan (\frac{\mathrm{\φ}}{\mathrm{\ϵ}}\left)\right)---\left(9\right)$

ε为Heavisde函数参数，通常取0.1～10；ω(d)为权重函数，其中d＝x-y：

$\mathrm{\ω}\left(d\right)=\left\{\begin{array}{cc}\frac{1}{a}\exp (-\frac{|d{|}^2}{2{\mathrm{\σ}}^2})& if\left|d\right|\≤\mathrm{\ρ}\\ 0& if\left|d\right|>\mathrm{\ρ}\end{array}\right.---\left(10\right)$

a和σ为权重函数参数，ρ为距离阈值；和分别为演化曲线内部和外部的均值：

$L_{p-in}^i\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)L_p^i\left(y\right)H\left(\mathrm{\φ}\right(y\left)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\right(y\left)\right)dy}---\left(11\right)$

$L_{p-out}^i\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)L_p^i\left(y\right)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}{\∫\mathrm{\ω}(y-x)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}---\left(12\right)$

则纹理能量项E^LSS如式(13)：

$E^{LSS}={\∫}_{\mathrm{\Ω}}{E}_x^{LSS}dx---\left(13\right).$

5.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤4中，LGDF能量项记为E^LGDF，通过式(14)计算：

$\begin{array}{c}{E}^{LGDF}=-{\mathrm{\λ}}_1\∫\∫\mathrm{\ω}\left(x-y\right)\ln p_{1,x}\left(I\left(y\right)H\left(\mathrm{\φ}\left(y\right)\right)\right)dydx\\ -{\mathrm{\λ}}_2\∫\∫\mathrm{\ω}\left(x-y\right)\ln p_{2,x}\left(I\left(y\right)\left(1-H\left(\mathrm{\φ}\left(y\right)\right)\right)\right)dydx\end{array}---\left(14\right)$

p_j,x(I(y))(其中j∈{1,2})为概率分布：

$p_{j,x}\left(I\right(y\left)\right)=\frac{1}{\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_j\left(x\right)}\exp (-\frac{{\left(u_j\right(x)-I(y\left)\right)}^2}{2{\mathrm{\σ}}_j{\left(x\right)}^2})---\left(15\right)$

其中u₁、u₂、σ₁(x)²和σ₂(x)²分别为：

$u_1\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)I\left(y\right)H\left(\mathrm{\φ}\right(y\left)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\right(y\left)\right)dy}---\left(16\right)$

$u_2\left(x\right)=\frac{\∫\mathrm{\ω}(y-x)I\left(y\right)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}{\∫\mathrm{\ω}(y-x)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}---\left(17\right)$

${\mathrm{\σ}}_1{\left(x\right)}^2=\frac{\∫\mathrm{\ω}(y-x){\left(u_1\right(x)-I(y\left)\right)}^{2}H\left(\mathrm{\φ}\right(y\left)\right)dy}{\∫\mathrm{\ω}(y-x)H\left(\mathrm{\φ}\right(y\left)\right)dy}---\left(18\right)$

${\mathrm{\σ}}_2{\left(x\right)}^2=\frac{\∫\mathrm{\ω}(y-x){\left(u_2\right(x)-I(y\left)\right)}^2(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}{\∫\mathrm{\ω}(y-x)(1-H(\mathrm{\φ}\left(y\right)\left)\right)dy}---\left(19\right).$

6.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤5中，记本方法的能量函数为E^LL，则：

E^LL＝ξ₁E^LSS+ξ₂E^LGDF+ν·Length(C)+μ·P(x) (20)

ξ₁和ξ₂分别为LSS权重参数和LGDF权重参数，通常取1～40，为曲线长度，为防止重新初始化项，避免了曲线演化过程中的重新初始化过程，ν和μ分别为曲线长度和防止重新初始化项的系数，ν通常取0.0001*255²～0.01*255²，μ通常取0.1～1.5。

7.如权利要求1所述的一种基于水平集模型的纹理图像分割方法，其特征还在于：

步骤6，具体为：

步骤6.1，通过步骤5所得能量函数获取梯度下降流公式。其中，式(20)能量函数对应的梯度下降流公式为：

$\begin{array}{c}\frac{\∂\mathrm{\φ}}{\∂t}=-{\mathrm{\ξ}}_1\mathrm{\δ}\left(\mathrm{\φ}\right)(e_{11}-e_{12})-{\mathrm{\ξ}}_2\mathrm{\δ}\left(\mathrm{\φ}\right)\left(e_{21}-e_{22}\right)\\ +\mathrm{\ν}\mathrm{\δ}\left(\mathrm{\φ}\right)div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\mathrm{\μ}\left({\▿}^2\mathrm{\φ}-div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)\right)\end{array}---\left(21\right)$

δ为Dirac函数，即式(9)中Heaviside函数的导数，e₁₁、e₁₂、e₂₁和e₂₂分别为：

$e_{11}\left(x\right)=\∫\mathrm{\ω}(x-y)\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-in}^i\left(x\right)|dy---\left(22\right)$

$e_{11}\left(x\right)=\∫\mathrm{\ω}(x-y)\underset{i=1}{\overset{M}{\Σ}}|L_p^i\left(y\right)-L_{p-out}^i\left(x\right)|dy---\left(23\right)$

$e_{21}\left(x\right)=\∫\mathrm{\ω}(x-y)\[ln\left({\mathrm{\σ}}_1\right(y\left)\right)+\frac{{\left(u_1\right(y)-I(x\left)\right)}^2}{{\mathrm{\σ}}_1{\left(y\right)}^2}\]---\left(24\right)$

$e_{22}\left(x\right)=\∫\mathrm{\ω}(x-y)\[ln\left({\mathrm{\σ}}_2\right(y\left)\right)+\frac{{\left(u_2\right(y)-I(x\left)\right)}^2}{{\mathrm{\σ}}_2{\left(y\right)}^2}\]---\left(25\right)$

步骤6.2，将梯度下降流公式转化为LBM形式，以便更新水平集函数；

其中，根据LBM，梯度下降流公式的流体扩散方程形式为：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\γ}div\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+F---\left(26\right)$

$\mathrm{\γ}=\frac{2}{9}(2\mathrm{\τ}-1)---\left(27\right)$

联立梯度下降流公式(15)与LBM公式(20)、式(21)，解得LBM的松弛时间τ和内在引导力F为：

$\mathrm{\τ}=\frac{9}{4}\left(\mathrm{\δ}\right(\mathrm{\φ})-\mathrm{\μ})+\frac{1}{2}---\left(28\right)$

$F=-{\mathrm{\ξ}}_1\mathrm{\δ}\left(\mathrm{\φ}\right)(e_{11}-e_{12})-{\mathrm{\ξ}}_2\mathrm{\δ}\left(\mathrm{\φ}\right)(e_{21}-e_{22})+\mathrm{\μ}{\▿}^2\mathrm{\φ}---\left(29\right)$

对于每个LBM中的晶格(在本方法中即像素)，它们的碰撞、对流和反弹结果分别为：

$\begin{array}{c}{f}_c(x^{*},t+\mathrm{\Δ}t)=f_c(x,t)+\frac{1}{\mathrm{\τ}}(A_c\mathrm{\φ}-f_c(x,t\left)\right)\\ +3\×\frac{2\mathrm{\τ}-1}{2\mathrm{\τ}}{A}_c\left(F\·v_c\right)\end{array}---\left(30\right)$

f_c(x+v_cΔt,t+Δt)＝f_c(x^*,t+Δt) (31)

f_c(x+v_cΔt,t+Δt)＝f_c(x+v_-cΔt,t+Δt) (32)

其中A_c为LBM中的方向参数，v_c为晶格的速度向量，v_-c为与v_c方向相反的速度向量。叠加以上三种结果，利用式(33)更新水平集函数：

$\mathrm{\φ}=\underset{c}{\Σ}{f}_c---\left(33\right)$

步骤6.3，重复步骤6.2的水平集函数更新过程，直到水平集函数收敛或误差达到预设的阈值内，即可得到曲线演化结果，即纹理图像分割曲线为：

C(t+Δt):＝{x|φ(x)＝0} (34)。