在数字图像处理的研究和应用中, 图像分割问题是一类非常重要的问题.图像分割就是利用相关数学模型及计算机实现方法将图像中的目标和背景区分开来, 以便进一步的研究和分析.它广泛地应用于工农业生产自动化、字符识别、医学以及遥感图像分析、故障处理等方面.至今, 在各种文献中提出的图像分割算法也多种多样, 归结一下主要有以下几种类型: 1) 基于阈值的分割方法, 比较常用的有大律法(Ostu, 1978)、最小误差法(Kittler, 1986) 等. 2) 基于边缘的分割方法, 基于边缘的分割方法指的是基于灰度值的边缘检测, 它是建立在边缘灰度值会呈现出阶跃型或屋顶型变化这一观测基础上的方法. 3) 基于区域的分割方法, 此类方法是将图像按照相似性准则分成不同的区域, 主要包括种子区域生长法、区域分裂合并法和分水岭法等几种类型. 4) 基于能量泛函的分割方法, 该方法的基本思想是用连续曲线来表达目标边缘, 并定义一个能量泛函使得其自变量包括边缘曲线, 因此分割过程就转变为求解能量泛函的最小值过程, 一般可以通过求解函数对应的欧拉方程来实现, 能量达到最小时的曲线位置就是目标的轮廓所在.水平集方法就是一种基于能量泛函的图像分割方法, 由于其鲁棒性和效率高等特点, 使其受到越来越多的国内外研究者重视.该方法的核心思想是将图像空间中的二维曲线嵌入到三维曲面中作为其零水平集.当曲面演变时, 嵌入到三维曲面中的零水平集曲线也随之改变, 最终只要确定零水平集即可确定移动曲面演化的结果. C-V模型^[1]是由Chan和Vese在2001年提出的一种区域性的水平集方法, 它利用曲线的内外灰度均值来促进零水平集的演化.该模型对灰度均匀的图像很有效, 且分割不受噪声影响^[2].但该模型在每次迭代时都要检验水平集函数是否为符号距离函数, 从而计算量很大^[3]; 其次C-V模型由于没有考虑图像的边缘信息, 从而边缘控制能力较弱, 不利于分割一些弱边缘或边缘模糊的图像^[4].于是很多作者提出了一些改进的C-V模型.文献[5]在C-V模型里加了一项内部能量项来约束水平集函数逼近符号距离函数, 从而避免了水平集函数的重新初始化过程, 文献[6-8]在C-V模型里加入了图像边缘信息, 提高了图像的边缘捕捉能力, 文献[9]利用窄带法进行迭代实验, 提高了分割效率.小波变换由于其紧支撑性和良好的去噪性能等特点使其在图像分割也发挥着越来越重要的作用.本文在C-V模型和小波变换的基础上, 提出了一种添加图像边缘信息和无需重新初始化的水平集分割模型, 并在数值实验时, 利用改进的窄带法来实现, 数值实验表明, 该方法比传统的水平集方法具有更高的分割效率和更好的分割效果.

水平集方法的基本思想是将图像空间中的二维曲线嵌入到三维曲面中作为其零水平集.当曲面演变时, 嵌入到三维曲面中的零水平集曲线也随之改变, 最终只要确定零水平集即可确定移动曲面演化的结果. C-V模型是由Chan和Vese在2001年提出的一种区域性的水平集方法, 它利用曲线的内外灰度均值来促进水平集的演化. C-V模型的能量泛函如下所示:

$\begin{array}{l} E({c_1},{c_2},C) = \alpha \int\limits_{{\rm{inside}}(C)} {{{\left| {{u_0}(x,y) - {c_1}} \right|}^2}} dxdy + \beta \\ \int\limits_{{\rm{outside}}(C)} {{{\left| {{u_0}(x,y) - {c_2}} \right|}^2}} dxdy + \mu \cdot {\rm{Length}}(C), \end{array}$

其中$C$表示光滑闭合的轮廓曲线, $u_0 $表示源图像, $c_1 $和$c_2 $分别是演化曲线$C$内部和外部的图像灰度均值, 一般情形下取$\alpha =\beta =1$, $\mu $为正常数.引入水平集函数$\varphi (x, y)$来代替演化曲线$C$, 并且规定点$(x, y)$在$C$的内部时, $\varphi (x, y)>0$; 点$(x, y)$在$C$的外部时, $\varphi (x, y)＜0$; 点$(x, y)$在$C$上时, $\varphi (x, y)=0$.从而上述能量泛函可写为如下形式:

$\begin{eqnarray*} \displaystyle E(c_1, c_2, C)&=&\mu \displaystyle\int\limits_\Omega {\delta _\varepsilon (\varphi (x, y))\left| {\nabla \varphi (x, y)} \right|} dxdy+\alpha \int\limits_\Omega {\left| {u_0 (x, y)-c_1 } \right|} ^2H_\varepsilon (\varphi (x, y))dxdy \\ && +\beta \displaystyle\int\limits_\Omega {\left| {u_0 (x, y)-c_2 } \right|} ^2(1-H_\varepsilon (\varphi (x, y)))dxdy, \end{eqnarray*}$

其中$H_\varepsilon (z)$和$\delta _\varepsilon (z)$分别是海氏函数$H(z)=\left\{\begin{array}{*{20}c} {1, z\ge 0}, \hfill \\ {0, z\le 0} \hfill \\ \end{array}\right. $和狄拉克函数$\delta (z)=\dfrac{d}{dz}H(z)$的正则化形式, $\Omega $表示图像域.上述能量泛函最小化问题可以通过求解能量泛函对应的欧拉方程来实现, 从而得到如下的水平集演化方程

$\dfrac{\partial \varphi }{\partial t}=\delta _\varepsilon (\varphi )[\mu \cdot \mbox{div}(\dfrac{\nabla \varphi }{\left| {\nabla \varphi } \right|})-\alpha (u_0-c_1 )^2+\beta (u_0-c_2 )^2], $

其中$\varphi (0, x, y)=\varphi _0 (x, y) \in \Omega $.

上式中的灰度均值$c_1 $和$c_2 $可分别在每次迭代中采用如下的方式进行更新

$\begin{equation} \label{eq1} c_1 (\varphi )=\frac{\displaystyle\int\limits_\Omega {u_0 } (x, y)H_\varepsilon (\varphi (x, y))dxdy}{\displaystyle\int\limits_\Omega {H_\varepsilon (\varphi (x, y))dxdy} }, \quad c_2 (\varphi )=\frac{\displaystyle\int\limits_\Omega {u_0 } (x, y)(1-H_\varepsilon (\varphi (x, y)))dxdy}{\displaystyle\int\limits_\Omega {(1-H_\varepsilon (\varphi (x, y)))dxdy} }. \end{equation}$

(2.1)

小波分析由于具有良好的时频局部性和多分辨率分析的功能, 使其在图像处理中得到广泛的应用.利用小波分析可以实现对图像进行编码压缩、边缘提取、滤波、添加数字水印等等.

首先利用二维离散小波对原图像进行$n$级MALLAT分解, 然后计算第$i$层的高频分量$H_i $, 它由水平、竖直和对角三部分组成, 对$H_i $做归一化处理, 使之数据更加规范, 记处理后的数据为$H_i ^\prime $, 将$H_i ^\prime $放大到原图像大小, 记为$cH{ }_i^\prime $, 从而可以得到原图像的平均高频信息量: $cH=\frac{1}{n}\sum\limits_i {cH_i } ^\prime $.边界的模糊程度能够从高频分量的幅值体现, 如果边界比较清晰, 高频分量幅值较大; 反之, 边界较模糊, 高频分量幅值较小, 甚至在当前级别的小波分解中不明显, 需要进行更高级别的分解, 这就是小波变换的多分辨率分析的优势.

小波变换的模极大值点对应于图像的边缘点, 所以可以利用检测模极大值点来获得图像的边缘点.首先对原图像上的每一个点进行二维离散小波变换并计算出其模值和幅角, 找出模极大值点对应于原图像的突变点, 记下其位置, 得到原图像可能的轮廓边缘点和轮廓边缘的感兴趣区域, 从而在利用C-V模型分割图像时, 初始曲线的设置可以通过连接图像的边缘点及一些图像中点得到, 这样设置的好处在于初始曲线已经位于图像边缘附近, 从而只需较少的迭代次数就可以得到分割结果.此处只需得到可能的边缘点和图像边缘轮廓的感兴趣区域, 从而去指导C-V模型中初始轮廓线的设置以及后续窄带的设置; 如果要得到精确的边缘点, 可对原图像进行多尺度分解, 然后在不同的尺度上分别检测边缘点, 最后通过插值得到原图像的边缘点; 因为在不同的尺度上边缘点和噪声点的表现特性不同, 从而可以消去噪声对边缘检测的影响^[7].

C-V模型对噪声不敏感, 并且该模型对边缘模糊的图像以及含有内部轮廓的图像很有效^[2].但C-V模型也存在一些不足, 首先C-V模型没有利用图像的边缘信息, 从而对图像的局部控制能力较弱, 于是在分割一些灰度不均匀的图像时效果不好^[3]; 其次C-V模型在迭代过程中需要不断地重新初始化水平集函数, 计算量很大^[6].基于此, 本文在C-V模型的基础上, 提出了一种加入边缘信息、并且无需重新初始化的改进模型.首先由于水平集函数$\varphi $通常是取由初始轮廓线生成的符号距离函数, 而水平集函数在迭代时往往会发生退化, 使其不再保持为符号距离函数, 所以在进行若干次迭代之后, 必须水平集函数$\varphi $重新初始化, 使其恢复为对当前零水平集的符号距离函数.但重新初始化的计算量很大, 从而导致分割效率较低.为了避免重新初始化, 很多文献[5]提出了一种改进方案, 即在关于水平集函数$\varphi $的能量泛函中增加一项内部约束能量:

$E(\varphi )=\int\limits_\Omega {\frac{1}{2}} (\left| {\nabla \varphi } \right|-1)dxdy, $

它的梯度下降流为

$\frac{\partial \varphi }{\partial t}=\Delta \varphi-\mbox{div}(\frac{\nabla \varphi }{\left| {\nabla \varphi } \right|}).$

其次传统的C-V模型只考虑了图像的全局区域信息, 而没有考虑图像的边缘梯度信息, 从而水平集的演化速度和边缘信息检测会受到影响, 从而导致在演化时往往容易陷入到局部极小值, 于是考虑在能量泛函中引入边缘检测函数^[6]来替换上式中的$\delta _\varepsilon (\varphi )$.利用3.1中得到图像平均高频分量构造边缘检测函数$g(cH)=\frac{1}{1+\left| {\nabla cH} \right|^2}$.从而得到改进后的C-V水平集演化方程如下:

$\begin{equation} \label{eq2} \frac{\partial \varphi }{\partial t}=g(\left| {\nabla G_\sigma \ast u} \right|)[\mu \cdot \mbox{div}(\frac{\nabla \varphi }{\left| {\nabla \varphi } \right|})-\alpha (u_0-c_1 )^2+\beta (u_0-c_2 )^2]+\Delta \varphi -\mbox{div}(\frac{\nabla \varphi }{\left| {\nabla \varphi } \right|}). \end{equation}$

(3.1)

在使用水平集模型对图像进行分割时, 需要计算每个时刻全图像范围内的所有网格点到当前轮廓曲线的距离, 然后依次连接距离为零的点来获得新的演化曲线, 如此反复进行下去直到分割结束.该方法计算量非常大, 算法复杂度高, 设$N$为图像离散后的网格点数, 则该算法的复杂度为$O(N^3)$.为了提高模型的分割速度和效率, 文献[4]提窄带水平集算法.该方法将距离计算限制在以当前演化曲线为中心的一个窄带内, 每次演化只需要计算窄带内的点到当前演化曲线的距离; 当曲线演化到窄带边界时, 再以当前曲线为中心建立新的窄带, 直到分割结束; 该算法的复杂度仅为$O(KN_1 ^2)$(其中$k$为窄带的宽度, $N_1 $为窄带内的图像网格点数).本文在这种方法的基础上, 提出一种基于$3\times 3$模板的快速水平集窄带算法, 通过分析计算可得, 该算法的复杂度仅为$O(KN_1 )$, 从而大大提高了分割效率.

给定初始轮廓曲线之后, 在上述图像轮廓的感兴趣区域以初始轮廓曲线上的点作为$3\times 3$模板中心, 与模板中其他点相对应的图像中点标记为窄带点, 让模板遍历整个初始轮廓曲线以形成窄带(图 1).这种模板的选择能够全面体现水平集的搜索方向.因为窄带内任何非水平集上一点与其零水平集上对应节点的连线方向通常在零水平集轮廓曲线的法线上(图 1), 而该方向可以用以零水平集节点为中心, 并与其3$\times $3模板内的8个节点和节点之间连线方向予以确定; 通过确定第一步的搜索方向, 可以很方便地将下一步的零水平集节点用模板内的8个节点来加以控制(图 3), 而不是窄带内的所有像素点, 从而大大提高算法的速度. 图 1中, $Q(x, y)$为水平集函数$\varphi =c$上点$P(x, y)$到零水平集$\varphi =0$的最短距离节点; 这里用3$\times $3模板, 其中0为中心点, $V_1 -V_8 $为其相邻的8个节点. 图 3中的$V_1 $和$V_2 $为模板节点, 0为零水平集上的一个节点; 假设$V_1 $和$V_2 $的水平集函数值分别为$\varphi _1 $和$\varphi _2 $; $\alpha _1 $和$\alpha _2 $为参数, 且

$\alpha _1 =\frac{\varphi _2 }{\varphi _2 -\varphi _1 }, {\kern 1pt}\;\alpha _2 =1-\alpha _1, $

从而可得$\alpha _1 \varphi _1 +\alpha _2 \varphi _2 =0$.该算法的另一个优点是无需设定窄带重置次数, 因为每次对邻域节点的近似求解过程已包含有重置零水平集过程.

上述的水平集演化方程可以使用有限差分法^[7]来实现数值求解, 即将所有空间域上的的偏导数用中心差分来实现, 所有时域上的偏导数用向前差分来实现, 具体表示如下: $\dfrac{\partial \varphi }{\partial t}=\dfrac{\varphi _{i, j}^{n+1} -\varphi _{i, j}^n }{\Delta t}, $其中$\Delta t$是时间步长.演化方程中的曲率$\mbox{div}(\dfrac{\nabla \varphi }{\left| {\nabla \varphi } \right|})$可以使用二阶中心差分来进行求解

$\mbox{div}(\frac{\nabla \varphi }{\left| {\nabla \varphi } \right|})=\frac{\varphi _{xx} \varphi _y^2-2\varphi _{xy} \varphi _x \varphi _y +\varphi _{yy} \varphi _x^2 }{(\varphi _x^2 +\varphi _y^2 )^{\frac{3}{2}}},$

其中$\varphi _x \, , \varphi _{y}, \varphi _{xx}, \varphi _{xy}, \varphi _{yy} $可用如下所给公式进行计算

$\begin{eqnarray*}&&\varphi _x =\dfrac{1}{2h}(\varphi _{i+1, j} -\varphi _{i-1, j} ), \varphi _y =\dfrac{1}{2h}(\varphi _{i, j+1} -\varphi _{i, j-1} ), \varphi _{xx} =\dfrac{1}{h^2}(\varphi _{i+1, j} +\varphi _{i-1, j} -2\varphi _{i, j} ), \\ &&\varphi _{yy} =\dfrac{1}{h^2}(\varphi _{i, j+1} +\varphi _{i, j-1} -2\varphi _{i, j} ), \varphi _{xy} =\dfrac{1}{h^2}(\varphi _{i+1, j+1} -\varphi _{i-1, j+1} -\varphi _{i+1, j-1} +\varphi _{i-1, j-1}, \end{eqnarray*}$

其中$h$是离散网格的间隔, 演化方程中的$c_1, c_2 $可按上式(2.1) 求解.水平集演化终止准则:当曲线长度变化的绝对值在给定的迭代次数(n)里一直小于某个事先给定的阈值$\eta $, 水平集演化终止^[8].在本文的实验里, 固定$n$=8, $\eta $=4.

该算法流程如下:

1) 输入原始图像$u_0 $.

2) 对$u_0 $作二维离散小波变换提取其可能的边缘点, 将集中的可能边缘点进行链接, 作为初始轮廓曲线C, 并设置如下相应的参数值:时间步长$\Delta t$, 网格间隔$h$, 以及$\alpha, \beta, \mu $的值.

3) 利用$3\times 3$模板在初始轮廓曲线附近设置窄带并按(2) 式演化水平集函数$\varphi $.

4) 从水平集函数$\varphi $中提取零水平集即演化曲线.

5) 判断水平集演化是否终止, 如果是, 算法结束, 否则转第三步.

下面的实验均使用Matlab7进行编程仿真, 实验所用的计算机配置为: Intel Core4处理器, 10G内存, Windows7操作系统.模型的参数设置为: $\alpha =\beta =1$, $\mu =0.01\ast 255^2$, 时间步长$\Delta t=0.1$, 网格间隔$h=1$, 正则化参数$\varepsilon =1$.实验1为C-V模型和改进的C-V模型以及本文方法(改进的C-V模型和窄带法)分别对一副随机生成的图像分割实验对比. 图 4为原始图像及初始轮廓, 图 5为C-V模型分割结果, 图 6为改进后C-V模型分割结果, 图 7为本文方法分割结果; 其中C-V模型需要迭代95次, 用时7.634s完成分割; 改进的C-V模型需要迭代40次, 用时3.022s;而本文方法只需迭代32次, 用时1.037s即可收敛.实验2是对一幅SAR卫星云图像分别用C-V模型、改进的C-V模型以及本文方法进行分割. 图 8位原始图形和初始轮廓, 结果C-V模型迭代了98次, 用时8.524s收敛(图 9); 改进的C-V模型需要迭代60次, 用时5.022s (图 10), 而本文方法只迭代了51次, 用时3.426s就得到了分割结果(图 11).实验结果表明, 改进的C-V模型比传统的C-V模型更能有效的捕捉图像的边缘信息, 从而分割效率也更高; 而窄带法能降低算法的复杂度, 大大提高图像的分割速度.

本文依据传统的C-V模型在图像分割时只考虑图像的全局信息的基础上, 利用小波多分辨率分析提出了一种加入图像梯度边缘信息的改进模型, 并且利用小波变换检测出的图像边缘点知道初始分割曲线的设置, 从而大大减少了迭代次数, 在数值实现时提出一种基于$3\times 3$模板的快速窄带法, 通过数值实验验证, 该方法比传统的C-V模型有稳定性好、捕捉图像边缘能力强、分割速度快、算法复杂度低等优点.