CN101425184A - 基于第二代Bandelet域隐马尔科夫树模型的图像分割方法 - Google Patents

Abstract

本发明公开了一种基于第二代Bandelet域隐马尔科夫树模型的图像分割方法，主要解决图像分割质量差的问题。其实现过程为：输入待分割图像并提取训练图像块；对输入图像和训练图像进行第二代Bandelet变换；采用期望最大化EM算法对训练图像进行Bandelet域隐马尔科夫树模型训练，得到模型参数；利用模型参数计算输入图像的多尺度似然值；利用高斯混合模型计算输入图像的像素级似然值；利用最大似然准则估计得到图像的初始分割结果；采用两种上下文模型对初始分割结果进行多尺度融合，得到后发分割结果；选择后分割结果中的像素级分割结果作为图像最终分割结果输出。本发明具有分割图像边缘清晰和区域一致性好的优点，可用于合成纹理图像、SAR图像和航拍图像的分割。

Description

基于第二代Bandelet域隐马尔科夫树模型的图像分割方法

技术领域

本发明属于图像处理技术领域，涉及该技术在图像分割领域的应用，具体地说是一种基于第二代Bandelet域隐马尔科夫树模型的图像分割方法。该方法可应用于合成纹理图像、SAR图像和航拍图像的分割中。

背景技术

图像分割是数字图像处理技术中的关键技术之一。它是将图像中有意义的特征部分提取出来，其有意义的特征包括图像的边缘、纹理等，这是进一步进行图像识别、分析和理解的基础。图像分割是一种重要的图像技术，在理论研究和实际应用中都得到人们的广泛的重视。图像分割的方法和种类很多，没有标准的分割方法，许多不同种类的图像或景物都可以作为待分割的图像数据，同时，某些分割方法也只适用于对某些特殊类型的图像进行分割，而分割结果的好坏需要根据具体的场合及要求去评价。因此，图像分割仍是目前研究的热点之一。

图象分割算法一般是基于亮度值的两个基本特性之一：不连续性和相似性。早期的图像研究中，分割方法主要可分为两大类，即基于边界的分割方法和基于区域的分割方法，而后来兴起的基于统计域的图像分割方法得到了更多的关注。

基于统计域模型的图象分割方法按变换域进行划分，到目前为止最典型的有这三类：基于小波域的方法、基于复小波域的方法以及基于contourlet域的方法，不同域的特征表征不同的图象纹理信息。按训练模型的选择，则可分为独立高斯混合IMM模型、隐马尔科夫树HMT模型、子带绑定的隐马尔科夫树HMT-3S模型以及矢量隐马尔科夫模型，不同的模型选择决定变换域特征之间的某些固定关系，如选用HMT模型，则决定了多尺度变换域上每一层的特征应该对应相同，HMT-3S则不仅要求每一层上对应的特征相同而且决定要求不同子带特征之间应该具有某联系。相比较之下，IMM认为每一子带每一尺度之间均是独立不相关的，其要求最少，因此可塑性越强。对于该领域后融合方面的方法，首先韩国的H.choi等人提出一种基于背景的贝叶斯分割方法，参见H.Choi，R.G.Baraniuk.Multiscale Image Segmentation Using Wavelet-Domain Hidden Markov Models.IEEE Transactions on Image Processing，2001，10(9)：1309-1321。之后中国的Fan G.L等人总结各种背景模型并提出了一种JMCMS融合方法，融合分割，参见Fan G.L.，Xia X.G..Ajoint multi-context and multi-scale approach to Bayesian image segmentation.IEEETransactions on Geoscience and Remote Sensing，2001，39(12)：2680-2688。值得注意的是这些方法中所涉及到的上下文背景模型，均只利用较粗尺度以及本层分割尺度上的信息，而不曾考虑其下一较细尺度上更能表征边缘的信息。

基于小波域隐马尔科夫树模型的图像分割WD-HMTseg和基于contourlet域隐马尔科夫树模型的图像分割CHMTseg是目前的研究热点之一。但由于小波对非线性目标的弱稀疏性，使得WD-HMTseg的分割结果的边缘不太准确且有较多杂块产生；而contourlet变换加入了更多方向信息，反而使得CHMTseg结果的区域一致性并不是很令人满意。因此，需要找到一种在分割结果的边缘准确性和区域一致性上都有较好性能的图像分割方法。

发明内容

本发明的目的在于为了克服上述已有技术分割边缘不平滑和区域一致性较差的不足，提出了一种基于第二代Bandelet域隐马尔科夫树模型的图像分割BHMTseg方法，以提高图像分割的质量。

本发明的技术方案是：将第二代Bandelet变换与隐马尔科夫树模型结合对图像进行初始分割，然后采用上下文模型context-1和context-2进行后融合，实现对合成纹理图像、航拍图像和SAR图像的分割。其实现过程如下：

(1)输入待分割图像，并截取具有相同纹理区域的子块作为训练图像；

(2)采用第二代Bandelet变换对输入图像及训练图像进行变换，得到Bandelet系数；

(3)采用期望最大化EM算法对训练图像进行Bandelet域隐马尔科夫树模型训练，得到模型参数；

(4)利用模型参数采用隐马尔科夫树模型的似然算法；

(5)应用高斯混合模型计算输入图像的像素级似然值；

(6)利用最大似然法则估计输入图像每个子块的似然值和像素级似然值，得到输入图像的初始分割结果；

(7)对输入图像的初始分割结果先进行第一次多尺度融合，得到输入图像的第一次分割结果；再对输入图像的第一次分割结果进行第二次多尺度融合，得到输入图像各尺度的后分割结果，并将该分割结果中的像素级分割结果作为输入图像最终输出的分割结果。

本发明与现有的技术相比具有以下优点：

1、由于本发明运用了第二代Bandelet变换，因而在表示图像时能充分利用图像几何正则性的特点，提高了分割结果中边缘的准确性；

2、由于本发明充分考虑了第二代Bandelet系数的稀疏性，并对其运用两种不同的隐马尔科夫树上下文模型进行多尺度融合，提高了分割结果的区域一致性；

3、由于本发明充分利用了Bandelet系数的鲁棒性和隐马尔科夫树模型概率计算的稳定性，提高了分割结果的稳定性。

仿真实验表明，本发明提高了合成纹理图像、SAR图像和遥感图像分割结果的边缘准确性，增强了分割结果的区域一致性，并提高了分割结果的稳定性。

附图说明

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

图2是本发明多尺度融合采用的context-1和context-2两种上下文模型示意图；

图3(a)是本发明采用的一幅具有两类纹理的合成纹理图像；

图3(b)是对采用基于小波域隐马尔科夫树模型对图3(a)进行分割的最终结果图；

图3(c)是采用本发明对图3(a)进行分割的最终结果图；

图3(d)是本发明采用的一幅具有三类纹理的合成纹理图像；

图3(e)是采用基于小波域隐马尔科夫树模型对图3(d)进行分割的最终结果图；

图3(f)是采用本发明对图3(d)进行分割的最终结果图；

图3(g)是本发明采用的一幅具有四类纹理的合成纹理图像；

图3(h)是采用基于小波域隐马尔科夫树模型对图3(g)进行分割的最终结果图；

图3(i)是采用本发明对图3(d)进行分割的最终结果图；

图4(a)是本发明采用的一幅航拍图；

图4(b)是采用基于小波域隐马尔科夫树模型对图4(a)进行分割的最终结果图；

图4(c)是采用基于contourlet域隐马尔科夫树模型对图4(a)进行分割的最终结果图；

图4(d)是采用本发明对图4(a)进行分割的最终结果图；

图4(e)是本发明采用的一幅SAR图像；

图4(f)是采用基于小波域隐马尔科夫树模型对图4(e)进行分割的最终结果图；

图4(g)是采用基于contourlet域隐马尔科夫树模型对图4(e)进行分割的最终结果图；

图4(h)是采用本发明对图4(e)进行分割的最终结果图；

图4(i)是本发明采用的一幅SAR图像；

图4(j)是采用基于小波域隐马尔科夫树模型对图4(i)进行分割的最终结果图；

图4(k)是采用基于contourlet域隐马尔科夫树模型对图4(i)进行分割的最终结果图；

图4(l)是采用本发明对图4(i)进行分割的最终结果图。

具体实施方式

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

步骤一、输入待分割图像，并截取具有相同纹理区域的子块作为训练图像。

输入待分割图像图3(a)、图3(d)、图3(g)、图4(a)、图4(e)和图4(i)，分别对这些待分割图像截取具有相同纹理区域的子块作为训练图像块，截取子块的尺寸大小为64×64，截取的子块个数等于待分割图像要分割的纹理类别数。其中，图3(a)要截取的子块数为2，图3(d)要截取的子块数为3，图3(g)要截取的子块数为4，图4(a)、图4(e)和图4(i)要截取的子块数为2。

步骤二、采用Stéphane Mallat等提出的第二代Bandelet变换对输入图像进行第二代Bandelet变换，得到Bandelet系数。

实现该步骤的具体过程如下：

1)对输入图像进行二维小波变换；

2)二进剖分：在Bandelet变换的每个尺度上，连续地把一个区域剖分为四个相等大小的方形子区域，该方形区域的尺寸限制为宽度为L个像素且4≤L≤2^J/2，例如，对一幅2ⁿ×2ⁿ大小的图像而言，设定方形区域的最小尺寸为2^m×2^m，则二进剖分的结果是将小波分解后的系数图分成了2^2(n-m)个宽度为2^m×2^m的子块；

3)选择步骤2)中的每个子块，并对选择的子块进行方向采样。对N×N大小的子块可设定其采样方向数为4×N，每个方向对应于[0，π)区间的一个角度θ，步长为π/4N；

4)将子块区域沿每个采样方向θ进行正交投影。即将子块区域的每个点(i，j)这样处理：t＝-i*sinθ+j*cosθ；

5)对投影信号t进行重新排序，得到一个一维离散信号f_d；

6)对一维信号f_d进行一维离散小波变换，得到一维小波系数f_dR；

7)利用得到的一维离散信号f_d和一维小波系数f_dR建立Lagrangian目标函数，即

其中，T为门限，R_G是在使用熵编码器对几何特征参数d编码时所需的比特数，R_B表示对计算Bandelet系数的参数数目编码时所需的比特数，λ取3/28；

8)计算采样区域的每个采样方向对应的Lagrangian目标函数值，将使Lagrangian目标函数值最小的采样方向作为采样区域的最佳方向；

9)对有几何流的采样区域沿对应的最佳方向投影，并对投影信号重新排序得到一维离散信号，对该离散信号进行一维离散小波变换得到的一维小波系数作为Bandelet系数输出，对无几何流的区域则将二维小波变换得到的系数作为Bandelet系数输出。

步骤三、对训练图像在Bandelet域采用期望最大化EM算法训练隐马尔科夫树模型参数。

在已知分类数的图像中选取相同数目的对应于不同纹理的均一区域作为训练图像进行隐马尔科夫模型的训练，可得到模型参数向量Θ(，

μ_i，m，

)，其中，m＝S，L，表示节点的状态，S代表小状态，L代表大状态；

表示根节点S_i的初始状态概率； ${\mathrm{\ϵ}}_{i,\mathrm{\ρ}\left(i\right)}^{m,n}=p_{S_i|S_{\mathrm{\ρ}\left(i\right)}}$ (m|S_ρ(i)＝n)，表示父节点状态ρ(i)为n时，子节点状态i为m的状态转移概率；μ_i，m和

分别表示节点状态为m时的均值和方差。

步骤四、利用模型参数采用隐马尔科夫树模型的似然算法，计算在不同尺度下输入图像的每个子块的似然函数。

利用得到的模型参数计算不同尺度下每个图像子块的似然值，具体计算过程是：

给定一组二维Bandelet域隐马尔科夫模型参数M和测试图的第二代Bandelet系数b；计算Bandelet系数在根节点b的一个子带T的似然函数，即先假定这些子带的BHMT参数为Θ，由期望最大化EM算法中的一个向上的四叉树扫描的步骤得到条件似然函数为：

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

其中，p(S_i＝m|Θ)是由模型Θ得到的状态转移概率，

表示节点s_i的初始状态概率。

由于假设Bandelet系数各子带独立，则二值图像块d_i的Bandelet系数包括系数即每个子块有三个Bandelet子带，由独立的假设可以得到二值图像块在模型Θ下的概率f(d_i|M)，即每个子块的似然函数：

$f\left(d_i\right|M)=f\left(T_i^{\mathrm{LH}}\right|{\mathrm{\Θ}}^{\mathrm{LH}})f\left(T_i^{\mathrm{HH}}\right|{\mathrm{\Θ}}^{\mathrm{HH}})f\left(T_i^{\mathrm{HL}}\right|{\mathrm{\Θ}}^{\mathrm{HL}})---\left(3\right)$

式中，

表示LH子带在模型Θ下的似然函数，

表示HH子带在模型Θ下的似然函数，

表示HL子带在模型Θ下的似然函数。

步骤五、应用高斯混合模型计算输入图像的像素级似然值。

概率密度函数可以通过对应于每一纹理的像素直方图来建模，纹理可被近似为一个高斯混合模型。对于训练纹理的像素，给其一个高斯混合模型来建模像素灰度值，并通过该高斯混合模型计算每一像素的相似度，再通过比较相似度大小来获得图像像素级似然值。

步骤六、利用最大似然法则估计输入图像每个子块的似然值和像素级似然值，得到输入图像的初始分割结果。

选择(3)式中的最大似然值中的最大值对应的纹理类别作为子带的初始类别分割结果。其中，最大似然值的计算式是：

$c_i^{\mathrm{ML}}=\arg \underset{c\∈\{\mathrm{1,2},...,N_c\}}{\max }f\left(d_i\right|M_c)---\left(4\right)$

式中c∈{1，2，...，N_c}，是每个像素的类别标号，N_c代表输入图像总的分割类别数。

步骤七、对输入图像的初始分割结果进行多尺度融合，得到输入图像的最终分割结果。

首先，对输入图像的初始分割结果利用context-1模型进行第一次多尺度融合，得到输入图像的第一次分割结果。该context-1模型如图2(a)所示，它是由父类标及其邻域的主类标和父类标组成的上下文向量；

然后，对输入图像的第一次分割结果利用context-2模型进行多尺度融合，得到输入图像各尺度的后分割结果。该context-2模型如图2(b)所示，它是由父类标、父类标及其邻域的主类标和子邻域的主类标组成的一个上下文向量。

最后，将后分割结果中的像素级分割结果作为输入图像的最终分割结果输出。具体实现过程如下：

1)计算在context-1模型下，上下文矢量为

时，尺度j处第i个图像块

属于类别c_i的概率 $p\left(c_i\right|d_i^j,v_i^j),$ 即

$p(c_i=m|d_i^j,v_i^j)=\frac{p_{c_i}\left(m\right)p(v_i={\stackrel{\‾}{v}}_k|c_i=m)f\left(d_i^j\right|c_i=m)}{\underset{l=1}{\overset{N_c}{\mathrm{\Σ}}}{p}_{c_i}\left(l\right)p(v_i={\stackrel{\‾}{v}}_k|c_i=l)f\left(d_i^j\right|c_i=l)}$

式中，

的上下文信息用上下文矢量v_i来代表，表示尺度j处的上下文矢量，j＝0，...，log₂(N)-1，N表示图像块的宽度；

表示尺度j处的第i个图像块，为多尺度似然计算得出的类别为c_i时的似然值；m＝S，L，代表节点的状态，S代表小状态，L代表大状态；是节点类别为c_i时的初始概率，N_c代表总的类别数，p(v_i|c_i)表示类别为c_i时v_i的概率；

2)根据得到的 $p(c_i=m|d_i^j,v_i^j)$ 计算图像块

的最大后验概率MAP，得到各节点的初始类标，即输入图像的第一次分割结果；

$c_i^{\mathrm{MAP}}=\underset{c_i\∈\{1,...,N_c\}}{\arg \max }p\left(c_i\right|d_i^j,v_i^j);$

3)利用context-2模型对第一次分割结果再次计算在上下文向量下二值图像块

量；

是节点类别为c_i时的初始概率；m＝S，L，代表节点的状态，S代表小状态，L代表大状态；

表示节点类别为c_i时

的概率；

代表类别为c_i时像素块

的似然值；

4)根据得到的

再次计算图像块

的最大后验概率，即 $c_i^{\′\mathrm{MAP}}=\underset{c_{i\′}\∈\{1,...,N_c\}}{\arg \max }p\′\left(c_{i\′}\right|d_i^{\′j},{v_i^{\′}}^j),$ 得到输入图像的各尺度后分割结果；

5)将后分割结果中的像素级分割结果作为输入图像的最终分割结果输出，如图3(c)、图3(f)、图3(i)、图4(d)、图4(h)和图4(1)所示。

图3(b)、图3(e)和图3(h)分别为图3(a)、图3(d)和图3(g)的基于小波域隐马尔科夫树模型分割方法的的分割结果，通过将图3(c)与3(b)、图3(f)与图3(e)，以及图3(i)与图3(h)比较，看出，本发明的分割方法使得输入图像的分割边界更准确，分割的区域一致性更好。表1给出了基于小波域隐马尔科夫树模型的分割WD-HMTseg方法和本发明分割BHMTseg方法的错分率数据，这些错分率数据也表明，本发明明显优于基于小波域隐马尔科夫树模型的分割方法。

表1 WD-HMTseg和BHMTseg分割结果的错分率Pe

图4(b)，图4(f)和图4(j)分别为图4(a)，图4(e)和图4(i)的基于小波域隐马尔科夫树模型的仿真分割结果；图4(c)，图4(g)和图4(k)分别为图4(a)，图4(e)和图4(i)的基于contourlet域隐马尔科夫树模型的仿真分割结果；图4(d)，图4(h)和图4(1)分别为对图4(a)，图4(e)和图4(i)采用本发明进行分割的仿真结果。通过分别比较4(d)与4(b)和图4(c)、图4(h)与图4(f)和图4(g)以及图4(1)与图4(j)和图4(k)，可以看出本发明的分割结果边缘更加准确，区域一致性更强。

仿真实验表明，本发明能比现有的基于小波域隐马尔科夫树模型的分割方法和基于contourlet域隐马尔科夫树模型的分割方法得出的分割结果更稳定，分割质量更高。

Claims (3)

1、一种基于第二代Bandelet域隐马尔科夫树模型的图像分割方法，包括如下过程：

(1)输入待分割图像，并截取具有相同纹理区域的子块作为训练图像；

(2)采用第二代Bandelet变换对输入图像及训练图像进行变换，得到Bandelet系数；

(3)采用期望最大化EM算法对训练图像进行Bandelet域隐马尔科夫树模型训练，得到模型参数；

(4)利用模型参数采用隐马尔科夫树模型的似然算法，计算在不同尺度下输入图像的每个子块的似然函数值；

(5)采用高斯混合模型计算输入图像的像素级似然值；

(6)对输入图像每个子块的似然值和像素级似然值应用最大似然准则，得到输入图像的初始分割结果；

(7)对输入图像的初始分割结果先进行第一次多尺度融合，得到输入图像的第一次分割结果；再对输入图像的第一次分割结果进行第二次多尺度融合，得到输入图像各尺度的后分割结果，并将该分割结果中的像素级分割结果作为输入图像最终输出的分割结果。

2、根据权利要求1所述的图像分割方法，其中过程(2)所述的用第二代Bandelet变换对输入图像及训练图像进行变换，按如下过程进行：

(2a)对输入图像进行二维小波变换；

(2b)对二维小波变换后的输入图像进行二进剖分，即在Bandelet变换的每个尺度上，连续地把一个区域剖分为四个大小相等的方形子块区域，并将每个方形子块区域的尺寸限制为宽度为L个像素，且4≤L≤2^j/2，j为输入图像的宽度；

(2c)对每个方形子块区域进行方向采样，并将子块区域沿每个采样方向进行正交投影，得到投影信号；

(2d)对投影信号进行重新排序，得到一个一维离散信号f_d；

(2e)对一维离散信号进行一维离散小波变换，得到一维小波系数f_dR；

(2f)利用得到的一维离散信号f_d和一维小波系数f_dR建立Lagrangian目标函数，即

其中，T为门限，R_G是在使用熵编码器对几何特征参数d编码时所需的比特数，R_B表示对计算Bandelet系数的参数数目编码时所需的比特数，λ取3/28；

(2g)计算采样区域的每个采样方向对应的Lagrangian目标函数值，将使Lagrangian目标函数值最小的采样方向作为采样区域的最佳方向；

(2h)对有几何流的采样区域沿对应的最佳方向投影，并对投影信号重新排序得到一维离散信号，对该离散信号进行一维离散小波变换得到的一维小波系数作为Bandelet系数输出，对无几何流的区域则将二维小波变换得到的系数作为Bandelet系数输出。

3、根据权利要求1所述的图像分割方法，其中步骤(7)按如下过程进行：

(3a)计算在context-1模型下，上下文矢量为

时，尺度j处第i个图像块

属于类别c_i的概率

，即

$p(c_i=m|d_i^j,v_i^j)=\frac{p_{c_i}\left(m\right)p(v_i={\stackrel{\‾}{v}}_k|c_i=m)f\left(d_i^j\right|c_i=m)}{\underset{l=1}{\overset{N_c}{\mathrm{\Σ}}}{p}_{c_i}\left(l\right)p(v_i={\stackrel{\‾}{v}}_k|c_i=l)f\left(d_i^j\right|c_i=l)}$

式中，的上下文信息用上下文矢量v_i来代表，

表示尺度j处的上下文矢量，j＝0，...，log₂(N)-1，N表示图像块的宽度；

为多尺度似然计算得出的类别为c_i时

的似然值；m＝S，L，代表节点的状态，S代表小状态，L代表大状态；是节点类别为c_i时的初始概率，N_c代表总的类别数，p(v_i|c_i)表示类别为c_i时v_i的概率；

(3b)根据得到的 $p=\left(c_i\right|d_i^j,v_i^j)$ 计算图像块

的最大后验概率MAP，得到各节点的初始类标，即输入图像的第一次分割结果；

$c_i^{\mathrm{MAP}}=\underset{c_i\∈\{1,...,N_c\}}{\arg \max }p\left(c_i\right|d_i^j,v_i^j);$

(3c)利用context-2模型对第一次分割结果再次计算在上下文向量

下二值图像块属于类别c_i的概率

，即：

$p\′(c_{i\′}=m|d_i^{\′j},v_i^{\′j})=\frac{{p\′}_{c_{i\′}}\left(m\right)p(v_i^{\′}={\stackrel{\‾}{v}}_k|c_{i\′}=m)f\left(d_i^{\′j}\right|c_{i\′}=m)}{\underset{l=1}{\overset{N_c}{\mathrm{\Σ}}}{p\′}_{c_{i\′}}\left(l\right)p(v_i^{\′}={\stackrel{\‾}{v}}_k|c_{i\′}=l)f\left(d_i^{\′j}\right|c_{i\′}=l)};$

式中，

表示第一次分割结果尺度j处的第i个图像块，

代表

对应的上下文矢量；是节点类别为c_i时的初始概率；m＝S，L，代表节点的状态，S代表小状态，L代表大状态；

表示节点类别为c_i时

的概率；代表类别为c_i时像素块

的似然值；

(3d)根据得到的

再次计算图像块

的最大后验概率，即 $c_i^{\′\mathrm{MAP}}=\underset{c_{i\′}\∈\{1,...,N_c\}}{\arg \max }p\′\left(c_{i\′}\right|d_i^{\′j},v_i^{\′j}),$ 得到输入图像的各尺度后分割结果；

(3e)将后分割结果中的像素级分割结果作为输入图像的最终分割结果输出。

Images