一种基于局部灰度聚类特征的图像水平集分割方法与流程

本发明属于计算机视觉、模式识别、图像处理与分析的技术领域，具体涉及一种基于局部灰度聚类特征的图像水平集分割方法。

背景技术：

图像分割在计算机视觉、模式识别、图像处理领域起着重要作用。其目的是将图像划分为多个感兴趣区域，每个区域内的像素组成一个具有特定形状和结构特征的感兴趣对象。虽然被研究多年且有多种分割方法被提出，从存在噪声和灰度偏移场的结构复杂的图像中分割感兴趣对象仍然是一个挑战。水平集分割方法将图像分割转化为定义在水平集函数上的能量泛函的能量最小化问题。水平集分割方法能够提供亚像素级的平滑封闭的分割曲线，且可以将形状或灰度分布的先验知识包含在能量公式中，从而增强分割方法的健壮性。因此，这类方法在图像分割领域具有明显优势。一般地，水平集分割方法用比水平集函数低一维的零水平集曲线表示分割对象的边界。零水平集曲线随着水平集函数进行演化，直到曲线正好运动到分割对象边界时，停止演化。此时，以水平集函数为自变量的能量泛函取最小值，反之亦然。但是，传统的水平集方法不具备应对影像灰度不一致性问题的能力。

技术实现要素：

针对现有技术的不足，本发明提出一种基于局部灰度聚类特征的图像水平集分割方法。

本发明的技术方案是：

一种基于局部灰度聚类特征的图像水平集分割方法，包括以下步骤：

步骤1：读取待分割图像I；

步骤2：采用M个正交基函数的线性加权和拟合偏移场b，并初始化各基函数的权重值；

所述偏移场b的公式如下所示：

b＝w^TG；

其中，w＝(w₁，w₂，...，w_M)^T为基函数权重列向量，G＝(g₁，g₂，...，g_M)^T为M个基函数构成的列向量，g₁，g₂，...，g_M为两两正交的基函数，w₁，w₂，...，w_M为初始化的各基函数的权重值；

步骤3：初始化图像的水平集函数集合：根据待分割图像I待划分的区域数N，设定初始的k个零水平集分割曲线，根据设定的零水平集分割曲线生成高一维的水平集函数集合Φ＝(φ₁，φ₂，...，φ_k)；

步骤4：建立图像水平集分割的能量泛函F(Φ，w，c)，根据待分割图像设定水平集分割控制参数：图像数据项权重系数λ_n、水平集函数正则项系数μ、水平集分割曲线平滑项系数v；

所述能量泛函如下所示：

$F(\mathrm{\Φ},w,c)=D(\mathrm{\Φ},w,c)+\mathrm{\μ}\underset{i=1}{\overset{k}{\Σ}}R\left({\mathrm{\φ}}_i\right)+v\underset{i=1}{\overset{k}{\Σ}}L\left({\mathrm{\φ}}_i\right);$

其中，为图像数据项，n∈{1，2，...，N}，1为函数值全为1的常函数，*为卷积操作符号，K为归一化的非负偶函数，c＝(c₁，c₂，...，c_N)为聚簇中心集合，M_n(Φ(x))为第n个子区域Ω_n的隶属度函数，Φ(x)为x像素点处k个水平集函数值构成的集合；为水平集函数正则项，为φ_i(x)的梯度，为水平集分割曲线平滑项，δ(φ_i(x))为φ_i(x)的Dirac函数，为φ_i(x)的Heaviside函数；

步骤5：设定最大迭代次数J和图像水平集分割的能量泛函F(Φ，w，c)的收敛阈值；

步骤6：分别更新聚簇中心集合c^j、图像水平集函数集合Φ^j、基函数权重列向量w^j，得到第j次迭代的能量泛函F^j(Φ，w，c)；

步骤7：判断是否达到迭代终止条件：若当前迭代次数达到最大迭代次数或者F^j(Φ，w，c)与F^j-1(Φ，w，c)的差值小于图像水平集分割的能量泛函F(Φ，w，c)的收敛阈值，则执行步骤8，否则，迭代次数加1，返回步骤6；

步骤8：根据当前更新后的图像水平集函数集合构造图像的隶属度函数即待分割图像I的分割结果，并根据更新后的基函数权重列向量w^j与基函数列向量G得到待分割图像I的偏移场估计。

可选地，所述设定的初始零水平集分割曲线可以是指定形状的闭合曲线或是随机生成的闭合曲线，所述根据设定的零水平集分割曲线生成高一维的水平集函数φ₁，φ₂，...，φ_k为阶梯函数或符号距离函数。

可选地，所述更新聚簇中心集合c^j、图像水平集函数集合Φ^j、基函数权重列向量w^j的具体方法为：

所述更新聚簇中心集合c^j的公式如下所示：

$c_n^j=\frac{\∫I\left(x\right)M_n\left({\mathrm{\Φ}}^{j-1}\right(x\left)\right)(K*({\left(w^{j-1}\right)}^{T}G\left)\right)\left(x\right)dx}{\∫M_n\left({\mathrm{\Φ}}^{j-1}\right(x\left)\right)(K*{\left({\left(w^{j-1}\right)}^{T}G\right)}^2)\left(x\right)dx};$

所述更新图像水平集函数集合Φ^j的公式如下所示：

$\frac{\∂{\mathrm{\φ}}_n^j}{\∂t}=-\underset{n=1}{\overset{N}{\Σ}}\frac{\∂M_n\left({\mathrm{\Φ}}^{j-1}\right)}{\∂{\mathrm{\φ}}_n^{j-1}}{\mathrm{\λ}}_n{e}_n^{j-1}+\mathrm{\μ}\left({\▿}^2{\mathrm{\φ}}_n^{j-1}-div\left(\frac{\▿{\mathrm{\φ}}_n^{j-1}}{|\▿{\mathrm{\φ}}_n^{j-1}|}\right)\right)+v\mathrm{\δ}\left({\mathrm{\φ}}_n^{j-1}\right)div\left(\frac{\▿{\mathrm{\φ}}_n^{j-1}}{|\▿{\mathrm{\φ}}_n^{j-1}|}\right);$

其中，y为像素点x周围一个半径为ρ的足够小的圆形邻域O_x内的像素点，G(y)＝(g₁(y)，g₂(y)，...，g_M(y))^T；

所述更新基函数权重列向量w^j的公式如下所示：

$w^j={\{\∫\[K*\left(\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_n{\left(c_n^{j-1}\right)}^{2}M\left({\mathrm{\Φ}}^{j-1}\right)\right)\]\left(y\right)G\left(y\right)G^T\left(y\right)dy\}}^{-1}\∫\left(K*\left(I\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_n{c}_n^{j-1}M\left({\mathrm{\Φ}}^{j-1}\right)\right)\right)\left(y\right)G\left(y\right)dy.$

可选地，所述根据当前更新后的图像水平集函数集合构造图像的隶属度函数的具体方法为：

由水平集函数的Heaviside变换的逻辑运算得出，且若待分割图像I中的像素点x属于待划分的N个区域中的第n个子区域Ω_n中，则隶属度函数否则

本发明的有益效果：

本发明提出一种基于局部灰度聚类特征的图像水平集分割方法，该方法克服了弱边界、图像噪声、以及灰度不一致性(偏移场)对图像分割精度的负面影响，得到精确分割结果的同时，能够估计原图像的灰度不一致性，进而起到校正图像灰度的作用。

附图说明

图1为实施例1中基于局部灰度聚类特征的图像水平集分割方法的流程图；

图2为实施例1中设定的2个初始零水平集分割函数及初始分割结果；

其中，(a)为实施例1中设定的2个初始零水平集分割函数；

(b)为实施例1中初始分割结果；

图3为实施例1中迭代2次后的水平集函数及分割结果；

其中，(a)为实施例1中迭代2次后的水平集函数；

(b)为实施例1迭代2次后的分割结果；

图4为本实施例1中迭代5次后的水平集函数及分割结果；

其中，(a)为实施例1中迭代5次后的水平集函数；

(b)为实施例1中迭代5次后的分割结果；

图5为实施例1中迭代10次后的水平集函数及分割结果；

其中，(a)为实施例1中迭代10次后的水平集函数；

(b)为实施例1迭代10次后的分割结果；

图6为实施例1中迭代50次后的水平集函数及分割结果；

其中，(a)为实施例1中迭代50次后的水平集函数；

(b)为实施例1中迭代50次后的分割结果；

图7为实施例1中得到待分割图像I的偏移场估计、灰度不一致性得到校正的图像、以及校正后图像的灰度直方图；

其中，(a)为实施例1中得到待分割图像I的偏移场估计；

(b)为实施例1中得到的灰度不一致校正后的图像；

(c)为本实施例1中校正后图像的灰度直方图；

图8为实施例2中k＝1的四幅待分割图像的初始零水平集分割曲线及其分割结果；

其中，(a)为实施例2中k＝1的四幅待分割图像及初始零水平集分割曲线；

(b)为实施例2中k＝1的四幅待分割图像的分割结果；

图9为实施例3中存在噪声、偏移场、弱边界信息的真实图像进行分割后的初始零水平集分割曲线、得到的偏移场估计、偏移场校正后的图像和分割结果；

其中，(a)为实施例3中四幅待分割图像及初始零水平集分割曲线；

(b)为实施例3中四幅待分割图像得到的偏移场估计；

(c)为实施例3中四幅偏移场校正后的图像；

(d)为实施例3中四幅待分割图像的分割结果。

具体实施方式

下面结合附图对本发明具体实施方式加以详细说明。

实施例1：

一种基于局部灰度聚类特征的图像水平集分割方法，如图1所示，包括以下步骤：

步骤1：读取待分割图像I。

步骤2：采用M个正交基函数的线性加权和拟合偏移场b，并初始化各基函数的权重值。

本实施方式中，M为15，偏移场b(x)的公式如式(1)所示：

b＝w^TG (1)

其中，w＝(w₁，w₂，...，w_M)^T为基函数权重列向量，T为转置操作符号，G＝(g₁，g₂，...，g_M)^T为M个基函数构成的列向量，g₁，g₂，...，g_M为两两正交的4阶Legendre多项式函数，w₁＝w₂＝...＝w_M＝1为初始化的各基函数的权重值。

步骤3：初始化图像的水平集函数集合：根据待分割图像I待划分的区域数N，设定初始的k个零水平集分割曲线，根据设定的零水平集分割曲线生成高一维的水平集函数集合Φ＝(φ₁，φ₂，...，φ_k)。

本实施方式中，根据待分割图像I待划分的区域数N＝3，设定初始零水平集分割曲线的个数k的计算公式如式(2)所示：

k＝upper(log₂N) (2)

其中，upper为上整数界符号，本实施方式中，由式(2)可得k＝2。

本实施方式中，设定的初始零水平集分割曲线为闭合圆，根据设定的零水平集分割曲线生成高一维的水平集函数φ₁，φ₂分别为闭合圆为边界的阶梯函数，如图2所示。

步骤4：建立图像水平集分割的能量泛函F(Φ，w，c)，根据待分割图像设定图像水平集分割的控制参数：图像数据项权重系数λ_n、水平集函数正则项系数μ、水平集分割曲线平滑项系数v。

本实施方式中，能量泛函的公式如式(3)所示：

$F(\mathrm{\Φ},w,c)=D(\mathrm{\Φ},w,c)+\mathrm{\μ}\underset{i=1}{\overset{k}{\Σ}}R\left({\mathrm{\φ}}_i\right)+v\underset{i=1}{\overset{k}{\Σ}}L\left({\mathrm{\φ}}_i\right)---\left(3\right)$

其中，为图像数据项，n∈{1，2，3}，1为函数值全为1的常函数，*为卷积操作符号，K为归一化的非负偶函数，c＝(c₁，c₂，c₃)为聚簇中心集合，M_n(Φ(x))为第n个子区域Ω_n的隶属度函数，Φ(x)为x像素点处k个水平集函数值构成的集合；为水平集函数正则项，为φ_i(x)的梯度，为水平集分割曲线平滑项，δ(φ_i(x))为φ_i(x)的Dirac函数，为φ_i(x)的Heaviside函数。

本实施方式中，图像数据项D(Φ，w，c)用于将零水平集分割曲线演化吸引到待分割对象边界位置，水平集函数正则项R(φ_i)用于确保水平集函数在演化中保持特定的平滑属性，水平集分割曲线平滑项L(φ_i)用于确保演化过程中分割曲线的平滑性。

本实施方式中，图像数据项系数越大，零水平集分割曲线越切近等分割对象在图像中表现出的边界，水平集函数正则项系数越大水平集函数越接近于符号距离函数，水平集分割曲线平滑项系数越大得到的分割曲线越光滑，因此根据待分割图像设定图像水平集分割的控制参数：图像数据项权重系数λ₁＝1.2，λ₂＝1.0，λ₃＝0.7，水平集函数正则项系数μ＝1.0，水平集分割曲线平滑项系数v＝0.005*255*255。

步骤5：设定最大迭代次数J和图像水平集分割的能量泛函F(Φ，w，c)的收敛阈值。

本实施方式中，设定最大迭代次数J为50，图像水平集分割的能量泛函F(Φ，w，c)的收敛阈值为0.01。

步骤6：分别更新聚簇中心集合c^j、图像水平集函数集合Φ^j、基函数权重列向量w^j，得到第j次迭代的能量泛函F^j(Φ，w，c)。

本实施方式中，当零水平集分割曲线正好位于待分割对象的边界时，公式(3)给出的能量泛函取得最小值，因此，为使得取得最小值，保持Φ与w不变，解即更新聚簇中心集合c^j的公式如式(4)所示：

$c_n^j=\frac{\∫I\left(x\right)M_n\left({\mathrm{\Φ}}^{j-1}\right(x\left)\right)(K*({\left(w^{j-1}\right)}^{T}G\left)\right)\left(x\right)dx}{\∫M_n\left({\mathrm{\Φ}}^{j-1}\right(x\left)\right)(K*{\left({\left(w^{j-1}\right)}^{T}G\right)}^2)\left(x\right)dx}---\left(4\right)$

本实施方式中，保持w与c不变，采用最大梯度降方法求得更新后的图像水平集函数集合Φ^j的公式如式(5)所示：

$\frac{\∂{\mathrm{\φ}}_n^j}{\∂t}=-\underset{n=1}{\overset{N}{\Σ}}\frac{\∂M_n\left({\mathrm{\Φ}}^{j-1}\right)}{\∂{\mathrm{\φ}}_n^{j-1}}{\mathrm{\λ}}_n{e}_n^{j-1}+\mathrm{\μ}\left({\▿}^2{\mathrm{\φ}}_n^{j-1}-div\left(\frac{\▿{\mathrm{\φ}}_n^{j-1}}{|\▿{\mathrm{\φ}}_n^{j-1}|}\right)\right)+v\mathrm{\δ}\left({\mathrm{\φ}}_n^{j-1}\right)div\left(\frac{\▿{\mathrm{\φ}}_n^{j-1}}{|\▿{\mathrm{\φ}}_n^{j-1}|}\right)---\left(5\right)$

其中，y为像素点x周围一个半径为ρ的足够小的圆形邻域O_x内的像素点，G(y)＝(g₁(y)，g₂(y)，...，g_M(y))^T。

本实施方式中，保持Φ与c不变，解得更新后的基函数权重列向量w^j的公式如式(6)所示：

$w^j={\{\∫\[K*\left(\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_n{\left(c_n^{j-1}\right)}^{2}M\left({\mathrm{\Φ}}^{j-1}\right)\right)\]\left(y\right)G\left(y\right)G^T\left(y\right)dy\}}^{-1}\∫\left(K*\left(I\underset{n=1}{\overset{N}{\Σ}}{\mathrm{\λ}}_n{c}_n^{j-1}M\left({\mathrm{\Φ}}^{j-1}\right)\right)\right)\left(y\right)G\left(y\right)dy---\left(6\right)$

步骤7：判断是否达到迭代终止条件：若当前迭代次数达到最大迭代次数或者F^j(Φ，w，c)与F^j-1(Φ，w，c)的差值小于图像水平集分割的能量泛函F(Φ，w，c)的收敛阈值，则执行步骤8，否则，迭代次数加1，返回步骤6。

本实施方式中，迭代2次、5次、10次和50次后的水平集函数及分割结果如图3、图4、图5、图6所示。

步骤8：根据当前更新后的图像水平集函数集合构造图像的隶属度函数即待分割图像I的分割结果，并根据更新后的基函数权重列向量w^j与基函数列向量G得到待分割图像I的偏移场估计。

本实施方式中，由水平集函数的Heaviside变换的逻辑运算得出，且待分割图像I中的像素点x属于待划分的N个区域中的第n个子区域Ω_n中，则隶属度函数否则

本实施方式中，得到待分割图像I的偏移场估计、灰度不一致性得到校正的图像、以及校正后图像的灰度直方图如图7所示。

实施例2：

本发明最简单的情况，即将图像分割为两个区域N＝2，根据公式(2)计算得到k＝1。

如图8所示情况，(a)为四幅待分割图像及初始零水平集分割曲线，(b)为四幅待分割图像的分割结果，其中，四幅待分割图像的水平集分割曲线平滑项系数v依次取值为0.1×255×255，0.5×255×255，0.3×255×255，0.02×255×255。

实施例3：

本发明方法对存在噪声、偏移场、弱边界信息的真实图像的进行分割，将图像分割为两个区域N＝2，根据公式(2)计算得到k＝1。如图9所示，(a)为四幅待分割图像及初始零水平集分割曲线，(b)为本发明方法得到的四幅待分割图像的偏移场估计，(c)为四幅偏移场校正后的图像，(d)为四幅待分割图像的分割结果，其中，四幅待分割图像的水平集分割曲线平滑项系数v依次取值0.0045×255×255，0.003×255×255，0.005×255×255，0.01×255×255。