CN109816660A - 一种图像分割方法、终端设备及存储介质 - Google Patents

Abstract

本发明涉及一种图像分割方法、终端设备及存储介质，在该方法中利用马尔可夫随机场能够有效描述图像的空间相关性这一特性，将局部空域信息结合至欧式距离测度中，构造一个具有邻域约束的约束项，同时设计了一个具有全局特征的非局部正则项，而后将两者有机结合，建立新的模糊C均值算法的目标函数。本发明所述方法结合图像分割问题的随机性和模糊性特点，合理获取图像的先验知识，一定程度上消除了图像中低对比度的影响，同时算法也能够对图像中存在的亮度不均匀问题加以考虑，有效地提高了图像的分割精度。

Description

一种图像分割方法、终端设备及存储介质

技术领域

本发明涉及图像处理技术领域，尤其涉及一种图像分割方法、终端设备及存储介质。

背景技术

在临床诊断过程中，有时为了检测生理或病理状态的脑组织容量变化，往往需要对医学影像中的脑磁共振(Magnetic Resonance，MR)图像进行精确的组织分割，这样可以定量测量脑组织或病灶的截面面积或体积，为进一步诊疗方案的制定提供参考依据。但核磁成像过程中由于受到成像设备自身性能、射频场强的不均匀性以及患者体位变化等因素的影响，往往在影像中存在着噪声、亮度不均匀以及对比度低等问题，这给精确的脑组织分割带来了干扰。

模糊C均值(Fuzzy C-means，FCM)聚类是基于模糊集理论处理数据划分问题的一种方法，近些年来，FCM算法已成为图像分割的有效手段之一，其成功原因要归功于它所具有的模糊性质，能够从原始图像中保留更多的信息。但FCM也存在一定缺陷，它通常对未被噪声破坏的图像表现良好，而当图像遭受噪声、异常值或者伪影攻击的时候，其分割结果不尽人意，对真实的脑MR图像分割效果较差。分析其原因，在FCM的目标函数中仅仅使用了非鲁棒的欧式距离作为相似测度，它忽略了图像邻域内相邻像素间相互关系。为此，学者们提出了许多的改进方案，2001年，Pham等人提出了一种模糊聚类的空域模型，在目标函数中引入局部空间正则项，用以消除图像中像素不太可能的分类结果，该方法被称为鲁棒模糊C均值聚类算法(Robust Fuzzy C-Means，R-FCM)。2013年，Cui等提出一种模糊能量最小化方法，首先定义了围绕每个局部邻域中心点的局部FCM聚类的目标函数，然后，该目标函数关于整个图像中的邻域中心积分以形成全局模糊能量，通过最小化全局模糊能量以实现图像的分割。这些方法虽然在一定程度上克服了噪声的影响，但是图像的细节信息和边缘处理效果不好。

发明内容

针对上述问题，本发明旨在提供一种图像分割方法、终端设备及存储介质，结合马尔科夫随机场模型，充分考虑到空域各像素点间相互关系，对图像邻域关系进行有效描述，从而更好的分割图像。

具体方案如下：

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

S1：输入待分割的图像；

S2：初始化模型参数，设定聚类个数c、隶属度模糊指数m、调节系数β、马尔可夫随机场的邻域尺寸、基团势能参数γ、搜索窗的半径u，以及聚类停止迭代条件ε；

S3：在图像中，选定一个邻域空间，根据马尔可夫随机场的Potts模型，计算该邻域空间的先验概率

其中，

设定{z₁,z₂,…,z_N}表示一组随机变量，i∈[1,N]，Λ表示该组随机变量的配置空间，k表示类别，k∈[1,c]，z_i＝k表示随机变量z_i属于第k类，

在图像中，选取一个邻域空间，则下标i表示邻域中的中心像素，N_i表示以第i个像素为中心像素的邻域，下标j表示在邻域N_i内的第j个像素，表示邻域中所有像素的类别；

S4：根据先验概率和像素与聚类中心的欧式距离d_ki，结合空间邻域信息，计算像素与聚类中心的距离测度D_ki为：

其中，Sqrt函数表示平方根。

S5：设定非局部距离测度加权测度S_ij的计算公式为：

其中，σ表示噪声的方差；

h表示与噪声的方差σ成比例的常数，它控制着非局部距离测度加权测度S_ij的衰减程度，并满足0＜S_ij＜1且

表示以第i个像素为中心像素，半径为u的搜索窗；

x(N_i)表示以x_i为中心的邻域矢量；

E_i表示归一化常数，通过以下公式计算获得：

S6：融合马尔可夫随机场模型和模糊C均值算法，设定目标函数为：

其中，c表示聚类个数，N表示图像中的像素个数，u_ki表示模糊隶属度，u_jl表示隶属度变量，为邻域内第l个像素属于第j类的概率，L_k＝{1,2,…,c}\{k}，k∈[1,c]，表示不包括第k类的其他类；

S7：根据拉格朗日最优准则，计算目标函数的模糊隶属度u_ki和聚类中心v_k：

其中，t表示迭代次数；

S8：利用k均值聚类算法计算初始聚类中心同时设定迭代次数t的初始值t＝0，并初始化模糊隶属度

S9：计算并更新模糊隶属度和聚类中心

S10：判断是否满足max||v^(t+1)-v^(t)||＜ε，如果满足，进入S11，否则，设定t＝t+1，返回S9；

S11：输出分割后的图像。

进一步的，所述模糊隶属度u_ki的计算过程如下：

S911：由于且m＞1，根据拉格朗日最优准则，设：

其中，λ_i表示拉格朗日乘子法中的约束系数；

S912：计算L关于模糊隶属度u_ki的偏导数，并设得

S913：根据步骤S902中公式得：

S914：计算L关于约束系数λ_i的偏导数，并设得

S915：根据步骤S903和S904中公式得：

S916：根据步骤S903和S905中公式得模糊隶属度u_ki为：

进一步的，所述聚类中心v_k的计算过程如下：

S921：计算L关于聚类中心v_k的偏导数，并设得：

S922：根据步骤S921得聚类中心v_k为：

一种图像分割终端设备，包括处理器、存储器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本发明实施例上述的方法的步骤。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现本发明实施例上述的方法的步骤。

本发明采用如上技术方案，利用马尔可夫随机场能够有效描述图像的空间相关性这一特性，将局部空域信息结合至欧式距离测度中，构造一个具有邻域约束的约束项，同时设计了一个具有全局特征的非局部正则项，而后将两者有机结合，建立新的FCM聚类的目标函数。该方法结合图像分割问题的随机性和模糊性特点，合理获取图像的先验知识，一定程度上消除了图像中低对比度的影响，同时算法也能够对图像中存在的亮度不均匀问题加以考虑，有效地提高了图像的分割精度。

附图说明

图1所示为本发明实施例一的流程示意图。

图2所示为本发明实施例一中各种算法的图像分割结果比较示意图。

图3所示为本发明实施例一中各种算法的误分类率的结果比较示意图。

图4所示为本发明实施例一中存在灰度不均匀性和噪声的各种算法的结果比较示意图。

具体实施方式

为进一步说明各实施例，本发明提供有附图。这些附图为本发明揭露内容的一部分，其主要用以说明实施例，并可配合说明书的相关描述来解释实施例的运作原理。配合参考这些内容，本领域普通技术人员应能理解其他可能的实施方式以及本发明的优点。

现结合附图和具体实施方式对本发明进一步说明。

实施例一：

参考图1所示，本发明提供了一种图像分割方法，包括以下步骤：

S1：输入待分割图像。

该实施例中选用的待分割图像为脑MR图像。

S2：初始化模型参数：设定聚类个数c、隶属度模糊指数m、调节系数β、马尔可夫随机场的邻域尺寸、基团势能参数γ、搜索窗的半径u，以及聚类停止迭代条件ε。

S3：在图像中，选定一个邻域空间，根据马尔可夫随机场的Potts模型，计算该邻域空间的先验概率

其中，

设定{z₁,z₂,…,z_N}表示一组随机变量，i∈[1,N]，Λ表示该组随机变量的配置空间，k表示类别，k∈[1,c]，z_i＝k表示随机变量z_i属于第k类。

在图像中，选取一个邻域空间，则下标i表示邻域中的中心像素，N_i表示以第i个像素为中心的邻域，下标j表示在邻域N_i内的第j个像素，表示邻域中所有像素的类别。则表示一个已知邻域内所有像素的类别的条件下该邻域内中心像素属于第k类的概率。

所述马尔科夫随机场的Potts模型表示邻域中的中心像素的所属类别仅受到邻域中的各个像素所归属类别的影响。

S4：根据像素与聚类中心之间的欧式(Euclidean)距离d_ki，结合Potts模型将空间邻域信息结合进来，得到新的距离测度D_ki：

其中，Sqrt函数表示平方根。：

S5：计算非局部距离测度加权测度S_ij：

其中，σ表示噪声的方差；

h表示与噪声的方差σ成比例的常数，它控制着非局部距离测度加权测度S_ij的衰减程度，并满足0＜S_ij＜1且

表示以第i个像素为中心像素，半径为u的搜索窗；

x(N_i)表示以x_i为中心的邻域矢量；

E_i表示归一化常数，通过以下公式计算获得：

S6：融合了马尔可夫随机模型和非局部约束的FCM(简记为NL-FCM-MRF)算法的目标函数设定如下：

其中，调节系数β用来调节目标函数中的前一项与非局部正则项之间的平衡，它的取值可采用交叉验证的方法来进行设定；c表示聚类个数；N表示图像中的像素个数；u_ki表示模糊隶属度；u_jl表示隶属度变量，表示第l个像素属于第j类的概率，这里是在非局部处理中为了区分局部邻域和搜索窗中的邻域而设定的隶属度函数；L_k＝{1,2,…,c}\{k}，k∈[1,c]，表示不包括第k类的其他类；J_m为改进的模糊c均值聚类(FCM)目标函数的一种通用表示方法。

S7：根据拉格朗日最优准则，计算目标函数的模糊隶属度u_ki和聚类中心v_k：

其中，t表示迭代次数，每计算1次模糊隶属度和聚类中心，t增加1。

S8：利用k均值聚类(K-means)算法设定初始聚类中心同时设定迭代次数的初始值t＝0，并初始化模糊隶属度

S9：分别计算并更新模糊隶属度和聚类中心

所述模糊隶属度u_ki和聚类中心v_k的计算过程为：

由于且m＞1，根据拉格朗日最优准则，设

其中，λ_i表示拉格朗日乘子法中的约束系数(也称为乘子)，将约束条件函数与原目标函数联系到一起，用以求出得到原目标函数极值的各个变量的解。

分别计算L关于模糊隶属度u_ki和约束系数λ_i的偏导数，并设得

其中，u_ki和其邻域的乘积项加上相应的邻域和u_ki的反向乘积项通过偏导运算后使得调节系数β的系数消失，将该偏导数置零后，即根据式(2)可得到：

将式(4)带入式(3)得到

再将式(5)带入式(4)得到

相同的方法，求L关于聚类中心v_k的偏导数，并设得到

从式(7)可得到

S10：判断是否满足max||v^(t+1)-v^(t)||＜ε，如果满足，进入S11，否则，设定t＝t+1，返回S9。

S11：输出分割后图像。

本实施例利用马尔可夫随机场能够有效描述图像的空间相关性这一特性，将局部空域信息结合至欧式距离测度中，构造一个具有邻域约束的约束项，同时设计了一个具有全局特征的非局部正则项，而后将两者有机结合，建立新的FCM聚类的目标函数。该方法结合图像分割问题的随机性和模糊性特点，合理获取图像的先验知识，一定程度上消除了图像中低对比度的影响，同时算法也能够对图像中存在的亮度不均匀问题加以考虑，有效地提高了图像的分割精度。

在相同实验条件下，与标准FCM算法、R-FCM算法、MRF算法及所提出的模型的分割结果进行比较如下：

1、图2所示为四组受到强度为9％Rician噪声破坏的2D脑MR切片图像，这四组图像提取的是索引为45、60、93和113的轴向切片图像，图像已被预处理，去除了头盖骨、肌肉和血管等组织，对比结果见图2，其中，图2(b)为FCM算法分割结果，图2(c)为R-FCM算法分割结果，图2(d)为MRF算法分割结果，图2(e)为本实施例中模型的分割结果。

2、为了定量地分析分割效果，采用误分类率(rate，MCR)这一性能指标来进行对比测试，误分类率MCR的计算公式如下式所示。选取八组含有不同噪声等级(从5％-35％)的脑MR图像，分别采用FCM、R-FCM、MRF和本实施例中模型对脑组织(脑白质WM、脑灰质GM和脑脊液CSF)进行分割，对这四种算法的误分类率MCR平均值进行统计，统计结果如图3所示。

3、第三组实验选取的脑MR图像中既存在噪声问题，同时也受到亮度不均匀的影响。这些图像分别从轴、冠和矢位三个不同方向对3D脑MR图像进行切片选取，含有强度为70％的亮度不均匀性并被强度为7％噪声所破坏，分割结果如图4所示。

实施例二：

本发明还提供一种图像分割终端设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现本发明实施例一的上述方法实施例中的步骤。

进一步地，作为一个可执行方案，所述图像分割终端设备可以是桌上型计算机、笔记本、掌上电脑及云端服务器等计算设备。所述图像分割终端设备可包括，但不仅限于，处理器、存储器。本领域技术人员可以理解，上述图像分割终端设备的组成结构仅仅是图像分割终端设备的示例，并不构成对图像分割终端设备的限定，可以包括比上述更多或更少的部件，或者组合某些部件，或者不同的部件，例如所述图像分割终端设备还可以包括输入输出设备、网络接入设备、总线等，本发明实施例对此不做限定。

进一步地，作为一个可执行方案，所称处理器可以是中央处理单元(CentralProcessing Unit，CPU)，还可以是其他通用处理器、数字信号处理器(Digital SignalProcessor，DSP)、专用集成电路(Application Specific Integrated Circuit，ASIC)、现成可编程门阵列(Field-Programmable Gate Array，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。通用处理器可以是微处理器或者该处理器也可以是任何常规的处理器等，所述处理器是所述图像分割终端设备的控制中心，利用各种接口和线路连接整个图像分割终端设备的各个部分。

所述存储器可用于存储所述计算机程序和/或模块，所述处理器通过运行或执行存储在所述存储器内的计算机程序和/或模块，以及调用存储在存储器内的数据，实现所述图像分割终端设备的各种功能。所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序；存储数据区可存储根据手机的使用所创建的数据等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器，例如硬盘、内存、插接式硬盘，智能存储卡(Smart Media Card,SMC)，安全数字(Secure Digital,SD)卡，闪存卡(Flash Card)、至少一个磁盘存储器件、闪存器件、或其他易失性固态存储器件。

本发明还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现本发明实施例上述方法的步骤。

所述图像分割终端设备集成的模块/单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。所述计算机可读介质可以包括：能够携带所述计算机程序代码的任何实体或装置、记录介质、U盘、移动硬盘、磁碟、光盘、计算机存储器、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)以及软件分发介质等。

尽管结合优选实施方案具体展示和介绍了本发明，但所属领域的技术人员应该明白，在不脱离所附权利要求书所限定的本发明的精神和范围内，在形式上和细节上可以对本发明做出各种变化，均为本发明的保护范围。

Claims (5)

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

S1：输入待分割的图像；

S2：初始化模型参数，设定聚类个数c、隶属度模糊指数m、调节系数β、马尔可夫随机场的邻域尺寸、基团势能参数γ、搜索窗的半径u，以及聚类停止迭代条件ε；

S3：在图像中，选定一个邻域空间，根据马尔可夫随机场的Potts模型，计算该邻域空间的先验概率

其中，

设定{z₁,z₂,…,z_N}表示一组随机变量，i∈[1,N]，Λ表示该组随机变量的配置空间，k表示类别，k∈[1,c]，z_i＝k表示随机变量z_i属于第k类，

在图像中，选取一个邻域空间，则下标i表示邻域中的中心像素，N_i表示以第i个像素为中心像素的邻域，下标j表示在邻域N_i内的第j个像素，表示邻域中所有像素的类别；

S4：根据先验概率和像素与聚类中心的欧式距离d_ki，结合空间邻域信息，计算像素与聚类中心的距离测度D_ki为：

S5：设定非局部距离测度加权测度S_ij的计算公式为：

其中，σ表示噪声的方差；

h表示与噪声的方差σ成比例的常数，它控制着非局部距离测度加权测度S_ij的衰减程度，并满足0＜S_ij＜1且

表示以第i个像素为中心像素，半径为u的搜索窗；

x(N_i)表示以x_i为中心的邻域矢量；

E_i表示归一化常数，通过以下公式计算获得：

S6：融合马尔可夫随机场模型和模糊C均值算法，设定目标函数为：

其中，c表示聚类个数，N表示图像中的像素个数，u_ki表示模糊隶属度，u_jl表示隶属度变量，为邻域内第l个像素属于第j类的概率，L_k＝{1,2,…,c}\{k}，k∈[1,c]，表示不包括第k类的其他类；

S7：根据拉格朗日最优准则，计算目标函数的模糊隶属度u_ki和聚类中心v_k：

其中，t表示迭代次数；

S8：利用k均值聚类算法计算初始聚类中心同时设定迭代次数t的初始值t＝0，并初始化模糊隶属度

S9：计算并更新模糊隶属度和聚类中心

S10：判断是否满足max||v^(t+1)-v^(t)||＜ε，如果满足，进入S11，否则，设定t＝t+1，返回S9；

S11：输出分割后的图像。

2.根据权利要求1所述的图像分割方法，其特征在于：所述模糊隶属度u_ki的计算过程如下：

S911：由于且m＞1，根据拉格朗日最优准则，设：

其中，λ_i表示拉格朗日乘子法中的约束系数；

S912：计算L关于模糊隶属度u_ki的偏导数，并设得

S913：根据步骤S902中公式得：

S914：计算L关于约束系数λ_i的偏导数，并设得

S915：根据步骤S903和S904中公式得：

S916：根据步骤S903和S905中公式得模糊隶属度u_ki为：

3.根据权利要求2所述的图像分割方法，其特征在于：所述聚类中心v_k的计算过程如下：

S921：计算L关于聚类中心v_k的偏导数，并设得：

S922：根据步骤S921得聚类中心v_k为：

4.一种图像分割终端设备，其特征在于：包括处理器、存储器以及存储在所述存储器中并在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如权利要求1～3中任一所述方法的步骤。

5.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1～3中任一所述方法的步骤。