计算机系统应用  2019, Vol. 28 Issue (11): 10-18   PDF    
LATE水平集图像分割模型的矩形窄带法
曾笑云1, 杨晟院1, 潘园园1, 刘洋1, 左国才2     
1. 湘潭大学 信息工程学院, 湘潭 411105;
2. 湖南软件职业学院 软件与信息工程学院, 湘潭 411100
摘要:窄带法是水平集图像分割的一种常见的加速方法. 传统窄带仍然存在冗余的计算区域; 传统窄带法与LATE (Local Approximation of Taylor Expansion)水平集模型结合时, 图像分割效率反而可能下降. 针对这些问题, 本文提出了一种基于LATE水平集图像分割模型的矩形窄带法. 在每次LATE水平集迭代之前, 对水平集做如下窄带处理. 首先找出水平集的所有过零点; 然后对过零点做活动约束, 剔除不活动的过零点, 有效缩小窄带范围; 再对活动约束的过零点生成矩形窄带; 对重叠的矩形窄带进行合并优化, 使得矩形窄带总面积尽可能小. 最后, 在矩形窄带范围内求解水平集微分方程, 更新水平集, 完成本次迭代. 在水平集演化的不同阶段, 对传统窄带法的窄带面积与本文矩形窄带面积进行了比较. 随着迭代次数增加, 矩形窄带面积与传统窄带法的窄带面积之比逐渐减小到0, 说明矩形窄带法有效地减少了冗余计算量. 针对不同程度的灰度不均匀图像, 本文方法与LATE方法、结合LATE模型的直接窄带法、以及结合LATE模型的DTM窄带法进行了比较. 直接窄带法和DTM窄带法的分割速度反而慢于LATE方法. 对灰度严重不均匀的图像, 直接窄带法和DTM窄带法的分割质量受到了较大影响. 本文方法在保持较好分割效果的条件下, 分割速度快于LATE方法. 本文的矩形窄带方法有效地降低了算法复杂度, 提高了图像分割效率.
关键词: 活动约束    矩形窄带    LATE水平集模型    灰度不均匀    图像分割    
Rectangular Narrow-Band Method for Image Segmentation of LATE Level Set Model
ZENG Xiao-Yun1, YANG Sheng-Yuan1, PAN Yuan-Yuan1, LIU Yang1, ZUO Guo-Cai2     
1. College of Information Engineering, Xiangtan University, Xiangtan 411105, China;
2. School of Software and Information Engineering, Hunan Software Vocational College, Xiangtan 411100, China
Foundation item: National Natural Science Foundation of China (11571293)
Abstract: The narrow-band method is a common acceleration method for level set image segmentation. The traditional narrow-band still has redundant computational regions; When the traditional narrow-band method is combined with the LATE (Local Approximation of Taylor Expansion) level set model, the image segmentation efficiency may be reduced. In order to solve these problems, a rectangular narrow-band method based on LATE level set image segmentation model is proposed in this study. The level set is subjected to the following narrow-band processing before each LATE level set iteration. First, find out all the points of zero crossings of the level set; second constrict the points of zero crossings by the activity constraints, eliminate the inactive points of zero crossings, and effectively reduce the area of the narrow-band, then generate a rectangular narrow-band for the points of zero crossings by the active constraints, optimize the overlapping rectangular narrow-band so that the total area of the rectangular narrow-band is as small as possible. Finally, the level set differential equation is solved in the narrow-band of the rectangle, and the level set is updated to complete this iteration. In the different stages of the level set evolution, the area of the traditional narrow-band and the rectangular narrow-band of this study are compared. As the number of iterations increases, the ratio of the area of rectangular narrow-band to the area of traditional narrow-band is gradually reduced to zero, indicating that the rectangular narrow-band method effectively reduces the amount of redundancy calculation. For images with different degrees of intensity inhomogeneity, the proposed method is compared with the LATE method, the direct narrow-band method, and the DTM narrow-band method. The direct narrow-band method and the DTM narrow-band method have lower segmentation efficiency than the LATE method, and the segmentation quality is greatly affected for some images with severe intensity inhomogeneity. Under the condition of maintaining good segmentation effect, the segmentation speed of the proposed method is faster than that of LATE method. The rectangular narrow-band method in this study effectively reduces the complexity of the algorithm and improves the efficiency of image segmentation.
Key words: active constraint     rectangular narrow-band method     LATE level set method     intensity inhomogeneity     image segmentation    

水平集最先由Osher[1]提出, 是一种将低维问题嵌入高维问题求解的方法. 水平集广泛应用于图像分割[2,3], 它可获得亚像素精度的封闭轮廓和区域. 以2D灰度图像二相水平集为例, 水平集取值的正负将图像域划分为目标和背景两个区域. 水平集的零等高线被视为目标和背景的分界线, 即二维的活动轮廓曲线. 活动轮廓随着水平集的演化不断逼近真实轮廓, 最终完成分割.

早期的水平集模型, 如MS模型[4], CV模型[5], SLGS模型[6], 对灰度不均匀的图像分割效果不理想. 为了有效分割灰度不均匀图像, 水平集模型朝着越来越复杂的方向发展, 如RSF模型[7], LIC模型[8], LSACM模型[9], LATE模型[10]. 其中, LATE模型利用泰勒展开式对拟合函数进行非线性逼近, 极大的提高了分割灰度严重不均匀图像的能力.

水平集方法虽然能获得较好分割效果, 但也提高了计算复杂度. 为了提高运算效率, Chopp[11]提出窄带法, 并由Adalsteinsson等[12]给出了详细的实现方法. 窄带法的核心思想就是把计算区域约束到活动轮廓附近的带状区域, 避免了对整个图像域进行计算, 以此来提高水平集方法的计算效率.

窄带方法能够提高水平集的效率, 首先要求窄带生成过程要尽可能快, 生成窄带新增的计算代价应小于由于窄带减少计算面积而节约的计算代价. 为了提高窄带生成速度, 目前的方法主要有快速进行法[13,14], 快速扫描法[15,16], DTM窄带法[17]等.

文献[13]认为可设置较宽的窄带. 当活动轮廓未达到窄带边缘时, 无需更新窄带, 当活动轮廓达到窄带边缘, 但不发生波动时, 说明活动轮廓已经达到图像边缘, 也无需更新窄带. 只有当窄带边界点上有活动轮廓波动变化时, 才需更新窄带.

文献[18,19]抛弃了窄带更新策略, 在每次水平集迭代中均更新窄带. 这样可规避繁琐的窄带更新条件判断. 在不损害分割质量的前提下, 设置尽可能窄的窄带, 能更好地发挥窄带缩小计算区域的作用.

文献[20]提出窄带压缩数据结构. 为了保留邻域信息, 分别做行方向和列方向两次压缩. 压缩的窄带结构虽然可以规避对非窄带区域的范围判断, 但使得邻域结构变得复杂.

虽然传统窄带法有效地减少了计算范围, 提高了计算效率, 但传统窄仍然存在冗余的计算区域. 在水平集演化的过程中, 一部分活动轮廓先到达图像边缘不再运动, 另一部分活动轮廓还需继续演化逼近图像边缘. 先到达图像边缘的活动轮廓在后续演化过程中的计算属于冗余计算, 还未达到图像边缘的活动轮廓才真正需要进行演化计算. 为此, 我们提出活动约束策略, 将窄带的范围进一步约束在未达到图像边缘的活动轮廓的区域.

约束的活动轮廓区域形状不规则, 可采用最小矩形覆盖不规则的窄带区域, 从而构造矩形窄带. 为了保证矩形窄带的总面积尽可能小, 对窄带区域进行了合并优化. 相比传统不规则窄带, 矩形窄带结构更简单, 更有利于演化计算.

为此, 本文提出了一种基于LATE水平集图像分割模型的矩形窄带法. 本文剩余部分结构安排如下: 第1节介绍水平集和传统窄带等相关知识; 第2节介绍本文的矩形窄带法; 第3节为数值实验; 第4节为结论.

1 相关知识 1.1 水平集

$\Omega \in {{{R}}^{\rm{2}}}$ 为图像域, $I:\Omega \to {{R}}$ 为给定的灰度图像. $\phi :\Omega \to {{R}}$ 为水平集函数[21]. 活动轮廓曲线 ${{C}}$ $\phi $ 的零等高线, 如图1所示.

图 1 水平集的原理示意图

水平集方法可概括为: 首先确立合理的水平集微分方程 $\partial \phi /\partial t$ 并给出恰当的初始水平集. 然后利用水平集微分方程不断更新水平集, 驱动活动轮廓 ${{C}}$ 向图像的真实轮廓运动, 直到完成分割.

以LATE水平集模型[10]为例, 它的水平集微分方程为:

$\begin{split} \dfrac{{\partial \phi }}{{\partial t}}\! =&\! - \delta \left(\! \phi \! \right)\! \cdot \!\displaystyle\sum\limits_{i = 1}^2\! {\bigg(\int_\Omega\!\! {{k_\sigma }\left( {x\! - \!y} \right)} I\left( x \right)} - LI{M_i}\left( x \right) - b'\left( y \right){C_i}{^2}dy\bigg) \\ &+ \mu \delta \left( \phi \right) \cdot div\left(\dfrac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}\right) + \nu \left({\nabla ^2}\phi - div\left(\dfrac{{\nabla \phi }}{{\left| {\nabla \phi } \right|}}\right)\right) \end{split}\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! $ (1)

其中, $\mu $ , $\nu $ 均为实数常量, $\delta \left( \cdot \right)$ 为单位冲击函数, ${k_\sigma }$ 为高斯内核函数. $LI{M_i}$ , ${C_i}$ , $b$ 的表达式分别为:

$ LI{M_i}\left( x \right) = mean\left( {I(y):\;y \in {\Omega _x} \cap {R_i}} \right) $ (2)

其中, $y$ 为以 $x$ 为中心的邻域, ${R_i}\left( {i \in {\rm{\{ }}1,2{\rm{\} }}} \right)$ 分别表示目标区域和背景区域.

${C_i} = \dfrac{{\iint_\Omega {{k_\sigma }(x - y)\left( {I(x) - LI{M_i}(x)} \right){M_i}(\phi )dydx}}}{{\iint_\Omega {{k_\sigma }(x - y)b'(y){M_i}(\phi )dydx}}}$ (3)
$ b'(y) = \dfrac{{\int_\Omega {{k_\sigma }(x - y)} \displaystyle\sum\limits_{i = 1}^2 {\left( {\left( {I(x) - LI{M_i}(x)} \right){C_i}{M_i}(\phi )} \right)} dx}}{{\int_\Omega {{k_\sigma }(x - y)} \displaystyle\sum\limits_{i = 1}^2 {\left( {{C_i}^2{M_i}^2(\phi )} \right)} dx}} $ (4)

其中, ${M_1}(\phi ) = H(\phi )$ , ${M_2}(\phi ) = 1 - H(\phi )$ , $H( \cdot )$ 为单位阶跃函数.

水平集方法中, 称每次迭代的计算复杂度为水平集复杂度, 记为 ${C_{ls}}$ . LATE模型中, 复杂度最高的运算是卷积运算. 假设图像宽度和长度均为N, 卷积模板宽度为D, 其中, D远远小于N. 则LATE的水平集复杂度为: $C_{ls}^{late} = {\rm{O}}({N^2}{D^2})$ .

由于光照磁场或成像装置缺陷等因素[7], 可能导致获得的图像是灰度不均匀的. 灰度不均匀给图像分割带来很大的挑战. 在文献[8]中, 灰度不均匀图像被看作灰度均匀图像与偏置光源场的混合. 蝴蝶图像与严重不均匀偏置场的混合如图2所示.

图 2 灰度严重不均匀图像的合成

不同水平集模型对灰度严重不均匀蝴蝶图像的分割结果如图3所示. 可见, LATE模型具有很强的分割灰度不均匀图像的能力. 但是, LATE模型的计算代价较高, 分割缓慢.

图 3 不同水平集模型针对灰度严重不均匀图像的分割结果. (a) CV模型; (b) SLGS模型; (c) RSF模型; (d) LIC模型; (e) LSACM模型; (f) LATE模型

1.2 传统窄带及其复杂度分析

虽然水平集能够获得较好的分割结果, 但它将问题提高了一个维度, 增加了计算复杂度. 直观地看, 活动轮廓 $C$ 受过零点附近的水平集 $\phi $ 值的变化影响较大, 受远离过零点的水平集 $\phi $ 值影响较小. 传统窄带法通过把整个图像域的水平集的计算区域限制在过零点附近的窄带上来减小计算量, 如图4所示.

图 4 窄带的原理示意图

假设传统窄带的宽度与卷积模板宽度具有相同规模, 传统窄带长度规模与图像宽度相当. 则传统窄带的面积规模可以记为: ${S_{nb}} = O(ND)$ . 本文中传统窄带法泛指窄带面积规模为 ${S_{nb}} = O(ND)$ 的窄带法, 如直接生成窄带法, 快速进行法[13,14], 快速扫描法[15,16]和DTM[17]窄带法等.

窄带水平集方法每次迭代的复杂度称为总复杂度, 记为 ${C_{all}}$ . 总复杂度可以分成两部分: 生成窄带的计算复杂度称为窄带复杂度, 记为 ${C_{nb}}$ ; 对应窄带区域的水平集复杂度 ${S_{nb}}/{N^2} \cdot {C_{ls}}$ . 总复杂度可表示为:

${C_{all}} = {C_{nb}} + \dfrac{{{S_{nb}}}}{{{N^2}}}{C_{ls}}$ (5)

一般地, 生成窄带的复杂度 ${C_{nb}}$ 不应高于水平集复杂度 ${C_{ls}}$ , 否则总复杂度反而上升. 因此, 窄带水平集方法的总复杂度主要取决于窄带面积的规模 ${S_{nb}}$ .

${S_{nb}}$ 固定时, 生成窄带复杂度 ${C_{nb}}$ 越小越好. 直接生成窄带法对每个像素点进行扫描, 判断其与过零点的距离是否小于半个窄带宽度, 如果是则将其加入窄带. 因此, 直接生成窄带法的生成窄带复杂度为: $C_{nb}^{drct} = O({N^2}{D^2})$ . DTM窄带法则遍历每个过零点, 在半个窄带宽度范围内标记窄带. 因此, DTM的生成窄带复杂度为: $C_{nb}^{dtm} = O(N{D^2})$ .

快速进行法, 快速扫描法以及DTM方法等窄带法都只是降低了生成窄带复杂度 ${C_{nb}}$ . 由于窄带水平集方法的总复杂度主要取决于窄带面积的规模 ${S_{nb}}$ . 因此, 传统窄带法对提高水平集分割图像的效率有限.

2 本文的矩形窄带法 2.1 寻找过零点

$i,j$ 分别表示2D图像域上 $x$ 轴和 $y$ 轴上的坐标值, $I(i,j)$ 为图像, $\phi (i,j)$ 为水平集, $C$ 为活动轮廓, 则 $C$ $\phi (i,j)$ 的零等高线.

$CRS(i,j)$ 为二值矩阵, $CRS(i,j) = 1$ 表示像素点 $(i,j)$ 为过零点, $CRS(i,j) = 0$ 表示像素点 $(i,j)$ 为非过零点. $VIS(i,j)$ 为访问矩阵, 用来标记某像素点 $(i,j)$ 被标记为过零点的次数. 初始时 $CRS$ $VIS$ 均为零矩阵, 每次更新水平集 $CRS$ 重新置零. 可在 $x$ 轴和 $y$ 轴正负4个方向寻找过零点. 对每个像素点 $(i,j)$ 做以下操作:

$\begin{array}{l} {\rm{if}}\;\;\;\;\;\;\;\;(\phi (i - 1,j) > 0\;\;{\rm{and}}\;\;(\phi (i + 1,j) < 0) \\ \;\;\;\;\;\;\;\;\;\;{\rm{or}}\;(\phi (i - 1,j) < 0\;\;{\rm{and}}\;\;(\phi (i + 1,j > 0) \\ \;\;\;\;\;\;\;\;\;\;{\rm{or}}\;(\phi (i,j - 1) > 0\;\;{\rm{and}}\;\;(\phi (i,j + 1 < 0) \\ \end{array} $
$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;{\rm{or}}\;(\phi (i,j - 1) > 0\;\;{\rm{and}}\;\;(\phi (i,j + 1 < 0) \\ {\rm{then}}\;\;\;\;CRS(i,j) = 1,\;\;VIS(i,j) = VIS(i,j) + 1 \end{array}$

则可找到当前所有的过零点 $CRS(i,j) = 1$ , 以及统计当前像素点被访问的次数 $VIS$ .

在水平集演化的过程中, 一部分活动轮廓先到达图像边缘不再运动, 另一部分活动轮廓还需继续演化逼近图像边缘. 先到达图像边缘的活动轮廓在后续演化过程中的计算属于冗余计算, 还未达到图像边缘的活动轮廓才真正需要进行演化计算. 如果在静止的过零点附近生成窄带, 对活动轮廓的演化并没有作用. 可以设置延时参数 $delay$ , 若某像素点被访问的次数超过 $delay$ 次, 则认为它是静止不动的, 将不再运动的过零点从 $CRS(i,j)$ 中剔除. 称这样的处理方法为活动约束. 对每个像素点 $(i,j)$ 做活动约束, 具体操作如下:

$\begin{array}{l} {\rm{if}}\;\;\;\;\;\;VIS(i,j) >{{ delay}} \\ {\rm{then}}\;\;\;CRS(i,j) = 0 \\ \end{array} $

在活动约束之前, 过零点长度的规模为N. 假设 $delay$ 的规模与窄带宽相同, 则活动约束处理之后, 可认为过零点的长度规模为D.

2.2 生成矩形窄带

传统窄带的生成, 是在过零点基础上通过偏移生成带状区域的过程. 传统方法获得的过零点曲线是封闭的, 而经过活动约束, 本文方法获得的过零点集合只是传统过零点曲线的一部分, 是多段开放的链, 称为过零点链. 原则上, 应该在每条过零点链的邻域内生成不规则的窄带. 由于活动约束下每条过零点链长度较短, 窄带区域接近矩形区域, 不妨用规则的矩形区域代替不规则的窄带区域. 传统的窄带区域是不规则的, 需要单独建立窄带中每个点与原图像域之间的映射, 才能将窄带的计算结果返回到图像域. 而矩形区域则只需给出矩形的两个对角端点便可确定一个区域, 且矩形域与图像域的像素点位置对应非常简单, 这更有利于快速卷积运算以及边界处理.

可利用扫描法获得每条过零点链对应的矩形窄带区域, 步骤如下:

1)分别对 $VIS$ 中每个值为1的像素点作为扫描的起点, 逐一进行扫描. 将该点坐标插入链码 $q$ , 将该像素点的 $VIS$ 值置零.

2)若在8邻域内扫描到 $VIS$ 为1的点, 将该点作为新的扫描起点, 将其坐标点插入链码 $q$ , 对应 $VIS$ 值置零.

3)重复步骤2), 直到在8邻域内找不到下一个 $VIS$ 为1的点.

4)确定矩形域. 取链码 $q$ 中最小的横坐标和纵坐标作为矩形框的第一个端点, 取最大的横坐标和纵坐标作为矩形框的第二个端点. 得到最小的覆盖矩形域. 将矩形区域向四周扩宽大D/2个宽度, 防止过零点暴露在窄带最外层.

活动约束后, 剩余过零点的长度规模为 ${\rm{O}}(D)$ . 生成矩形窄带的过程既是对活动约束后剩余过零点的扫描过程, 该过程不再对半个窄带宽度范围内的点进行判断. 因此, 本文的方法生成窄带的计算复杂度仅为 $C_{nb}^{rec} = {\rm{O}}(D)$ . 矩形窄带的长宽规模均为 ${\rm{O}}(D)$ , 因此, 窄带的面积规模仅为 $S_{nb}^{rec} = {\rm{O}}({D^2})$ .

寻找过零点方法具有简单高效的优点, 但可能存在个别过零点的遗落, 导致过零点曲线不连续, 从而产生面积重叠的矩形区域. 如图5中的左图所示, 矩形A与矩形B产生重叠. 记矩形A的面积为 ${S_A}$ , 矩形B的面积为 ${S_B}$ . 对于重叠的矩形A和矩形B, 将右上和左下的区域补充形成一个大的矩形C, 如图5所示, 记矩形C的面积为 ${S_C}$ . 如果s, 则矩形A和矩形B不需合并成矩形C. 反之, 则采用矩形C替换矩形A和矩形B. 这样处理可以减少计算区域, 从而提高计算效率.

图 5 矩形区域合并优化

具体操作如下: 设 $\left( {s{x_1},s{y_1}} \right)$ , $\left( {e{x_1},e{y_1}} \right)$ 分别表示矩形框A的两个端点, $\left( {s{x_2},s{y_2}} \right)$ , $\left( {e{x_2},e{y_2}} \right)$ 分别表示矩形框B的两个端点, 如图5所示.

则判断它们相交的条件为:

$\begin{array}{l} s{x_1} \le e{x_2}\;{\rm and}\;e{x_1} \le s{x_2} \\ {\rm and}\;s{y_1} \le e{y_2}\;{\rm and}\;e{y_1} \le s{y_2} \\ \end{array} $

${x_1} = \min \left( {s{x_1},e{x_1}} \right)$ , ${y_1} = \min \left( {s{y_1},e{y_1}} \right)$ , ${x_1} = \min $ $\left( {s{x_1},e{x_1}} \right)$ , ${y_2} = \min \left( {s{y_2},e{y_2}} \right)$ . 判断合并后面积更小的条件为:

$\begin{array}{l} ({x_2} - {x_1})*({y_2} - {y_1}) < (s{x_2} - s{x_1})*(s{y_2} - s{y_1}) \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\; + (e{x_2} - e{x_1})*(e{y_2} - e{y_1}) \\ \end{array} $

每个矩形其实可以看做一个已知两个端点的二维闭区间, 面积优化只需对闭区间的端点进行操作, 计算代价远远小于寻找的遗落过零点的代价.

在每次水平集迭代中, 我们生成窄带的算法总结如算法1所示.

算法1. 生成窄带算法

输入: 原矩阵图像

输出: 若干矩形窄带

1. 初始化 $\scriptstyle CRS$ $\scriptstyle VIS$ , 根据当前水平集 $\scriptstyle \varphi $ , 获得所有过零点. 标记过零点矩阵 $\scriptstyle CRS$ , 更新访问矩阵 $\scriptstyle VIS$ 的访问次数;

2. 对访问次数超过delay的像素点进行限制, 将过零点( $\scriptstyle CRS=1$ )约束在满足 $\scriptstyle VIS(i,j) \le delay$ 的区域内;

3. 利用扫描法生成每段过零点链的最小覆盖矩阵框;

4. 将最小矩形框两两之间进行判断, 如果它们的区域相交且合并后面积更小, 则对它们进行合并优化.

由于图像也是矩形数组, 矩形框的窄带结果和原图像能保持一致, 直接将每个矩形框代入原水平集函数计算即可. 返回结果时只需做坐标偏移处理, 便可实现窄带到原图像的映射, 从而在本次迭代中实现水平集的更新. 可见, 本文的方法与水平集方法非常容易实现对接.

2.3 与传统窄带法的比较

传统窄带示意图如图6(a)所示. 传统窄带存在很大部分的无效计算区域, 实际的活动窄带区域只占总窄带的一小部分. 以DTM方法为例, DTM方法生成窄带的复杂度为 $C_{nb}^{dtm} = {\rm{O}}(N{D^2})$ , 窄带面积规模为 $S_{nb}^{dtm} = {\rm{O}}(ND)$ . 根据式(5), DTM方法的总复杂度为:

$\begin{split} C_{all}^{dtm}& = C_{nb}^{dtm} + \dfrac{{S_{nb}^{dtm}}}{{{N^2}}}C_{ls}^{late} \\ & = {\rm{O}}(N{D^2}) + \dfrac{{{\rm{O}}(ND)}}{{{N^2}}}{\rm{O}}({N^2}{D^2}) \\ & = {\rm{O}}(N{D^3}) \end{split} $

本文的矩形窄带示意图如图6(b)所示. 通过活动约束, 对长期禁止不动的点进行了屏蔽, 大幅减少了窄带的生成范围. 扫描生成矩形窄带的复杂度为: $C_{nb}^{rec} = {\rm{O}}(D)$ , 窄带面积规模为 $S_{nb}^{rec} = {\rm{O}}({D^2})$ . 根据式(5), 矩形窄带法的总复杂度为:

$\begin{split} C_{all}^{rec} &= C_{nb}^{rec} + \dfrac{{S_{nb}^{rec}}}{{{N^2}}}C_{ls}^{late} \\ &= {\rm{O}}(D) + \dfrac{{{\rm{O}}({D^2})}}{{{N^2}}}{\rm{O}}({N^2}{D^2}) \\ &= {\rm{O}}({D^4}) \end{split} $

由于D远远小于N, 所以 $C_{all}^{rec} < C_{all}^{dtm}$ . 可见, 本文的矩形窄带方法不仅生成窄带复杂度低于DTM方法, 而且结合LATE水平集的总复杂度也低于DTM方法. 其中, 窄带面积规模 $S_{nb}^{rec}$ 的降低, 对进一步提高窄带效率起到关键作用.

图 6 矩形窄带和传统窄带示意图

3 数值实验 3.1 窄带面积比较

本文的实验结果均在MATLAB R2016a上实现, 操作系统为Win10. 图7(a)(b)分别为针对不同灰度不均匀情形的两组实验. 实验采用矩形窄带法结合LATE模型[10]生成水平集, 并利用水平集过零点分别生成传统窄带(TradNb)、仅添加活动约束的不规则窄带(AcNb)以及矩形窄带(RecNb).

图 7 窄带演化过程. (a)从左至右迭代次数分别为1、44、88、132、176、220、266; (b)从左至右迭代次数分别为1、25、50、75、100、125、150. 第1行: 图像的活动轮廓; 第2行: 传统窄带(TradNb); 第3行: 仅添加活动约束的不规则窄带(AcNb,本文提出的过渡方案); 第4行: 活动约束矩形窄带(RecNb,本文最终采用方案).

本文中传统窄带法(TradNb)泛指窄带面积规模为 ${S_{nb}} = {\rm{O}}(ND)$ 的一类窄带法, 如直接生成窄带法, 快速进行法, 快速扫描法和DTMP窄带法等. AcNb为本文提出的过渡方案, RecNb为本文最终采用方案. 根据前文分析, 不规则窄带矩形化的过程可能略微增加窄带区域. 因此, RecNb相比AcNb面积可能略微增加, 但AcNb的计算和实现更为简单. RecNb和AcNb的面积规模均为 ${S_{nb}} = {\rm{O}}({D^2})$ , 其中D<N.

图7(a)为灰度不均匀程度一般的图像, 窄带半径为5. 图7(b)为灰度严重不均匀的图像, 窄带半径为10. LATE利用泰勒展开式对灰度不均匀进行调节. 当灰度变化平缓时, 各像素点的水平集受其邻域的影响较小. 当灰度不均匀程度较严重时, 较远邻域的泰勒展开权值变大, LATE模型能够利用较远邻域信息对灰度不均匀进行修正. 因此, 在灰度严重不均匀区域, LATE水平集分割缓慢且更容易受到窄带的影响. 可知, 灰度不均匀程度越高, 所需要设置的最小窄带半径越大.

图7中第3行的窄带面积明显小于第2行, 表明活动约束能够很大程度地减少窄带范围. RecNb为AcNb的最小矩形区域, RecNb与AcNb的面积规模相当, 但矩形区域更方便计算机处理. 图7中第3行与第4行表明RecNb与AcNb的窄带位置和面积相差不大.

图7中TradNb、AcNb和RecNb 对应的窄带面积如表1表2所示. 可见, AcNb的窄带不大于TradNb的窄带面积. 理论上, RecNb的面积略大于AcNb的面积. 但在实际扫描生成矩形窄带的过程中, 一些不必要的过零点被抛弃, RecNb的面积也可能略小于AcNb. 从表1表2的数据来看, RecNb与AcNb的窄带面积相差不大. 图8表明随着水平集的演化, 矩形窄带的面积与传统窄带面积之比逐渐减少到0. 可见, 本文的矩形窄带法能有效地减少窄带面积, 从而提高计算效率.

表 1 图7(a)对应窄带面积

表 2 图7(b)对应的窄带面积

图 8 图7中的矩形窄带面积与传统窄带面积的比值

3.2 运行效率比较

为了表述方便, 本节将直接窄带与LATE模型结合的窄带水平集记为DRCTLS; DTM窄带与LATE模型结合的窄带水平集记为DTMLS; 矩形窄带与LATE模型结合的窄带水平集记为RECLS. LATE方法, DRCTLS方法, DTMLS方法以及RECLS方法对图像的分割结果对比如图9所示. 在对比实验中, 保证各方法的共有的参数完全相同.

当灰度不均匀程度一般时(图9(a)–图9(e)), 各方法都能得到较好的分割结果. 当灰度严重不均匀时, 灰度不均匀区域分割缓慢, 且容易受到窄带范围的影响. 在相同的窄带宽度下, DRCTLS方法和DTMLS方法的分割结果可能受到损坏, 如图9(f)图9(h)图9(j)所示. 而RECLS方法对图9中不同程度灰度不均匀图像均能保持稳定的分割结果, 且分割效率高于LATE水平集方法以及其它窄带LATE方法.

图9中LATE、DRCTLS、DTMLS和RECLS方法对应的迭代次数, 运行时间, 以及平均每次迭代所需时间如表3所示. DRCTLS和DTMLS具有相同的窄带规模, 且DRCTLS的生成窄带复杂度高于DTMLS, 除图9(B)的极端情形外, DTMLS的分割效率总体上高于DRCTLS.

由于DRCTLS和DTMLS减少计算区域的收益不足以弥补生成窄带增加的额外计算开支, 反而可能导致总计算效率的下降. RECLS矩形窄带生成窄带复杂度和窄带面积规模均小于DTMLS方法, 特别是窄带面积规模的下降, 使得窄带计算效率明显提升. 可见, RECLS方法的计算效率明显优于DRCTLS窄带法, DTMLS窄带法以及未使用窄带的原始LATE方法.

3.3 分割精度分析

图像分割的准确性可以用Jaccard Similarity Coefficient (JSC)[9,10]标准来衡量.

$JSC({O_m},{O_t}) = \frac{{A({O_m} \cap {O_t})}}{{A({O_m} \cup {O_t})}}$ (6)
图 9 分割结果对比. 第1行: 原图像; 第2行: LATE模型分割结果; 第3行: LATE结合DRCT窄带法(DRCTLS)的分割结果; 第4行: LATE结合DTM窄带法(DTMLS)的分割结果; 第5行: 我们方法(RECLS)的分割结果. (a)为灰度均匀图像; (b)–(e)为灰度不均匀图像; (f)–(j)为灰度严重不均匀图像.

表 3 图9中对应的迭代次数以及运行时间

其中, ${O_t}$ 为标准分割区域, ${O_m}$ 为实际分割区域, 算符 $A( \cdot )$ 表示求对应区域的面积. JSC的取值范围在0到1, JSC的值越大, 分割结果越准确.

图10(a)为二值图像, 图10(b)图10(e)灰度不均匀程度依次递增. 以二值图10(a)的黑色区域作为标准分割区域 ${O_t}$ , CV模型, SLGS模型, RSF模型, LIC模型, LSACM模型, LATE模型以及本文的RECLS方法对图10(a)图10(e)的分割结果对应的JSC值如表4所示. 实验中, 图片尺寸均为100×100. 图10(a)图10(c) RECLS方法的窄带半径为5, 图10(d)图10(e) RECLS方法的窄带半径为10.

图 10 灰度严重不均匀图像的合成. (a)灰度均匀的二值图像; (b)(c)灰度不均匀程度一般的合成图像; (d)(e)灰度严重不均匀的合成图像.

可见, LATE方法对灰度严重不均匀图像具有较高的分割精度. RECLS方法对应的JSC值与LATE几乎一致. 结合上一节的结论, 本文提出的矩形窄带方法能在不影响LATE模型分割精度的条件下, 提高对灰度严重不均匀图像的分割效率.

4 结论与展望

本文提出一种新的矩形窄带方法. 通过活动约束进一步缩小了窄带的范围. 利用矩形窄带代替不规则窄带, 使其更容易与水平集方法相结合, 减少了更新水平集的计算量. 实验表明, 本文的方法即使在灰度严重不均匀情形下也能够保持稳定的分割结果.

本文的方法在窄带演化的过程中, 可能存在多个矩形窄带, 而这些窄带没有实现并行运算. 如何让多个矩形窄带区域实现并行运算, 进一步提高计算的效率, 是我们下一步研究的内容.

表 4 7种方法对图10(a)图10(e)的分割结果对应的JSC值

参考文献
[1]
Osher S, Sethian JA. Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations. Journal of Computational Physics, 1988, 79(1): 12-49. DOI:10.1016/0021-9991(88)90002-2
[2]
江少锋, 杨素华, 陈震, 等. 水平集中符号距离函数并行降维计算. 中国图象图形学报, 2018, 23(2): 174-181.
[3]
周力, 闵海. 基于局部连接度和差异度算子的水平集纹理图像分割. 中国图象图形学报, 2019, 24(1): 39-49.
[4]
Mumford D, Shah J. Optimal approximations by piecewise smooth functions and associated variational problems. Communications on Pure and Applied Mathematics, 1989, 42(5): 577-685. DOI:10.1002/cpa.3160420503
[5]
Chan TF, Vese LA. Active contours without edges. IEEE Transactions on Image Processing, 2001, 10(2): 266-277. DOI:10.1109/83.902291
[6]
Zhang KH, Zhang L, Song HH, et al. Active contours with selective local or global segmentation: A new formulation and level set method. Image and Vision Computing, 2010, 28(4): 668-676. DOI:10.1016/j.imavis.2009.10.009
[7]
Li CM, Kao CY, Gore JC, et al. 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
[8]
Li CM, Huang R, Ding ZH, et al. 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
[9]
Zhang KH, Zhang L, Lam KM, et al. A level set approach to image segmentation with intensity inhomogeneity. IEEE Transactions on Cybernetics, 2016, 46(2): 546-557. DOI:10.1109/TCYB.2015.2409119
[10]
Min H, Jia W, Zhao Y, et al. LATE: A level-set method based on local approximation of Taylor expansion for segmenting intensity inhomogeneous images. IEEE Transactions on Image Processing, 2018, 27(10): 5016-5031. DOI:10.1109/TIP.2018.2848471
[11]
Chopp DL. Computing minimal surfaces via level set curvature flow. Journal of Computational Physics, 1993, 106(1): 77-91. DOI:10.1006/jcph.1993.1092
[12]
Adalsteinsson D, Sethian JA. A fast level set method for propagating interfaces. Journal of Computational Physics, 1995, 118(2): 269-277. DOI:10.1006/jcph.1995.1098
[13]
Adalsteinsson D, Sethian JA. The fast construction of extension velocities in level set methods. Journal of Computational Physics, 1999, 148(1): 2-22. DOI:10.1006/jcph.1998.6090
[14]
Herrmann M. A domain decomposition parallelization of the fast marching method. Bonn, Germany: German Research Foundation (DFG), 2003.
[15]
Tsai YHR. Rapid and accurate computation of the distance function using grids. Journal of Computational Physics, 2002, 178(1): 175-195. DOI:10.1006/jcph.2002.7028
[16]
Zhao HK. A fast sweeping method for eikonal equations. Mathematics of Computation, 2005, 74(250): 603-627.
[17]
周则明, 项杰, 王洪元, 等. 水平集方法中窄带构造技术. 系统工程与电子技术, 2007, 29(7): 1201-1204. DOI:10.3321/j.issn:1001-506X.2007.07.042
[18]
Li CM, Xu CY, Konwar K M, et al. Fast distance preserving level set evolution for medical image segmentation. Proceedings of the 9th International Conference on Control, Automation, Robotics and Vision. Singapore. 2006. 1–7.
[19]
Li CM, Xu CY, Gui CF, et al. 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
[20]
柳周, 李宏伟. 窄带水平集方法. 计算机工程与设计, 2009, 30(14): 3348-3351.
[21]
Mitiche A, Ayed IB. Variational and Level Set Methods in Image Segmentation. Berlin, Heidelberg: Springer, 2010.