CN102831608A - 基于非稳态测量算法的改进规则距离水平集的图像分割方法 - Google Patents

Abstract

基于非稳态测量算法的改进规则距离水平集的图像分割方法，它属于数字图像处理领域，本发明的目的是为了更精确的分割高斯噪声和椒盐噪声干扰的图像。方法：一、构造均值算子，用构造的算子设置停止函数；二、手动设置初始轮廓，根据该轮廓初始化水平集函数；三、将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果。本发明的优点是比传统几何Snake模型更精确的分割高斯噪声和椒盐噪声干扰的图像。

Description

基于非稳态测量算法的改进规则距离水平集的图像分割方法

技术领域

本发明涉及基于Snake模型的图像分割方法，属于数字图像处理领域。

背景技术

Snake模型，又称为活动轮廓模型、蛇模型。活动轮廓的含义是首先在待分割目标周围手动或自动设置初始轮廓，并赋予初始轮廓能量。轮廓在模型自身的拓扑内力和图像的灰度数据产生的外力的作用下发生形变，像蛇一样活动，达到轮廓能量最小时，轮廓的所在位置即是待分割目标的边缘。

Snake模型是将图像分割问题转化为求能量泛函最小化的过程。活动轮廓突出优点是一旦设置初始轮廓，后续的轮廓演化不需要人为参与，自动化程度高，而采用力学原理，对轮廓本身定义了弹性力和刚性力，这样，无论分割目标的结构如何复杂，分割的最终轮廓都是光滑封闭的，这符合一般情况下的自然图像，尤其是医学图像的客观规律。Snake模型因计算效率高，适合建模，而被广泛应用于图像分割、目标跟踪、模式识别等应用。

Snake模型根据轮廓的表达方式的不同分为两种：参数Snake模型和几何Snake模型。几何Snake模型的轮廓曲线是由高一维曲面函数的零水平集表示的。几何Snake模型基于曲线演化理论和水平集方法(Level Set)采用三维曲面的零水平集代表二维轮廓曲线，通过Hamilton-Jacobi方程求解三维曲面演化过程来获得最终轮廓。水平集方法基于波前熵守恒理论，通过高一维的处理，避免了曲线演化过程中出现的奇异性，能够很容易实现曲线的融合、断裂等拓扑变化。

V.Casellse等人首先将水平集方法应用于活动轮廓模型中，在传统的Snake模型中，用零水平集函数代替活动轮廓，并将力平衡方程改写为曲线的演化方程，提出了著名的测地线活动轮廓模型(GAC)。

设活动轮廓C(s)，曲线的能量表达式

$E\left(C\left(s\right)\right)={\∫}_0^1\mathrm{\α}{\left|C{\left(s\right)}^{\′}\right|}^2+\mathrm{\β}{\left|C{\left(s\right)}^{\′\′}\right|}^2\mathrm{ds}-\mathrm{\λ}{\∫}_0^1|\▿I\left(C\left(s\right)\right)|\mathrm{ds}---\left(1\right)$

I表示为图像，α、β分别为轮廓的弹性系数和刚性系数，控制轮廓的光滑性，一般为常数。λ为常数，控制图像外能的权重。为梯度算子。由于水平集函数曲线本身具有较为光滑的特性，且曲线演化主要在曲线法线方向上的伸缩，令β＝0，同时替换外能为g(I)，满足

g(I)→0，保证在边缘处获得最小值。g为停止函数，控制水平集函数的演化和停止，一般采用如下定义：

$g\left(I\right)=\frac{1}{1+{|{\▿G}_{\mathrm{\σ}}*I|}^2}$

带入式(1)得到曲线的能量泛函，

$E\left(C\left(s\right)\right)={\∫}_0^1\mathrm{\α}{\left|C{\left(s\right)}^{\′}\right|}^2\mathrm{ds}-\mathrm{\λ}{\∫}_0^{1}g{\left(\right|\▿I\left(C\left(s\right)\right)\left|\right)}^2\mathrm{ds}$

根据Maupertuis’定律，最小化能量泛函，转化为求

$\mathrm{Min}\underset{0}{\overset{1}{\∫}}g\left(\right|\▿I\left(C\left(s\right)\right)\left|\right)\left|C^{\′}\left(s\right)\right|\mathrm{ds}$

曲线的弧长为可以发现上式表达了对曲线进行加权后求和。求最小问题即求轮廓的位置使g(I)的加权和最小，显然是轮廓在目标的边缘的时候值最小。

上述最小化问题是黎曼空间的测地线问题，根据最速下降法，得到曲线的演化方程：

$\frac{\∂C(s,t)}{\∂t}=(g\left(I\right)\mathrm{\κ}-\▿g\·N-c)N$

其中

κ为轮廓上各点的曲率，N为轮廓上各点的法向量，c为常数，g为定义的停止函数，在边缘处的值最小。

定义水平集函数φ(x，y，t)，(x，y)表示为轮廓上点的坐标位置，则轮廓表示为水平集函数的零水平集{C：φ(x，y，t)＝0}，将之代入曲线的演化方程，得到水平集函数的演化方程：

$\frac{\∂\mathrm{\φ}}{\∂t}=(g\·\mathrm{div}\left(\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}\right)+\▿g\·\frac{\▿\mathrm{\φ}}{|\▿\mathrm{\φ}|}+c)\▿\mathrm{\φ}$

式中div(·)为散度算子。

GAC模型存在的主要问题是，随着水平集演化的进行，水平集函数逐渐偏离初始的符号距离函数，导致

水平集演化不稳定，出现错误结果。传统的解决办法是在演化过程中随时对水平集函数进行重新初始化。

重新初始化难以确定需要重新初始化的时间，且计算量大。C.Li等人针对这一问题，提出了规则距离水平集模型(DRLSE)，通过在轮廓能量中引入规则距离项，自动调整水平集函数与符号距离的偏移，保持规则距离项如下：

是关于水平集函数变化率的势能函数，

为补偿项，在Snake模型的能量函数中充当内能的角色，调节水平集函数始终保持规则距离即

因此称为规则距离项。零水平集到达目标边缘时，轮廓能量达到最小值，这要求内能在

获得最小值才能稳定。根据极值定理，要求势函数

的选择为当变量满足

的时候达到最小值，并且为0。为使

满足条件，定义了如下的函数：

$p\left(s\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\π}\right)}^2}(1-\cos \left(2\mathrm{\πs}\right))& s\≤1\\ \frac{1}{2}{(s-1)}^2& s\≥1\end{array}\right.$

DRLSE模型是基于水平集方法的一种新型图像分割算法，它相对于GAC模型，具有无须重新初始化的优点，但还存在对噪声敏感的问题。对于受高斯和椒盐噪声干扰的图像，DRLSE模型的分割准确度低，不能满足分割要求。

发明内容

本发明的目的是为了更精确的分割高斯噪声和椒盐噪声干扰的图像，而提供了一种基于非稳态测量算法的改进规则距离水平集图像分割方法。

本发明基于非稳态测量算法的改进规则距离水平集图像分割方法，实现该方法的步骤如下：

一、构造均值算子，用构造的算子设置停止函数；

二、手动设置初始轮廓，根据该轮廓初始化水平集函数；

三、将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法求解能量函数最小化问题，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果。

其中步骤一所述构造均值算子，用构造的算子设置停止函数的方法为：：

a、构造两个边长为L的方形矩阵，且矩阵内的元素值均为1/L²，即均值滤波器h₁和h₂，其中L一般取3，5或9。

b、用构造的滤波器h₁，h₂，采用公式 $g=\frac{1}{1+k|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2|}$ 设置停止函数g，其中k为常数，采用公式 $k=\max \left({|{\▿G}_{\mathrm{\σ}}*I|}^2\right)\frac{{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2}{\max \left[\right|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2\left|\right]}$ 获得，I为图像，max(·)表示求最大值，*表示卷积，G_σ是均值为0标准差为σ的高斯滤波器，其中σ的取值范围一般为0.1～0.3。

其中步骤二手动设置初始轮廓R₀，根据该轮廓位置采用公式

初始化水平集函数，c₀为常数，一般取c₀＝5。

其中步骤三将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果的方法为：

a、将步骤一的停止函数g带入到DRLSE的能量方程，得到基于非稳态测量算法的改进规则距离水平集的能量方程：

${\mathrm{\ϵ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)=\mathrm{\μ}\underset{\mathrm{\Ω}}{\∫}p\left(\right|\▿\mathrm{\φ}\left|\right)\mathrm{dx}+\mathrm{\λ}\underset{\mathrm{\Ω}}{\∫}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\δ}}_{\mathrm{\ω}}\left(\mathrm{\φ}\right)|\▿\mathrm{\φ}|\mathrm{dx}$

$+\mathrm{\α}\underset{\mathrm{\Ω}}{\∫}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{H}_{\mathrm{\ω}}(-\mathrm{\φ})\mathrm{dx}$

式中μ、λ、α为常数，一般μ选择为0.2，λ选择为3～5，α控制轮廓演化的方向，α＜0轮廓膨胀，α＞0轮廓收缩，α值根据边缘的强弱一般选取1～5之间的某个值，弱边缘则选较小的值，强边缘选较大的值，式中势能函数p定义为：

$p\left(s\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\π}\right)}^2}(1-\cos \left(2\mathrm{\πs}\right))& s\≤1\\ \frac{1}{2}{(s-1)}^2& s\≥1\end{array}\right.$

Heaviside函数H_ω(x)和Dirac函数δ_ω(x)在实际中分别采用如下近似表达，其中ω的取值范围一般为0.1～2.5：

$H_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}(1+\frac{x}{\mathrm{\&omega;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right))& \left|x\right|\&le;\mathrm{\&omega;}\\ 1& x>\mathrm{\&omega;}\\ 0& x<-\mathrm{\&omega;}\end{array}\right.$

${\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&omega;}}[1+\cos \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right)]& \left|x\right|\&le;\mathrm{\&omega;}\\ 0& \left|x\right|>\mathrm{\&omega;}\end{array}\right.$

b、采用梯度下降法，得到水平集函数的演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}\&dtri;\mathrm{\&phi;})+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)$

c、采用中心差分法得到水平集的迭代方程：

$\frac{{{\mathrm{\&phi;}}_{i,j}}^{n+1}-{{\mathrm{\&phi;}}_{i,j}}^n}{\mathrm{\&Delta;t}}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n\left|\right)}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n)+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)$

其中Δt为迭代步长，一般选择Δt＝1。

以步骤二得到的初始水平集作为初始条件：

φ(x，y，0)＝φ₀(x，y)

φ₀(x，y)为设置的初始轮廓，迭代的稳态解的零水平集即为分割结果。

本发明是在DRLSE模型的基础上，采用非稳态测量算法重新定义停止函数，构建了新的能量模型，提出了基于非稳态测量算法的改进规则距离水平集的图像分割方法。该方法分割过程无须重新初始化水平集，并且由于非稳态测量算法对噪声不敏感，该方法可以准确分割含有高斯噪声和椒盐噪声的图像。

附图说明

图1为具体实施方式一中带噪声的米粒及初始轮廓示意图；

图2为具体实施方式一中采用本发明的米粒分割结果示意图；

图3为具体实施方式一中采用本发明和DRLSE模型的分割结果示意图；

图4为具体实施方式二中采用本发明对原图像的分割结果示意图；

图5为具体实施方式二中采用本发明对原图像混有标准差为0.1的高斯噪声的分割结果示意图；

图6为具体实施方式二中采用本发明对原图像混有标准差为0.2的高斯噪声的分割结果示意图；

图7为具体实施方式二中采用本发明对原图像混有标准差为0.5的高斯噪声的分割结果示意图；

图8为具体实施方式二中采用DRLSE模型对原图像的分割结果示意图；

图9为具体实施方式二中采用DRLSE模型对原图像混有标准差为0.1的高斯噪声的分割结果示意图；

图10为具体实施方式二中采用DRLSE模型对原图像混有标准差为0.2的高斯噪声的分割结果示意图；

图11为具体实施方式二中采用DRLSE模型对原图像混有标准差为0.5的高斯噪声的分割结果示意图。

具体实施方式

具体实施方式一：本实施方式采用基于非稳态测量算法的改进规则距离水平集的图像分割方法，实验用图为混有标准差为0.1的高斯噪声，实现该方法的步骤如下：

一、构造均值算子，用构造的算子设置停止函数；构造方法为：

a、构造两个边长为L的方形矩阵，且矩阵内的元素值均为1/L²，即均值滤波器h₁和h₂，其中L一般取3，5或9。

本实施例中参数选择L＝3。

b、用构造的滤波器h₁，h₂，采用公式 $g=\frac{1}{1+k|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2|}$ 设置停止函数g，其中k为常数，采用公式 $k=\max \left({|{\&dtri;G}_{\mathrm{\&sigma;}}*I|}^2\right)\frac{{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2}{\max \left[\right|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2\left|\right]}$ 获得，I为图像，max(·)表示求最大值，*表示卷积，G_σ是均值为0标准差为σ的高斯滤波器，其中σ的取值范围一般为0.1～0.3。

本实施例中参数选择σ＝0.2。

二、步骤二手动设置初始轮廓R₀，根据该轮廓位置采用公式

初始化水平集函数。

本实施例中，选取参数c₀＝2。

三、步骤三将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果的方法为：

a、将步骤一的停止函数g带入到DRLSE的能量方程，得到基于非稳态测量算法的改进规则距离水平集的能量方程：

${\mathrm{\&epsiv;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dx}+\mathrm{\&lambda;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dx}$

$+\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{H}_{\mathrm{\&omega;}}(-\mathrm{\&phi;})\mathrm{dx}$

式中μ、λ、α为常数，一般μ选择为0.2，λ选择为3～5，α控制轮廓演化的方向，α＜0轮廓膨胀，α＞0轮廓收缩，α值根据边缘的强弱一般选取1～5之间的某个值，弱边缘则选较小的值，强边缘选较大的值，式中势能函数p定义为：

$p\left(s\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;s}\right))& s\&le;1\\ \frac{1}{2}{(s-1)}^2& s\&GreaterEqual;1\end{array}\right.$

Heaviside函数H_ω(x)和Dirac函数δ_ω(x)在实际中分别采用如下近似表达，其中ω的取值范围一般为0.1～2.5：

$H_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}(1+\frac{x}{\mathrm{\&omega;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right))& \left|x\right|\&le;\mathrm{\&omega;}\\ 1& x>\mathrm{\&omega;}\\ 0& x<-\mathrm{\&omega;}\end{array}\right.$

${\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&omega;}}[1+\cos \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right)]& \left|x\right|\&le;\mathrm{\&omega;}\\ 0& \left|x\right|>\mathrm{\&omega;}\end{array}\right.$

b、采用梯度下降法，得到水平集函数的演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}\&dtri;\mathrm{\&phi;})+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)$

c、采用中心差分法得到水平集的迭代方程：

$\frac{{{\mathrm{\&phi;}}_{i,j}}^{n+1}-{{\mathrm{\&phi;}}_{i,j}}^n}{\mathrm{\&Delta;t}}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n\left|\right)}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n)+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)$

其中Δt为迭代步长，一般选择Δt＝1。

以步骤二得到的初始水平集作为初始条件：

φ(x，y，0)＝φ₀(x，y)

φ₀(x，y)为设置的初始轮廓，迭代的稳态解的零水平集即为分割结果。

本实施例中参数选择μ＝0.2、λ＝5、α＝3、ω＝0.2。

本实施例最终得到的含有噪声的图像的分割结果如图2中的曲线轮廓所示，DRLSE模型的分割结果如图3中箭头所指曲线所示，可以看出基于非稳态测量算法的改进规则距离模型的分割结果准确度更高。

具体实施方式二：本实施例采用仿真的双目标图，然后对该图分别加入方差为0.0、0.1、0.2、0.5的高斯噪声，从而形成四组具有不同噪声强度的仿真数据；

实现该方法的步骤如下：

一、构造均值算子，用构造的算子设置停止函数；构造方法为：

a、构造两个边长为L的方形矩阵，且矩阵内的元素值均为1/L²，即均值滤波器h₁和h₂，其中L一般取3，5或9。

本实施例中参数选择L＝3。

b、用构造的滤波器h₁，h₂，采用公式 $g=\frac{1}{1+k|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2|}$ 设置停止函数g，其中k为常数，采用公式 $k=\max \left({|{\&dtri;G}_{\mathrm{\&sigma;}}*I|}^2\right)\frac{{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2}{\max \left[\right|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2\left|\right]}$ 获得，I为图像，max(·)表示求最大值，*表示卷积，G_σ是均值为0标准差为σ的高斯滤波器，其中σ的取值范围一般为0.1～0.3。

本实施例中参数选择σ＝0.2。

二、步骤二手动设置初始轮廓R₀，根据该轮廓位置采用公式

初始化水平集函数。

本实施例中，参数c₀＝2。

三、步骤三将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果的方法为：

a、将步骤一的停止函数g带入到DRLSE的能量方程，得到基于非稳态测量算法的改进规则距离水平集的能量方程：

${\mathrm{\&epsiv;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dx}+\mathrm{\&lambda;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dx}$

$+\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{H}_{\mathrm{\&omega;}}(-\mathrm{\&phi;})\mathrm{dx}$

式中μ、λ、α为常数，一般μ选择为0.2，λ选择为3～5，α控制轮廓演化的方向，α＜0轮廓膨胀，α＞0轮廓收缩，α值根据边缘的强弱一般选取1～5之间的某个值，弱边缘则选较小的值，强边缘选较大的值，式中势能函数p定义为：

$p\left(s\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;s}\right))& s\&le;1\\ \frac{1}{2}{(s-1)}^2& s\&GreaterEqual;1\end{array}\right.$

Heaviside函数H_ω(x)和Dirac函数δ_ω(x)在实际中分别采用如下近似表达，其中ω的取值范围一般为0.1～2.5：

$H_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}(1+\frac{x}{\mathrm{\&omega;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right))& \left|x\right|\&le;\mathrm{\&omega;}\\ 1& x>\mathrm{\&omega;}\\ 0& x<-\mathrm{\&omega;}\end{array}\right.$

${\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&omega;}}[1+\cos \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right)]& \left|x\right|\&le;\mathrm{\&omega;}\\ 0& \left|x\right|>\mathrm{\&omega;}\end{array}\right.$

b、采用梯度下降法，得到水平集函数的演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}\&dtri;\mathrm{\&phi;})+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)$

c、采用中心差分法得到水平集的迭代方程：

$\frac{{{\mathrm{\&phi;}}_{i,j}}^{n+1}-{{\mathrm{\&phi;}}_{i,j}}^n}{\mathrm{\&Delta;t}}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n\left|\right)}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n)+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)$

其中Δt为迭代步长，一般选择Δt＝1。

以步骤二得到的初始水平集作为初始条件：

φ(x，y，0)＝φ₀(x，y)

φ₀(x，y)为设置的初始轮廓，迭代的稳态解的零水平集即为分割结果。

本实施例中参数选择μ＝0.2、λ＝5、α＝3、ω＝0.2。

本实施例最终得到的对含有不同标准差的高斯噪声的图像的分割结果如图4、5、6、7所示，图中的曲线轮廓即为最终分割轮廓。DRLSE模型的分割结果如图8、9、10、11所示。可以看出基于非稳态测量算法的改进规则距离模型的分割结果准确度更高。

本实施例中应用平均距离误差和Jaccard相似性这两项指标来定量化评价分割结果，其中平均距离误差越小、Jaccard相似性越高表示分割结果精度越高，将采用本发明方法计算得到的分割结果与DRLSE模型的分割结果进行比较，结果如表1所示，可以看出本实施例的方法计算得到分割结果精度更高。

表1

Claims (4)

1.一种基于非稳态测量算法的改进规则距离水平集的图像分割方法，其特征在于，实现该方法的步骤如下：

一、构造均值算子，用构造的算子设置停止函数；

二、手动设置初始轮廓，根据该轮廓初始化水平集函数；

三、将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到的稳态解的零水平集即最终的分割结果。

2.根据权利要求1所述的基于非稳态测量算法的改进规则距离水平集的图像分割方法，其特征在于，步骤一所述构造均值算子，用构造的算子设置停止函数的方法为：

a、构造两个边长为L的方形矩阵，且矩阵内的元素值均为1/L²，即均值滤波器h₁和h₂，其中L取3，5或9；

b、用构造的滤波器h₁，h₂，采用公式 $g=\frac{1}{1+k|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2|}$ 设置停止函数g，其中k为常数，采用公式 $k=\max \left({|{\&dtri;G}_{\mathrm{\&sigma;}}*I|}^2\right)\frac{{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2}{\max \left[\right|{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2\left|\right]}$ 获得，I为图像，max(·)表示求最大值，*表示卷积，G_σ是均值为0标准差为σ的高斯滤波器，其中σ的取值范围为0.1～0.3。

3.根据权利要求2所述的基于非稳态测量算法的改进规则距离水平集的图像分割方法，其特征在于，步骤二手动设置初始轮廓R₀，根据该轮廓位置采用公式

初始化水平集函数，c₀为常数，取c₀＝5。

4.根据权利要求1、2或3所述的基于自适应扩散基函数分解的组织纤维束结构信息提取方法，其特征在于，步骤三将步骤一设置的停止函数带入到DRLSE模型的能量方程，采用中心差分法对能量方程最小化，以步骤二得到的初始化水平集作为初始条件进行迭代，得到稳态解的具体过程为：

a、将步骤一的停止函数g带入到DRLSE的能量方程，得到基于非稳态测量算法的改进规则距离水平集的能量方程：

${\mathrm{\&epsiv;}}_{\mathrm{\&epsiv;}}\left(\mathrm{\&phi;}\right)=\mathrm{\&mu;}\underset{\mathrm{\&Omega;}}{\&Integral;}p\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)\mathrm{dx}+\mathrm{\&lambda;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)|\&dtri;\mathrm{\&phi;}|\mathrm{dx}$

$+\mathrm{\&alpha;}\underset{\mathrm{\&Omega;}}{\&Integral;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{H}_{\mathrm{\&omega;}}(-\mathrm{\&phi;})\mathrm{dx}$

式中μ、λ、α为常数，一般μ选择为0.2，λ选择为3～5，α控制轮廓演化的方向，α＜0轮廓膨胀，α＞0轮廓收缩，α值根据边缘的强弱一般选取1～5之间的某个值，弱边缘则选较小的值，强边缘选较大的值，式中势能函数p定义为：

$p\left(s\right)=\left\{\begin{array}{cc}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}(1-\cos \left(2\mathrm{\&pi;s}\right))& s\&le;1\\ \frac{1}{2}{(s-1)}^2& s\&GreaterEqual;1\end{array}\right.$

Heaviside函数H_ω(x)和Dirac函数δ_ω(x)分别采用如下近似表达，其中ω的取值范围为0.1～2.5：

$H_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2}(1+\frac{x}{\mathrm{\&omega;}}+\frac{1}{\mathrm{\&pi;}}\sin \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right))& \left|x\right|\&le;\mathrm{\&omega;}\\ 1& x>\mathrm{\&omega;}\\ 0& x<-\mathrm{\&omega;}\end{array}\right.$

${\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(x\right)=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\&omega;}}[1+\cos \left(\frac{\mathrm{\&pi;x}}{\mathrm{\&omega;}}\right)]& \left|x\right|\&le;\mathrm{\&omega;}\\ 0& \left|x\right|>\mathrm{\&omega;}\end{array}\right.$

b、采用梯度下降法，得到水平集函数的演化方程：

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;\mathrm{\&phi;}\left|\right)}{|\&dtri;\mathrm{\&phi;}|}\&dtri;\mathrm{\&phi;})+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;\mathrm{\&phi;}}{|\&dtri;\mathrm{\&phi;}|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left(\mathrm{\&phi;}\right)$

c、采用中心差分法得到水平集的迭代方程：

$\frac{{{\mathrm{\&phi;}}_{i,j}}^{n+1}-{{\mathrm{\&phi;}}_{i,j}}^n}{\mathrm{\&Delta;t}}=\mathrm{\&mu;}\mathrm{div}(\frac{p^{\&prime;}\left(\right|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n\left|\right)}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n)+{\mathrm{\&lambda;\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)\mathrm{div}\left(\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}\frac{\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n}{|\&dtri;{{\mathrm{\&phi;}}_{i,j}}^n|}\right)$

$+\mathrm{\&alpha;}\frac{1}{1+k[{(I*h_1)}^2*h_2-{(I*h_1*h_2)}^2]}{\mathrm{\&delta;}}_{\mathrm{\&omega;}}\left({{\mathrm{\&phi;}}_{i,j}}^n\right)$

其中Δt为迭代步长，选择Δt＝1；

以步骤二得到的初始水平集作为初始条件：

φ(x，y，0)＝φ₀(x，y)

φ₀(x，y)为设置的初始轮廓，迭代的稳态解的零水平集即为分割结果。

Images