CN106709928A - 一种快速的含噪图像二维最大类间方差阈值法 - Google Patents

Abstract

本发明涉及一种快速的含噪图像二维最大类间方差阈值法，首先求出噪音图像的灰度均值和灰度标准差；对图像每个像素采用3×3邻域的平均灰度值进行平滑，得到平滑图像；然后用最大类间方差阈值法计算平滑图像的类间方差，通过灰度均值和标准差缩小类间方差解的搜索空间，遍历搜索空间，记录使类间方差最大时的解为最优一维阈值T0；用二维最大类间方差法计算目标类与背景类的类间方差离散度矩阵的迹，通过T0和噪音图像的灰度标准差缩小迹的解的搜索空间，最后遍历解的搜索空间，记使离散度矩阵的迹最大时的灰度值二元组为最优二维切割阈值。本发明能够避免遍历全部灰度级，大大缩小计算量的同时也能得到准确的解。

Description

一种快速的含噪图像二维最大类间方差阈值法

技术领域

本发明属于数字图像处理领域，具体涉及图像阈值分割中的一种快速的含噪图像二维最大类间方差阈值法。

背景技术

图像分割是由图像处理进到图像分析的关键步骤，它指把图像分成各具特性的区域并提取出感兴趣的目标的技术和过程。通过把感兴趣的目标提取出来，更高层的分析和理解就有了基础。

现有的图像分割方法主要分以下几类：基于阈值的分割方法，基于区域的分割方法，基于边缘的分割方法以及基于特定理论的分割方法等。其中，基于阈值的图像分割方法是一种计算相对简单，运算效率相对较高的一种广泛应用的分割技术，利用图像中要提取的目标区域与其背景在灰度特性上的差异，把图像看作具有不同灰度级的两类区域(目标区域和背景区域)的组合，选取一个合理的阈值，以确定图像中每个像素点应该属于目标区域还是背景区域，从而产生相应的二值图像。阈值分割法实际上就是输入图像f到输出图像g的变换：

公式(1)中，T为阈值，f(m,n)表示输入图像f中像素点(m,n)处的灰度值，g(m,n)表示输出图像g中像素点(m,n)处的灰度值，对于目标的图像元素g(m,n)＝1，对于背景的图像元素g(m,n)＝0。

在众多基于阈值的分割方法中，最大类间方差法(简称Otsu法)以其分割效果好、性能稳定且适用范围广得到了广泛的应用。最大类间方差法是一种全局的自动非参数无监督的阈值选取算法，该方法以图像的一维直方图为依据，以目标和背景的类间方差最大为阈值选取准则，在很多情况下都能取得很好的阈值。但是，由于图像在采集或传输过程中容易受到噪声干扰，导致图像直方图没有明显的双峰或图像信噪比较低时，传统的最大类间方差法或其迭代算法很难获得满意的分割效果，其原因在于传统最大类间方差法是基于一维直方图的分割方法，它仅考虑图像的灰度信息而没有考虑像素的空间邻域信息。

针对传统最大类间方差法抗噪性差的不足,二维最大类间方差法(简称二维Otsu法)被提出，该方法建立了既能反映像素点的灰度分布又能体现像素点与其邻域空间相关信息的灰度均值二维直方图，最佳的阈值是在一个二维的类间方差测度准则取最大值时得到的一个二维矢量，并以此二维矢量作为分割门限进行图像分割，对噪音图像能获得满意的分割效果。但是，二维最大类间方差法在提高抗噪性的同时也大大增加了算法的计算搜索复杂性，不利于实时场合的广泛应用。

事实上，我们在对大量图片分别进行一维最大类间方差法和二维最大类间方差法求阈值时，二者得到的结果在灰度级上差距并不大，特别是在图像灰度方差较小时，二者仅相差几个灰度级。因为二维直方图仅仅是在一维直方图的基础上，将‘灰度’变成‘灰度-邻域灰度均值’的二元组，加入了邻域灰度均值后，二维最大类间方差法能有效消除噪音的影响。本发明提出一种快速的含噪图像二维最大类间方差法，利用图像的均值和方差特征降低最佳二维阈值的搜索范围，在大大降低计算复杂度的同时分割效果不受影响，提高图像分割的效率。

发明内容

本发明的目的在于提出一种快速的含噪图像二维最大类间方差阈值法，解决了原始二维最大类间方差阈值法计算复杂度过高的问题。

本发明的技术方案为：

一种快速的含噪图像二维最大类间方差阈值法，其特征在于，包含如下步骤：

步骤1，输入待分割含噪音图像F，用用f(i,j)表示原始图像在像素点(i,j)处的灰度值，M和N分别表示图像的高度和宽度，单位是像素，读取原始图像各个像素点的灰度值，统计图像每个灰度级L_k像素出现的频数H_k，这里k的取值范围是0到255；

步骤2，求出含噪图像F的灰度均值μ和灰度标准差D，公式为：

式中，M，N，L_k,H_k的定义同步骤1，k是循环变量；

步骤3，对图像每个像素采用设定邻域的平均灰度值进行平滑，得到平滑图像G；

步骤4,利用最大类间方差阈值法，将灰度级从到开始依次遍历获取平滑图像G的阈值T_G，具体包括：

步骤4.1,读入平滑图像G，统计图像灰度级a像素出现的频数f_a，设分割阈值为T_Gx，目标类像素点总数为n₀，背景类的像素点总数为n₁，目标类和背景类的点占图像比例分别为w₀,w₁，目标类和背景类的灰度均值分别为u₀，u₁，总的灰度均值为u_G，求解公式如下：

u_G＝w₀×u₀+w₁×u₁

上式中，a表示灰度级变量，f_a表示灰度值为a的像素点的个数，目标类与背景类的类间方差为G(a)：

G(a)＝w₀(u₀-u_G)²+w₁(u₁-u_G)²

步骤4.2,将a从开始依次遍历到记录当G(a)的值最大时对应的a为平滑图像G的最优阈值T_G，μ和D定义同步骤2；

步骤5,利用二维直方图P_ij和二维最大类间方差阈值法，在二维空间搜索最优阈值(S^*,T^*)，D为步骤2中求得的含噪图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，使得二维最大类间方差最大时对应的阈值为最优二维阈值(S^*,T^*)，具体包括；

步骤5.1统计噪音图像F的二维直方图P_ij，公式如下，

f_ij为像素灰度值为i且3×3邻域平均灰度值为j的像素点的个数，M和N的定义同步骤1；

步骤5.2,设分割阈值为(s,t)，将灰度值小于s且邻域灰度均值小于t的像素点划分为背景类；将灰度值大于s且邻域灰度均值大于t的像素点划分为目标类；统计目标类和背景类的像素点个数，计算目标类和背景类的灰度均值矢量μ₀，μ₁，全局的灰度均值矢量μ_T，求解公式为，

上式中，P_ij定义同步骤5.1,T表示转置矩阵，

计算目标类与背景类的离散度矩阵G_(s,t)的公式为：

P₀，P₁分别为目标类像素和背景类像素占总像素数的比例，如下式，

P₀＝∑_{(i＞s且j＞t)}P_ij P₁＝Σ_{(0＜i＜s且0＜j＜t)}P_ij

计算两类的离散度矩阵的迹tr(G_(s,t))公式为：

tr(G_(s,t))＝P₀[(μ_0i-μ_Ti)²+(μ_0j-μ_Tj)²]+P₁[(μ_1i-μ_Ti)²+(μ_1j-μ_Tj)²]

μ_0i，μ_0j，μ_1i，μ_1j，μ_Ti，μ_Tj是步骤5.2中灰度均值向量μ₀，μ₁，μ_T的分量；

定义tr(G_(s,t))的解的搜索空间为D为步骤2中求得的噪音图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，且s和t的取值范围均为0到255，将(s,t)从开始依次遍历到当tr(G_(s,t))的值最大时对应的(s,t)为最优二维阈值(S^*,T^*)；

步骤6，利用最佳阈值(S^*,T^*)，对原始图像进行分割并输出；

在上述的一种快速的含噪图像二维最大类间方差阈值法,所述步骤2求出含噪图像F的灰度均值μ和灰度标准差D的公式为，

式中，M，N，L_k,H_k的定义同步骤1，k是循环变量。

在上述的一种快速的含噪图像二维最大类间方差阈值法,所述步骤6利用最佳阈值(S^*,T^*)，对原始图像进行分割，按如下方式进行：将图像F灰度值大于S^*且设定邻域灰度均值大于T^*的像素点的灰度值重置为n，将其余像素点的灰度值重置为m，且n和m为0至255的任意整数，且n不等于m。

本发明具有如下优点：原始二维最大类间方差法需要计算L²个离散度矩阵的迹，每一次计算迹时，需要对s×t+(L-s)×(L-t)个点做累加运算，因此，原始二维最大类间方差法的计算复杂度为O(L⁴)，当L＝256时，需要进行约2³²次基本运算。本发明提出的快速方法，只需要计算个离散度矩阵的迹，且每一次计算迹时，只需要对个点做累加运算。大大缩小了二维最大类间方差法的计算工作量，提高了图像分割的效率。

附图说明

图1为本发明实施例的流程图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，需要说明的是最大类间方差法和二维最大类间方差法为已有方法，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明所采用的技术方案是：一种快速的含噪图像二维最大类间方差阈值法,其特征在于，包括以下步骤：

步骤1，输入待分割含噪音图像F，用用f(i,j)表示原始图像在像素点(i,j)处的灰度值，M和N分别表示图像的高度和宽度，单位是像素，读取原始图像各个像素点的灰度值，统计图像每个灰度级L_k像素出现的频数H_k，这里k的取值范围是0到255；

步骤2，用如下公式求出含噪图像F的灰度均值μ和灰度标准差D；

式中，M，N，L_k,H_k的定义同步骤1，k是循环变量。

步骤3，对图像每个像素采用3×3邻域的平均灰度值进行平滑，得到平滑图像G；步骤4,利用最大类间方差阈值法，将灰度级从到开始依次遍历获取平滑图像G的阈值T_G；

步骤4.1,读入平滑图像G，统计图像灰度级a像素出现的频数f_a，设分割阈值为T_Gx，目标类像素点总数为n₀，背景类的像素点总数为n₁，目标类和背景类的点占图像比例分别为w₀,w₁，目标类和背景类的灰度均值分别为u₀，u₁，总的灰度均值为u_G，求解公式如下：

u_G＝w₀×u₀+w₁×u₁

上式中，a表示灰度级变量，f_a表示灰度值为a的像素点的个数。目标类与背景类的类间方差为G(a)：

G(a)＝w₀(u₀-u_G)²+w₁(u₁-u_G)²

步骤4.2,将a从开始依次遍历到记录当G(a)的值最大时对应的a为平滑图像G的最优阈值T_G，μ和D定义同步骤2。

步骤5,利用二维直方图P_ij和二维最大类间方差阈值法，在二维空间搜索最优阈值(S^*,T^*)，D为步骤2中求得的含噪图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，使得二维最大类间方差最大时对应的阈值为最优二维阈值(S^*,T^*)；

步骤5.1统计噪音图像F的二维直方图P_ij，公式如下，

f_ij为像素灰度值为i且3×3邻域平均灰度值为j的像素点的个数，M和N的定义同步骤1；

步骤5.2,设分割阈值为(s,t)，将灰度值小于s且邻域灰度均值小于t的像素点划分为背景类；将灰度值大于s且邻域灰度均值大于t的像素点划分为目标类。统计目标类和背景类的像素点个数，计算目标类和背景类的灰度均值矢量μ₀，μ₁，全局的灰度均值矢量μ_T，求解公式为

上式中，P_ij定义同步骤5.1,T表示转置矩阵，

计算目标类与背景类的离散度矩阵G_(s,t)的公式为：

P₀，P₁分别为目标类像素和背景类像素占总像素数的比例，如下式，

P₀＝Σ_{(i＞s且j＞t)}P_ij P₁＝Σ_{(0＜i＜s且0＜j＜t)}P_ij

计算两类的离散度矩阵的迹tr(G_(s,t))公式为：

tr(G_(s,t))＝P₀[(μ_0i-μ_Ti)²+(μ_0j-μ_Tj)²]+P₁[(μ_1i-μ_Ti)²+(μ_1j-μ_Tj)²]

μ_0i，μ_0j，μ_1i，μ_1j，μ_Ti，μ_Tj是灰度均值向量μ₀，μ₁，μ_T的分量；

定义tr(G_(s,t))的解的搜索空间为D为步骤2中求得的噪音图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，且s和t的取值范围均为0到255，将(s,t)从开始依次遍历到当tr(G_(s,t))的值最大时对应的(s,t)为最优二维阈值(S^*,T^*)。

步骤6，利用最佳阈值(S^*,T^*)，对原始图像进行分割并输出，为了简化分割操作，分割方式如下：将图像F灰度值大于S^*且临域灰度均值大于T^*的像素点(即目标类中的像素点)的灰度值重置为n，将其余像素点的灰度值重置为m。其中，n＝255，m＝0，需要注意的是，如果是把图像分成两类通常是0和255，当然也可以不是0和255，而是其他的指定值都可以，只要两个值有较大的区分度都可以，一般只需n-m>＝100即可。

本发明提出的一种快速的含噪图像二维最大类间方差法，在求解最优阈值效率上比原始二维最大类间方差法有明显优越性，且阈值准确度也与原始二维最大类间方差法的相同，以下通过一组实验来说明，随机选用12幅标准测试图像，均为512*512像素，并添加10％的高斯噪音处理,选用测试图像分别是lena,cameraman,woman,bubbles,pirate,livingroom,mandril,peppers,walkbridge,jetpla ne,blonde,airfield，格式为.GIF,仿真测试环境为Matlab2014b，3.2GHz CPU，4G RAM，操作系统为Windows 7旗舰版。本发明提出方法和原始二维最大类间方差法的最优阈值和运行时间对比表如表1所示：

表1：原始二维Otsu和优化二维Otsu对比表

见表1，12幅图像用原始二维最大类间方差法求得的阈值均落在区域内，说明优化二维最大类间方差法求得的阈值与原始二维最大类间方差法求得的阈值一样，而其搜索区间为远远小原始算法的搜索区间和搜索次数，说明本发明的高效性和正确性。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (3)

1.一种快速的含噪图像二维最大类间方差阈值法，其特征在于，包含如下步骤：

步骤1，输入待分割含噪音图像F，用用f(i,j)表示原始图像在像素点(i,j)处的灰度值，M和N分别表示图像的高度和宽度，单位是像素，读取原始图像各个像素点的灰度值，统计图像每个灰度级L_k像素出现的频数H_k，这里k的取值范围是0到255；

步骤2，求出含噪图像F的灰度均值μ和灰度标准差D，公式为：

$\mathrm{\μ}=\frac{1}{M\×N}\underset{k=0}{\overset{255}{\Σ}}{L}_k\×H_k$

$D=\sqrt{\frac{1}{M\×N}\underset{k=0}{\overset{255}{\Σ}}{(L_k-\mathrm{\μ})}^2\×H_k}$

式中，M，N，L_k,H_k的定义同步骤1，k是循环变量；

步骤3，对图像每个像素采用设定邻域的平均灰度值进行平滑，得到平滑图像G；

步骤4,利用最大类间方差阈值法，将灰度级从到开始依次遍历获取平滑图像G的阈值T_G，具体包括：

步骤4.1,读入平滑图像G，统计图像灰度级a像素出现的频数f_a，设分割阈值为T_Gx，目标类像素点总数为n₀，背景类的像素点总数为n₁，目标类和背景类的点占图像比例分别为w₀,w₁，目标类和背景类的灰度均值分别为u₀，u₁，总的灰度均值为u_G，求解公式如下：

$w_0=\frac{n_0}{M\×N}$

$w_1=\frac{n_1}{M\×N}$

$u_0=\frac{1}{n_0}\underset{a=T_{Gx}+1}{\overset{255}{\Σ}}a\×f_a$

$u_1=\frac{1}{n_1}\underset{a=0}{\overset{T_{Gx}}{\Σ}}a\×f_a$

u_G＝w₀×u₀+w₁×u₁

上式中，a表示灰度级变量，f_a表示灰度值为a的像素点的个数，目标类与背景类的类间方差为G(a)：

G(a)＝w₀(u₀-u_G)²+w₁(u₁-u_G)²

步骤4.2,将a从开始依次遍历到记录当G(a)的值最大时对应的a为平滑图像G的最优阈值T_G，μ和D定义同步骤2；

步骤5,利用二维直方图P_ij和二维最大类间方差阈值法，在二维空间]搜索最优阈值(S^*,T^*)，D为步骤2中求得的含噪图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，使得二维最大类间方差最大时对应的阈值为最优二维阈值(S^*,T^*)，具体包括；

步骤5.1统计噪音图像F的二维直方图P_ij，公式如下，

$P_{ij}=\frac{f_{ij}}{M\×N},$

f_ij为像素灰度值为i且3×3邻域平均灰度值为j的像素点的个数，M和N的定义同步骤1；

步骤5.2,设分割阈值为(s,t)，将灰度值小于s且邻域灰度均值小于t的像素点划分为背景类；将灰度值大于s且邻域灰度均值大于t的像素点划分为目标类；统计目标类和背景类的像素点个数，计算目标类和背景类的灰度均值矢量μ₀，μ₁，全局的灰度均值矢量μ_T，求解公式为，

${\mathrm{\μ}}_0={({\mathrm{\μ}}_{0i},{\mathrm{\μ}}_{0j})}^T={\[\underset{i=s+1}{\overset{255}{\Σ}}\underset{j=i+1}{\overset{255}{\Σ}}{\mathrm{iP}}_{ij},\underset{i=s+1}{\overset{255}{\Σ}}\underset{j=i+1}{\overset{255}{\Σ}}{\mathrm{jP}}_{ij}\]}^T$

${\mathrm{\μ}}_1={({\mathrm{\μ}}_{1i},{\mathrm{\μ}}_{1j})}^T={\[\underset{i=1}{\overset{s}{\Σ}}\underset{j=1}{\overset{t}{\Σ}}{\mathrm{iP}}_{ij},\underset{i=1}{\overset{s}{\Σ}}\underset{j=1}{\overset{t}{\Σ}}{\mathrm{jP}}_{ij}\]}^T$

${\mathrm{\μ}}_T={({\mathrm{\μ}}_{Ti},{\mathrm{\μ}}_{Tj})}^T={\[\underset{i=0}{\overset{255}{\Σ}}\underset{j=0}{\overset{255}{\Σ}}{\mathrm{iP}}_{ij},\underset{i=0}{\overset{255}{\Σ}}\underset{j=0}{\overset{255}{\Σ}}{\mathrm{jP}}_{ij}\]}^T$

上式中，P_ij定义同步骤5.1,T表示转置矩阵，

计算目标类与背景类的离散度矩阵G_(s,t)的公式为：

$G_{(s,t)}=\underset{b=0}{\overset{1}{\Σ}}{P}_b\[({\mathrm{\μ}}_b-{\mathrm{\μ}}_T){({\mathrm{\μ}}_b-{\mathrm{\μ}}_T)}^T\]$

P₀，P₁分别为目标类像素和背景类像素占总像素数的比例，如下式，

P₀＝∑_{(i＞s且j＞t)}P_ijP₁＝∑_{(0＜i＜s且0＜j＜t)}P_ij

计算两类的离散度矩阵的迹tr(G_(s,t))公式为：

tr(G_(s,t))＝P₀[(μ_0i-μ_Ti)²+(μ_0j-μ_Tj)²]+P₁[(μ_1i-μ_Ti)²+(μ_1j-μ_Tj)²]

${\mathrm{\μ}}_{0i}=\underset{i=0}{\overset{s}{\Σ}}\underset{j=0}{\overset{t}{\Σ}}{\mathrm{iP}}_{ij},{\mathrm{\μ}}_{0j}=\underset{i=0}{\overset{s}{\Σ}}\underset{j=0}{\overset{t}{\Σ}}{\mathrm{jP}}_{ij},{\mathrm{\μ}}_{1i}=\underset{i=s+1}{\overset{255}{\Σ}}\underset{j=t+1}{\overset{255}{\Σ}}{\mathrm{iP}}_{ij},$

${\mathrm{\μ}}_{1j}=\underset{i=s+1}{\overset{255}{\Σ}}\underset{j=t+1}{\overset{255}{\Σ}}{\mathrm{jP}}_{ij},{\mathrm{\μ}}_{Ti}=\underset{i=0}{\overset{255}{\Σ}}\underset{j=0}{\overset{255}{\Σ}}{\mathrm{iP}}_{ij},{\mathrm{\μ}}_{Tj}=\underset{i=0}{\overset{255}{\Σ}}\underset{j=0}{\overset{255}{\Σ}}{\mathrm{jP}}_{ij}$

μ_0i，μ_0j，μ_1i，μ_1j，μ_Ti，μ_Tj是步骤5.2中灰度均值向量μ₀，μ₁，μ_T的分量；

定义tr(G_(s,t))的解的搜索空间为D为步骤2中求得的噪音图像F的灰度标准差，T_G为步骤4中求得的平滑图像G的一维阈值，且s和t的取值范围均为0到255，将(s,t)从开始依次遍历到当tr(G_(s,t))的值最大时对应的(s,t)为最优二维阈值(S^*,T^*)；

步骤6，利用最佳阈值(S^*,T^*)，对原始图像进行分割并输出。

2.如权利要求1所述的一种快速的含噪图像二维最大类间方差阈值法，其特征在于，所述步骤2求出含噪图像F的灰度均值μ和灰度标准差D的公式为，

$\mathrm{\μ}=\frac{1}{M\×N}\underset{k=0}{\overset{255}{\Σ}}{L}_k\×H_k$

$D=\sqrt{\frac{1}{M\×N}\underset{k=0}{\overset{255}{\Σ}}{(L_k-\mathrm{\μ})}^2\×H_k}$

式中，M，N，L_k,H_k的定义同步骤1，k是循环变量。

3.如权利要求1所述的一种快速的含噪图像二维最大类间方差阈值法，其特征在于，所述步骤6利用最佳阈值(S^*,T^*)，对原始图像进行分割，按如下方式进行：将图像F灰度值大于S^*且设定邻域灰度均值大于T^*的像素点的灰度值重置为n，将其余像素点的灰度值重置为m，且n和m为0至255的任意整数，且n不等于m。