主动轮廓模型

自述：该方法的核心是将像素转换成距离场，当目标区域与背景区域分割效果较好时，能量达到最小，即能量最不活跃。
由于时间的关系，没有将此文献完全看懂并且与代码对应；代码中的kappa关于曲率或者散度的计算与该论文给定的散度div计算公式不能实现匹配，这是导致没有继续研究的技术问题；
文献：
1、TFT-LCD Mura缺陷机器视觉检测方法研究-卢小鹏（本文主要参考）
2、https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2796112/局部区域分割
3、http://www.shawnlankton.com/?s=chan-vese局部区域分割的相关代码
4、http://mathworld.wolfram.com/Curvature.html曲率的计算
5、https://www.mathworks.com/matlabcentral/fileexchange/23445-chan-vese-active-contours-without-edges Chan-vese分割代码
6、http://www.mathworks.com/matlabcentral/fileexchange/30284-active-contours-implementation—test-platform-gui/content/Activeontours/localized_seg.m代码
7、A level set approach for computing solutions to incompressible twophase flow
停止迭代的条件一

1、主动轮廓模型

1988年，Kass等人提出了主动轮廓模型，将图像分割问题转换为求解能量泛函最小值问题，为图像分割提供一种全新的思路，称为研究的重点和热点。主动轮廓模型的主要原理通过构造能量泛函，在能量函数最小值驱动下，轮廓曲线逐渐向待检测物体的边缘逼近，最终分割出目标。由于主动轮廓模型利用曲线演化定位目标的边缘，因此也称为Snake模型。主动轮廓模型是当前应用最多的利用变分思想求解的图像分割方法。其最大优点是在高噪声的情况下，也能得到连续、光滑的闭合分割边界。按照能量函数构造方式的不同，可以将主动轮廓模型主要分为基于边缘和基于区域两类，同时也有一些研究人员提出了基于边缘和区域相结合的主动轮廓模型。

2、基于边缘的主动轮廓模型

目标的边缘是表征目标的重要信息，也是分割目标的重要依据。Kass等人提出的Snake模型是一种典型的基于边缘的主动轮廓模型，其主要是根据目标边缘的梯度跳变来检测目标的。
对于图像<script type="math/tex" id="MathJax-Element-1">I(x,y),C(S)=C(x(s),y(s))</script>为图像内的演化曲线，则Snake模型的能量函数定义为:

3、基于区域的主动轮廓模型

C-V是一个典型的区域主动轮廓模型，它以图像的像素灰度信息作为能量，巧妙的构造能量函数，然后通过求取能量函数的最小值，最终把目标分割出来。
假设待分割的图像为<script type="math/tex" id="MathJax-Element-15">I</script>,<script type="math/tex" id="MathJax-Element-16">I_0</script>表示带分割的目标，<script type="math/tex" id="MathJax-Element-17">I_b</script>表示图像背景，<script type="math/tex" id="MathJax-Element-18">C</script>表示图像内的一条闭合曲线。则闭合曲线<script type="math/tex" id="MathJax-Element-19">C</script>把图像<script type="math/tex" id="MathJax-Element-20">I</script>分为曲线<script type="math/tex" id="MathJax-Element-21">C</script>内部区域和外部区域两个部分，其能量函数如下:

4、水平集方法(Level Set)

文摘1：基于水平集的图像分割算法是一种进化版的Snake算法，需要给定初始的轮廓曲线，然后根据泛函能量最小化，进行曲线演化。水平集方法用的是一种隐式函数方法，跟传统的snake算法相比在思想上差异很大，snake算法曲线演化的时候，是曲线上离散点显示坐标的位置更新移动，只要懂得能量最小化的曲线演化规则，就可以很快理解算法，并写出代码；然而水平集方法，更新的不是曲线离散点的坐标，而是更新整张图片像素点到曲线的有向距离场。因此算法最关键的是理解这个距离场的更新规则。
我们先了解显示二维曲线与隐式曲线(水平集函数)的相互转换公式。给定初始的轮廓曲线C，我们将它转换成水平集函数，水平集函数定义为:

{ d[( x, y), C]− d[( x, y), C](x,y)在封闭曲线C的外部(x,y)在封闭曲线C的内部

</script>也就是说，如果给你一条初始封闭轮廓曲线C，进行水平集图像分割，我们需要些的第一个函数就是计算图像的每个像素点 <script type="math/tex" id="MathJax-Element-55">p(x,y)</script>到曲线的最短距离 <script type="math/tex" id="MathJax-Element-56">d</script>，如果该像素点 <script type="math/tex" id="MathJax-Element-57">p</script>位于曲线C的内部，那么有向距离为 <script type="math/tex" id="MathJax-Element-58">-d</script>，反之为 <script type="math/tex" id="MathJax-Element-59">d</script>。这样遍历图像每个像素点，每个像素点都可以求得对应的有向距离 <script type="math/tex" id="MathJax-Element-60">u(x,y)</script>。
隐式距离场向显示离散曲线转换:
只要求出满足 <script type="math/tex" id="MathJax-Element-61">u(x,y)=0</script>的点，就是曲线上的点，因为如果该像素点到曲线C的最短距离为0，那么这个像素点肯定在这条曲线上，因此只要把满足 <script type="math/tex" id="MathJax-Element-62">u(x,y)=0</script>的像素点全部提取出来，获得这些点的坐标 <script type="math/tex" id="MathJax-Element-63">p(x,y)</script>而这些点便是曲线 <script type="math/tex" id="MathJax-Element-64">C</script>的离散点，这样就完成了从隐式距离场到显式离散场曲线的转换。
据此我们可以得到算法的大体流程：
输入：给定离散的初始轮廓曲线C，待分割图像T
输出：分割结果曲线C’
步骤：
1、根据公式，计算每个像素点到离散曲线C的最短有向距离u(x,y)
2、根据图像梯度等信息，对u(x,y)进行演化，使得其沿着能量最小化的方向演化，这个过程说的简单一点，就是更新每个像素点的u(x,y)值。
3、根据第2步的演化结果，遍历每个像素点(x,y)，判断其水平集函数值是否为零。 If u(x,y)==0:保存像素点坐标(x,y)(因为这个点就是曲线C’上的点)
得到所有u(x,y)==0的点，就是最后我们想要的图像分割结果曲线C’
文献引用2：在二维平面内，显性表示的函数曲线在演化过程中无法随着目标的合并或者分裂进行拓扑变化，因此在表示目标拓扑变化方面，显性表示的二维函数曲线就无能为力。为了解决这个问题，1988年，Osher和Sethian用隐式参数函数表示演化曲线，进而提出了水平集方法。水平集方法的提出，最早是为了解决火苗在热力学方程下的膨胀、收缩、分裂等形状变化过程。由于水平集函数能够灵活地表示目标的拓扑变化，在图像分割邻域得到了广泛的应用。
水平集方法的基本原理是把目标曲线或者曲面作为零水平集嵌入到更高一维的水平集函数中，也即用零平面截取水平集函数得到的封闭曲线或者曲面来代替演化曲线或者曲面，随着水平集函数的变化，演化曲线或者曲面也随之改变，并能适应拓扑的变化。
在二维平面，如果隐式表示的闭合曲线为： <script type="math/tex" id="MathJax-Element-65">C(x,y)=0</script>，根据水平集原理，把此闭合曲线嵌入到三维水平集函数 <script type="math/tex" id="MathJax-Element-66">z=\varphi(x,y)</script>中，则用 <script type="math/tex" id="MathJax-Element-67">z=0</script>平面截此水平集函数的曲面，即可得到闭合曲线 <script type="math/tex" id="MathJax-Element-68">C(x,y)=0</script>。如下图所示。

主动轮廓模型
Figure 1. 二维曲线的水平集表示

当三维水平集曲面 <script type="math/tex" id="MathJax-Element-69">z=\varphi(x,y)</script>在驱动力的作用下发生变化时，其零水平集平面截得的闭合曲线 <script type="math/tex" id="MathJax-Element-70">C(x,y)=0</script>也随之变化。
在Mura缺陷分割中，设C-V模型的轮廓曲线为 <script type="math/tex" id="MathJax-Element-71">C</script>，轮廓曲线在驱动力的作用下，随着时间的推移进行演化，可得到演化曲线集 <script type="math/tex" id="MathJax-Element-72">C(t)</script>。根据水平集原理，把轮廓曲线嵌入到三维水平集函数 <script type="math/tex" id="MathJax-Element-73">\varphi(x,y,t)</script>，其零水平集可以表示为 <script type="math/tex" id="MathJax-Element-74">\{\varphi(x,y,t)=0\}</script>。轮廓曲线随水平集函数演化过程如下图所示:

主动轮廓模型
Figure 1. 基于水平集方法的演化曲线变化过程

在上图中， <script type="math/tex" id="MathJax-Element-75">C(t)</script>表示轮廓曲线， <script type="math/tex" id="MathJax-Element-76">z=\varphi(x,y,t)</script>表示三维水平集函数， <script type="math/tex" id="MathJax-Element-77">C(t_0)</script>表示在 <script type="math/tex" id="MathJax-Element-78">t_0</script>时刻获得的轮廓曲线，此时的分割结果是一个单联通区域。随着时间的推移，在 <script type="math/tex" id="MathJax-Element-79">t_1</script>时刻获得的轮廓曲线曲线变为 <script type="math/tex" id="MathJax-Element-80">C(t_1)</script>，此时轮廓曲线分裂为两个单联通区域。也即随着轮廓曲线的演化，可以检测出两个Mura曲线区域，甚至更多，这是显性的二维函数无法实现的，这正是本文采用水平集方法的原因。
三维水平集函数 <script type="math/tex" id="MathJax-Element-81">z=\varphi(x,y,t)</script>的零水平为:

通过水平集方法推导，我们可以得到轮廓曲线内部的平均灰度<script type="math/tex" id="MathJax-Element-93">c_0</script>和轮廓曲线外部的平均灰度<script type="math/tex" id="MathJax-Element-94">c_b</script>的计算公式如下:

{10 z≥0 z<0

</script>
<script type="math/tex" id="MathJax-Element-100">\delta(z)</script>为Dirac delta函数:

2、改进的Chan-Vese模型

为了提高C-V模型的分割能力，把目标的平均灰度与背景的平均灰度的差的平方作为能量项加入模型中，该能量项表达式如下:

F( c0, c b, C)== E( C) μ L( C)+ v∗ A r e a( i n s i d e( C))∗( c0− c b)2+ λ1∫ i n s i d e( C)| I( x, y)− c0|2d xd y+ λ2∫ o u t s i d e( C)| I( x, y)− c b|2d xd y

</script>
其中:

3、半隐差分格式的水平集解法

由于水平集方法在图像处理领域的广泛应用，人们对水平集方法的数值实现也做了大量的研究，并且提出了许多求解方法。
<script type="math/tex" id="MathJax-Element-106">\blacksquare </script>有限差分法：为了使轮廓曲线能够收敛至目标的边缘，一般采用较小的时间步长，导致曲线演化速度很慢。
在曲线演化的过程中，零水平集表示轮廓曲线，因此可以在不更新整个水平集函数的情况下，只更新零水平集即可达到曲线演化的目的。基于此，人们又提出了许多水平集快速解法，如窄带法，快速行进法，AOS法等。
<script type="math/tex" id="MathJax-Element-107">\blacksquare </script>窄带法求解:由于C-V模型的能量函数涉及了图像中所有像素点的灰度值，因此并不适合利用窄带法求解。
<script type="math/tex" id="MathJax-Element-108">\blacksquare </script>快速行进法:快速行进法求解效率高，但演化曲线只能向单一的方向收敛，对于位置不定的mura缺陷分割，快速行进法应用受到了限制。
<script type="math/tex" id="MathJax-Element-109">\blacksquare </script>AOS格式求解方法:AOS格式求解方法是无条件稳定，并且时间步长较大，但是求解过程中需要矩阵的求逆运算，不仅增加了运算量，还会影响分割结果的准确性。
<script type="math/tex" id="MathJax-Element-110">\blacksquare </script>半隐差分格式水平集解法：本文构造具有无条件稳定特性的半隐差分格式对改进的C-V模型进行求解，在不降低模型分割能力的同时，提高分割的速度。
<script type="math/tex" id="MathJax-Element-111">h</script>:空域的步长(space step)对于图像我们取该值为1
1）半隐差分格式构造
对水平函数集<script type="math/tex" id="MathJax-Element-112">\varphi</script>在X轴和Y轴的偏微分：

1、Our level set function is denoted as <script type="math/tex" id="MathJax-Element-128">\phi</script> and it is taken positive outside the bubble and negative inside the bubble.Therefore,the bubble interface is the zero level set of <script type="math/tex" id="MathJax-Element-129">\phi</script>.(We shall initialize <script type="math/tex" id="MathJax-Element-130">\phi</script> to be the signed normal distance from the interface)
2、It should be noted that while <script type="math/tex" id="MathJax-Element-131">\phi</script> is initially a distance function it will not remain so,we shall present a novel way of reinitializing <script type="math/tex" id="MathJax-Element-132">\phi</script>so that it remains a distance function.
3、One possible discretization of (17) is as follows. we define

⎧⎩⎨⎪⎪ m a x(( a+)2,( b−)2)+ m a x(( c+)2,( d−)2)−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√−1 m a x(( a−)2,( b+)2)+ m a x(( c−)2,( d+)2)−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−√−10 ϕ0 i, j>0 ϕ0 i, j<0 o t h e r w i s e

</script>
where the + superscript denotes the positive part and the - superscript denotes the negative part.Equation(17) is then updataed using

参考：

1、TFT-LCD Mura缺陷机器视觉检测方法研究-卢小鹏。
2、A level set approach for computing solutions to incompressible twophase flow《曲线重新初始化方式》