CN106296675A - 一种高噪声灰度不均匀图像的分割方法 - Google Patents

Abstract

本发明公开了一种高噪声灰度不均图像的分割方法，它是基于距离保持水平集方法和基于标记分水岭方法相结合的图像分割方法，该方法利用二者的优点弥补了相互的缺陷，实现了对同时含有高噪声和灰度不均匀目标的分割，能较准确地分割出感兴趣物体的边缘，能克服传统分水岭方法中存在的图像过度分割问题，且与传统水平集方法相比，运行速度较快，算法的运行时间不随图像尺寸大小的变化而变化。

Description

一种高噪声灰度不均匀图像的分割方法

技术领域

本发明涉及医学图像处理，尤其针对同时存在高噪声和灰度不均匀的医学图像进行分割的技术。

背景技术

图像分割和目标轮廓的提取对于图像理解、图像分析、模式识别、计算机视觉等具有重要意义。图像分割是图像处理到图像分析的关键，目的在于把图像分成各具特征的区域并提取感兴趣的目标过程。医学图像包括CT、MRI及其它医学影像设备所获得的图像。医学图像分割技术是医学图像处理和分析中的关键技术，它是一个根据区域内相似性和差异性把图像分割成若干区域的过程。

根据图像分割的不同特点可以将图像分割方法分成很多类。通常从分割技术的角度来进行分类，根据像素与其邻域间的相似性和不连续性可分成基于边缘的分割方法、基于区域的方法、基于特征空间的方法、基于能量泛函方法和其他方法。

基于边缘的分割方法通过检测不同区域的边缘来解决图像分割问题。通常所指的基于边缘的分割方法是基于灰度值的边缘检测，其边缘灰度值呈现出阶跃型变化的基础上，使用微分算子进行边缘检测，具体实现时可以使用图像与模板进行卷积来完成。

基于区域的分割方法是将图像按照相似性准则分成不同的区域，主要包括区域分裂合并方法、种子区域生长法。区域分裂合并方法从一个初始的包含许多非同质区域的分割开始进行分裂，然后再对分裂结果进行合并，直到获得符合某种要求的同质区域为止，分水岭变换是常用的一种分裂操作。区域生长法是将具有相似性质的像素点合并起来构成同质区域。

基于特征空间分割方法基于图像分类思想，将图像分割问题看成样本在特征空间的聚类问题，具有类似分布的样本将被划分为同一类别，从而实现图像的分割。基于特征空间的分割方法可划分为两大类：一类是直接采用各种聚类算法作用于特征空间，另一类是对各类样本点进行统计分析，再根据统计直方图来设置相应的分割阈值。

基于能量泛函的分割方法是指在活动轮廓模型基础上发展出来的算法，首先定义一个能量泛函使得其自变量包括边缘曲线和区域数据项，然后通过求解函数对应的欧拉-拉格朗日方程来实现，能量达到最小时的曲线位置就是目标的轮廓所在。活动轮廓模型可以分为参数活动轮廓模型和几何活动轮廓模型。

其他分割方法。现有的分割方法还包括基于物理模型的方法，基于模糊技术的方法，以及基于学习的分割方法等。

与本申请有关的技术是结合基于区域分割方法的分水岭变换和基于能量泛函的水平集方法两种分割方法，使用了改进的距离保持水平集演化方法。

分水岭变换的图像分割方法通常是对待分割图像的形态梯度信号进行分水岭变换，即将梯度图像看成是假想的地形表面，每个像素的梯度值代表该点的海拔高度，一般待分割图像的边缘像素通常有较大梯度值，对应于地表的“分水岭线”，而每个区域的内部像素通常有较小梯度值，对应于地表的“集水盆地”，“分水岭线”将梯度图像分割成若干个不同的“集水盆地”，得到待分割图像中的均匀一致区域。分水岭变换的图像分割方法虽然具有计算负担轻、分割精度高的优点，但是，存在的图像过度分割问题。

水平集方法基本思想是将二维闭合曲线演化问题转化为三维曲面问题来求解，在二维固定坐标系中不断更新水平集函数，从而达到演化隐含在水平集函数中闭合曲线的目的。即使隐含在水平集函数中的闭合曲线发生了拓扑结构的变化，水平集函数依旧保持有效。距离保持水平集方法在传统的能量泛函中添加了一项内部能量泛函，主要用来纠正水平集函数与符号距离函数的偏差，从而实现在演化过程中不需要再次初始化的目标，引入核函数来定义局部二值拟合能量，即寻找最优化轮廓曲线C使曲线两侧局部区域的图像灰度拟合函数近似为两个常数。水平集方法虽然分割效果高，但是分割速度慢，效率低。

本发明技术方案中涉及到现有距离保持水平集分割方法，为了便于理解，先详述现有的距离保持水平集分割方法。

首先需要定义一个能量泛函，使得其自变量包括边缘曲线和区域数据项，然后求解函数对应的Euler-Lagrange方程，当能量泛函取得最小值时的曲线位置就是目标的轮廓所在。

完整的距离保持水平集模型定义如下：

E(φ,f₁,f₂)＝ε(φ,f₁,f₂)+μP(φ)+vL(φ)

式中：

为弧长规则项，用来保证轮廓曲线C的光滑性；

为惩罚项，用来约束水平集函数，P(φ)是水平集函数的内部能量泛函，定量地表示了水平集函数偏离符号距离函数的程度，用来纠正水平集函数与符号距离函数的偏差；

ε(φ,f₁,f₂)+vL(φ)是外部能量泛函，记作：E_m(φ)＝ε(φ,f₁,f₂)+vL(φ),用来驱使零水平集向着图像中的目标边界运动；

μ和v为正的权重系数；

为了能够最小化能量函数E，采用梯度下降法和变分方法得到如下曲线演化方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\μ}(\▿{\mathrm{\φ}}^2-div(\frac{V\mathrm{\φ}}{\left|V\mathrm{\φ}\right|}\left)\right)+{\mathrm{v\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)div\left(\frac{V\mathrm{\φ}}{\left|V\mathrm{\φ}\right|}\right)-{\mathrm{\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)({\mathrm{\λ}}_1{e}_1-{\mathrm{\λ}}_2{e}_2)$

式中δ_ε(φ)是Dirac函数，e₁和e₂定义如下：

e₁＝∫_ΩK_σ(x-y)|I(y)-f₁(x)|²dy，e₂＝∫_ΩK_σ(x-y)|I(y)-f₂(x)|²dy

拟合函数f₁(x)和f₂(x)在迭代过程中按如下的方程式进行更新，其中*代表卷积运算：

$f_1\left(x\right)=\frac{K_{\mathrm{\σ}}{\left(x\right)}^{*}\[H\left(\mathrm{\φ}\right(x\left)\right)I\left(x\right)\]}{K_{\mathrm{\σ}}{\left(x\right)}^{*}H\left(\mathrm{\φ}\right(x\left)\right)}$

$f_2\left(x\right)=\frac{K_{\mathrm{\σ}}{\left(x\right)}^{*}\[(1-H(\mathrm{\φ}\left(x\right)\left)\right)I\left(x\right)\]}{K_{\mathrm{\σ}}{\left(x\right)}^{*}(1-H(\mathrm{\φ}\left(x\right)\left)\right)}$

由于权重系数v是一个常数，使得模型缺乏方向和大小的自适应性。

中国知网2014年1月公开了一篇《基于水平集和分水岭的图像分割算法研究》(作者：西安电子科技大学余敏锋)，该文献公开的技术是一篇硕士论文，申请人用该文献提到的技术对图像进行分割，根本就达不到文中提到的效果，在实际上它是根本不可能实现的，比如它公开的CT图像分割迭代达到1500次时时间仅为5.4s，分割迭代达到5000次时时间仅为9.6s，这在现实中是不可能的。另外，该文献对同时存在灰度不均匀和高噪声的医学图像分割效果不好，正如该文献最后一段记载的内容“分割效果不理想，对含有高噪声的目标分割效果不好”，正好印证了这一点。从它提供的实验结果图片也可以看出：它针对的主要是灰度均匀的图像。

发明内容

为了克服现有技术对同时含有高噪声和灰度不均匀的目标分割效果不理想的技术缺陷，满足医学图像分割的质量和效率要求，本发明提供一种基于水平集方法和基于标记分水岭方法相结合的一种高噪声灰度不均图像的分割方法。该方法利用二者的优点弥补了相互的缺陷，实现了对同时含有高噪声和灰度不均匀目标的分割。

为达到上述目的，本发明采取的技术方案是：

一种高噪声灰度不均图像的分割方法。其特征在于，它是基于距离保持水平集方法和基于标记分水岭方法相结合的图像分割方法，能较准确地分割出感兴趣物体的边缘，能克服传统分水岭方法中存在的图像过度分割问题，且与传统水平集方法相比，运行速度较快，算法的运行时间不随图像尺寸大小的变化而变化。具体描述如下：

第一步：采用自适应距离保持水平集方法对目标图像进行边缘的检测，得到初始目标图像的轮廓边缘。

本发明的自适应距离保持水平集方法是在现有的距离保持水平集分割方法引入了可变的权重系数v，可变权重系数v定义为：

式(1)中I(x,y)为图像函数，sgn(·)是符号函数，G_σ是标准差为σ的高斯滤波器，是图像经过高斯滤波后的梯度，ΔG_σ×I为Laplace算子作用于高斯滤波后的图像的结果，c>0是常数。

定义能量泛函如下：

式(2)中，λ>0和v为常数，δ(x)和H(x)分别是一维Dirac函数和Heaviside函数，g(x)为停止速度函数，g(x)＝(1+x²)^-1。

最小化能量泛函E(φ)，应用Euler-Lagrange方程，并引入人工变量t，运用最速下降法，即可得到偏微分方程：

其水平集形式为：

其中，为演化曲线的曲率。

在数值实现时，采用了如下的正则化Dirac函数来替代Dirac函数δ(x)：

式(3)作用于整个水平集函数，从而获得能量泛函的全局最小，进一步提高了零水平集检测多层轮廓的能力以及对深度凹陷区域与多目标物体边界的捕获能力。

本方法选取ε＝1.6，σ＝2.0，μ＝0.04，λ＝7.2，时间步长为0.5；

第二步：将第一步得到的目标图像轮廓边缘采用形态学灰度膨胀、腐蚀后，得到形态学梯度图像，将得到的形态学梯度图像作为基于标记的分水岭算法的内部标记。

所述的基于标记的分水岭算法中，形态学梯度是一种通过检测目标图像中某点的梯度值大小来确定该点是否存在轮廓边缘的方法。本发明需要利用结构元素b对目标图像f(x,y)分别做膨胀和腐蚀，求出f的局部极大值和局部极小值，其形态学梯度图像g(x,y)可表示为：

其中表示形态学的膨胀，表示形态学的腐蚀，b(x,y)表示圆盘状结构元素。

第三步：将第二步得到的结果使用MATLAB函数im2bw和graythresh将目标图像转变为二值图像，然后求补，求出距离变换，然后利用函数watershed来计算距离变换的负分水岭分割，将得到的结果作为标记分水岭算法的外部标记，最后得到待处理的图像；

使用距离变换的分水岭变换分割二值图像的伪代码如下：

>>g＝im2bw(f,graythresh(f))；//把图像f变换为二值图像

L＝watershed(A,conn)；//计算距离变换的负分水岭变换，A是输入数组，conn指定了连通性(对于二维数组默认值为8)

>>gc＝～g；//求补

>>D＝bwdist(gc)；//工具箱函数bwdist，计算距离变换

>>L＝watershed(-D)；

>>w＝L＝＝0；//L的零值像素指出分水岭的脊线像素

>>g2＝g&-w；//使用原始的二值图像和图像w的补，通过逻辑&操作完成分割；

第四步：利用形态学重建，对待处理的图像进行处理，详细的处理过程是：

对于形态学梯度图像来说，不能彻底去除噪声，需要利用形态学重建进一步消除目标图像中含有的噪声以及那些无法包含结构元素的像素点。形态学开重建和闭重建运算是通过结合形态学膨胀和腐蚀两个运算方法形成的。对于形态学梯度图像g(x，y)、比对图像p(x，y)与结构元素b，形态学膨胀定义为：

其中b为圆盘状结构元素，∧表示逐点求最小值，形态学膨胀为迭代性运算。当迭代的次数达到预设值或者时，中止。因此，形态学开重建定义为：

式(8)中表示测地学膨胀收敛的结果。

形态学开、闭重建是互为对偶的。因此，形态学腐蚀及闭重建可定义为：

其中，V表示的是逐点求最大值，表示形态学腐蚀时的结果。

因此，形态学重建定义为：

第五步：对形态学重建后的结果使用Watershed函数进行基于标记的分水岭分割，输出图像分割结果：详细分割过程是：

第5.1步：利用扩展最小变换(H-minima)技术，首先给出一个图像阈值H，通过与之对比消去局部区域小于H的极小值点，以此减少过分割区域的数目。本方法采用所有极小值的平均值来修正H值的大小。假设中极小值为则

其中mean(·)为平均值函数；为图像梯度值；

第5.2步：使用Gaussian滤波器获得低频图像

其中，Hmin为扩展最小变换，为二值标记图像，也就是经扩展最小变换后的图像。

第5.3步：使用极小值强制标记运算修改图像梯度值即可得到可表示为：

第5.4步：最后对进行分水岭分割，得到的分割结果用I_w表示为：

其中，Watershed表示分水岭算法的MATLAB函数。

本发明使用基于标记的分水岭分割图像算法的伪代码如下：

>>rm＝imregionalmin(g)；//工具箱函数imregionalmin计算图像中所有的局部小区域位置，其中g是灰度图像，rm是二值图像，rm的前景像素标记出局部小区域的位置

>>im＝imextendedmin(f,h)；//工具箱函数imextendedmin消除无关小区域，其中f是灰度图像，h是高度阈值，im是二值图像，im的前景像素标记了深的局部小区域的位置，用函数imextendedmin得到内部标记符集合

>>fim＝f；

>>Lim＝watershed(bwdist(im))；//

>>em＝Lim＝＝0；//

>>g2＝imimposemin(g,im|em)；//em的前景像素标记了输出图像中局部最小区域的期望位置。通过在内部和外部标记符的位置覆盖局部最小区域，改进梯度图像

>>L2＝watershed(g2)；

>>f2＝f；

>>f2(L2＝＝0)＝255；//计算改进后基于标记的梯度图像的分水岭变换

本方法使用MATLAB 2014b进行实验，实验所用的计算机配置为：在OS X EICapitan操作系统下CPU为2.7GHz Intel Core i5，内存为8GB。

本发明的优点是：

1、本发明采用的自适应距离保持水平集方法，能使得模型具有了方向和大小的自适应性，提高了图像分割的效率，克服了传统水平集分割方法速度慢的缺陷。

2、本发明将形态学重建与扩展最小变化技术结合，避免了图像的过分割和欠分割问题，提高了图像的分割质量。

3、本发明将自适应距离保持水平集方法和形态学重建的膨胀和腐蚀结合，既消除了图像中噪声，又明确了细节，减少了因目标图像的噪声、明暗细节等引起的分水脊线偏移，从而避免出现过分割的问题，使得最后的分割效果好。

附图说明

图1是本发明流程图。

图2.1是待分割的人工合成图像；图2.2是图2.1的分割实验结果图。

图3.1是待分割的MRI脑部图像；图3.2是图3.1的分割实验结果图。

图4.1是待分割的血管图；图4.2是图4.1的分割实验结果图。

下面结合附图进一步说明本发明的技术方案

具体实施方式

图1是本发明对于高噪声灰度不均匀图像的分割方法流程图，具体的分割过程在发明内容已经详述，在此不再详述。

为了说明本发明的积极效果，下面分别以高噪声且灰度不均匀的人工合成图像、MRI脑部图像和血管图为例，在同样配置的电脑上分别进行控制标记符的分水岭图像分割实验、距离保持水平集演化图像分割实验。为了定量分析和评价各种方法的分割质量，本发明采用Jaccard相似度(Jaccard Similarity，JS)系数来度量图像分割方法的精度。Jaccard相似度系数定义如下：

JS(A,B)＝|A∩B|/|A∪B|

其中A为标准的分割结果，B为待评价的分割结果。JS指标值越高则表明算法的性能越好。三种分割方法的迭代次数和运行时间的实验数据对比见表1。三种分割方法的精度比较见表2。对比实验前后图分别见图2-4

表1：迭代次数和运行时间对比表

表2：精度比较

分割精度(JS系数)	人工合成图像	MRI脑部图像	血管图像
				距离保持水平集方法	0.9092	0.7538	0.9142
分水岭方法	0.6342	0.7868	0.8017
				本发明方法	0.9113	0.9128	0.9387

从表一和表二可以看出，本发明提出的基于水平集和分水岭方法的图像分割算法分割速度更快，分割效果更好。

从图2-4比较图可以看出：分割前，目标图像的灰度不均匀，噪声严重，但是经过本发明分割后，图像轮廓光滑，尤其是分割后的血管图，主血管附近的小血管都能看得很清晰。

Claims (1)

1.一种高噪声灰度不均匀图像的分割方法，其特征在于，它是基于距离保持水平集方法和基于标记分水岭方法相结合的图像分割方法，具体步骤如下：

第一步：采用自适应距离保持水平集方法对目标图像进行边缘的检测，得到初始目标图像的轮廓边缘；

所述的自适应距离保持水平集方法是在距离保持水平集分割方法引入了可变的权重系数v，可变权重系数v定义为：

式(1)中I(x,y)为图像函数，sgn(·)是符号函数，G_σ是标准差为σ的高斯滤波器，是图像经过高斯滤波后的梯度，ΔG_σ×I为Laplace算子作用于高斯滤波后的图像的结果，c>0是常数；

定义能量泛函如下：

式(2)中，λ>0和v为常数，δ(x)和H(x)分别是一维Dirac函数和Heaviside函数，g(x)为停止速度函数，g(x)＝(1+x²)^-1；

最小化能量泛函E(φ)，应用Euler-Lagrange方程，并引入人工变量t，运用最速下降法，即可得到偏微分方程：

其水平集形式为：

其中，为演化曲线的曲率；

在数值实现时，采用了如下的正则化Dirac函数来替代Dirac函数δ(x)：

式(3)作用于整个水平集函数，从而获得能量泛函的全局最小；

上述ε＝1.6；σ＝2.0，μ＝0.04，λ＝7.2，时间步长为0.5；

第二步：将第一步得到的目标图像轮廓边缘采用形态学灰度膨胀、腐蚀后，得到形态学梯度图像，将得到的形态学梯度图像作为基于标记的分水岭算法的内部标记；具体为：

利用结构元素b对目标图像f(x,y)分别做膨胀和腐蚀，求出f的局部极大值和局部极小值，其形态学梯度图像g(x,y)可表示为：

其中⊕表示形态学的膨胀，表示形态学的腐蚀，b(x,y)表示圆盘状结构元素；

第三步：将第二步得到的结果使用MATLAB函数im2bw和graythresh将目标图像转变为二值图像，然后求补，求出距离变换，然后利用函数watershed来计算距离变换的负分水岭分割，将得到的结果作为标记分水岭算法的外部标记，最后得到待处理的图像；

第四步：利用形态学重建，对待处理的图像进行处理，详细的处理过程是：对于形态学梯度图像g(x,y)、比对图像p(x,y)与结构元素b，形态学膨胀定义为：

其中b为圆盘状结构元素，∧表示逐点求最小值，形态学膨胀为迭代性运算。当迭代的次数达到预设值或者时，中止。因此，形态学开重建定义为：

式(8)中表示测地学膨胀收敛的结果；

形态学开、闭重建是互为对偶的，因此，形态学腐蚀及闭重建表示为：

其中，∨表示的是逐点求最大值，表示形态学腐蚀时的结果；

因此，形态学重建表示为：

第五步：对形态学重建后的结果使用Watershed函数进行基于标记的分水岭分割，输出图像分割结果：详细分割过程是：

第5.1步：利用扩展最小变换H-minima技术，首先给出一个图像阈值H，通过与之对比消去局部区域小于H的极小值点，以此减少过分割区域的数目。本方法采用所有极小值的平均值来修正H值的大小。假设中极小值为则

其中mean(·)为平均值函数；为图像梯度值；

第5.2步：使用Gaussian滤波器获得低频图像

其中，H min为扩展最小变换，为二值标记图像，也就是经扩展最小变换后的图像；

第5.3步：使用极小值强制标记运算修改图像梯度值即可得到表示为：

第5.4步：最后对进行分水岭分割，得到的分割结果用I_w表示为：

其中，Watershed表示分水岭算法的MATLAB函数。