CN101447080B - 基于非下采样Contourlet变换的HMT图像分割方法 - Google Patents

Abstract

本发明公开了一种基于非下采样Contourlet变换的HMT图像分割方法，主要解决已有分割方法区域一致性和边缘保持差的缺点。其步骤为：(1)取待分割图像和每类的训练图像，分别对其进行非下采样Contourlet变换得到多尺度变换系数；(2)根据训练图像的非下采样Contourlet系数，按照一一对应的父子状态关系的隐马尔科夫树模型，估计模型参数；(3)计算待分割图像在各尺度系数子带对应的似然值，使用标记树融合多尺度似然函数得到多尺度最大后验概率分类；(4)对每一尺度依次根据上下文信息context-5模型，更新类别标记；(5)结合隐马尔科夫随机场模型和待分割图像相邻像素空间相关性信息，更新类别标记，得到最终分割结果。本发明具有区域一致性和边缘保持性好的优点，可用于合成纹理图像的分割。

Description

基于非下采样Contourlet变换的HMT图像分割方法 

技术领域

本发明属于图像处理技术领域，涉及一种多尺度几何分析技术在图像分割领域的应用，具体地说是一种图像分割方法，可用于纹理图像的分割处理。 

背景技术

图像分割是一种重要的图像处理技术。在图像的研究和应用中，人们往往对图像中的某些部分感兴趣，它们一般对应图像中特定的、具有独特性质的区域。为了辨识和分析目标，需要将它们分离提取出来，在此基础上才有可能对目标进一步处理。图像分割就是把图像分成各具特性的区域并提取出感兴趣的目标的技术和过程。这里特性可以是像素的灰度、颜色、纹理等，对应的目标可以对应单个区域或多个区域。 

近年来，小波分析已广泛应用于图像处理和分析中。小波对含点状奇异的目标函数而言是最优的基，但对具有直线或曲线状奇异的函数，小波系数则不是最稀疏的函数表示方法，所以小波在图像处理应用中并不能很好地表示图像中的轮廓边缘信息。多尺度几何分析是为分析二维或更高维奇异性产生的新工具。它致力于构建最优逼近意义下的高维函数表示方法，能够更优的表示图像的边缘信息，因此被广泛的应用在图像处理中。Arthur LCunha、Jianping Zhou和Minh N Do于2005年提出非下采样的Contourlet变换，并研究了相关应用。非下采样Contourlet变换NSCT是基于非下采样塔型结构和非下采样方向滤波带分解，经由àtrous algorithm实现并生成一种灵活的多尺度，多方向和平移不变的的图像分解方法。 

近年来，随着多尺度几何分析理论和方法的提出与发展，人们也提出了许多结合特定理论、方法和工具的图像分割技术。隐马尔科夫树模型HMT是近年来发展起来并应用于图像处理的一种参数统计模型，该模型可以充分挖掘多尺度变换域系数之间的统计特性，从而将系数分布的求解问题转化为HMT中对应的隐状态确定问题。2001年，韩国的Choi提出小波域隐马尔科夫树模型的图像分割方法。该方法中，HMT模型利用了小波的的多尺度特性，使系数之间的关系转变成隐状态的归属问题，并进行了在不同尺度上结合最大后验概率对图像进行分割。该方法的缺点是没有良好的区域一致性，且边缘不够光滑。2005年，中国的孙强提出一种改进的小波域HMT模型图像分割方法。其采用了一种新的上下文模型，结合了父节点邻域和子节点邻域的信息进行上下文的融合分割。该方法减少了小波产 生的杂块信息，增强了边缘信息，但该方法仅限于遥感图的分割，有一定的局限性。2005年，中国的沙宇恒等人提出了基于Contourlet域HMT模型的图像分割方法。该方法对纹理图，航拍图和SAR图像都有较好的分割结果，其图像的边缘信息保持较好，但是对均一区域分割效果欠佳。 

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于非下采样Contourlet变换的HMM图像分割方法，以得到具有较好的均一区域的纹理图像分割结果。 

实现本发明目的的技术方案是：利用非下采样Contourlet变换对图像处理的优点，将非下采样Contourlet变换系数通过高斯混合模型估计参数，对其按照一一对应的父子状态关系的隐马尔可夫树建模，利用训练得到的参数逐步对待分割图像进行最大似然函数的分类和最大后验概率的分类，最后结合马尔科夫随机场模型，充分考虑图像的空间相关性信息处理，得到最终的分割结果。具体实现过程如下： 

(1)输入待分割图像，从输入图像中截取N类具有均一区域的训练图像块，其中N为待分割图像的类别数； 

(2)对每类训练图像块分别进行非下采样Contourlet变换，得到多尺度的非下采样Contourlet变换系数； 

(3)采用期望最大化算法对每类训练图像块的非下采样Contourlet变换系数，按照一一对应的父子状态关系的隐马尔可夫树模型训练，得到隐马尔可夫模型参数Θ； 

(4)输入待分割的图像，对其进行非下采样Contourlet变换，根据模型参数Θ，计算待分割图像在各尺度下图像子块对应的似然值； 

(5)根据各个尺度似然值，采用最大化似然函数算法，得到图像各尺度上的初分割结果； 

(6)利用上下文模型对各个尺度的初分割结果进行修正，即对各尺度采用上下文context-5模型，通过最大化后验概率方法对该尺度初分割结果进行类别标记更新，得到各尺度的修正结果。 

(7)对步骤(6)得到的结果，结合隐马尔科夫随机场模型，基于待分割图像的相邻像素的空间相关性信息，不断更新分割的类别，取该分割结果为待分割图像最终的分割结果。 

本发明由于利用非下采样Contourlet变换的多尺度，多方向和平移不变的特点，避免引起吉布斯现象，弥补了小波域隐马尔科夫树分割方法对线性奇异的非稀疏的缺点；同时由于本发明基于非下采样Contourlet变换按照一一对应的父子状态关系的HMT模型建模，采用上下文模型对各个尺度的结果进行融合，并且结合马尔科夫随机场模型，充分挖掘相邻 像素的空间相关信息，克服了Contourlet域隐马尔可夫树分割方法分割结果有较差的边缘保持性和区域一致性的缺点。仿真结果表明，本发明方法较小波域HMT模型分割方法和基于Contourlet域HMT模型分割方法，对于纹理图像分割效果有明显改进。 

附图说明

图1是本发明的流程示意图； 

图2是本发明使用的一一对应的父子状态关系的隐马尔可夫树模型示意图； 

图3是本发明使用的上下文模型context-5示意图； 

图4是本发明用于第一种两类合成纹理图像的仿真分割结果图 

图5是本发明用于第二种两类合成纹理图像的仿真分割结果图； 

图6是本发明用于第一种三类合成纹理图像的仿真分割结果图； 

图7本发明用于第二种三类合成纹理图像的仿真分割结果图； 

图8是本发明用于第三种三类合成纹理图像的仿真分割结果图。 

具体实施方式

参照图1，本发明的具体实现过程如下： 

步骤1，输入待分割图像，从输入图像中选取N类具有均一区域的训练图像块，其中N为待分割图像的类别数，每幅训练图像块大小为128×128。 

步骤2，对每类训练图像分别进行非下采样Contourlet变换，采用‘maxflat’塔型滤波器和‘diamond maxflat’方向滤波器进行3层变换，每层8个方向子带，得到多尺度的非下采样Contourlet变换系数C_i。 

步骤3，采用期望最大化算法，对每类训练图像的非下采样Contourlet变换系数按照一一对应的父子状态关系的隐马尔可夫树模型进行训练，得到隐马尔可夫模型参数Θ；如图2所示，非下采样Contourle变换隐马尔可夫树模型中，两个节点场分别描述了在尺度J和J—1的系数，每个父隐状态变量连接一个子隐状态变量，相邻尺度之间的隐状态变量是一一对应的关系。 

步骤4，输入待分割的图像，对其进行非下采样Contourlet变换，根据隐马尔可夫模型参数Θ，求出图像在各尺度下图像子块对应的似然值，具体方法如下： 

假设图像的非下采样Contourlet变换的系数子带独立，用Θ¹，Θ²，...，Θ^k分别表示k子带的参数向量，有M＝{Θ¹，Θ²，...，Θ^k}，计算： 

$f\left(C_i\right|M)=f\left(T_i^1\right|{\mathrm{\Θ}}^1)f\left(T_i^2\right|{\mathrm{\Θ}}^2)...f\left(T_i^k\right|{\mathrm{\Θ}}^k)$

其中，相似度

的计算式为： 

$f\left(T_i\right|\mathrm{\Θ})=\underset{m=S,L}{\mathrm{\Σ}}{\mathrm{\β}}_{s_i}\left(m\right)p(S_i=m|\mathrm{\Θ})$

其中，β_i(m)＝f(T_i|S_i＝m，Θ)是由一个向上的树扫描算法得来的条件似然函数，p(S_i＝m|Θ)是模型训练时得到的状态概率，T_i表示系数子带对应的一个子树，S_i为非下采 样Contourlet变换系数C_i的高斯混合模型中的一个离散随机状态变量，它有两个状态  $p_{S_i}\left(m\right),m=S,L,$ 其中S，L分布表示小状态和大状态。 

步骤5，根据各个尺度似然值，采用最大化似然值算法，计算最大似然值  $c_i^{\mathrm{ML}}=\arg \underset{c\∈\{\mathrm{1,2},...,N\}}{\max }f\left(C_i\right|M),$ 由

得到图像的各尺度上的初分割结果。其中，f(C_i|M)为每个系数的似然值。 

步骤6，对初分割的结果的各尺度，根据上下文context-5模型进行类别标记更新，得到修正结果，具体步骤如下： 

(6a)选择以子节点为中心的3×3邻域内的主类标为背景向量V₁，组成上下文模型context-5，如图3所示，计算子节点的背景向量V； 

(6b)根据节点的背景向量V和初分割时得到的似然值，计算条件后验概率： 

$p(c_i=m|d_i^j,v_i^j)=\frac{e_{j,m}{\mathrm{\α}}_{j,v_i,m}f\left(d_i^j\right|c_i=m)}{\underset{l=1}{\overset{N_c}{\mathrm{\Σ}}}{e}_{j,l}{\mathrm{\α}}_{j,v_i,l}f\left(d_i^j\right|c_i=l)}$

其中，

表示尺度j上位置i处的特征系数，

表示尺度j上位置i处的背景向量，c_i为位置i上的像素点所属的类标，  $e_{j,m}=p_{c_i}\left(m\right)$ 表示尺度j上取类标为m的概率，α_j，vi，m＝p(v_i＝v_k|c_i＝m)表示尺度j上背景向量为v_i下取类标为m的概率； 

(6c)利用下式更新参数e_j，m和

的值： 

$e_{j,m}=\frac{1}{2^{2j}}\underset{i}{\mathrm{\Σ}}p(c_i=m|v_i^j,d_i^j)$

${\mathrm{\α}}_{j,v_k,m}=\frac{1}{2^{2j}\·e_{j,m}}\underset{\mathrm{iwith}{v}_i^j=v_k}{\mathrm{\Σ}}p(c_i=m|v_i^j,d_i^j),$ 对每个v_k，k∈{1，...，N_v}。 

其中，v_k为背景向量V的确切取值，对于N类纹理v_k有N²种取值，

为6(b)中公式计算的条件后验概率； 

(6d)重复步骤(6a)～(6b)，直到满足迭代停止条件，即达到允许误差为止，得到修正的分割结果。 

步骤7，对步骤(6)得到的结果，结合隐马尔科夫随机场模型，基于待分割图像的相邻像素的空间灰度相关信息，不断更新分割的类别，取该分割结果为待分割图像最终的分割结果，具体方法如下：

(7a)根据待分割图像和步骤(6)的修正结果建立马尔可夫随机场模型，即假设待分割图像Y中每一个像素点的强度y_i，都对应于标记场X中的一个标号x_i，根据贝叶斯定理，给定待分割图像Y时，标记场X的后验概率分布表示为： 

$P\left(X\right|Y)=\frac{P\left(Y\right|X)P\left(X\right)}{P\left(Y\right)}$

其中P(X)是先验吉布斯分布，P(Y|X)是给定标记场X时图像Y的条件概率，P(Y)是一个未知常数； 

(7b)建立待分割图像的能量函数，即假设每个像素灰度服从高斯分布，图像中每个像素的能量函数的表达式为： 

$E=\ln \left|{\mathrm{\Σ}}_k\right|+\frac{{(y_i-{\mathrm{\μ}}_k)}^2}{2{\mathrm{\Σ}}_k}+2\mathrm{m\β}$

其中，μ_k和∑_k分别表示图像中属于第k类的灰度均值向量和协方差矩阵，m表示像素y_i的邻域中属于不同于k类的像素的个数，β是一个权值，表示邻域中相邻像素间相互关系的重要性； 

(7c)按照(7b)步骤中所述的能量函数表达式，使用迭代条件模式算法寻找能量函数的最小解，当满足允许误差时停止，将标记场X的结果作为最终的图像分割结果。 

本发明的效果可以通过以下仿真结果进一步说明。 

1.仿真内容 

应用本发明方法与基于小波的隐马尔可夫树模型的图像分割方法W-HMT与基于Contourlet的隐马尔可夫树图像分割方法C-HMT，分别对五幅合成纹理图像，包括两幅两类纹理及三幅三类纹理的图像进行分割，并给出合成纹理图像的错分率作为客观评价标准，错分率定义为图像总错分像素点个数与图像总像素点个数的百分比。小波变换采用‘9-7’滤波器，Contourlet变换采用‘9-7’塔型和方向滤波器组。 

2.仿真结果 

图4和图5为本发明方法应用于含有两类合成纹理图像的仿真分割结果。其中，图4(a)和图5(a)为两类纹理的合成纹理图像，图4(b)和图5(b)分别为已有W-HMT方法对图4(a)和图5(b)的分割结果，图4(c)和图5(c)分别为已有C-HMT方法对图4(a)和图5(a)的分割结果，图4(d)和图5(d)分别为本发明方法对图4(a)和图5(a)的分割结果。图6、图7和图8为本发明方法应用于含有三类纹理的合成纹理图像的仿真分割结果。其中，图6(a)、图7(a)和图 8(a)分别为三类纹理的合成纹理图像，图6(b)、图7(b)和图8(b)分别为已有W-HMT方法对图6(a)、图7(a)和图8(a)的分割结果，图6(c)、图7(c)和图8(c)分别为已有C-HMT方法对图6(a)、图7(a)和图8(a)的分割结果，图6(d)、图7(d)和图8(d)分别为本发明方法对图6(a)、图7(a)和图8(a)的分割结果。可以看出，图4(c)、图5(c)、图6(b)和图6(c)分割结果中有明显的锯齿状边缘，而图4(d)、图5(d)和图6(d)中的边缘较光滑。图6(b)、图6(c)、图7(b)、图7(c)、图8(b)和图8(c)中存在大小不一的错分割杂块，而图6(d)、图7(d)和图8(d)中只有很少的错分块，分割的区域比较一致。本发明方法由于利用了非下采样Contourlet变换良好的几何逼近特性，结合了隐马尔可夫模型的优点，在合成纹理图像的分割上，较小波域HMT模型和基于Contourlet域HMT模型的分割方法，使分割结果区域一致性有明显改善，同时分割结果的边缘也有所改善。 

本发明方法与W-HMT，C-HMT方法对合成纹理图像错分率的对比如表1所示。 

表1　W-HMT，C-HMT方法和本发明方法分割结果的错分率比较 

从表1的错分率可以看出，本发明方法的错分率小于现有的W-HMT与C-HMT方法。 

综上，本发明利用非下采样Contoulet变换的特点，弥补了小波域隐马尔科夫树分割方法对线性奇异的非稀疏的缺点，以及Contourlet域隐马尔可夫树模型分割方法区域一致性差的弊端。利用其良好的高维奇异性的几何逼近特性以及根据非下采样Contoulet变换系数的分布特性，选择适当的上下文信息模型，以及结合马尔可夫随机场模型，充分考虑多尺度间、同一尺度相邻子带间系数，以及相邻像素的关系，较好的保留了边缘信息并减少分割结果中的杂块效应，本发明的方法对于纹理图像的处理效果有明显改进。

Claims (3)

1.一种基于非下采样Contourlet变换的HMT图像分割方法，包括如下步骤：

(1)输入待分割图像Y，从输入的待分割图像Y中截取N类具有均一区域的训练图像块，其中N为待分割图像Y的类别数；

(2)对每类训练图像块分别进行非下采样Contourlet变换，得到多尺度的非下采样Contourlet变换系数；

(3)采用期望最大化算法对每类训练图像块的非下采样Contourlet变换系数，按照一一对应的父子状态关系的隐马尔可夫树模型进行训练，得到隐马尔可夫模型参数Θ；

(4)根据隐马尔可夫模型参数Θ，计算待分割图像Y在各尺度下图像子块对应的似然值：

假设非下采样Contourlet变换的系数子带独立，用Θ¹，Θ²，...，Θ^k分别表示k个子带的参数向量，有M＝{Θ¹，Θ²，...，Θ^k}，计算：

$f\left(C_i\right|M)=f\left({T_i}^1\right|{\mathrm{\Θ}}^1)f\left({T_i}^2\right|{\mathrm{\Θ}}^2)...f\left({T_i}^k\right|{\mathrm{\Θ}}^k),$

其中，相似度f(T_i|Θ)的计算式为：

$f\left(T_i\right|\mathrm{\Θ})=\underset{m=\mathrm{Sma},\mathrm{Lar}}{\mathrm{\Σ}}{\mathrm{\β}}_i\left(m\right)p(S_i=m|\mathrm{\Θ}),$

其中，β_i(m)＝f(T_i|p(S_i＝m|Θ)是由一个向上的树扫描算法得来的条件似然函数，p(S_i＝m|Θ)是模型训练时得到的状态概率，T_i表示系数子带对应的一个子树，S_i为非下采样Contourlet变换系数C_i的高斯混合模型中的一个离散随机状态变量，它有两个状态概率p(S_i＝m|Θ)，m＝Sma，Lar，其中Sma，Lar分布表示小状态和大状态；

(5)根据各个尺度下图像子块对应的似然值，采用最大化似然值的算法，得到各尺度的初分割结果；

(6)利用上下文模型context-5对各尺度的初分割结果进行修正，即对各尺度采用上下文模型context-5，通过最大化后验概率法对各尺度初分割结果更新类别标记，得到各尺度的修正结果：

(6a)选择以子节点为中心的3×3邻域内的主类别标记为背景向量V₁，组成上下文模型context-5，计算子节点的背景向量V；

(6b)根据子节点的背景向量V和初分割时得到的似然值，计算条件后验概率；

$p(c_i=n|v_i^j,d_i^j)=\frac{e_{j,n}{\mathrm{\α}}_{j,v_i,n}f\left(d_i^j\right|c_i=n)}{\underset{l=1}{\overset{N_c}{\mathrm{\Σ}}}{e}_{j,l}{\mathrm{\α}}_{j,v_i,l}f\left(d_i^j\right|c_i=l)},$

其中，

表示尺度j上位置i处的特征系数，表示尺度j上位置i处的背景向量，c_i为位置i上的像素点所属的类别标记，e_j，n表示尺度j上取类别标记为n的概率，表示尺度j上背景向量为v_i下取类别标记为n的概率；

(6c)利用下式更新参数e_j，n和

的值：

$e_{j,n}=\frac{1}{2^{2j}}\underset{i}{\mathrm{\Σ}}p(c_i=n|d_i^j,v_i^j),$

${\mathrm{\α}}_{j,v_k,n}=\frac{1}{2^{2j}\·e_{j,n}}\underset{\mathrm{iwith}{v}_i^j=v_k}{\mathrm{\Σ}}p(c_i=n|d_i^j,v_i^j),$

其中，v_k为背景向量V的确切取值，对于N类纹理，v_k有N²种取值，

为(6b)中的条件后验概率；

(6d)重复步骤(6a)～(6b)，直到满足迭代停止条件，即达到允许误差为止，得到修正的分割结果；

(7)对步骤(6)得到的各尺度的修正结果，结合隐马尔可夫随机场模型，基于待分割图像Y的相邻像素的空间灰度相关信息，不断更新分割的类别，取分割结果为待分割图像Y最终的分割结果：

(7a)根据待分割图像Y和步骤(6)的修正结果建立隐马尔可夫随机场模型，即假设待分割图像Y中每个像素的强度y_i，都对应于标记场X中的一个标号x_i，根据贝叶斯定理，给定待分割图像Y时，标记场X的后验概率分布表示为：

$P\left(X\right|Y)=\frac{P\left(Y\right|X)P\left(X\right)}{P\left(Y\right)},$

其中P(X)是先验吉布斯分布，P(Y|X)是给定标记场X时待分割图像Y的条件概率，P(Y)是一个未知常数；

(7b)建立待分割图像Y的能量函数，即假设每个像素灰度服从高斯分布，待分割图像Y中每个像素的能量函数的表达式为：

$E=\ln \left|{\mathrm{\Σ}}_n\right|+\frac{{(y_i-{\mathrm{\μ}}_n)}^2}{2{\mathrm{\Σ}}_n}+2m^{\′}\mathrm{\γ},$

其中，μ_n和∑_n分别表示待分割图像Y中属于第n类的灰度均值向量和协方差矩阵，m’表示像素y_i的邻域中属于不同于第n类的像素的个数，γ是一个权值，表示邻域中相邻像素间相互关系的重要性；

(7c)按照(7b)步骤中所述的能量函数表达式，使用迭代条件模式算法寻找能量函数的最小解，当满足允许误差时停止，将标记场X的结果作为最终的图像分割结果。

2.根据权利要求1所述的基于非下采样Contourlet变换的HMT图像分割方法，其中步骤(2)所述的对每类训练图像块分别进行非下采样Contourlet变换，是采用‘maxflat’塔型滤波器和‘diamond maxflat’方向滤波器对每类训练图像块分别进行3层变换，得到每层8个方向子带的非下采样Contourlet变换系数。

3.根据权利要求1所述的基于非下采样Contourlet变换的HMT图像分割方法，其中步骤(3)所述的一一对应的父子状态关系隐马尔可夫树模型，是非下采样Contourlet变换系数的每个父隐状态变量连接一个子隐状态变量，该子隐状态变量在相邻尺度之间为一一对应的父子关系。

Images