一种基于直方图的EM多阈值图像分割方法及装置与流程

本公开涉及图像处理领域，具体涉及一种基于直方图的em多阈值图像分割方法及装置。

背景技术：

图像分割就是将一幅图像划分成若干聚类区域，每个区域都有相似的特性。在众多的图像分割方法中，直方图阈值法因其简单高效备受注意。它的基本思想是目标和背景分布在直方图灰度轴的两端，只要找到这个分割点，就可以进行图像分割，最大类间方差法(大津otsu)、最大熵法、模糊聚类、期望最大(em)值法等，都是典型的阈值分割方法。随着研究的深入、应用领域的扩大，阈值分割算法也在不断进化，出现了多阈值分割、融合目标特征的分割、以及利用二维直方图的多维度分割算法等。

由于图像像素绝大部分都分布在二维直方图的主对角线附近，这是因为除非边缘点，否则像素灰度与其邻域的平均灰度非常接近。利用这个特点，现在的图像分割技术主要有两种分割方法，第一种图像分割方法为，以(s,t)为分割点，将图像根据对角线区域划分为目标和背景，其余是边缘或噪声。第一种图像分割方法采用双斜线分割，靠近主对角附近的区域是正常像素，远离对角线的点则为边缘或噪声点，第一种方法假设阈值附近像素点出现概率为零，是一种近似划分，造成区域错分而第二种分割方法虽然将概率不为零的点都包含了进来，但是斜线的截距却难以确定。

技术实现要素：

本公开提供一种基于直方图的em多阈值图像分割方法及装置，采用二维高斯混合模型(gaussianmixturemodel，gmm)，用二维gmm拟合二维直方图曲线，采用em算法获取gmm的分布参，在构建二维直方图的基础上，分析了二维混合模型参数估计推导方法，提出了一种自适应方法决定混合分量个数，利用em算法实现图像分割，获得像素点的统计分布，通过确定类属完成图像分割，通过获取像素点的分布类属，在斜分法的基础上提出利用置信度进行区域划分，在二维直方图的基础上利用em算法实现了多阈值分割，直方图阈值是图像分割的重要依据，二维直方图因为考虑了每个像素的邻域相关信息，对噪声图像的分割效果更好，有别于传统的最大类间方差及最大熵阈值法。

为了实现上述目的，根据本公开的一方面，提供一种基于直方图的em多阈值图像分割方法，所述方法包括以下步骤：

步骤1，根据信号点判断条件寻找图像中的信号点区域；

步骤2，根据贝叶斯最小错误率准则在信号点区域中对图像进行标记划分出图像的噪声点区域；

步骤3，在噪声点区域，计算邻域灰度值与各混合分量的欧拉距离，判断类属。

进一步地，在步骤1中，所述根据信号点判断条件寻找图像中的信号点区域的方法为，在图像中按照信号点判断条件的约束条件，所述约束条件为：

通过约束条件进行寻找图像中的信号点区域，即图中满足约束条件的区域，所述信号点区域为图像中的一部分区域，任意像素点p在坐标(x，y)处的灰度值为f(x，y)，取p的8-邻域平均灰度值g(x，y)，二者的值域都在[0,l]之间，其中，a1、a2、b1和b2指纵截距。

进一步地，在步骤2中，所述贝叶斯最小错误率准则的方法包括以下步骤：

步骤2.1，获取图像的二维直方图，输入一个的图像，所述图像的图像矩阵为m行n列，根据每个像素本身的灰度值，每个像素的邻域平均灰度值，由每个像素本身的灰度值和邻域平均灰度值构成了二维坐标系；

设任意像素点p在坐标(x，y)处的灰度值为f(x，y)，取p的8-邻域(n8)平均灰度值g(x，y)，二者的值域都在[0,l]之间；

根据公式计算邻域灰度，即获取到二维直方图，式中w是以p点为中心的滤波掩膜,

步骤2.2，对二维直方图建立二维高斯混合模型，图像有m*n个像素，样本rxy(i,j)表示第x行y列像素点的本体灰度和邻域灰度，且样本互相独立，建立每个rxy的二维混合概率密度函数p(rxy；θ)，构造似然函数l(θ；r)，其中，pk(rxy；θk)是rxy第k个二维高斯独立分布的密度函数，θk是其参数向量，θk＝{μ1k,μ2k,σ1k,σ2k}。ωk是第k簇的混合比例系数，满足ω1+...+ωk＝1的条件，k是混合模型的分量数目，θ是混合模型的参数集，θ＝{θ1,θ2,…,θk}；设r＝{r(i,j),i,j＝0,1,...,l}是图像二维直方图二元点集，h(i,j)是二维直方图，这样在em算法中就可以脱离图像，仅对二维直方图运算，则二维高斯混合模型为，

步骤2.3，计算混合分量的个数k与各个分布的权重ωk，引入隐含类别标签设是为示性函数，取值只有0和1，在给定类属标签后，每个r(i,j)的混合概率密度函数转化为混合权重ωk由类属标签的分布；设r(i,j)属于第k个类属的概率且在已知类属标签的情况下，每个r(i,j)的混合概率密度函数二维高斯混合模型转换为，

步骤2.4，根据类属标签，采用贝叶斯后验概率猜测隐含类别，借用em算法迭代更新，首先假定二维混合模型的参数已知，且类属的分布概率也已知，然后根据样本观测值，计算已知观测样本的类属分布，该分布就是关于已有参数的后验概率，记做观测样本数值引起变化，模型参数也发生变化，通过求取似然函数极大值，获得新的参数估计子对于二维高斯分布，设两个变量互相独立，那么单个二维高斯分布如公式

二维混合高斯模型的似然函数如公式，

将二维混合高斯模型的似然函数先求和再进行对数运算；

步骤2.5，根据jensen不等式，一个凹函数有f(ex)≥e[f(x)]，考虑到log(x)是凹函数，且就是的期望，那么由jensen不等式，二维混合高斯模型的似然函数转换为,

步骤2.6，获取每一个k分布，分开求取偏导数，得到下面的参数更新公式。

同时得到类属更新分布φk，

步骤2.7，重复执行步骤2.4，直到得到的各值都满足收敛条件为止。

进一步地，在步骤2中，所述对图像进行标记划分出图像的噪声点区域为取图像矩阵像素点的最大后验概率公式的k的区域。

进一步地，在步骤3中，在噪声点区域，计算邻域灰度值与各混合分量的欧拉距离，判断类属，

(i,j)∈k,ifd＝min(g(i,j)-μ2k),k＝1,2,,…,k，k为混合分量的个数，(i,j)为样本rxy(i,j)表示第i行j列像素点的本体灰度和邻域灰度，平均灰度值g(x，y)，μ2k为常数且取值范围为负无穷大到正无穷大之间，根据类属更新分布φk，判断图像中噪声点区域的各个像素点的类属。

优选地，采用em统计法，因为获得了像素点的统计分布，通过确定类属，即可完成图像分割，根据每个像素点的类属进行图像分割。

本发明还提供了一种基于直方图的em多阈值图像分割装置，所述装置包括：存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序运行在以下装置的单元中：

信号点区域寻找单元，用于根据信号点判断条件寻找图像中的信号点区域；

噪声点标记单元，用于根据贝叶斯最小错误率准则在信号点区域中对图像进行标记划分出图像的噪声点区域；

像素类属判断单元，用于在噪声点区域，计算邻域灰度值与各混合分量的欧拉距离，判断类属。

本公开的有益效果为：本发明提供一种基于直方图的em多阈值图像分割方法及装置，分量越多，迭代步数和训练时间越长，不会错分区域，方便确定领域灰度划分的斜线的截距，随着混合模型个数的增加，区域划分越来越精细，特征越来越明显。较好地自适应完成图像的多阈值分割。保证了图像分割的正确率。

附图说明

通过对结合附图所示出的实施方式进行详细说明，本公开的上述以及其他特征将更加明显，本公开附图中相同的参考标号表示相同或相似的元素，显而易见地，下面描述中的附图仅仅是本公开的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，在附图中：

图1为二维直方图的图例；

图2为现有的方法采取的图像分割方法；

图3为加入椒盐噪声后的图像和二维直方图；

图4为分割结果图；

图5所示为一种基于直方图的em多阈值图像分割装置图。

具体实施方式

以下将结合实施例和附图对本公开的构思、具体结构及产生的技术效果进行清楚、完整的描述，以充分地理解本公开的目的、方案和效果。需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。

本公开提出一种基于直方图的em多阈值图像分割方法，具体包括以下步骤：

获取二维直方图，二维直方图是在一维直方图的基础上考虑邻域后的二维分布。一个m行n列的图像，除了每个像素有自己的灰度值，还可以计算每个像素的邻域平均灰度值。这两个变量就构成了二维坐标系。设任意像素点p在坐标(x，y)处的灰度值为f(x，y)，取p的8-邻域(n8)平均灰度值g(x，y)，二者的值域都在[0,l]之间。邻域灰度的计算如公式(1)，式中w是以p点为中心的滤波掩膜。

定义二维直方图h(i,j)，它表示图像中当f(x,y)＝i且g(x,y)＝j时像素的个数或频次(i,j＝0,1,...l)，二元对(i,j)的联合概率为公式(3)。为了书写方便，二元对(i,j)简写为r(i,j)，图像中像素点(x,y)具有二元对(i,j)灰度记为rxy(i,j)，简记rxy。

以米粒图为例，它的二维直方图，如图1所示，图1为二维直方图的图例，其中，图1中的(a)为原图,(b)为一维直方图,(c)为二维直方图投影,(d)为二维直方图，，其形状表明可以用多个高斯分布叠加进行直方图分析。

采用em算法实现聚类分析，利用多个高斯分布的叠加对直方图进行逼近，得到聚类的高斯混合模型是图像分割的重要方法之一。在进行混合模型的参数估计过程中，观测数据到底属于哪个类属也就是类标签是未知的，是个不完全数据集，属于隐性问题，可以用期望最大化(expectationmaximization,em)算法进行参数估计。em算法充分挖掘不完全数据集的信息，通过引入丢失数据集，获得不完全数据集的最大似然函数得到有效的参数估计值。

对二维直方图建立二维高斯混合模型，图像有m×n个像素，样本rxy(i,j)表示第x行y列像素点的本体灰度和邻域灰度，且样本互相独立。建立每个rxy的二维混合概率密度函数p(rxy；θ)，构造似然函数l(θ；r)，如公式(4)(5)所示。

这里pk(rxy；θk)是rxy第k个二维高斯独立分布的密度函数，θk是其参数向量，θk＝{μ1k,μ2k,σ1k,σ2k}。ωk是第k簇的混合比例系数，满足ω1+...+ωk＝1的条件。k是混合模型的分量数目，θ是混合模型的参数集，θ＝{θ1,θ2,…,θk}。设r＝{r(i,j),i,j＝0,1,...,l}是图像二维直方图二元点集，h(i,j)是二维直方图，则公式(5)又可化为公式(6)，这样在em算法中就可以脱离图像，仅对二维直方图运算。

因为数据信息的不完备，存在两个问题，一是如何估计混合分量的个数k；二是怎样估计各个分布的权重ωk。针对这两个问题，引入隐含类别标签为了简化操作，设是一个示性函数，取值只有0和1，见公式(7)。在给定类属标签后，每个r(i,j)的混合概率密度函数就变化为混合权重ωk就由类属标签的分布代替。设r(i,j)属于第k个类属的概率且在已知类属标签的情况下，公式(4)(6)修改为公式(8)(9)。

上面是给定了类属标签，但是实际上类属依然是未知的，因此采用贝叶斯后验概率猜测隐含类别，借用em算法迭代更新。

e步骤：首先假定二维混合模型的参数已知，且类属的分布概率也已知，然后根据样本观测值，计算已知观测样本的类属分布，该分布就是关于已有参数的后验概率，记做

m步骤：在上面的e步骤里，观测样本数值引起变化，模型参数也发生变化，通过求取似然函数极大值，获得新的参数估计子对于二维高斯分布，设两个变量互相独立，那么单个二维高斯分布如公式(11)所示：

二维混合高斯模型的似然函数如公式(12)所示：

在上面的公式里，先求和再进行对数运算，如果将二者的顺序颠倒，将大大简化计算。根据jensen不等式，一个凹函数有f(ex)≥e[f(x)]，考虑到log(x)是凹函数，且就是的期望。那么由jensen不等式，公式(12)可以整理成先求对数再求和如公式(13)。

每一个k分布可以分开求取偏导数，得到下面的参数更新公式。

同时得到类属更新分布φk。

e步骤和m步骤形成迭代关系，它们重复执行，直到得到的各值都满足收敛条件为止。

进行图像分割，在图1(c)中，图像像素绝大部分都分布在二维直方图的主对角线附近，这是因为除非边缘点，否则像素灰度与其邻域的平均灰度非常接近。利用这个特点，现有的方法采取的图像分割方法如图2和图3所示，图2为现有的方法采取的图像分割方法，图3为加入椒盐噪声后的图像和二维直方图，以(s,t)为分割点，对角线区域1，2是目标和背景，其余是边缘或噪声。采取的图像分割如图所示，采用双斜线分割，靠近主对角附近的区域1，2是正常像素，远离对角线的点则为边缘或噪声点。第一种方法假设阈值附近像素点出现概率为零，是一种近似划分，造成区域错分。第二种分割方法虽然将概率不为零的点都包含了进来，但是斜线的截距却难以确定。

而采用em统计法，因为获得了像素点的统计分布，通过确定类属，即可完成图像分割。只要知道像素点的分布类属，即可获得本文在斜分法的基础上提出利用置信度进行区域划分。

图像分割的具体算法如下：

第一步：寻找信号点区域，判断条件：

第二步，在信号点区域，根据贝叶斯最小错误率准则对图像进行标记，就是最大后验概率公式(10)的k取值。

第三步，在噪声点区域，计算邻域灰度值与各混合分量的欧拉距离，判断类属。

(i,j)∈k,ifd＝min(g(i,j)-μ2k),k＝1,2,,…,k

在以上的步骤中，每个像素点都进行了类属判断。保证了图像分割的正确率。

实验结果与分析

为了验证基于二维直方图的em分割算法的分割效果，本文选取2张标准测试图像lenna图和cameraman，以及2张来自weizmannsegmentationdatabase(wsd)的自然景物图像作为实验对象，进行阈值分割实验。实验是在3.60ghzcpu和4.00gb内存的pc机、matlab2015b环境中进行的。

不同混合分量数量的分割实验

对四幅图像分别进行高斯混合模型数量是2、3、4三种em分割，分割后的图像如图4所示，图4为分割结果图，从上至下分别为lenna图、cameraman图、buggy图、pengium图，从左至右分别为原图、双混合模型、三混合模型和四混合模型em分割后的图像，四混合模型采用黑、白、红、绿四个颜色表示分割后图像。

采用em算法分割图像，只要给定混合模型的个数，算法可以自动收敛得到多阈值分割结果。图4可以看出，随着混合模型个数的增加，区域划分越来越精细，特征越来越明显。较好地自适应完成图像的多阈值分割。下表是混合模型建立的迭代步数以及训练时间表。

表1基于二维直方图em分割算法的速度

从表1中可以看出，分量越多，迭代步数和训练时间越长。

去噪性能比较

为了验证基于二维直方图em分割算法的优越性，将标准图像rice进行二值分割，并与一维em进行对比实验：

采用两种评价方法：pr(precision-recall)曲线以及f-measure评价法。

第一种方法将阈值从0逐步增加到255，用每个阈值图像分割，并计算出该阈值下的精确率(precision)和召回率(recall)。精确率的计算公式召回率的计算公式tp表示把正类预测为正类，fp表示把负类预测为正类，fn则表示把原来的正类预测为负类。精确率的分母是所有预测为正的样本数，精确率是针对预测结果而言的，表示预测为正的样本中有多少是真正的正样本，而召回率的分母原来样本中所有的正样本数，召回率是针对原来的样本而言的，表示样本中的正例有多少被预测正确了。

本公开的实施例提供的一种基于直方图的em多阈值图像分割装置，如图5所示为本公开的一种基于直方图的em多阈值图像分割装置图，该实施例的一种基于直方图的em多阈值图像分割装置包括：处理器、存储器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述一种基于直方图的em多阈值图像分割装置实施例中的步骤。

所述装置包括：存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序运行在以下装置的单元中：

信号点区域寻找单元，用于根据信号点判断条件寻找图像中的信号点区域；

噪声点标记单元，用于根据贝叶斯最小错误率准则在信号点区域中对图像进行标记划分出图像的噪声点区域；

像素类属判断单元，用于在噪声点区域，计算邻域灰度值与各混合分量的欧拉距离，判断类属。

所述一种基于直方图的em多阈值图像分割装置可以运行于桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备中。所述一种基于直方图的em多阈值图像分割装置，可运行的装置可包括，但不仅限于，处理器、存储器。本领域技术人员可以理解，所述例子仅仅是一种基于直方图的em多阈值图像分割装置的示例，并不构成对一种基于直方图的em多阈值图像分割装置的限定，可以包括比例子更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述一种基于直方图的em多阈值图像分割装置还可以包括输入输出设备、网络接入设备、总线等。

所称处理器可以是中央处理单元(centralprocessingunit，cpu)，还可以是其他通用处理器、数字信号处理器(digitalsignalprocessor，dsp)、专用集成电路(applicationspecificintegratedcircuit，asic)、现成可编程门阵列(field-programmablegatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等，所述处理器是所述一种基于直方图的em多阈值图像分割装置运行装置的控制中心，利用各种接口和线路连接整个一种基于直方图的em多阈值图像分割装置可运行装置的各个部分。

所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述一种基于直方图的em多阈值图像分割装置的各种功能。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、图像播放功能等)等；存储数据区可存储根据手机的使用所创建的数据(比如音频数据、电话本等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(smartmediacard,smc)，安全数字(securedigital,sd)卡，闪存卡(flashcard)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。

尽管本公开的描述已经相当详尽且特别对几个所述实施例进行了描述，但其并非旨在局限于任何这些细节或实施例或任何特殊实施例，而是应当将其视作是通过参考所附权利要求考虑到现有技术为这些权利要求提供广义的可能性解释，从而有效地涵盖本公开的预定范围。此外，上文以发明人可预见的实施例对本公开进行描述，其目的是为了提供有用的描述，而那些目前尚未预见的对本公开的非实质性改动仍可代表本公开的等效改动。