CN102819840B - 一种纹理图像的分割方法 - Google Patents

Abstract

本发明公开了一种纹理图像分割方法，提取图像的局部的各种特征信息表征图像中不同纹理的区域的特性得到特征图像，通过主成分分析等手段降低数据量再通过mean shift算法（中文名为：均值漂移算法）对特征向量进行聚类,从而完成较可靠的纹理图像分割。本发明相对于现有技术，图像分割可靠性高，人机交互少的有益效果。

Description

一种纹理图像的分割方法

技术领域

本发明涉及计算机图像处理领域，特别涉及一种纹理图像的分割方法。

背景技术

图像分割是图像分析的重要过程，图像分割一般认为是分割图像并找出图像中不同的物体以找到一种我们感兴趣的图像数据的紧凑描述。而传统的图像分割技术对纹理图像中的噪声敏感,易造成过分割现象。纹理作为图像的一个重要特征,在计算机视觉和图像处理中有重要应用例如医学早期的癌细胞识别诊断遥感图像中对军事及民用目标的识别。

对目前纹理图像的分割分为三大类:统计法,结构法和时域/频域方法。

发明内容

本发明所要解决的技术问题是，提供一种计算速度快、图像分割准确的纹理图像分割方法。

为便于理解，下文首先对本发明的基本原理阐述如下：

本方法基于时域和频域的方法使用近年来发展迅速的Gabor滤波器(中文名为:伽柏滤波器)提取图像的局部的各种特征信息表征图像中不同纹理的区域的特性得到特征图像，通过主成分分析等手段降低数据量再通过mean shift算法（中文名为：均值漂移算法）对特征向量进行聚类,从而完成较可靠的纹理图像分割。

Gabor滤波器由于其在时域和频域具有很好的分辨能力,被广泛用在计算机视觉和图像处理领域。作为在方向和频率上有选择性的带通滤波器，它的冲击响应是高斯包络的复数正弦函数。Jones和Palmer(1987)指出复高斯函数与猫的视觉接受场函数相似,所以通过对生物视觉的分析对Gabor滤波器的参数选择具有启发性。

mean shift这个概念最早是由Fukunaga等人在一篇关于概率密度梯度函数的估计中提出来的,其最初含义正如其名,就是偏移的均值向量,在这里mean shift是一个名词,它指代的是一个向量,但随着Mean Shift理论的发展,mean shift的含义也发生了变化,如果我们说mean shift算法,一般是指一个迭代的步骤,即先算出当前点的偏移均值,移动该点到其偏移均值,然后以此为新的起始点,继续移动,直到满足一定的条件结束。与其他常用的聚类算法相比本算法能自动确定聚类数目，广泛用在图像分割，图像跟踪等领域。然而对某些纹理丰富的图像直接使用mean shift算法容易导致过分割的现象。所以方法中提出先使用Gabor滤波器分解图像，再使用mean shift算法进行聚类,有效抑制了过分现象，并完成图像的自动分割。

本发明采用以下技术手段解决所述技术问题：

一种纹理图像分割方法，包括以下步骤：

1）对输入的图像I(x,y)做傅里叶变换,得到图像的频率域数据；

2）根据预定的旋转方向和滤波器频段构造一组Gabor滤波器：

{H_k(x,y)|k=1...N}，其中，H_k表示这组滤波器中第K个滤波器的核函数,N为滤波器的个数；

3）在频率域使用所述Gabor滤波器对图像滤波，得到一组滤波结果{I′_k(x,y)|k=1...N}，其中，I′_k(x,y)表示第k个滤波器对图像滤波的结果,k为滤波器个数,x、y为像素坐标；

4）对滤波结果I′_k(x,y)进行特征抽取，抽取的特征包括纹理的能量E_d(x,y)、反向灰度值M_d(x,y)和方向灰度方差V_d(x,y)，从而得到各频段和方向下的特征向量：

f_k(x,y)=[E_d(x,y),M_d(x,y),V_d(x,y)]|k=1..N，其中，f_k(x,y)为行向量,x,y是对应图像中的像素坐标；

5）使用meanshift算法对步骤4）得到的特征向量进行聚类，完成图像分割。

优选地：

所述步骤5）包括以下步骤：使用PCA算法（中文名为：主成分分析算法）对所述步骤4）得到的特征向量进行降维得到降维后的特征向量，然后再使用meanshift算法对该降维后的特征向量进行聚类完成图像分割。

所述预定的旋转方向为0°,45°,90°和135°，所述Gabor滤波器的频程为1~2,分为5个频段，其中，中频段的中心频率取0.25。

所述步骤4）中，

$E_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′2}(i,j),$

$M_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′}(i,j),$

$V_d(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[I_d^{\′}(i,j)-M_k(x,y)]}^2,$

其中，Ω为中心位于(x,y)大小为W_x×W_y的矩形窗口。

所述步骤3）包括以下步骤：3.1）将步骤1）所得到的频率域数据与步骤2）中构造的滤波器组分别相乘后得到滤波后的一组数据{I_k(x,y)|k=1...N}；3.2）对{I_k(x,y)|k=1...N}取傅里叶反变换，得到一组经过滤波的图像{I′_k(x,y)|k=1...N}。

与现有技术相比，本方法中在频率域使用Gabor滤波器组对图像滤波可提高运算速度,减少因截断和采样而造成的混叠；同时Gabor滤波器的使用有效抑制了因图像中含有复杂纹理而造成的mean shift过分割现象；而通过计算衡量纹理的多个特征，确保了对纹理分割的可靠性。

优选方案中，通过使用PCA算法可以对得到的特征向量进行降维,去除噪声和冗余信息的干扰并减少计算量；通过使用means shift算法对纹理特征向量聚类可以自动无监督的完成图像分割,降低人工的交互操作；对高频滤波器及低频滤波器频率极向中心的选择，有利于对纹理信息集中的中频段部分的处理。

附图说明

图1为Gabor滤波器的三维图像形状示意图；

图2a-2d分别是取旋转角度0°,45°,90°,135°在频域得到的滤波器平面形状示意图；

图3是分析得到的滤波器的性质示意图；

图4是具体实施例的输入图像；

图5是对图4的图像做傅里叶变换获得的频域图像的幅值；

图6是图4的图像在滤波器方向为为0°频率峰值0.25,频程1.2时的滤波结果。

图7是mean shift的原理示意图。

具体实施方式

下面对照附图并结合优选的实施方式对本发明作进一步说明。

一种纹理图像分割方法，其特征在于，包括以下步骤：

1）对输入的图像I(x,y)做傅里叶变换,得到图像的频率域数据；

2）根据预定的旋转方向和滤波器频段构造一组Gabor滤波器：

{H_k(x,y)|k=1...N}，其中，H_k表示这组滤波器中第K个滤波器的核函数,N为滤波器的个数；

3）在频率域使用所述Gabor滤波器对图像滤波，得到一组滤波结果{I′_k(x,y)|k=1...N}，其中，I′_k(x,y)表示第k个滤波器对图像滤波的结果,k为滤波器个数,x、y为像素坐标；

4）对滤波结果I′_k(x,y)进行特征抽取，抽取的特征包括纹理的能量E_d(x,y)、反向灰度值M_d(x,y)和方向灰度方差V_d(x,y)，从而得到各频段和方向下的特征向量：

f_k(x,y)=[E_d(x,y),M_d(x,y),V_d(x,y)]|k=1..N，其中，f_k(x,y)为行向量,x,y是对应图像中的像素坐标；

5）使用meanshift算法对步骤4）得到的特征向量进行聚类，完成图像分割。

为助于理解，下面结合更加具体的实施例，对上述方法各个步骤进行说明，但本发明并不限于下述实施例：

2.1设计Gabor滤波器

Gabor滤波器被广泛用于图像处理领域，但为便于公众理解，本节依然首先对Gabor滤波器的原理及公式演变予以详细说明：

二维Gabor函数是一个纵横比为λ,尺度参数为δ,主轴与x轴成角度φ的高斯调制的负指数正弦函数,其一般形式是

h(x,y)=g(x′,y′)exp[-2πj(u₀x+v₀y)]                        (1)

式中

(x′,y′)=(xcosφ+ysinφ,-xsinφ+ycosφ)                    (2)

为旋转坐标,u₀,v₀对应频域中的频率中心。

式(1)中的二维高斯函数为:

$g(x,y)=\frac{1}{2\mathrm{\π\λ}{\mathrm{\δ}}^2}\exp [-\frac{{(x/\mathrm{\λ})}^2+y^2}{2{\mathrm{\δ}}^2}]---\left(3\right)$

λ为纵横比

对公式(1)作二维傅里叶变换可得

H(u,v)=exp{-2π²δ²[(u'-u′₀)²λ²+(v'-v′₀)²]}              (4)式中

(u',v')=(ucosφ+vsinφ,-usinφ+vcosφ)

由式(4)可知H(u,v)是一个主轴与μ轴成φ角的带通高斯函数,其纵横比为1/λ,u₀,v₀为滤波器频率中心，极向中心频率为(单位为cycle/piexs)方向为θ=arctan(u₀/v₀)。

通常情况下,取Gabor函数中高斯函数部分的方向与复数调制部分函数的辅角相等(φ=θ),即滤波器在空间域和频率域有相同的取向,此时公式(1)和(4)为:

h(x,y)=g(x′,y′)exp[2πju₀x′]                           (6)

H(u,v)＝exp{-2π²δ²[(u′-F)²λ²+v'²]}                    (7)

将(6)展开得到

h(x,y)=g(x′,y′)[cos(2πFx′)+jsin(2πFx′)]             (8)

对(8)的实部做傅里叶变换得到滤波器的核函数公式：

$H_k(u,v)=\mathrm{Aexp}\{-\frac{1}{2}[\frac{{(u^{\′}-F)}^2}{{\mathrm{\δ}}_u^2}+\frac{{v^{\′}}^2}{{\mathrm{\δ}}_v^2}\left]\right\}+\mathrm{Aexp}\{-\frac{1}{2}[\frac{{(u^{\′}+F)}^2}{{\mathrm{\δ}}_u^2}+\frac{{v^{\′}}^2}{{\mathrm{\δ}}_v^2}\left]\right\}---\left(9\right)$

式中：

δu_＝1/2πδ_x;δ_v＝1/2πδ_y;A=2πδ_xδ_y，

(u',v')=(ucosφ+vsinφ,-usinφ+vcosφ)，

(u₀v₀)表示不同频段滤波器的的频率中心,δ_x,δ_y表示时域滤波器的对称轴长度，δ_u,δ_v表示频域滤波器对称轴的长度,根据选择的频程和旋转角度计算获得；φ为旋转角度。

如图1所示，为滤波器的三维图像形状示意图。选择合适的旋转角度和频段即可构建一组滤波器，例如：

分别取不同的旋转角度在频域得到的滤波器平面形状示意图如图2a-2d所示。图像中长宽比取δ_v/δ_u=0.9499，旋转角度φ取(0°,45°,90°,135°)。

对Gabor滤波器的参数选择已经有很长时间的研究,大部分基于生物心理学，其中的一个重要发现是人类视觉细胞的空间-频率域带宽在0.5到2.5频程。对本文使用的Gabor滤波器而言,

$\mathrm{\ΔF}={\log }_2\frac{F+{\mathrm{\δ}}_{u}C}{F-{\mathrm{\δ}}_{u}C}---\left(10\right)$

ΔF为滤波器频程,其中F为滤波器中心频率。

又根据滤波器的形状可近似得到

$\tan \left(\frac{1}{2}\mathrm{\Δ\φ}\right)=\frac{{\mathrm{\δ}}_{v}C}{F}---\left(11\right)$

由公式（10）（11）可知，确定相邻旋转角度差Δφ，滤波器组频程ΔF，和某一频段下的滤波器中心频率F,滤波器频段个数就可确定这个滤波器组。本算法中取频程为1.2，分为5个频段，相邻旋转角度差Δφ为45°。实际实施过程中，也可取其他的频程，但实验表明频程为1~2,分成5个频段，相邻旋转角度差Δφ取45°能够在较好的兼顾处理效果和运算速度。

由于纹理信息多集中在中频段部分,在本发明中选取的Gabor滤波器组的第3个频段的中心频率取0.25.第1,2和第4,5频段的频率中心根据(10)（11）获得。

此处根据分析得到的滤波器的形状如图3所示，由图示可见不同频段（图中取了三个频段）,不同角度的滤波器的混叠并不明显。

2.2频率域滤波

对输入的图像I(x,y)做傅里叶变换,得到图像的频率域数据。在频率域使用Gabor滤波器对图像滤波。例如，输入的图像如图4所示则得到如图5所示的频域图像的幅值，图5用对数函数对灰度做过拉伸处理。

将所述I(x,y)的频率域数据与2.1中设计的滤波器组分别相乘得到滤波后的一组数据{I_k(x,y)|k=1...N},N表示选定的滤波器个数，对{I_k(x,y)|k=1...N}取傅里叶反变换得到一组经过滤波的图像{I′_k(x,y)|k=1...N}，如图6所示，其所示图像是在滤波器方向为0°频率峰值0.25,频程1.2时的滤波结果。

2.3抽取特征

对得到的滤波图像,抽取纹理特征,抽取以下三种特征：

能量定义为

$E_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′2}(i,j)---\left(12\right)$

方向灰度均值定义为：

$M_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′}(i,j)---\left(13\right)$

方向灰度值方差定义为

$V_d(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[I_d^{\′}(i,j)-M_k(x,y)]}^2---\left(14\right)$

上述各式中，Ω为中心位于(x,y)大小为W_x×W_y的矩形窗口，W_x=W_y∈{3,5,7,9}。

通过对I′_k(x,y)的每个像素抽取的特征得到

{f_k(x,y)=[E_k(x,y),M_k(x,y),V_k(x,y)]|k=1...N}        (15)

f_k(x,y)为行向量,x,y是对应图像中的像素坐标。

2.4PCA特征降维处理

由上述过程得到的特征向量维数很大需要降维。以减少特征冗余，提高计算速度。本方法中使用PCA（Principal Component Analysis，主成分分析）算法对特征向量进行降维。PCA是寻找最小均方意义下最能够代表原始数据的交换方法，其具体方法如下：

将所有像素对应的特征向量组成一个矩阵假设像素的个数是M,每个像素对应的特征向量为X_i(i=1,2,...M)(行向量),则将X_i组成矩阵

X=[X₁,X₂,...,X_M]^T,                                 (16)

X的协方差矩阵定义为

$Q=\frac{X^{T}X}{M-1}$

求协方差矩阵Q的本征值λ_i及本征值对应的本征向量组E=[E₁,E₂,E₃,...,E_m],将E变换为单位正交矩阵,将本征矢量按照本征值的大小排序,取前p列本征向量即得到线性变换矩阵W.因此得到新的特征向量:

y_i=W^Tx′_i

其中y_i是p维的列向量,p的选取原则即保证取得的p个成分代表数据全部方差的比例是

$\mathrm{Ratio}=\frac{\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i}---\left(17\right)$

本发明中取该比例为90%来确定p.确定p之后,每个像素对应的特征向量变为X′_i=[x₁,x₂,..,x_p],(i=1，M为像素索引,p为使用PCA确定的特征向量维数)。

2.5使用Mean Shift聚类算法完成分割

给定d维空间R^d中的n个样本点x_i,i=1,…,n,在x点的Mean Shift向量的基本形式定义为:

$M_h\left(x\right)=\frac{1}{g}\underset{g\∈S_i}{\mathrm{\Σ}}(x_i-x)---\left(18\right)$

其中,S_h是一个半径为h的高维球区域,满足以下关系的y点的集合,

S_h(x)≡{y|(y-x)^T(y-x)≤h²}                     (19)

g表示在这n个样本点x_i中,有g个点落入S_k区域中。

我们可以看到(x_i-x)是样本点x_i相对于点x的偏移向量,式（18）定义的mean shift向量M_h(x)就是对落入区域S_h中的g个样本点相对于点x的偏移向量求和然后再平均.从直观上看,如果样本点x_i从一个概率密度函数f(x)中采样得到,由于非零的概率密度梯度指向概率密度增加最大的方向,因此从平均上来说,S_h区域内的样本点更多的落在沿着概率密度梯度的方向。因此,对应的,mean shift向量M_h(x)应该指向概率密度梯度的方向。

如图7所示,大圆圈所圈定的范围就是S_h,小圆圈代表落入S_h区域内的样本点x_i∈S_h,黑点就是Mean Shift的基准点x,箭头表示样本点相对于基准点x的偏移向量。很明显的可以看出,平均的偏移向量M_h(x)会指向样本分布最多的区域,也就是概率密度函数的梯度方向。为方便讨论,把2.4中得到的特征向量写为以下形式

X_i=[x₁,x₂,...,x_p]                              (21）

考虑到像素的坐标位置，特征向量变为

X′_i=[X^s,X_i]                                   (22)

其中为像素i在图像中的坐标。

我们分别用x_i和z_i,i=1,…,n表示原始和分割后的图像.用Mean Shift算法对每个像素对应的特征向量进行聚类,的具体步骤如下,

对每一个象素点,

1,随机选取像素的特征向量X_i。

2,将X_i带入(18)计算M_h得到新的特征向量,将特征向量再代入式(18)反复迭代。

3，直到收敛到中心M，且再次迭代中心M不再变化。将所有可以收敛到M的点标记为同一类，即赋值

4.反复执行1~3过程,直到所有的像素对应的特征都被聚到某一类。

5.把中心很接近的类合并为同一类并计算该新类的中心。

根据每个像素的标记,把类内元素较少的类舍弃，至此完成图像分割。以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干等同替代或明显变型，而且性能或用途相同，都应当视为属于本发明的保护范围。

Claims (5)

1.一种纹理图像分割方法，其特征在于，包括以下步骤：

1)对输入的图像I(x,y)做傅里叶变换,得到图像的频率域数据；

2)根据预定的旋转方向和滤波器频段构造一组Gabor滤波器：

{H_k(x,y)|k＝1...N}，其中，H_k表示这组滤波器中第K个滤波器的核函数,N为滤波器的个数；

3)在频率域使用所述Gabor滤波器对图像滤波，得到一组滤波结果{I'_k(x,y)|k＝1...N}，其中，I'_k(x,y)表示第k个滤波器对图像滤波的结果,k为滤波器个数,x、y为像素坐标；

4)对滤波结果I'_k(x,y)进行特征抽取，抽取的特征包括纹理的能量E_d(x,y)、方向灰度均值M_d(x,y)和方向灰度方差V_d(x,y)，从而得到各频段和方向下的特征向量：

fk(x,y)＝[E_d(x,y),M_d(x,y),V_d(x,y)]|k＝1...N，其中，f_k(x,y)为行向量,x,y是对应图像中的像素坐标；

5)使用meanshift算法对步骤4)得到的特征向量进行聚类，完成图像分割；

所述步骤5)包括以下步骤：使用PCA算法对所述步骤4)得到的特征向量进行降维得到降维后的特征向量，在降维之后的特征向量基础上添加像素坐标获得包含像素坐标的特征向量，然后再使用meanshift算法对包含像素坐标的特征向量进行聚类完成图像分割。

2.根据权利要求1所述的纹理图像分割方法，其特征在于：所述预定的旋转方向为0°,45°,90°和135°，所述Gabor滤波器的频程取为1～2。

3.根据权利要求1所述的纹理图像分割方法，其特征在于：所述步骤4)中，

$E_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′2}(i,j),$

$M_k(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{I}_k^{\′}(i,j),$

$V_d(x,y)=\frac{1}{W_x\×W_y}\underset{(i,j)\∈\mathrm{\Ω}}{\mathrm{\Σ}}{[I_d^{\′}(i,j)-M_k(x,y)]}^2,$

其中，Ω为中心位于(x,y)大小为W_x×W_y的矩形窗口，其中，W_x＝W_y∈{3,5,7,9}。

4.根据权利要求1、2或3所述的纹理图像分割方法，其特征在于：所述Gabor滤波器使用5个频段，中频段的滤波器中心频率选择为0.25。

5.根据权利要求1、2或3所述的纹理图像分割方法，其特征在于：所述步骤3)包括以下步骤：3.1)将步骤1)所得到的频率域数据与步骤2)中构造的滤波器组分别相乘后得到滤波后的一组数据{I_k(x,y)|k＝1...N}；3.2)对{I_k(x,y)|k＝1...N}取傅里叶反变换，得到一组经过滤波的图像{I'_k(x,y)|k＝1...N}。