基于区域的主动轮廓模型_从零到一用python实现Chan-Vese(CV)算法

1、Chan-Vese(CV) 算法原理

同样是构造能量函数，以snake为代表的基于边缘的主动轮廓模型,其目的在于使曲线演化至图像特定边缘，而以Chan-CV模型为代表的基于区域的主动轮廓模型目的则在于演化出以目标曲线为边界的特定区域。总的来说，snake等基于边缘的模型目的在于找目标曲线，而CV模型则在于找出特定区域，虽然本质任务相同，但实现方法上还是存在巨大差异的。
CV模型使用图像的像素灰度信息作为能量，构造能量，从而将曲线演化至目标区域,对于某一图像$(x,y)\in \Omega$, 用$I(x,y)$表示图像上点$(x,y)$的像素值,定义封闭曲线 $C$为图像上一闭合区域 $\omega$ 的边界,$inside(C)$表示曲线内闭合区域 $\omega$($C=\partial \omega$)，$outside(C)$则表示曲线外区域 $\overline{\omega }$，$C_1$和$C_2$分别表示图像上曲线内和曲线外像素平均值,并定义如下的能量项：$$E_1(C)+E_2(C)={\int }_{inside(C)}|I_0(x,y)-C_1{|}^2\mathrm{d}x\mathrm{d}y+{\int }_{outside(C)}|I_0(x,y)-C_2{|}^2\mathrm{d}x\mathrm{d}y$$ 曲线$C$和目标物体无非有以下四种关系：
可以分析当曲线C在目标图像中目标边界时，上述能量项取极小值，在补充一些正则化项目后，我们可以得出如下的能量泛函：$$E(C)=\mu Length(C)+\upsilon Area(C)+E_1(C)+E_2(C)$$

2、基于水平集的实现方法

2.1演化方程推导

具体的实现过程需要依赖水平集方法，我相信，即使看过网上很多的分析帖子，但是大家尤其是水平集初次接触者都会觉得这种方法很表面上看起来很有道理，但是仔细琢磨琢磨又说不出个所以然，心里可能存在一些困惑比如说，到底啥是水平集方法？为什么水平集函数的演化结果就能够实现和直接曲线演化同样的结果？这个水平集方法我回头会单独开一节进行讨论，不过一下几点心得方法可能会帮助大家理解水平集方法：

1.水平集是一种曲线演化的方法，本质上来讲就图像分割而言，和主动轮廓模型并没有什么关系，只不过水平集方法具有巨大的计算优势，能够极大方便模型后续的实现过程，因此很多时候大家都直接会把水平集函数嵌入到能量泛函，然后再将利用能量泛函极小值化(变分水平集方法)，获得演化方程；其实除了这种方法外，完全可以在能量函数中先不引进水平集函数，直接对能量泛函求极值获得演化方程，然后再引入水平集函数进行求解，这个我会专门写一篇文章进行讲解;
2. 回答一下上面我提出的问题：为什么用水平集函数的演化能够替代直接曲线演化,或者说基于同样的分割原理，为什么水平集函数的演化能够和直接用曲线演化产生同样的结果，并且效率更高呢？因为对于曲线演化来讲，任一点的演化速度都可以用曲线上的切线方向速度和法线速度构成，而切线方向速度分量并不会改变曲线形状(至于为什么以后我会用数学公式推导证明) 而曲线上任一点的法向速度方向恰好和三维水平集函数的梯度方向相同，那么你会法向，只要控制能量泛函的外力相同，那么你会发现控制曲线演化的速度本质上和控制水平集函数演化的函数本质上是一样的(以后再补图)。
3. 在具体的实现过程中，希望大家能够摆脱一种曲线演化的概念，在采用水平集方法CV算法实现时，我们的焦点已经不再局限于曲线以及曲线的演化之上，而在于基于整个图像区域的水平集函数演化,简言之，水平集是面(图像面)的函数，而不是演化曲线的函数,当水平集函数演化时，零水平集即目标曲线也会随之演化(从而可以自然地改变其拓扑结构)，当水平集函数的演化停止时，目标曲线也就自动勾勒出来了。
4. 使用水平集方法确实存在很多优势，但是这个优势是需要和曲线演化进行对比分析的，这个我以后写点东西让大家慢慢品。

有了上面这些概念，大家对水平集函数应该有一些具体的概念了，那么接下来我们就要引入水平集函数$\varphi (x,y)$对能量函数进行调整，首先对水平集函数做如下定义，
{ C = ∂ ω = { ( x , y ) ∈ Ω : ϕ ( x , y ) = 0 } i n s i d e ( C ) = ω = { ( x , y ) ∈ Ω : ϕ ( x , y ) > 0 } o u t s i d e ( C ) = Ω \ ω ‾ = { ( x , y ) ∈ Ω : ϕ ( x , y ) < 0 } \begin{cases} C = \partial {\omega} = \lbrace(x,y)\in \Omega :\phi(x,y) =0 \rbrace \\[2ex] inside(C) = \omega = \lbrace(x,y)\in \Omega :\phi(x,y) >0 \rbrace \\[2ex] outside(C) = \Omega \backslash \overline{\omega} = \lbrace(x,y)\in \Omega :\phi(x,y) < 0 \rbrace \end{cases} ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​C=∂ω={(x,y)∈Ω:ϕ(x,y)=0}inside(C)=ω={(x,y)∈Ω:ϕ(x,y)>0}outside(C)=Ω\ω={(x,y)∈Ω:ϕ(x,y)<0}​这样就完成了使用水平集函数对演化曲线的替代，进一步我们定义Heaviside Fucntion(阶跃函数)：
$$H(z)=\{\begin{array}{ll}1,& \text{if z>0}\\ 0,& \text{if z<0},\end{array}，{\delta }_0(z)=\frac{\mathrm{d}}{\mathrm{d}\mathrm{z}}H(z)$$
从而可得：$$Length(\varphi =0)={\int }_{\Omega }|\nabla H(\varphi (x,y))|dxdy$$ $$Area(\varphi \ge 0)={\int }_{\Omega }H(\varphi (x,y))dxdy$$ $${\int }_{inside(C)}|I_0(x,y)-C_1{|}^2\mathrm{d}x\mathrm{d}y={\int }_{\Omega }|I_0(x,y)-C_1{|}^{2}H(\varphi (x,y)\mathrm{d}x\mathrm{d}y$$ $${\int }_{outside(C)}|I_0(x,y)-C_2{|}^2(1-H(\varphi (x,y))\mathrm{d}x\mathrm{d}y$$总的能量函数也就可以写为：$$\begin{gathered} F(\varphi (x,y))=\mu {\int }_{\Omega }|\nabla H(\varphi (x,y))|dxdy+\nu {\int }_{\Omega }H(\varphi (x,y)dxdy \\ +{\lambda }_1{\int }_{\Omega }|I_0(x,y)-C_1{|}^{2}H(\varphi (x,y)\mathrm{d}x\mathrm{d}y+ \\ {\lambda }_2{\int }_{outside(C)}|I_0(x,y)-C_2{|}^2(1-H(\varphi (x,y))\mathrm{d}x\mathrm{d}y \end{gathered}$$其中 $$C_1=\frac{{\int }_{\Omega }I(x,y)H(\varphi (x,y)dxdy}{{\int }_{\Omega }H(\varphi (x,y)dxdy}$$ $$C_2=\frac{{\int }_{\Omega }I(x,y)(1-H(\varphi (x,y))dxdy}{{\int }_{\Omega }(1-H(\varphi (x,y))dxdy}$$,因为理论上来讲，Heaviside Fucntion(阶跃函数)是不存在的，所以我们选择了一个替代函数：$$H_{ϵ}(z)=\frac{1}{2}(1+\frac{2}{\pi }arctan(\frac{z}{ϵ}))$$ $${\delta }_{ϵ}(z)=\frac{1}{\pi }\cdot \frac{ϵ}{{ϵ}^2+{\pi }^2}$$，当$ϵ$的取值越小，替代函数也就越接近理想中的的阶跃函数。
至此，有关能量泛函的所有细节已经交代完毕,然后利用变分法和梯度下降流(梯度下降法求解能量泛函极值)获得能量泛函的偏微分方程：$$\frac{\partial \varphi }{\partial t}={\delta }_{ϵ}(\varphi )[\mu div(\frac{\nabla \varphi }{|\nabla \varphi |})-\nu -{\lambda }_1(I-C_1{)}^2+{\lambda }_2(I-C_2{)}^2]$$

2.2离散化过程

在对上述偏微分方程进行离散化前，首先定义$h$表示空间补步长,$\Delta t$表示时间步长，那么我们就可以用${\varphi }_{i,j}^n=\varphi (n\Delta ,x_i,y_i)$作为$\varphi (t,x,y)$的近似，然后定义如下差分式子: $${\Delta }_{-}^x\varphi (x,y)={\varphi }_{i,j}-{\varphi }_{i-1,j},{\Delta }_{+}^x\varphi (x,y)={\varphi }_{i+1,j}-{\varphi }_{i,j}$$ $${\Delta }_{-}^y\varphi (x,y)={\varphi }_{i,j}-{\varphi }_{i,j-1},{\Delta }_{+}^x\varphi (x,y)={\varphi }_{i,j+1}-{\varphi }_{i,j}$$然后我们根据$C_1$和$C_2$的计算公式计算出对应值.
最后对应的离散化方程式就为(这可能是我打出的史上最长公式)：$$\begin{gathered} \frac{{\varphi }_{i,j}^{n+1}-{\varphi }_{i,j}^n}{\Delta t}={\delta }_n({\varphi }_{i,j}^n)[\frac{\mu }{h^2}{\Delta }_{-}^x\cdot (\frac{{\Delta }_{+}^x{\varphi }_{i,j}^{n+1}}{\sqrt{({\Delta }_{+}^x{\varphi }_{i,j}^n{)}^2/(h^2)+({\varphi }_{i,j+1}-{\varphi }_{i,j-1}{)}^2/(2h{)}^2}}) \\ +\frac{\mu }{h^2}{\Delta }_{-}^y\cdot (\frac{{\Delta }_{+}^y{\varphi }_{i,j}^{n+1}}{\sqrt{({\Delta }_{+}^y{\varphi }_{i,j}^n{)}^2/(h^2)+({\varphi }_{i+1,i}-{\varphi }_{j-1,i}{)}^2/(2h{)}^2}})-\upsilon -{\lambda }_1(I_{i,j}-C_1({\varphi }^n){)}^2+{\lambda }_2(I_{i,j}-C_2({\varphi }^n){)}^2] \end{gathered}$$

2.3 代码实现

那么接下来就是进行代码实现啦,代码实现主要分为一下主要部分：

1. 读入图片
   

Image = cv2.imread('1.bmp',1)
image = cv2.cvtColor(Image,cv2.COLOR_BGR2GRAY)
img = np.array(image,dtype='float64')
plt.imshow(img,cmap = 'gray')
print(img.shape)

2. 初始化水平集函数

IniLSF = np.ones([img.shape[0],img.shape[1]],img.dtype)
IniLSF[40:80,40:80]=-1
IniLSF = -IniLSF

3. 定义水平集算法

def CV(LSF,img,nu,mu,epison,step):
    
    Drc = (epison/math.pi)/(epison*epison+LSF*LSF)
#    Hea = 0.5*(1+(2/math.pi)*mat_math(LSF/epison,'atan'))
    Hea = 0.5*(1+(2/math.pi)*np.arctan(LSF/epison))
    Iy,Ix = np.gradient(LSF)##q4#
#    s = mat_math(Ix*Ix+Iy*Iy,"sqrt")
    s = np.sqrt(Ix*Ix+Iy*Iy)
    Nx = Ix/(s+0.000001)
    Ny = Iy/(s+0.000001)
    Mxx,Nxx = np.gradient(Nx)
    Nyy,Myy = np.gradient(Ny)
    cur = Nxx + Nyy
    
    Length = nu*Drc*cur
    
    Area = mu*Drc
    
    s1 = Hea*img
    s2 = (1-Hea)*img
    s3 = 1-Hea
    C1 = s1.sum()/Hea.sum()
    C2 = s2.sum()/s3.sum()
    
    CVterm = Drc*(-1*(img-C1)*(img-C1)+1*(img-C2)*(img-C2))
    LSF = LSF+step*(Length+Area+CVterm)
    return LSF

4. 定义关键参数，并迭代计算，

nu = 0.0001*255*255
mu = 1
num = 10
epison = 1
step = 0.1
LSF = IniLSF
for i in range(1,num):
    LSF = CV(LSF,img,nu,mu,epison,step)
    if i%2 == 0:
        plt.contour(LSF,[0],linewidths = 3.0,linestyles = 'dotted',colors='r')
#    plt.draw()
plt.show()
print("over")

2.4 讨论分析