LATE水平集图像分割模型的矩形窄带法①

曾笑云，杨晟院，潘园园，刘 洋，左国才

1(湘潭大学 信息工程学院，湘潭 411105)

2(湖南软件职业学院 软件与信息工程学院，湘潭 411100)

水平集最先由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 水平集

设 Ω ∈R2为图像域，I:Ω →R为给定的灰度图像.φ:Ω →R为 水平集函数[21].活动轮廓曲线C为φ 的零等高线，如图1所示.

图1 水平集的原理示意图

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

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

其中，µ，ν 均为实数常量，δ (·) 为 单位冲击函数，kσ为高斯内核函数.LIMi，Ci，b的表达式分别为：

其中，y为以x为 中心的邻域，Ri(i∈{1,2})分别表示目标区域和背景区域.

其中，M1(φ)=H(φ)，M2(φ)=1-H(φ)，H(·)为单位阶跃函数.

水平集方法中，称每次迭代的计算复杂度为水平集复杂度，记为Cls.LATE 模型中，复杂度最高的运算是卷积运算.假设图像宽度和长度均为N，卷积模板宽度为D，其中，D远远小于N.则LATE 的水平集复杂度为：

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

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

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

图3 不同水平集模型针对灰度严重不均匀图像的分割结果.

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

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

图4 窄带的原理示意图

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

窄带水平集方法每次迭代的复杂度称为总复杂度，记为Call.总复杂度可以分成两部分：生成窄带的计算复杂度称为窄带复杂度，记为Cnb；对应窄带区域的水平集复杂度Snb/N2·Cls.总复杂度可表示为：

一般地，生成窄带的复杂度Cnb不应高于水平集复杂度Cls，否则总复杂度反而上升.因此，窄带水平集方法的总复杂度主要取决于窄带面积的规模Snb.

当Snb固 定时，生成窄带复杂度Cnb越小越好.直接生成窄带法对每个像素点进行扫描，判断其与过零点的距离是否小于半个窄带宽度，如果是则将其加入窄带.因此，直接生成窄带法的生成窄带复杂度为：=O(N2D2).DTM 窄带法则遍历每个过零点，在半个窄带宽度范围内标记窄带.因此，DTM 的生成窄带复杂度为：=O(ND2).

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

2 本文的矩形窄带法

2.1 寻找过零点

设i,j分别表示2D 图像域上x轴和y轴上的坐标值，I(i,j) 为图像，φ (i,j) 为水平集，C为活动轮廓，则C为φ(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)做以下操作：

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

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

在活动约束之前，过零点长度的规模为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 个宽度，防止过零点暴露在窄带最外层.

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

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

图5 矩形区域合并优化

具体操作如下：设 (sx1,sy1)，(ex1,ey1)分别表示矩形框A的两个端点，(sx2,sy2)，(ex2,ey2)分别表示矩形框B的两个端点，如图5所示.

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

令x1=min(sx1,ex1)，y1=min(sy1,ey1)，x1=min(sx1,ex1)，y2=min(sy2,ey2).判断合并后面积更小的条件为：

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

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

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

2.3 与传统窄带法的比较

传统窄带示意图如图6(a)所示.传统窄带存在很大部分的无效计算区域，实际的活动窄带区域只占总窄带的一小部分.以DTM 方法为例，DTM 方法生成窄带的复杂度为=O(ND2)，窄带面积规模为=O(ND).根据式(5)，DTM 方法的总复杂度为：

本文的矩形窄带示意图如图6(b)所示.通过活动约束，对长期禁止不动的点进行了屏蔽，大幅减少了窄带的生成范围.扫描生成矩形窄带的复杂度为：=O(D)，窄带面积规模为=O(D2).根据式(5)，矩形窄带法的总复杂度为：

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

图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)泛指窄带面积规模为Snb=O(ND)的一类窄带法，如直接生成窄带法，快速进行法，快速扫描法和DTMP 窄带法等.AcNb 为本文提出的过渡方案，RecNb 为本文最终采用方案.根据前文分析，不规则窄带矩形化的过程可能略微增加窄带区域.因此，RecNb 相比AcNb 面积可能略微增加，但AcNb 的计算和实现更为简单.RecNb 和AcNb 的面积规模均为Snb=O(D2)，其中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]标准来衡量.

图9 分割结果对比.第1 行：原图像；第2 行：LATE 模型分割结果；第3 行：LATE 结合DRCT 窄带法(DRCTLS)的分割结果；第4 行：LATE 结合DTM 窄带法(DTMLS)的分割结果；第5 行：我们方法(RECLS)的分割结果.(a)为灰度均匀图像；(b)-(e)为灰度不均匀图像；(f)-(j)为灰度严重不均匀图像.

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

其中，Ot为标准分割区域，Om为实际分割区域，算符A(·)表示求对应区域的面积.JSC 的取值范围在0 到1，JSC 的值越大，分割结果越准确.

图10(a)为二值图像，图10(b)-图10(e)灰度不均匀程度依次递增.以二值图10(a)的黑色区域作为标准分割区域Ot，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 值