Distance Regularized Level Set Image Segmentation Algorithm by Means of Dislocation Theory
-
摘要: 把材料科学中的位错理论引入到水平集方法中.图像中水平集曲线的演化被看作刃位错中位错线的滑移过程，运用位错动力学机制推导出驱使水平集曲线演化的位错组态力.结合距离正则化水平集方法，把水平集方法的边缘检测函数替换为基于位错动力学理论的速度停止函数，并构建了新的距离正则化水平集函数演化方程.水平集曲线在位错组态力和速度停止函数的驱使下移动.位错组态力反映了单位长度曲线上的平均受力情况，不仅包括了图像梯度信息，也包括了位错组态力的作用范围等信息，因此可以有效地避免在局部图像梯度异常的情况下发生曲线停止演进的现象，或者避免在弱边缘处由于图像梯度较小发生局部边界泄漏的现象.实验结果表明，本文算法对弱边缘图像具有较好的分割效果.Abstract: Dislocation theory of materials science is introduced into the level set method. Curve evolution of the level set method is viewed as a slipping process of edge dislocation, and the evolution of zero level set is driven by the dislocation configuration force which is derived from the dislocation dynamics mechanism. Combined with distance regularized level set method, the edge indicator function is replaced by the speed stopping function, and a new evolution equation of distance regularized level set method is constructed with the dislocation theory. In the proposed algorithm, the dislocation configuration force and the speed stopping function drive the curve evolution of level set. The dislocation configuration force reflects the average force on the unit length curve, not only including the gradient information of image but also including the information of the effective range of the dislocation configuration force. The proposed algorithm can effectively avoid the phenomenon that the level set function stops evolution because of abnormal image gradient, and the phenomenon of boundary leakage because of the smaller image gradient. Experimental results show that the proposed algorithm has good segmentation performance for images with weak boundaries.
-
Key words:
- Image segmentation /
- level set method /
- dislocation /
- distance regularized evolution
-
水平集方法是在活动轮廓模型的基础上发展而来的一种分割算法.水平集方法将低维曲线嵌入到高一维的空间中并将嵌入位置当作零水平集位置, 使用连续曲线来表达目标边缘, 由闭合曲线演化理论得到曲线演化方程.定义一个与曲线参数有关的能量泛函, 因此分割过程就转变为求解能量泛函最小值的过程.在图像分割应用中, 最终零水平集截面上点的集合就是初始轮廓曲线的演化结果, 该点集为图像分割结果.
Caselles算法是最早的水平集模型, 水平集函数基于图像梯度和曲线曲率进行演化[1].该模型在本质上和Snake模型相似, 是基于能量泛函的优化问题.所以运用该模型解决图像分割问题时不可避免地存在对初始轮廓过分依赖的缺点.在对能量泛函求解的过程中, 由于局部极值问题, 导致初始轮廓的选取直接影响到是否能得到正确的分割结果. CV (Chan-Vese)模型是基于区域分割的水平集模型[2].该模型利用能量泛函最小值进行曲线演化得到分割结果. CV模型的缺点是需要不断初始化, 这个过程大大增加了计算量, 降低了演化效率. Bernard算法利用二维三阶B样条基函数把水平集表示成离散的形式[3], 能够有效抑制噪声. Shi算法采用双链表的方式来实现曲线的演化, 无需重新初始化、无需求解偏微分方程, 是一种快速水平集算法[4]. Li等提出了一种新的基于区域的水平集图像分割方法, 能够在图像强度分布不均匀的情况实现分割[5]. Min等针对复杂自然图像的分割, 提出一种新颖的基于灰度和纹理信息的模型, 主要融合了灰度信息和纹理信息两方面的特征.相比传统的水平集图像分割方法, 能够更好地捕获图像灰度信息[6]. Wang等提出了一种新的复杂图像分割水平集方法, 其中局部统计分析和全局相似性度量都纳入到能量泛函的构建中[7]. Zhao等应用二次区域生长法进行初步分割, 构造新的边缘指示函数和水平集能量函数, 有效地解决了肝脏分割问题[8].周则明等采用水平集方法进行图像预处理, 提出了一种基于能量最小化的星载SAR (Synthetic aperture radar)图像建筑物分割方法[9].张迎春等提出了一种基于粗糙集和新能量公式的水平集图像分割方法[10].
在图像分割中, 如果目标边界比较模糊且灰度分布不均匀, 使用传统的水平集方法进行图像分割难以准确地分割出目标轮廓.本文在水平集方法中引入位错理论, 把水平集函数的演化视为位错线的运动, 提出了一种基于位错理论的水平集方法.
1. 距离正则化水平集方法
假设由$y=f(x)$表示平面曲线, 也可写为隐函数形式, $y-f(x)=0$, 表示了$x$与$y$的对应关系.
$ \begin{equation} \phi (x, y) = y - f(x)=0 \end{equation} $
(1) 给定一个初始封闭曲线$C$, 该曲线可以不断沿其法线方向朝外或朝内以一定的速度进行演化, 引入时间变量$t$后, 则形成了随时间变化的曲线族$ C(t)$.可以把曲线族看作是更高维空间曲面函数$\phi$的零水平集$\{\phi= 0, t\}$, $\phi( x, y, t)$就是随时间变化的水平集函数, 而水平集就是曲面上具有相同函数值的点的集合, 零水平集就是曲面上函数值为0的所有点的集合.
给定水平集函数$\phi( C, t)$, 则$\phi(C (t), t) = 0$为$t $时刻演化曲线$C (t )$的零水平集, 根据复合函数求导对$\phi(C (t), t) = 0$进行求导可得(求关于时间的偏导数):
$ \begin{equation} \frac{{\partial \phi }}{{\partial t}} + \nabla \phi \cdot \frac{{\partial C}}{{\partial t}} = 0 \end{equation} $
(2) 其中, $\nabla$为汉密尔顿算子, $\nabla \phi$为$\phi$的梯度.
由平面内曲线性质可知, 初始轮廓曲线的单位法向矢量为
$ \begin{equation} \pmb N = - \frac{{\nabla \phi }}{{|\nabla \phi |}} \end{equation} $
(3) 曲线演化方程可简写为
$ \begin{equation} \frac{{\partial C}}{{\partial t}} = F\pmb{N} \end{equation} $
(4) 其中, $F$为法向速度函数.
将曲线演化理论中的式(3)与式(4)代入式(2)中, 可得到水平集的基本演化方程:
$ \begin{equation} \frac{{\partial \phi }}{{\partial t}} = F|\nabla \phi | \end{equation} $
(5) 水平集方法的实质就是求解一个随时间变化的偏微分方程, 该演化方程属于Hamilton-Jacobi方程, 可以通过分离变量法求解, 即按照对时间项$\triangle t$和空间项(将图像按像素进行网格化)分开处理的方法对此类方程进行数值求解.
为了保证数值计算的精度, 水平集函数应具备如下两个性质:
1) 水平集函数必须具有一定的光滑性.从水平集演化方程中可以得知, 演化时需要计算水平集函数的梯度$|\nabla \phi |$, 而梯度对数值变化是非常敏感的.且大部分方法中的速度函数$F$都会包含有常数项和曲率项, 曲率项在数值求解时需要计算水平集函数的二阶导数, 而常数项也要计算一阶导数.因此水平集函数的计算精度由水平集函数的光滑性直接决定.
2) 水平集函数应保持为符号距离函数(Signed distance function, SDF).当函数$\phi$满足$|\nabla \phi |=1$时, $\phi$即为符号距离函数.要求水平集函数保持为符号距离函数的原因有两方面, 一方面因为曲线的隐式表达与符号距离函数是等价的, 另一方面从数值计算的角度来说, 由于符号距离函数的梯度$|\nabla \phi |=1$, 因此可以保证离散网格的大小为1, 使得数值计算具有较高的精度.
所以, 运用水平集方法进行图像分割时需要把零水平集函数初始化为符号距离函数[11].当水平集函数在二维空间中演化时, 设$C(t=0)$为初始轮廓曲线, 则零时刻的水平集函数为
$ \begin{align} &f(x, y, t = 0) = \nonumber\\&\quad {\rm{sgn}}(x, y, C(t = 0)) \times {\rm{dist}}(x, y, C(t = 0)) \end{align} $
(6) 其中, ${\rm{sgn}}(\cdot)$是符号函数, ${\rm{sgn}}(x, y, C(t=0))$的取值为$+1$或$-1$.当点$(x, y)$在$C(t=0)$的外部时, ${\rm{sgn}}(x, y, C(t=0))=+1$, 当点$(x, y)$在$C(t=0)$时, ${\rm{sgn}}(x, y, C(t=0))=-1$. ${\rm{dist}}(x, y, C(t=0))$为点$(x, y)$到初始曲线$C(t=0)$的最短距离.所以, 计算SDF函数分为两步: 1)判断点在初始轮廓曲线的内部还是外部; 2)计算点到曲线的最短距离.
在基于水平集方法处理图像分割问题时, 初始轮廓会随演化而发生震荡, 导致其逐渐失去光滑性及符号距离函数的特征, 演化过程中的误差积累将直接使演化结果严重偏离目标轮廓.为了解决这一问题, 需要定时为零水平集重新初始化.标准的重新初始化方法是通过解以下的Hamilton-Jacobi方程实现:
$ \begin{equation} \frac{{\partial \phi }}{{\partial t}} = {\rm{sgn}}({\phi})(1 - \left| {\nabla \phi } \right|) \end{equation} $
(7) 零水平集重新初始化虽然保证了演化的准确性, 但计算量非常大, 牺牲了效率.为了解决这一问题, Li等提出了距离正则化水平集模型(Distance regularized level set, DRLSE)[12-13].在该模型中引入了能量惩罚项使水平集函数自我约束为符号距离函数, 避免了需要不断初始化零水平集函数的问题, 并构造一个距离正则化项, 与外部能量相结合驱动初始轮廓到达目标轮廓位置.该模型能量泛函可定义为
$ \begin{equation} E\left( \phi \right) = \mu {R_p}\left( \phi \right) + {E_{\rm ext}}(\phi ) \end{equation} $
(8) 其中, $R_{p}(\phi)$为水平集函数$\phi$的正则化项, 为常数, $E_{\rm ext}(\phi)$为外部能量泛函, 取决于具体应用.在图像处理应用中, $E_{\rm ext}(\phi)$为图像区域或曲线的能量.
正则化项$R_{p}(\phi)$定义为
$ \begin{equation} {R_p}\left( \phi \right) := \int_\Omega {p\left( {|\nabla \phi |} \right)} \text{d}x \end{equation} $
(9) 其中, $p$为正则化项$R_{p}(\phi)$的势能(能量密度)函数, 定义为
$ \begin{equation} p\left( s \right) = \left\{ \begin{array}{ll} \dfrac{1}{2\pi}\left( 1 - \cos \left( 2\pi s \right) \right), &s < 1\\ \dfrac{1}{2}\left(s - 1\right)^2, &s \ge 1 \end{array} \right. \end{equation} $
(10) 外部能量泛函$E_{\rm ext}(\phi)$定义为
$ \begin{align} &{E_{\rm ext}}\left( \phi \right) =\lambda \int_\Omega {g{\delta _\varepsilon }} \left( \phi \right)|\nabla \phi |\text{d}x + \nonumber\\ &\qquad \alpha \int_\Omega {g{H_\varepsilon }\left( { - \phi } \right)} \text{d}x \end{align} $
(11) 其中, $\lambda$和$\alpha$为常数, $\delta _\varepsilon(x)$为狄拉克(Dirac)函数, $H_\varepsilon(x)$为海氏(Heaviside)函数. $\delta _\varepsilon(x)$和$H_\varepsilon(x)$分别定义为
$ \begin{equation} \delta _\varepsilon (x) = \left\{\begin{array}{ll} \dfrac{1}{{2\varepsilon }}\left[ {1 + \cos \left( {\dfrac{{\pi x}}{\varepsilon }} \right)} \right], &|x| \le \varepsilon\\[2mm] 0, &|x| > \varepsilon \end{array} \right. \end{equation} $
(12) $ \begin{equation} {H_\varepsilon }\left( x \right) = \left\{ \begin{array}{ll} \dfrac{1}{2}\left( 1 + \dfrac{x}{\varepsilon } + \dfrac{1}{\pi }\sin \left( \dfrac{{\pi x}}{\varepsilon}\right)\right), &|x| \le \varepsilon\\ 1, &x > \varepsilon\\ 0, &x < - \varepsilon \end{array} \right. \end{equation} $
(13) 正则化水平集模型的目标是最小化能量函数(式(8)).合理地选择正则项$R_{p}(\phi)$, 也就是合理选择势能函数$p$, 可以使水平集函数$\phi( x, y)$保持符号距离函数性质, 即$|\nabla \phi |=1$.
令$|\nabla \phi |=s$, $p=(s-1)^{2}/{2}$, 则正则化项$R_{p}(\phi)$可表示为
$ \begin{equation} {R_p}\left( \phi \right) = \frac{1}{2}{\int_\Omega {\left( {\left| {\nabla \phi } \right| - 1} \right)} ^2}\text{d}x \end{equation} $
(14) $ \begin{equation} \frac{{\partial {R_p}}}{{\partial \phi }} = - {\rm{div}}\left( {{d_p}\left( {\left| {\nabla \phi } \right|} \right)\nabla \phi } \right) \end{equation} $
(15) 其中, $\rm{div}$是散度.
$ \begin{equation} {d_p}\left( s \right) = \frac{{p'\left( s \right)}}{s} \end{equation} $
(16) 在变分法中, 最小化一个能量函数$E(\phi)$的方法是找到梯度下降流方程的稳定解:
$ \begin{equation} \frac{{\partial \phi }}{{\partial t}} = - \frac{{\partial E}}{{\partial \phi }} \end{equation} $
(17) 式中, 负号表示水平集函数$\phi( C, t)$的方向与能量函数$E(\phi)$的偏导数相反, 是能量函数$E(\phi)$的最陡下降方向(最陡梯度下降流).
对能量函数式(8)求偏导:
$ \begin{equation} \frac{{\partial E}}{{\partial \phi }} = \mu \frac{{\partial {R_p}}}{{\partial \phi }} + \frac{{\partial {E_{\rm ext}}}}{{\partial \phi }} \end{equation} $
(18) 根据式(15)与式(17),
$ \begin{align} \frac{{\partial \phi }}{{\partial t}} = & - \mu \frac{{\partial {R_p}}}{{\partial \phi }} - \frac{{\partial {E_{\rm ext}}}}{{\partial \phi }} = \nonumber\\ &\mu {\rm{div}}\left( {{d_p}\left( {\left| {\nabla \phi } \right|} \right)\nabla \phi } \right) - \frac{{\partial {E_{\rm ext}}}}{{\partial \phi }} \end{align} $
(19) 式(19)就是距离正则水平集演化方程.
如果取$p=(s-1)^{2}/{2}$, 则$d_{p}(s)=1-(1/s)$, 距离正则水平集演化方程可以表示为
$ \begin{equation} \frac{{\partial \phi }}{{\partial t}} = \mu \left[ {{\nabla ^2}\phi - {\rm{div}}\left( {\frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}} \right)} \right] - \frac{{\partial {E_{\rm ext}}}}{{\partial \phi }} \end{equation} $
(20) 结合式(11), 距离正则化水平集演化方程可以表示为
$ \begin{align} \frac{{\partial \phi }}{{\partial t}} = &\mu {\rm{div}}\left( {{d_p}\left( {|\nabla \phi |} \right)\nabla \phi } \right) + \nonumber\\ &\lambda {\delta _\varepsilon }\left( \phi \right){\rm{div}}\left( {g\frac{{\nabla \phi }}{{|\nabla \phi |}}} \right) + \alpha g{\delta _\varepsilon }\left( \phi \right) \end{align} $
(21) 式中, 等号右边第一项$ \mu {\rm{div}}\left( {{d_p}\left( {|\nabla \phi |} \right)\nabla \phi } \right)$具有距离正则化的作用, 避免了一般水平集方法在演化过程中需要不断对零水平集函数进行初始化的问题.后两项为外部能量项, 其中第二项$\lambda {\delta _\varepsilon }\left( \phi \right){\rm{div}}\left( g\nabla \phi /|\nabla \phi | \right)$为零水平集曲线的能量项, 当它最小时就说明初始轮廓已演化至目标轮廓位置.第三项$\alpha g{\delta _\varepsilon }\left( \phi \right)$可以确定演化过程中曲线的方向. $\alpha\in {\bf R}$, 当$\alpha < 0$, 曲线将向外演化, 当$\alpha>0$时, 曲线将向内演化.
在水平集演化方程后两项中都用到了边缘检测函数$g$.在图像分割应用中, 边缘检测函数主要由图像梯度控制, 使曲线在图像梯度大的地方缓慢演化, 而在梯度小的地方快速演化. $g$的表达式为
$ \begin{equation} g := \frac{1}{{1 + |\nabla {G_\sigma }*I|^{2}}} \end{equation} $
(22) 其中, $G_\sigma$是方差为$\sigma$的高斯函数, $I$表示图像, $*$表示卷积运算.
2. 位错理论
位错(Dislocation)是原子的一种特殊组态.在材料科学中, 位错是指原子间局部不规则排序(晶体学缺陷).由于外力或内部应力的影响, 使晶体内部质点排列变形、原子行列间相互滑移, 不再符合理想晶格的有序排列而形成线状的缺陷, 称为位错.从几何角度分析, 可把位错视为材料内部已发生滑移部分和未发生滑移部分的分界线, 即线缺陷.从数学理论方面分析, 位错就是材料拓扑结构上的缺陷[14-15].
位错的理想形式主要有两种, 分别是刃位错(Edge dislocations)和螺位错(Screw dislocations)[16-17].刃位错是指晶体在大于屈服值的切应力$τ$的作用下发生滑移.已滑移部分和未滑移部分的交线叫位错线.位错线犹如砍入晶体中一把刀的刀刃, 所以这种位错称为刃位错(或棱位错).螺位错是指一个晶体的某一部分相对于其余部分发生滑移, 原子平面沿着一根轴线盘旋上升, 每绕轴线一周, 原子面上升一个晶面间距, 在中央轴线处形成一个螺位错.
位错在滑移面上受到垂直于位错线的作用力, 当此力足以克服运动阻力时, 位错便可以沿着作用力方向移动, 这种沿着滑移面移动的位错运动称为滑移.通常用柏格斯矢量(Burgers vector, 简称柏氏矢量)表示晶体滑移的方向和大小.柏氏矢量是一个反映由位错引起的点阵畸变大小的物理量.该矢量的模表示畸变的程度, 称为位错的强度.位错运动导致晶体滑移时, 滑移量大小即柏氏矢量$\pmb{b}$, 滑移方向即为柏氏矢量的方向.
位错在切应力作用下发生运动.在刃位错中, 由于位错移动的方向总是与位错线垂直, 故可设想有一个作用在位错线上且垂直于位错线的力造成了位错的移动.如果把位错这个原子组态抽象为一个实体, 本文把作用在这个实体上并使之运动的力称为位错组态力, 以$\pmb{F}_{\rm d}$表示.
假设作用在滑移面上的应力为$τ$, 当使长度为$l$的位错线移动距离$\text{d}s$之后, 晶体正好位移了位错的一个柏氏矢量$\pmb{b}$, 设此面上切应力使晶体滑移$b$所做的功为$W_1$:
$ \begin{equation} {W_1} = τ \cdot \left( {l\text{d}s} \right) \cdot b \end{equation} $
(23) 位错组态力$\pmb{F}_{\rm d}$所做的功为$W_2$:
$ \begin{equation} {W_2} = \pmb{F}_{\rm d} \cdot \left( {l\text{d}s} \right) \end{equation} $
(24) 根据虚功原理:
$ \begin{equation} W_1 =W_2 \end{equation} $
(25) 因此在外力场中单位长度位错的受力$\pmb{F}_{\rm d}$大小为
$ \begin{equation} \pmb{F}_{\rm d} =τ\cdot b \end{equation} $
(26) 由此可见, 对任意位错线, 在滑移面上沿柏氏矢量方向施加一个均匀应力$τ$, 则单位长度位错线上要受到一个作用力$\pmb{F}_{\rm d}$, 其大小为$τb$, 方向是沿该点的法线方向.位错线在这个作用力$\pmb{F}_{\rm d}$的驱使下移动.
在主动轮廓模型中, 活动曲线在速度函数$F$的驱动下沿着曲线法线方向演化, 导致曲线发生变形.在水平集方法中, 零水平集函数在图像梯度和边缘检测函数的驱使下向物体边界靠拢.它们与位错理论中位错线的移动有相似之处, 在刃位错中, 位错线在位错组态力$\pmb{F}_{\rm d}$的驱使下移动.
3. 基于位错理论的水平集方法
距离正则化水平集方法大大减小了对初始轮廓选取的依赖性, 选择简单的初始轮廓形式便可得到较好的分割结果.引入了惩罚项的距离正则化水平集方法, 避免了需要不断初始化零水平集函数而带来的计算量增加的问题.
在距离正则化水平集方法的图像分割中, 耦合了图像信息的边缘检测函数$g$, 在驱动初始轮廓向目标轮廓演化时起到了重要作用.边缘检测函数$g$主要依赖于图像的梯度信息, 对于图像梯度信息不明显的图像, 即包含弱边缘的图像, 采用距离正则化水平集方法无法得到令人满意的分割结果.
为了解决水平集方法对梯度信息不明显的弱边缘图像的分割问题, 本文引入材料物理学中的位错理论, 将水平集算法中曲线的移动看作刃位错中位错线的滑移过程, 并运用位错机制中的弹性性质、能量关系、位错动力学等理论, 推导出驱使初始轮廓曲线移动的位错组态力.
首先定义一个位错能量函数, 根据水平集方法的思想, 将二维图像平面放在一个三维空间中, 然后在这个三维空间中定义曲线的能量[18-19].将三维空间中的曲线$\gamma(s)$参数化, 与$\gamma(s)$相关的能量为
$ \begin{equation} E(\gamma ) = {\rm min}\int {\frac{1}{2}} {\left\| \pmb{w}(x, y, z) \right\|^2}\text{d}x\text{d}y\text{d}z \end{equation} $
(27) 约束条件为
$ \begin{equation} \nabla \times \mathit{\boldsymbol{w}} = \delta (\gamma )\tau \end{equation} $
(28) 其中, $\nabla$为哈密尔顿算子, $\pmb{w}=(w_1, w_2, w_3)$是在三维空间中作用在曲线上的应力, $\nabla \times \pmb{w}$是$\pmb{w}$在三维空间中的旋度, $\delta(\gamma)$是曲线的$\delta$函数, $τ$是$\delta(\gamma)$上的切向量.
应力$\pmb{w}$可如下求解:
$ \begin{equation} \mathit{\boldsymbol{w}}(x, y, z) = - \frac{1}{{4\pi }}\int_\gamma {\frac{{\mathit{\boldsymbol{r}} \times \text{d}\mathit{\boldsymbol{l}}}}{{{r^3}}}} \end{equation} $
(29) 其中, $\text{d}\pmb{l}$是曲线$\gamma$上的线元, $\pmb{r}=(x-x(s), y=y(s), z-z(s))$是点$(x, y, z)$与点$(x(s), y(s), z(s))$之间的向量, 这两点间的距离为
$ \begin{gather} r = \sqrt {{{(x - x(s))}^2} + {{(y - y(s))}^2} + {{(z - z(s))}^2}} \end{gather} $
(30) 位错线可以看作是零水平集$\phi$表示的曲线, 曲线的单位切向量$τ$为
$ \begin{equation} τ = \frac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}} \end{equation} $
(31) 曲线$\gamma$的$\delta$函数可以表示为
$ \begin{equation} \delta (\gamma ) = \delta (\phi )\left| {\nabla \phi } \right| \end{equation} $
(32) 曲线的线元$\text{d}\pmb{l}$可以表示为
$ \begin{equation} {\rm d}{\pmb{l}} = τ\delta (\gamma )\text{d}x\text{d}y\text{d}z = \nabla \phi \delta (\phi )\text{d}x\text{d}y\text{d}z \end{equation} $
(33) 把式(33)代入到式(29)可得:
$ \begin{equation} {\pmb{w}}(x, y, z) = - \frac{1}{{4\pi }}\int_\gamma {\frac{{{\pmb{r}} \times \nabla \phi \delta (\phi )}}{{{r^3}}}} \text{d}x\text{d}y\text{d}z \end{equation} $
(34) 位错线在位错组态力$\pmb{F}_{\rm d}$的驱使下沿着曲线法线方向移动.设$I$是图像灰度函数, 图像中目标边界的方向可以定义为
$ \begin{equation} {\pmb{n}} = \frac{{\nabla ({G_\sigma } * I)}}{{\left| {\nabla ({G_\sigma } * I)} \right|}} \end{equation} $
(35) 虚拟的位错组态力$\pmb{F}_{\rm d}$的大小可以表示为
$ \begin{align} &\pmb{F}_{\rm d} = {\pmb{w}}(x, y, z) \times {\pmb{n}} = \nonumber\\&\qquad - \frac{1}{{4\pi }}\int_\gamma {\frac{{{\pmb{r}} \times \nabla \phi \delta (\phi )\nabla ({G_\sigma } * I)}}{{{r^3}\left| {\nabla ({G_\sigma } * I)} \right|}}} \text{d}x\text{d}y\text{d}z \end{align} $
(36) 位错线或零水平集曲线$\phi$在位错组态力$\pmb{F}_{\rm d}$的驱使下沿着曲线法线方向移动.当位错线接近目标边界时, 图像梯度较大, 位错组态力$\pmb{F}_{\rm d}$也较大, 远离目标边界时, 图像梯度较小, 位错组态力$\pmb{F}_{\rm d}$也较小.本文方法的目标是使曲线在接近目标边界时缓慢演化, 而在远离目标边界时快速演化, 即当活动曲线靠近目标边界时, 演化速度逐渐减小直到停止.
所以本文建立如下的速度停止函数(边界检测函数):
$ \begin{equation} g_{\rm d} = \dfrac{1}{{1 - \dfrac{1}{{4\pi }}\displaystyle\int_\gamma {\dfrac{{{\pmb{r}} \times \nabla \phi \delta (\phi )\nabla ({G_\sigma } * I)}}{{{r^3}\left| {\nabla ({G_\sigma } * I)} \right|}}} \text{d}x\text{d}y\text{d}z} } \end{equation} $
(37) 本文用速度停止函数$g_{\rm d}$取代边缘检测函数$g$, 则引入位错理论后新的距离正则化水平集函数演化方程如式(38)所示.
$ \frac{\partial\phi}{\partial t} = \mu {\rm{div}}\left( {{d_p}\left( {|\nabla \phi |} \right)\nabla \phi } \right)+\\\lambda {\delta _\varepsilon } \left( \phi \right){\rm{div}}\left( {g_{\rm d}\dfrac{{\nabla \phi }} {{|\nabla \phi |}}} \right)+\alpha g_{\rm d}{\delta _\varepsilon }\left( \phi \right) =\mu {\rm{div}}\left( {{d_p}\left( {|\nabla \phi |} \right)\nabla \phi } \right) + \lambda {\delta _\varepsilon }( \phi )\times\nonumber\\ ~~~~~~~{\rm{div}}\left\{ \left[ \dfrac{1}{ 1 - \dfrac{1}{{4\pi }}\displaystyle\int_\gamma {\dfrac{{{\pmb{r}} \times \nabla \phi \delta (\phi )\nabla ({G_\sigma } * I)}}{{{r^3}\left| {\nabla ({G_\sigma } * I)} \right|}}} \text{d}x\text{d}y\text{d}z } \right]\dfrac{{\nabla \phi }} {{|\nabla \phi |}}\right\}+ \nonumber\\ ~~~~~~~\alpha \left[ \dfrac{1}{ 1 - \dfrac{1}{{4\pi }}\times\displaystyle\int_\gamma {\dfrac{{{\pmb{r}} \times \nabla \phi \delta (\phi )\nabla ({G_\sigma } * I)}}{{{r^3}\left| {\nabla ({G_\sigma } * I)} \right|}}} \text{d}x\text{d}y\text{d}z } \right]{\delta _\varepsilon }( \phi ) $
(38) 在距离正则化水平集方法中, 水平集函数的演化主要取决于图像梯度.在图像分割中, 如果目标边界比较模糊且灰度分布不均匀, 使用传统的水平集方法难以在弱边缘处准确地分割出目标轮廓.另外, 图像梯度对噪声非常敏感, 因此经常会造成边界泄漏.本文提出的基于位错理论的水平集方法中, 水平集函数的演化主要取决于位错组态力$\pmb{F}_{\rm d}$和速度停止函数$g_{\rm d}$.位错组态力包含了更多的图像信息.其中, $\nabla ({G_\sigma } * I)$表示了图像的梯度信息. $r$表示了活动曲线之间的距离信息.
首先, 位错组态力包含了图像梯度信息, 使得本文方法可以像传统水平集方法一样具有正确分辨边缘的能力.另外, 位错组态力中比传统水平集方法多了一个参数$r$, 反映了位错组态力的作用范围.靠近位错线(活动曲线)的区域位错组态力较大, 反之较小.在图像弱边缘处或灰度不均匀区域, 位错组态力由于包含了更多的图像信息, 与传统水平集方法相比可以更精确地控制水平集函数的演进.
其次, 位错组态力的表达式是一个积分公式, 所以实际上是表示了单位长度曲线上的平均受力情况.这种平均力可以在曲线受力不均的情况下抑制局部受力突变对整体的影响, 并对曲线的形状进行约束, 使得曲线在较窄的图像弱边缘处不至于发生局部边界泄漏, 或者在由于噪声引起的局部图像梯度异常的情况下避免曲线停止演进.
4. 实验结果
为了验证基于位错理论的水平集方法的有效性, 本文进行了图像分割实验, 并与其他经典水平集算法的分割结果进行了比较.实验环境如下: CPU型号为AMD Phenom Ⅱ X4 B97, CPU频率3.2 GHz, 8 GB内存.操作系统为Windows 7, 编程环境为Matlab 2014a.部分实验测试图像来自加州大学伯克利分校的公开图像数据库.
本文提出的基于位错理论的水平集方法中, 水平集函数演化方程(式(38))含有三个控制参数, 分别是$\mu$、$\lambda$和$\alpha$.其中参数$\mu$用于控制水平集函数的正则化, 本文实验中设置为0.005.参数$\lambda$用于控制水平集曲线的演化, 本文实验中设置为1.非零参数$\alpha$用于控制水平集曲线的收缩与膨胀.本文实验中, 如果初始轮廓在目标对象内部, 设置$\alpha=-1$; 否则, $\alpha=1$.本文实验中, 高斯核函数标准差$\sigma=9$, 时间步长$\Delta t =0.2$.初始轮廓曲线按式(38)进行演化.在距离正则化项(内部能量项)的约束下整个迭代过程无需重新初始化.当曲线靠近目标边界时, 速度停止函数(式(37))快速减小, 曲线演进缓慢直到停止.当达到最大预设迭代次数或曲线已收敛时, 迭代停止.判断曲线收敛的条件是前后间隔10次迭代中曲线内部面积之差小于一个阈值, 本文中阈值取为0.02.
图 1是本文算法的图像分割结果.第一列是原始图像.第二列、第三列和第四列分别是迭代20次、50次和100次后图像的分割结果.实验结果表明, 本文算法可以较好地跟踪出目标轮廓.
本文算法也适用于含有多个目标主体的图像分割. 图 2是本文算法的多目标主体图像分割结果.实验结果表明, 即使在目标主体较多、目标主体之间距离较远的情况下本文算法也可以取得较好的分割效果.
图 3是在加入噪声的情况下本文算法的图像分割结果.第一列是无噪声图像, 第二列的图像中加入了方差为4的高斯白噪声, 第三列的图像中加入了椒盐噪声, 第四列的图像中加入了脉冲噪声, 第五列的图像中同时加入了高斯白噪声和椒盐噪声.实验结果表明, 在加入上述各种噪声的情况下, 本文算法仍可以较好地跟踪出目标轮廓.
本文对在不同参数情况下的迭代次数和运算时间进行了比较.实验图像为模拟树叶图像, 默认参数为$\lambda=1$、$\mu=0.005$、$\alpha=-1$、$\sigma=9$、$\Delta t=0.2$.每次实验中修改其中一个参数, 而其他参数保持默认值.实验结果如图 4和图 5所示.
在其他参数不变的情况下, $\lambda$增大时, 迭代次数减少.水平集曲线或位错线能量较大, 图像分割比较精确.但是$\lambda$较大时, 会导致迭代次数急剧增加, 迭代无法收敛.在其他参数不变的情况下, $\mu$增大时, 迭代次数和运算时间出现波动.当$\mu$较大时, 水平集函数无法正确正则化, 迭代缓慢.当高斯核函数标准差$\sigma$较小时, 容易发生边界泄漏.当$\sigma$增大时, 迭代变慢, 运算时间比较长.
本文也把采用本文算法的图像分割实验结果与其他经典水平集算法的分割结果进行了比较.这些经典水平集算法包括Caselles算法[1]、CV (Chan-Vese)模型[2]、Bernard算法[3]、Shi算法[4]、Min算法[6]以及距离正则化水平集方法(DRLSE)[11].
当图像的边缘有比较好的对比度, 边缘梯度信息明显时, 经典水平集算法都可以有较好的分割效果.但是在实际应用中, 经常要处理边缘梯度信息不明显的图像分割问题, 即弱边缘图像分割.下面通过实验来验证本文算法对弱边缘的图像分割效果. 图 6和表 1是对具有弱边缘的模拟图像分割实验的对比结果. 图 6中, (a)初始轮廓; (b) Caselles算法[1], 收敛阈值设置为2, 推动曲线运动的气球力设置为1; (c) CV算法[2], 收敛阈值设置为2, 控制水平集函数正则化的曲率项设置为0.2; (d) Bernard算法[3], 收敛阈值设置为2, 控制水平集曲线演化平滑程度的尺度因子设置为1; (e) Shi算法[4], 曲线演化迭代次数设置为300, 正则化迭代次数设置为30, 高斯滤波器方差设置为4; (f) DRLSE算法[11], $\lambda=1$、$\mu=0.003$、$\alpha=-3$、高斯核函数方差$\sigma$设置为7; (g) Min算法[6]; (h)本文算法, $\lambda=1$、$\mu=0.005$、$\alpha=-1$、$\sigma=9$、$\Delta t=0.2$. 图 7所示实验对比实验中用到的各个经典水平集算法的参数设置均与此相同.实验结果表明, 只有本文算法正确地跟踪出了主体的边界, 其他算法均在弱边缘处停止演化或发生了边界泄漏.
图 6 弱边缘图像分割对比实验结果Fig. 6 The comparison of experimental results of weak boundary image segmentation经典的水平集算法中, 除了少数算法(如CV算法)采用曲线能量泛函最小化进行曲线演化以外, 多数算法都依赖于图像梯度信息进行曲线演化, 如Caselles算法和DRLSE算法等.在具有弱边缘的图像中, 弱边缘处的图像梯度较小, 采用基于图像梯度信息的水平集算法进行图像分割难以取得良好效果.这些算法在弱边缘处会出现不能完整地跟踪出边界以及边界泄漏现象.本文提出的基于位错理论的距离正则化水平集方法, 位错线或零水平集曲线$\phi$在位错组态力$\pmb{F}_{\rm d}$和速度停止函数$g_{\rm d}$的驱使下沿着曲线法线方向移动.因为位错组态力$\pmb{F}_{\rm d}$和速度停止函数$g_{\rm d}$综合了图像梯度、曲线能量、距离等信息, 因此本文算法对弱边缘图像的分割更加有效.从图 6中可以看出, 采用本文方法对模拟的弱边缘图像进行分割时, 没有在弱边缘处出现边界泄漏现象.
为了能够定量地对本文算法的图像分割性能进行评价并与其他经典算法进行更有效的比较, 本文参考了Jaccard similarity (JS)等图像分割性能指标[20-22], 提出如下两个图像分割性能评价指标.
1) 面积重叠误差
$ \begin{eqnarray} E_{\rm overlap}=\left( 1- \frac{A\bigcap M}{A\bigcup M} \right)\times 100 \% \end{eqnarray} $
(39) 其中, $A$和$M$分别代表水平集算法分割区域和手动分割区域的像素数量.
这个面积重叠误差指标取值为零时表示完美分割, 取值为100 $\%$时表示算法分割区域与真实目标区域完全不重叠.
2) 边界平均距离
设$C$是算法分割出的目标边界, 对于边界$C$上的每一点$P_{i}~ ( i=1, 2, \cdots, N)$计算$P_{i}$与真实目标边界$S$的距离$\text{dist}(P_{i}, S)$, 然后求平均值.
$ \begin{eqnarray} D_{\rm mean}=\frac{1}{N}\sum\limits_{i=1}^N \text{dist}(P_{i}, S) \end{eqnarray} $
(40) 基于上述两个图像分割性能评价指标, 本文算法与其他经典算法的对比实验结果如图 7和表 2所示.
表 2 图 7所示实验中各算法的迭代次数、运算时间、面积重叠误差以及边界平均距离的对比Table 2 The comparison of the iteration number, the operation time, area overlap error, and average boundary distance for the experiments shown in Fig. 7算法 迭代次数 运算时间(秒) $E_{\rm overlap}$ $D_{\rm mean}$ Caselles[1] 5 000 159.536302 - - CV[2] 2 915 103.690744 21.19 $\%$ 13.62 Bernard[3] 18 121.678605 14.68 $\%$ 9.54 Shi[4] 798 184.167387 18.81 $\%$ 11.71 DRLSE[11] 2 320 383.120423 0.42 $\%$ 0.31 Min[6] 500 71.325241 - - Wang[7] 2 032 258.625021 0.65 $\%$ 0.48 本文算法 1 865 285.869702 0.36 $\%$ 0.27 图 7所示实验中各算法的迭代次数、运算时间、面积重叠误差以及边界平均距离如表 2所示.实验结果表明, DRLSE算法、Wang算法和本文算法完整地跟踪出了血管轮廓.本文算法在面积重叠误差、边界平均距离性能指标上与DRLSE算法和Wang算法基本相当, 在迭代次数、运算时间上优于DRLSE算法.
5. 结论
1) 本文把材料科学中的位错理论引入到水平集方法中, 水平集函数演化被看作刃位错中位错线的滑移过程, 构建了新的速度停止函数(边界检测函数)和基于位错理论的距离正则化水平集演化方程.
2) 运用位错机制中的弹性性质、能量关系、位错动力学等理论, 推导出驱使初始轮廓曲线移动的位错组态力.
3) 采用基于距离正则化的算法, 使水平集函数自我约束为符号距离函数, 不需要不断重新初始化零水平集函数, 提高了运算效率.
4) 本文算法中零水平集曲线在位错组态力的驱使下移动, 与多数水平集算法依赖图像梯度信息进行曲线演化不同, 本文算法中的位错组态力不仅包括图像梯度信息, 也包括曲线能量信息、位错组态力的作用范围等信息, 因此采用本文算法可以有效地避免在局部图像梯度异常的情况下发生曲线停止演进的现象, 或者在弱边缘处由于图像梯度较小发生局部边界泄漏的现象.
-
图 6 弱边缘图像分割对比实验结果
Fig. 6 The comparison of experimental results of weak boundary image segmentation
表 1 图 6所示实验中各算法的迭代次数和运算时间的对比
Table 1 The comparison of the iteration number and the operation time for the experiments shown in Fig. 6
下载: 导出CSV
表 2 图 7所示实验中各算法的迭代次数、运算时间、面积重叠误差以及边界平均距离的对比
Table 2 The comparison of the iteration number, the operation time, area overlap error, and average boundary distance for the experiments shown in Fig. 7
算法 迭代次数 运算时间(秒) $E_{\rm overlap}$ $D_{\rm mean}$ Caselles[1] 5 000 159.536302 - - CV[2] 2 915 103.690744 21.19 $\%$ 13.62 Bernard[3] 18 121.678605 14.68 $\%$ 9.54 Shi[4] 798 184.167387 18.81 $\%$ 11.71 DRLSE[11] 2 320 383.120423 0.42 $\%$ 0.31 Min[6] 500 71.325241 - - Wang[7] 2 032 258.625021 0.65 $\%$ 0.48 本文算法 1 865 285.869702 0.36 $\%$ 0.27
下载: 导出CSV
-
[1] Caselles V, Kimmel R, Sapiro G. Geodesic active contours. International Journal of Computer Vision, 1997, 22(1):61-79 doi: 10.1023/A:1007979827043 [2] Chan T F, Vese L A. Active contours without edges. IEEE Transactions on Image Processing, 2001, 10(2):266-277 http://dl.acm.org/citation.cfm?id=2320071 [3] Bernard O, Friboulet D, Thévenaz P, Unser M. Variational b-spline level-set:a linear filtering approach for fast deformable model evolution. IEEE Transactions on Image Processing, 2009, 18(6):1179-1191 doi: 10.1109/TIP.2009.2017343 [4] Shi Y G, Karl W C. A real-time algorithm for the approximation of level-set-based curve evolution. IEEE Transactions on Image Processing, 2008, 17(5):645-656 doi: 10.1109/TIP.2008.920737 [5] Li C M, Huang R, Ding Z H, Gatenby J C, Metaxas D N, Gore J C. A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI. IEEE Transactions on Image Processing, 2011, 20(7):2007-2016 doi: 10.1109/TIP.2011.2146190 [6] Min H, Jia W, Wang X F, Zhao Y, Hu R X, Luo Y T, Xue F, Lu J T. An intensity-texture model based level set method for image segmentation. Pattern Recognition, 2015, 48(4):1547-1562 doi: 10.1016/j.patcog.2014.10.018 [7] Wang X F, Min H, Zou L, Zhang Y G. A novel level set method for image segmentation by incorporating local statistical analysis and global similarity measurement. Pattern Recognition, 2015, 48(1):189-204 doi: 10.1016/j.patcog.2014.07.008 [8] Zhao Y Q, Wang X H, Wang X F, Shih F Y. Retinal vessels segmentation based on level set and region growing. Pattern Recognition, 2014, 47(7):2437-2446 doi: 10.1016/j.patcog.2014.01.006 [9] 周则明, 孟勇, 黄思训, 胡宝鹏.基于能量最小化的星载SAR图像建筑物分割方法.自动化学报, 2016, 42(2):279-289 http://www.aas.net.cn/CN/abstract/abstract18817.shtmlZhou Ze-Ming, Meng Yong, Huang Si-Xun, Hu Bao-Peng. Building segmentation of spaceborne SAR images based on energy minimization. Acta Automatica Sinica, 2016, 42(2):279-289 http://www.aas.net.cn/CN/abstract/abstract18817.shtml [10] 张迎春, 郭禾.基于粗糙集和新能量公式的水平集图像分割.自动化学报, 2015, 41(11):1913-1925 http://www.aas.net.cn/CN/abstract/abstract18766.shtmlZhang Ying-Chun, Guo He. Level set image segmentation based on rough set and new energy formula. Acta Automatica Sinica, 2015, 41(11):1913-1925 http://www.aas.net.cn/CN/abstract/abstract18766.shtml [11] Adalsteinsson D, Sethian J A. A fast level set method for propagating interfaces. Journal of Computational Physics, 1995, 118(2):269-277 doi: 10.1006/jcph.1995.1098 [12] Li C M, Kao C Y, Gore J C, Ding Z H. Minimization of region-scalable fitting energy for image segmentation. IEEE Transactions on Image Processing, 2008, 17(10):1940-1949 doi: 10.1109/TIP.2008.2002304 [13] Li C M, Xu C Y, Gui C F, Fox M D. Distance regularized level set evolution and its application to image segmentation. IEEE Transactions on Image Processing, 2010, 19(12):3243-3254 doi: 10.1109/TIP.2010.2069690 [14] Lemoine G, Delannay L, Idrissi H, Colla M S, Pardoen T. Dislocation and back stress dominated viscoplasticity in freestanding sub-micron Pd films. Acta Materialia, 2016, 111:10-21 doi: 10.1016/j.actamat.2016.03.038 [15] Nagasako N, Asahi R, Isheim D, Seidman D N, Kuramoto S, Furuta T. Microscopic study of gum-metal alloys:a role of trace oxygen for dislocation-free deformation. Acta Materialia, 2016, 105:347-354 doi: 10.1016/j.actamat.2015.12.011 [16] Guiu F. The stress dependence of the dislocation velocity in molybdenum. Physica Status Solidi, 2016, 25(1):203-207 http://www.researchgate.net/publication/243143363_The_stress_dependence_of_the_dislocation_velocity_in_molybdenum [17] Yuan R, Beyerlein I J, Zhou C Z. Statistical dislocation activation from grain boundaries and its role in the plastic anisotropy of nanotwinned copper. Acta Materialia, 2016, 110:8-18 doi: 10.1016/j.actamat.2016.02.064 [18] Luo Y S, Chung A C S. Nonrigid image registration with crystal dislocation energy. IEEE Transactions on Image Processing, 2013, 22(1):229-243 http://ieeexplore.ieee.org/document/6220249/ [19] Xiang Y, Chung A C S, Ye J. An active contour model for image segmentation based on elastic interaction. Journal of Computational Physics, 2006, 219(1):455-476 doi: 10.1016/j.jcp.2006.03.026 [20] Badakhshannoory H, Saeedi P. A model-based validation scheme for organ segmentation in CT scan volumes. IEEE Transactions on Biomedical Engineering, 2011, 58(9):2681-2693 doi: 10.1109/TBME.2011.2161987 [21] Heimann T, van Ginneken B, Styner M A, Arzhaeva Y, Aurich V, Bauer C, et al. Comparison and evaluation of methods for liver segmentation from CT datasets. IEEE Transactions on Medical Imaging, 2009, 28(8):1251-1265 doi: 10.1109/TMI.2009.2013851 [22] Zhang K H, Zhang L, Yang M H. Fast compressive tracking. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2014, 36(10):2002-2015 doi: 10.1109/TPAMI.2014.2315808 期刊类型引用(10)
1. 孟小燕. 融入自学习与多领导者策略的改进鲸鱼优化算法及多阈值图像分割. 计算机应用与软件. 2025(01): 205-212+233 . 
百度学术2. 孙超. 伪对立学习与差分进化的改进鲸鱼优化算法及图像分割应用. 计算机应用与软件. 2024(09): 265-272+287 . 百度学术3. 柳长安，冯雪菱，孙长浩，赵丽娟. 基于改进麻雀算法的最大2维熵分割方法. 激光技术. 2022(02): 274-282 . 百度学术4. 李建华，王艺衡，杨来侠，李素丽，张聚良. 基于Canny算子和距离正则化水平集的乳腺植入物图像分割算法. 计算机应用与软件. 2022(10): 212-216+245 . 百度学术5. 邢致恺，贾鹤鸣，宋文龙. 基于莱维飞行樽海鞘群优化算法的多阈值图像分割. 自动化学报. 2021(02): 363-377 . 本站查看6. 宋杰，肖亮，练智超. 级联稀疏卷积与决策树集成的病理图像细胞核分割方法. 自动化学报. 2021(02): 378-390 . 本站查看7. 龚进慧，张贵仓. 自适应尺度参数的非均匀图像快速水平集分割. 计算机应用与软件. 2021(04): 231-238+304 . 百度学术8. 侯木舟，陈英皞. 基于改进的水平集模型工业零件图像分割方法. 中北大学学报(自然科学版). 2019(02): 155-160+166 . 百度学术9. 符艳军，孙开锋. 基于自动阈值分割的管内火焰大小计算. 自动化与仪器仪表. 2019(04): 118-121 . 百度学术10. 姚志文. 基于频域特征的机械焊缝缺陷图像分割仿真. 计算机仿真. 2019(12): 159-162+463 . 百度学术其他类型引用(8)
-
下载:
计量
- 文章访问数: 2129
- HTML全文浏览量: 237
- PDF下载量: 926
- 被引次数: 18
