2. 湖南软件职业学院 软件与信息工程学院, 湘潭 411100
2. School of Software and Information Engineering, Hunan Software Vocational College, Xiangtan 411100, China
水平集最先由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 水平集设
水平集方法可概括为: 首先确立合理的水平集微分方程
以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) |
其中,
$ LI{M_i}\left( x \right) = mean\left( {I(y):\;y \in {\Omega _x} \cap {R_i}} \right) $ | (2) |
其中,
${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) |
其中,
水平集方法中, 称每次迭代的计算复杂度为水平集复杂度, 记为
由于光照磁场或成像装置缺陷等因素[7], 可能导致获得的图像是灰度不均匀的. 灰度不均匀给图像分割带来很大的挑战. 在文献[8]中, 灰度不均匀图像被看作灰度均匀图像与偏置光源场的混合. 蝴蝶图像与严重不均匀偏置场的混合如图2所示.
不同水平集模型对灰度严重不均匀蝴蝶图像的分割结果如图3所示. 可见, LATE模型具有很强的分割灰度不均匀图像的能力. 但是, LATE模型的计算代价较高, 分割缓慢.
1.2 传统窄带及其复杂度分析
虽然水平集能够获得较好的分割结果, 但它将问题提高了一个维度, 增加了计算复杂度. 直观地看, 活动轮廓
假设传统窄带的宽度与卷积模板宽度具有相同规模, 传统窄带长度规模与图像宽度相当. 则传统窄带的面积规模可以记为:
窄带水平集方法每次迭代的复杂度称为总复杂度, 记为
${C_{all}} = {C_{nb}} + \dfrac{{{S_{nb}}}}{{{N^2}}}{C_{ls}}$ | (5) |
一般地, 生成窄带的复杂度
当
快速进行法, 快速扫描法以及DTM方法等窄带法都只是降低了生成窄带复杂度
设
设
$\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}$ |
则可找到当前所有的过零点
在水平集演化的过程中, 一部分活动轮廓先到达图像边缘不再运动, 另一部分活动轮廓还需继续演化逼近图像边缘. 先到达图像边缘的活动轮廓在后续演化过程中的计算属于冗余计算, 还未达到图像边缘的活动轮廓才真正需要进行演化计算. 如果在静止的过零点附近生成窄带, 对活动轮廓的演化并没有作用. 可以设置延时参数
$\begin{array}{l} {\rm{if}}\;\;\;\;\;\;VIS(i,j) >{{ delay}} \\ {\rm{then}}\;\;\;CRS(i,j) = 0 \\ \end{array} $ |
在活动约束之前, 过零点长度的规模为N. 假设
传统窄带的生成, 是在过零点基础上通过偏移生成带状区域的过程. 传统方法获得的过零点曲线是封闭的, 而经过活动约束, 本文方法获得的过零点集合只是传统过零点曲线的一部分, 是多段开放的链, 称为过零点链. 原则上, 应该在每条过零点链的邻域内生成不规则的窄带. 由于活动约束下每条过零点链长度较短, 窄带区域接近矩形区域, 不妨用规则的矩形区域代替不规则的窄带区域. 传统的窄带区域是不规则的, 需要单独建立窄带中每个点与原图像域之间的映射, 才能将窄带的计算结果返回到图像域. 而矩形区域则只需给出矩形的两个对角端点便可确定一个区域, 且矩形域与图像域的像素点位置对应非常简单, 这更有利于快速卷积运算以及边界处理.
可利用扫描法获得每条过零点链对应的矩形窄带区域, 步骤如下:
1)分别对
2)若在8邻域内扫描到
3)重复步骤2), 直到在8邻域内找不到下一个
4)确定矩形域. 取链码
活动约束后, 剩余过零点的长度规模为
寻找过零点方法具有简单高效的优点, 但可能存在个别过零点的遗落, 导致过零点曲线不连续, 从而产生面积重叠的矩形区域. 如图5中的左图所示, 矩形A与矩形B产生重叠. 记矩形A的面积为
具体操作如下: 设
则判断它们相交的条件为:
$\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} $ |
令
$\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. 初始化
2. 对访问次数超过delay的像素点进行限制, 将过零点(
3. 利用扫描法生成每段过零点链的最小覆盖矩阵框;
4. 将最小矩形框两两之间进行判断, 如果它们的区域相交且合并后面积更小, 则对它们进行合并优化.
由于图像也是矩形数组, 矩形框的窄带结果和原图像能保持一致, 直接将每个矩形框代入原水平集函数计算即可. 返回结果时只需做坐标偏移处理, 便可实现窄带到原图像的映射, 从而在本次迭代中实现水平集的更新. 可见, 本文的方法与水平集方法非常容易实现对接.
2.3 与传统窄带法的比较传统窄带示意图如图6(a)所示. 传统窄带存在很大部分的无效计算区域, 实际的活动窄带区域只占总窄带的一小部分. 以DTM方法为例, 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)所示. 通过活动约束, 对长期禁止不动的点进行了屏蔽, 大幅减少了窄带的生成范围. 扫描生成矩形窄带的复杂度为:
$\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, 所以
3 数值实验 3.1 窄带面积比较
本文的实验结果均在MATLAB R2016a上实现, 操作系统为Win10. 图7(a)(b)分别为针对不同灰度不均匀情形的两组实验. 实验采用矩形窄带法结合LATE模型[10]生成水平集, 并利用水平集过零点分别生成传统窄带(TradNb)、仅添加活动约束的不规则窄带(AcNb)以及矩形窄带(RecNb).
本文中传统窄带法(TradNb)泛指窄带面积规模为
图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. 可见, 本文的矩形窄带法能有效地减少窄带面积, 从而提高计算效率.
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) |
其中,
图10(a)为二值图像, 图10(b)–图10(e)灰度不均匀程度依次递增. 以二值图10(a)的黑色区域作为标准分割区域
可见, LATE方法对灰度严重不均匀图像具有较高的分割精度. RECLS方法对应的JSC值与LATE几乎一致. 结合上一节的结论, 本文提出的矩形窄带方法能在不影响LATE模型分割精度的条件下, 提高对灰度严重不均匀图像的分割效率.
4 结论与展望本文提出一种新的矩形窄带方法. 通过活动约束进一步缩小了窄带的范围. 利用矩形窄带代替不规则窄带, 使其更容易与水平集方法相结合, 减少了更新水平集的计算量. 实验表明, 本文的方法即使在灰度严重不均匀情形下也能够保持稳定的分割结果.
本文的方法在窄带演化的过程中, 可能存在多个矩形窄带, 而这些窄带没有实现并行运算. 如何让多个矩形窄带区域实现并行运算, 进一步提高计算的效率, 是我们下一步研究的内容.
[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.
|