图像分割^[1]是将图像划分成不同类型区域，在同一区域呈现特征相似性，而在不同区域呈现特征差异性。它是图像模式识别和分析处理中的基础性工作。近年来，具有良好理论基础和较好分割结果的活动轮廓模型方法得到研究学者的广泛关注。

活动轮廓模型方法的基本原理是用3维函数 ${\textit z}{\rm{ = }}\phi (x(t),y(t))$ 的零水平集 $\phi (x(t),y(t)) = 0$ 表示2维演化曲线，它是一种基于曲线演化的水平集方法，其利用变分理论极小化基于图像信息水平集函数 $\phi (x)$ 的能量泛函，从而得到水平集函数的演化方程。利用演化方程更新水平集函数，驱动零水平集 $\phi (x) = 0$ 向目标边界运动并最终收敛停留在目标边界上得到最终的图像分割结果。几何活动轮廓模型在能量泛函极小化框架下可以引入图像信息，使得此种图像分割方法有很强的适应性。

活动轮廓模型一般分为基于边缘和基于区域的两种类型，分别利用边缘信息和区域的统计特性驱动活动轮廓向目标边界靠近。作者研究的模型属于后一种类型。Chan等^[2]利用图像的全局区域信息，使用全局二值分段常数驱动活动轮廓向目标边缘逼近CV模型，但不能有效分割灰度分布不均匀图像。为此研究人员提出基于混合区域活动轮廓模型^[3–4]。Li等^[5]利用图像像素的局部邻域区域信息拟合局部能量，提出了局部二值拟合LBF模型；Zhang等^[6]提出LIF模型，利用局部图像拟合能量提取局部图像信息以达到分割图像目的。但这两种模型对椒盐噪声污染图像分割效果不理想。为此研究人员提出基于局部熵LCK模型^[7]、GLCK模型^[8]。刘瑞娟等^[9]提出了一种改进CV模型和LIF模型线性组合的活动轮廓模型。戚世乐等^[10]提出先使用简化CV模型得到粗分割结果，再利用LBF模型进行细分割的两阶段活动轮廓模型。

LBF模型中的拟合函数对活动轮廓内外部区域的灰度局部加权平均，所以LBF模型对高斯噪声有一定的抗噪性。但是，其对椒盐噪声不能取得较好的分割结果，作者在分析LBF模型拟合项基础上提出一种新的局部区域的活动轮廓模型，并且提出了一种新的边缘检测算子，提高模型对噪声图像分割的抗噪性。

1 LBF模型

给定图像I，在图像区域Ω中C表示一闭合曲线，将图像区域Ω分成两部分Ω₁=outside（C），Ω₂=inside（C）。则曲线C可以看成目标图像的轮廓，LBF模型能量函数定义为：

$\begin{aligned}[b]\!\!\!\! \,{E_{\rm LBF}} \! \! = \, & E(\phi ,{f_1}(x),{f_2}(x)) = \\[-2pt]& \! {\lambda _1} \! \int \! {( \! \int \! {{K_\sigma }(x \! - \! y)|I(y) \! - \! {f_1}(x){|^2}} H(\phi (y)){\rm{d}}y){\rm{d}}x + } \\[-2pt]& \! {\lambda _2} \! \int \! {( \! \int \! {{K_\sigma }(x \! - \! y)|I(y) \! - \! {f_2}(x){|^2}} (1 \! - \! H(\phi (y))){\rm{d}}y){\rm{d}} \! x} \!\!\!\!\!\!\!\!\!\!\!\!\end{aligned}$

(1)

利用变分方法和梯度下降法，式（1）的梯度下降流为:

$\frac{{\partial \phi }}{{\partial t}} = - \delta (\phi )({\lambda _1}{e_1} - {\lambda _2}{e_2})$

(2)

式中， ${\lambda _1}$ 和 ${\lambda _2}$ 是2个正的常数， ${e_i}(x) = \int {{K_\sigma }(y - x)|I(x)} - $ ${f_i}(y){|^2}{\rm{d}}y (i = 1,2)$ 。式（1）中f₁（x）和f₂（x）是x邻域的加权均值拟合函数，定义为：

$\begin{aligned}& {f_1}(x) = \frac{{{K_\sigma }(x) \cdot (H(\phi (x))I(x))}}{{{K_\sigma }(x) \cdot H(\phi (x))}},\\[3pt]& {f_2}(x) = \frac{{{K_\sigma }(x) \cdot ((1 - H(\phi (x)))I(x))}}{{{K_\sigma }(x) \cdot (1 - H(\phi (x)))}}{\text{。}}\end{aligned}$

2 本文模型 2.1 提出的拟合项

在LBF模型中，从拟合函数f₁（x）和f₂（x）公式可知，f₁（x）和f₂（x）分别表示为以x为中心的活动轮廓内外区域的局部加权均值，这相当于对图像进行了局部均值滤波，所以LBF模型对高斯噪声有着一定的抗噪性，但其对椒盐噪声不能得到较好的分割结果。中值滤波^[11]对椒盐噪声具有较好的抑制作用。

根据文献[12]研究可知：设数据集 $\{ {x_i}\} _{i = 1}^n$ ，求解 $\min \displaystyle\sum\limits_{i = 1}^n {{{({x_i} - c)}^2}} $ ，c为极小值的必要条件为 $\displaystyle\sum\limits_{i = 1}^n {({x_i} - c) \! = \! 0} $ ，则解得 $c =\displaystyle \frac{1}{n}\sum\limits_{i = 1}^n {{x_i}} $ ，显然c为数据集 $\{ {x_i}\} _{i = 1}^n$ 的均值。求解 $\min \displaystyle\sum\limits_{i = 1}^n {|{x_i} - m|} $ ，m为极小值的必要条件为 $\displaystyle\sum\limits_{i = 1}^n {\frac{{{x_i} - m}}{{|{x_i} - m|}}} \! = \! 0$ 。假定数据集 $\{ {x_i}\} _{i = 1}^n$ 从小到大排序，并且假定1） ${x_i} = m$ 时， $\displaystyle\frac{{{x_i} - m}}{{|{x_i} - m|}} = 0$ ，又因为 ${x_i} > m$ ，则 $\displaystyle\frac{{{x_i} - m}}{{|{x_i} - m|}} \! = \! 1$ 。2） ${x_i} < m$ 时，则 $\displaystyle\frac{{{x_i} - m}}{{|{x_i} - m|}} = - 1$ ，又有：如果n为偶数，可以解得m = ${x_{(n + 1)/2}}$ ；如果n为奇数，可以解得 $m = ({x_{n/2}} + {x_{n/2 + 1}})/2$ ，即m为数据集 $\{ {x_i}\} _{i = 1}^n$ 的中值。

根据以上讨论，本文定义：

$\begin{aligned}{m_i}\left( x \right) = med\left( {\frac{{\displaystyle\int {{K_\sigma }(x - y)\left( {{M_i}(x)I(x)} \right){\rm{d}}y} }}{{\displaystyle\int {{K_\sigma }(x - y){M_i}(x){\rm{d}}y} }}} \right), \ i = 1,2 \!\!\!\!\!\!\end{aligned}$

(3)

式中， ${M_1}(x) = H(\phi (x))$ ， ${M_2}(x) = 1 - H(\phi (x))$ ，med是中值算子。从而定义新的拟合能量项：

$\begin{aligned}[b]& {E_n} \left( {\phi , {m_1}(x),{m_2}(x)} \right) = {\int_{outside\left( C \right)}\left| {I - {m_1}\left( x \right)} \right.\left| {{\rm d}x{\rm d}y} \right.} + \\& \qquad\quad\quad\quad\quad {\int _{inside(C)}\left| {I - {m_2}\left( x \right)} \right.\left| {{\rm d}x{\rm d}y} \right.} \end{aligned}$

(4)

采用水平集方法，式（4）可以改写成：

$\begin{aligned}[b]& {E_n}(\phi , {m_1}(x),{m_2}(x)) = \int {|I - {m_1}(x)|H(\phi ){\rm d}x{\rm d}y} + \\& \qquad \int {|I - {m_2}(x)|(1 - H(\phi )){\rm d}x{\rm d}y} \end{aligned}$

(5)

2.2 测地弧长正则项

在水平集方程的演化过程中为了保持活动轮廓的平滑性，通常将轮廓C的弧长作为正则项加入到能量函数中，活动轮廓弧长能量函数定义为：

${E_{l'}} = \int {|\nabla H(\phi )| {\rm{d}}x} = \int {\delta (\phi )|\nabla \phi |} {\rm{d}}x$

(6)

式（6）中，关于水平集函数 $\phi $ 的梯度下降流为：

$\frac{{\partial \phi }}{{\partial t}} = {\rm{div}}(\frac{{\nabla \phi }}{{|\nabla \phi |}})$

(7)

轮廓弧长能量函数梯度下降流方程为平均曲率运动方程，会使得水平集函数过度平滑，根据文献[13]研究可知，此正则项可以使用活动轮廓的测地弧长来替代。测地弧长能量函数定义为：

${E_l} = \int {g\delta (\phi )|\nabla \phi |} {\rm d}x$

(8)

式中，g为边缘停止函数，定义 $g = 1/{(1 + |\nabla {G_\sigma }*I|)^2}$ 其中，| $ * $ |表示绝对值运算。显然边缘停止函数是基于图像梯度信息的，难以得到图像不同区域间模糊的弱边缘和灰度渐进图像目标，并且易受到噪声的影响，即其噪声能力弱。根据文献[14]研究可知2阶导数能够有效地区别边缘和渐变区域。文献[14]中，义 $ {u_{\eta \eta }} =$ $ (u_x^2{u_{xx}} \! + \! 2{u_x}{u_y}{u_{xy}} \! + \! u_y^2{u_{yy}})/(u_x^2 + u_y^2)$ ， ${u_{\xi \xi }} = (u_y^2{u_{xx}} - 2{u_x}{u_y}{u_{xy}} + $ $ u_x^2{u_{yy}})/(u_x^2 + u_y^2)$ ，并且提出了一种差分曲率的边缘检测算子，表示为 $D = ||{u_{\eta \eta }}| - |{u_{\xi \xi }}||$ 。经过试验，作者提出一种新的边缘检测算子，应用到边缘停止函数g中。参考 ${u_{\eta \eta }}$ 和 ${u_{\xi \xi }}$ 定义：

$ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! u_{\eta \eta }' = u_x^2{u_{xx}} + 2{u_x}{u_y}{u_{xy}} + u_y^2{u_{yy}} $

(9)

$u_{\xi \xi }' = u_y^2{u_{xx}} - 2{u_x}{u_y}{u_{xy}} + u_x^2{u_{yy}}$

(10)

省略原定义形式的归一化操作，定义：

$M_1= |u_{\eta \eta }'| + |u_{\xi \xi }'|$

(11)

从而提出了一种新的边缘检测算子：

$N = {M_1}/2$

(12)

通过一个实例说明本文定义的边缘检测算子的有效性。图1（a）是具有渐变区域无噪声图像，图1（b）、（c）是分别添加方差为0.003、0.010高斯噪声的噪声图像，图1（d）、（e）、（f）是分别添加密度为0.001、0.005和0.010的椒盐噪声。图2（a）～（c）、（d）～（f）、（j）～（l）、（m）～（o）、（p）～（r）分别对应使用梯度算子、差分算子、本文定义算子对图1中每幅子图的边缘检测的比较，图1（b）、（c）噪声方差为0.003和0.010的高斯噪声的噪声图像，图1（d）、（e）、（f）是密度为0.001、0.005和0.010的椒盐噪声的边缘检测的比较。图2中数值都被归一化成[0， $40 \cdot {{e/{\rm max}}}\left( {{e}} \right)$ ]，其中， $e \in \{ \nabla ,D,N\} $ ， $\nabla $ 为梯度算子，D为差分算子，N为本文定义边缘检测算子。从图2（a）、（b）、（c）可以看出：在无噪声情况下，3种算子都能较为有效地检测出渐变图像的边缘，相较于梯度算子和差分算子，本文定义算子能很好地检测出边缘所在，不受图像其他区域信号影响。在图像受到较为轻微高斯噪声影响下，如图2（d）、（e）、（f）所示，本文定义算子几乎不受噪声影响，与无噪声情况下的边缘检测效果一致，而梯度算子和差分算子只有在人工目视情况下可见，受到噪声较为严重的影响。在图像受到较为严重噪声影响下各种边缘检测算子比较如图2（g）、（h）、（i）所示，本文定义算子依然可以检测出图像边缘，而此时的梯度算子和差分算子在人工目视情况下都已不可见。在受到轻微椒盐噪声影响下，3种算子都能检测出边缘所在，从目视情况来看，本文提出算子受到噪声影响最小，如图2（j）、（k）、（l）所示。随着椒盐噪声增强，梯度算子和差分算子受到噪声影响越来越严重，而本文提出算子在边缘检测几乎不受噪声强弱影响，如图2（m）、（n）、（o）和（p）、（q）、（r）所示。通过上述实验可知，本文提出的新的边缘检测算子能有效地检测出图像的边缘，特别是在图像受到噪声的影响下，相较于传统的边缘检测算子具有更好的鲁棒性。本文将新提出的边缘检测算子应用到水平集模型中的停止函数，重新定义边缘停止函数，表示为：

$g = \frac{1}{{1 + |N({G_\sigma } \cdot I){|^2}}}$

(13)

2.3 惩罚能量

为了保持活动轮廓在演化过程中的稳定性，初始水平集函数通常将曲线生成符号距离函数，并且在曲线演化迭代过程中必须周期地新初始化，使得水平集函数重新变成带符号的距离函数，重新初始化的计算量大，分割速度慢。文献[15]提出了一种符号距离惩罚能量函数，定义为：

${E_{\rm p}} = \frac{1}{2}\int {{{(|\nabla \phi - 1|)}^2}} {\rm{d}}x$

(14)

此能量函数只与水平集函数有关，在极小化E_p过程中， $|\nabla \phi |$ 逐渐逼近于1，约束水平集函数在演化过程中保持为1个符号距离函数，从而较好地解决了重新初始化的难题。

2.4 本文模型

结合LBF模型中局部拟合能量项E_LBF、测地弧长能量项E_l、惩罚能量项E_p和本文新提出的拟合能量项项E_n得到能量函数：

$E = \alpha {E_{\rm{LBF}}} + \beta {E_{\rm n}} + \kappa {E_{\rm l}} + \mu {E_{\rm p}}$

(15)

式中， $\alpha $ 、 $\beta $ 、 $\kappa $ 、 $\mu $ 为正的实数， $\alpha $ 和 $\beta $ 分别是调整局部拟合能量项E_LBF和本文新提出的拟合能量项E_n的权重。当图像中不含椒盐噪声时， $\alpha $ 取较大值，否则， $\alpha $ 取较小值。为了降低参数选取的难度，本文实验中固定 $\beta $ =0.01。局部拟合能量项E_LBF中 ${\lambda _1}$ 和 ${\lambda _2}$ 的值取1，则本文最终的能量函数为：

$\begin{aligned}[b]\! \! \! E =\, & \alpha \! \int \! {(\! \int \! {{K_\sigma }(x - y)|I(y) - {f_1}(x){|^2}} H(\phi (y)){\rm d}y){\rm{d}}x + } \\& \alpha \! \int \! {(\! \int \! {{K_\sigma }(x - y)|I(y) - {f_2}(x){|^2}} (1 - H(\phi (y))){\rm d}y){\rm{d}}x} + \! \! \! \! \! \\& \beta \! \int \! {|I - {m_1}(x)|H(\phi ){\rm{d}}x{\rm{d}}y} + \\& \beta \! \int \! {|I - {m_2}(x)|(1 - H(\phi )){\rm{d}}x{\rm{d}}y} + \\& \kappa \! \int \! {g\delta (\phi )|\nabla \phi |} {\rm{d}}x{\rm{d}}y + \mu \frac{1}{2}\int {{{(|\nabla \phi - 1|)}^2}} {\rm{d}}x{\rm{d}}y{\text{。}}\end{aligned}$

利用变分法和梯度下降法，关于水平集 $\phi $ 的梯度下降流为:

$\begin{aligned}[b]\frac{{\partial \phi }}{{\partial t}} = & - \alpha \delta (\phi )({e_1} - {e_2}) + \\[-2pt]& - \beta \delta (\phi )(|I - {m_1}| - |I - {m_2}|) + \\& \kappa \delta (\phi ){\rm{div}}\left(g\frac{{\nabla \phi }}{{|\nabla \phi |}}\right) + \mu \left({\nabla ^2}\phi - {\rm{div}}\left(\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}} \right)\right)\end{aligned}$

(16)

式中，

${e_i}(x) = \int {{K_\sigma }(y - x)|I(x) - {f_i}(y){|^2}{\rm{d}}y, \ \ i = 1,2}{\text{。}} $

e₁和e₂中，拟合函数f（x）为

${f_i}(x) = \frac{{{K_\sigma }(x)*({M_i}(x)I(x))}}{{{K_\sigma }(x)*{M_i}(x)}}, \ \ i = 1,2{\text{。}}$

式（16）中，

${m_i}(x) = med(\frac{{\int {{K_\sigma }(x - y)({M_i}(x)I(x)){\rm d}y} }}{{\int {{K_\sigma }(x - y){M_i}(x){\rm d}y} }}), \ \ i = 1,2{\text{。}}$

其中： ${M_1}(x) = H(\phi (x))$ ； ${M_2}(x) = 1 - H(\phi (x))$ ；med为中值算子，窗口尺寸大小以高斯核参数为标准，取2 $\sigma $ +1；K为核函数，本文选取高斯核函数，定义为：

${K_\sigma }(x) = \frac{1}{{{{(2\text{π} )}^{n/2}}{\sigma ^n}}}{{\rm{e}}^{ - |x{|^2}/2{\sigma ^2}}}$

(17)

式中， $ \sigma$ 为尺度参数，本文实验中取 $\sigma $ =3。

上述公式中采用正则化Heaviside和Dirac函数，分别定义为：

$H(x) = \frac{1}{2}\left( {1 + \frac{2}{\text{π}}{\rm{arctan}}\left( {\frac{x}{\varepsilon }} \right)} \right)$

(18)

$\delta (x) = H'(x) = \frac{1}{\text{π} }\frac{\varepsilon }{{{\varepsilon ^2} + {x^2}}}$

(19)

本文实验中选取 $\varepsilon $ =1。

3 实验结果及分析

在不同噪声污染的图像上验证本文提出模型的有效性，并且和CV模型^[7]、LBF模型^[10]、LIF模型^[11]进行对比试验，证明本文提出方法在椒盐图像分割上抗噪性优势。实验中，时间步长取Δt=0.1，空间步长Δx=Δy=1， $\beta $ =0.01，图像中有高斯噪声时取 $\alpha $ =5，含有椒盐噪声时取 $\alpha $ =1.0，参数 $\kappa $ 取65，参数 $\mu $ 取2.0，CV、LBF和LIF模型参数具体见相关文献[7, 10–11]。

图3（a）～（e）分别对应无噪声原始图像，添加方差0.001、0.002高斯噪声图像和添加密度为0.005、0.010的椒盐噪声。从图4（a）～（d）可以看出：LBF模型和本文模型都能对无噪声图像进行有效的分割；CV模型中拟合项考虑的是活动轮廓内外灰度平均值，不能很好地描述灰度不均匀图像特性，所以出现分割错误；LIF模型中拟合项使用局部拟合图像和原始图像差异程度描述图像特性，对灰度不均匀图像描述特性不够，所以出现大量分割错误。在受到高斯噪声污染的图像中，LBF模型出现分割错误，并且随着噪声水平的提高，错误率也随之提高，如图4（g）、（k）所示。而本文模型在高斯噪声影响下，依然能得到较为满意的分割结果，随着噪声水平的增高，本文模型方法分割结果几乎不受噪声影响，如图4（h）、（l）所示。这是因为本文模型结合了LBF模型局部拟合项；而LBF模型拟合项在高斯噪声图像分割中具有一定的抗噪性，但其对受噪声影响的弱边缘分割不理想，出现错误分割，如图4（g）、（k）所示；并且本文模型中提出一种新的边缘检测算子，增强了本文模型的抗噪性，如图4（h）、（l）所示。CV模型和LIF模型在受高斯噪声影响下，出现错误分割结果，如图4（e）、（f）、（i）、（j）所示。本文提出模型提高对椒盐噪声图像分割的鲁棒性，但是依然能对高斯噪声图像分割具有优良鲁棒性，本文后续试验将不再对受高斯噪声影响的图像进行测试。

图4（m）～（p）、（q）～（t）是各种模型对椒盐噪声图像的分割结果。从分割结果来看：CV模型和LIF模型在不同噪声水平上都不能得到正确的分割结果；LBF也出现了错误分割结果，不能收敛于弱边缘，并且大量收敛于图像中孤立的椒盐噪声点，如图4（o）、（s）所示，并且随着噪声水平的增强其分割错误率增大。本文模型在不同噪声水平上都能得到满意的分割结果，并且分割结果几乎不受噪声水平的影响，如图4（p）、（t）所示。其原因为：一方面，本文提出的新的边缘检测算子具有很强的抗噪性，第2.2节通过实验已经详细说明这一问题；另一方面，针对椒盐噪声，本文新提出的拟合项相当于是对噪声图像进行局部中值滤波计算，从而消弱椒盐噪声对分割的影响。

图5、6、7是3个实例的不同噪声图图像和不同方法分割的比较。

图5、6、7中（a）～（c）是不同图像添加不同密度水平的椒盐噪声，图5、6、7中（d）～（k）分别对应CV模型、LIF模型、LBF模型和本文模型对不同噪声水平图像的分割结果。从结果图像可以看出，CV模型、LIF模型、LBF模型在不同噪声水平上都出现了错误分割而不能得到较为满意的分割结果，本文提出的模型在不同噪声水平上都能取得较为满意的分割结果。通过不同图像的不同噪声水平的分割结果可以看出，本文模型相较于其他模型更具有优势。

4 结　论

为了提高活动轮廓模型对椒盐噪声的抗噪性，在分析局部二值拟合模型中局部拟合项基础上，提出了一种新的局部拟合项，新的能量函数中拟合项的极小值是局部区域的中值，中值的计算对椒盐噪声不敏感；提出一种新的边缘检测算子重新定义边缘停止函数，采用基于新的边缘停止函数的测地弧长作为正则项，进一步提高模型的抗噪性；在新的活动轮廓模型中引入惩罚能量，避免了水平集函数在演化过程中的重新初始化问题。实验结果表明，本文模型对噪声图像分割相较于其他模型具有更好的鲁棒性。