基于PDE的图像分割

基于偏微分方程的图像分割

比较常见的用于图像分割的 PDE 模型有蛇（Snake）模型、测地线活动轮廓模型（Geodesic activate contour, GAC）和 CV 模型。以下介绍这三种模型以及它们对应的梯度下降流。

一、蛇模型

蛇模型的方程为
$$E(C(p))={\int }_0^1\alpha |C_p(p,t)|+\beta |C_{pp}(p,t){|}^2-\lambda |\nabla I(C(p))|dp$$
第一项是曲线弧长积分，最小化能量泛函也就是使得找到的闭合曲线尽可能的短。第二项是曲率，此项要求曲线尽可能光滑。第三项是找到使得图像梯度达到极大值的曲线，即使曲线停在图像的边缘上。因此最小化整个能量泛函可找到分割曲线。

matlab 实现

% =========================================================================  
%                   Snakes：Active Contour Models  
% =========================================================================  
% By gujinjin 2012/12/10-12  Sunny  
% 基于KASS等的论文思想  
% 参考文献：  
% [1] KASS etc.Snakes：Active Contour Models  
% [2] CSDN 博客 - Author：乐不思蜀Tone  
% [3] Ritwik Kumar(Harvard University),D.Kroon（Twente University)的程序  
% [4] 《数学建模与数学实验》  
% ------  忘了参考的那篇博客找到的代码了。。。。。如有作者看到，可以标记上参考。
clc;clf;clear all;  
  
% =========================================================================  
%                      获取手动取点坐标  
% =========================================================================  
% 读取显示图像  
%I = imread('Coronary_MPR.jpg');  
I =  imread('t2.jpg');  
% 转化为双精度型  
%I = im2double(I);   
% 若为彩色，转化为灰度  
if(size(I,3)==3), I=rgb2gray(I); end  
figure(1),imshow(I);  
%---------------------------  
%        高斯滤波  
%---------------------------  
sigma=1.0;  
% 创建特定形式的二维高斯滤波器H  
H = fspecial('gaussian',ceil(3*sigma), sigma);  
% 对图像进行高斯滤波,返回和I等大小矩阵  
Igs = filter2(H,I,'same');  
%---------------------------  
%      获取Snake的点坐标  
%---------------------------  
figure(2),imshow(Igs);  
x=[];y=[];c=1;N=100; %定义取点个数c,上限N  
% 获取User手动取点的坐标  
% [x,y]=getpts  
while c<N  
    [xi,yi,button]=ginput(1);  
    % 获取坐标向量  
    x=[x xi];  
    y=[y yi];  
    hold on  % 用来保持要画的图坐标不变
    % text(xi,yi,'o','FontSize',10,'Color','red');  
    plot(xi,yi,'ro');  
    % 若为右击，则停止循环  
    if(button==3), break; end  
    c=c+1;  
end  
  
% 将第一个点复制到矩阵最后，构成Snake环  
xy = [x;y];  
c=c+1;  
xy(:,c)=xy(:,1);  
% 样条曲线差值  
t=1:c;  
ts = 1:0.1:c;  
xys = spline(t,xy,ts);  % 按后面ts为查询点进行插值
xs = xys(1,:);  
ys = xys(2,:);  
% 样条差值效果  
hold on  
temp=plot(x(1),y(1),'ro',xs,ys,'b.');  
legend(temp,'原点','插值点');  
  
% =========================================================================  
%                     Snakes算法实现部分  
% =========================================================================  
NIter =1000; % 迭代次数  
alpha=0.2; beta=0.2; gamma = 1; kappa = 0.1;  
wl = 0; we=0.4; wt=0;  
[row col] = size(Igs);  
  
% 图像力-线函数  
Eline = Igs;  
% 图像力-边函数  
[gx,gy]=gradient(Igs);  
Eedge = -1*sqrt((gx.*gx+gy.*gy));  % 边缘
% 图像力-终点函数  
% 卷积是为了求解偏导数，而离散点的偏导即差分求解  
m1 = [-1 1];   
m2 = [-1;1];  
m3 = [1 -2 1];   
m4 = [1;-2;1];  
m5 = [1 -1;-1 1];  
cx = conv2(Igs,m1,'same');  
cy = conv2(Igs,m2,'same');  
cxx = conv2(Igs,m3,'same');  
cyy = conv2(Igs,m4,'same');  
cxy = conv2(Igs,m5,'same');  
  
for i = 1:row  
    for j= 1:col  
        Eterm(i,j) = (cyy(i,j)*cx(i,j)*cx(i,j) -2 *cxy(i,j)*cx(i,j)*cy(i,j) + cxx(i,j)*cy(i,j)*cy(i,j))/((1+cx(i,j)*cx(i,j) + cy(i,j)*cy(i,j))^1.5);  
    end  
end  
  
%figure(3),imshow(Eterm);  
%figure(4),imshow(abs(Eedge));  
% 外部力 Eext = Eimage + Econ  
Eext = wl*Eline + we*Eedge + wt*Eterm;  
% 计算梯度  
[fx,fy]=gradient(Eext);  
  
xs=xs';  
ys=ys';  
[m n] = size(xs);  
[mm nn] = size(fx);  
  
% 计算五对角状矩阵  
% 附录: 公式（14） b(i)表示vi系数(i=i-2 到 i+2)  
b(1)=beta;  
b(2)=-(alpha + 4*beta);  
b(3)=(2*alpha + 6 *beta);  
b(4)=b(2);  
b(5)=b(1);  
  
A=b(1)*circshift(eye(m),2);  % circshift(A,K, m) 是循环移位函数，使得矩阵A移位K个，m决定行移位或列移位，m=1是列移，m=2是行移
A=A+b(2)*circshift(eye(m),1);  
A=A+b(3)*circshift(eye(m),0);  
A=A+b(4)*circshift(eye(m),-1);  
A=A+b(5)*circshift(eye(m),-2);  
% 以上操作构成了矩阵A
% 计算矩阵的逆  
[L U] = lu(A + gamma.* eye(m));  
Ainv = inv(U) * inv(L);   
  
% =========================================================================  
%                      画图部分  
% =========================================================================  
%text(10,10,'+','FontName','Time','FontSize',20,'Color','red');  
% 迭代计算  
figure(3)  
for i=1:NIter;  
    ssx = gamma*xs - kappa*interp2(fx,xs,ys);  % 二维网格数据插值，根据 xs ys对fx进行插值得到fx
    ssy = gamma*ys - kappa*interp2(fy,xs,ys);  
   
    % 计算snake的新位置  
    xs = Ainv * ssx;  
    ys = Ainv * ssy;  
      
    % 显示snake的新位置  
    imshow(I);   
    hold on;  
    plot([xs; xs(1)], [ys; ys(1)], 'r-');  
    hold off;  
    pause(0.001)      
end  

二、GAC模型

曲线演化理论的研究表明，曲线尽可能短的过程中就会逐渐光滑。因此蛇模型的第二项可以删去。另一方面，蛇模型不仅依赖曲线 $C$ 的形状，还依赖于自由参数 $p$，这提升了方程的求解难度。1997年，V. Caselles等人提出了不含自由参数的 GAC 模型。具体如下：

首先，将自由参数用弧长参数 $s$ 代替，有 $ds=|C_p|dp$，此时，蛇模型的第一项便成为常值。另外，为了避免能量泛函中使用符号，引入单减非负函数 $g$ 放入第三项。得到 GAC 模型
$$E(C(s))={\int }_0^{L(C)}g(\nabla I(C(s)))ds$$
最小化 $E(C(s))$ 得到使得曲线落在使得图像 $|\nabla I|$ 极大值的位置，也就是图像边缘，从而达到分割的目的。

1、GAC 的梯度下降流

假设模型的解可以从某一初始位置迭代到达局部极小值，则求解模型的梯度下降流是必不可少的一步。

梯度下降流类似随机梯度下降法中的梯度，计算能量随时间的变化量，决定曲线演化的下一步方向及大小。

我们仍旧利用参数法表示闭合曲线 $C,C(p):[0,1]\to R^2，C(0)=C(1)$，则 GAC 模型重新写为
$$E(C(p))={\int }_0^{1}g(C(p))|C_p|dp$$
为了表示曲线随时间演化，引入时间变量 $t$，
$$E(C(p,t))={\int }_0^{1}g(C(p,t))|C_p(p,t)|dp$$
求梯度
$$\begin{array}{rl}\frac{\partial E}{\partial t}& ={\int }_0^1\frac{d}{dt}g(C(p,t))|C_p(p,t)|+g(C(p,t))\cdot \frac{d}{dt}|C_p(p,t)|dp\\ & ={\int }_0^1\nabla g\cdot C_t\cdot |C_p|+g\cdot \frac{C_p\cdot C_{pt}}{|C_p|}dp\\ & ={\int }_0^1\nabla g\cdot C_t\cdot |C_p|+g\cdot \text{T}\cdot C_{pt}dp\end{array}$$
$\text{T}$ 表示曲线 $C$ 的单位切向量，对上式第二项做分布积分，
∫ 0 1 g ⋅ T ⋅ C p t d p = ∫ 0 1 − { g ⋅ T } p ⋅ C t d t = ∫ 0 1 − ∇ g ⋅ C p ⋅ T ⋅ C t − g ⋅ T p ⋅ C t d p \begin{aligned} \int_0^1g\cdot \textbf T\cdot C_{pt}dp&=\int_0^1-\{g\cdot \textbf T\}_p\cdot C_tdt \\ &=\int_0^1-\nabla g\cdot C_p\cdot \textbf T \cdot C_t-g\cdot \textbf T_p\cdot C_tdp \end{aligned} ∫01​g⋅T⋅Cpt​dp​=∫01​−{g⋅T}p​⋅Ct​dt=∫01​−∇g⋅Cp​⋅T⋅Ct​−g⋅Tp​⋅Ct​dp​
将此式带入上面梯度中，再有弧长参数 $s$ 代替 $p$，得
$$\frac{\partial E}{\partial t}={\int }_0^{L(C)}(\nabla g-[\nabla g\cdot \text{T}]\text{T}-g\cdot {\text{T}}_s)\cdot C_{t}ds$$
由于
$$\nabla g-[\nabla g\cdot \text{T}]\text{T}=[\nabla g\cdot \text{N}]\text{N},{\text{T}}_s=k\text{N}$$

于是
$$\frac{\partial E}{\partial t}={\int }_0^{L(C)}([\nabla g\cdot \text{N}]\text{N}-g\cdot k\text{N})\cdot C_{t}ds$$
当能量泛函达到稳定状态，即 $\frac{\partial E}{\partial t}=0$， $C_t$ 需要满足
$$\frac{\partial C}{\partial t}=g\cdot k\text{N}-[\nabla g\cdot \text{N}]\text{N}$$
也就是负梯度方向。这就是 GAC 模型的梯度下降流。采用水平集方法进行求解。具体如下：
按水平集方法，$\beta =gk-\nabla g\cdot \text{N},k=-div(\text{N})=-div(\frac{\nabla u}{|\nabla u|})$，因此嵌入函数的梯度为
$$\begin{array}{rl}\frac{\partial u}{\partial t}& =\beta |\nabla u|=(gk-\nabla g\cdot \frac{\nabla u}{|\nabla u|})\cdot |\nabla u|\\ & =(-gdiv(\frac{\nabla u}{|\nabla u|})-\nabla g\cdot \frac{\nabla u}{|\nabla u|})\cdot |\nabla u|\\ & =-div(g\cdot \frac{\nabla u}{|\nabla u|})\cdot |\nabla u|\end{array}$$
我这里的推导与书上差一个负号，不知道差在哪儿了。
以上方程即为水平集演化方程，接下来针对此方程进行数值离散，编码迭代计算即可求解。

2、matlab 实现

% 借鉴网址：https://www.math.hkust.edu.hk/~masyleung/Teaching/CAS/
function phi = GeodesicActiveContour( A, b, tfinal )

fileout = 5;    % number of intermediate output

dtout = tfinal/fileout;
time=0;
[m,n,p]=size(A);
xi=2:m+1;
yi=2:n+1;
phi=zeros(m+2,n+2);
B=zeros(m+2,n+2);
B(xi,yi)=A;

[x,y]=meshgrid(1:m,1:n);
xc=m/2;
yc=n/2;
r=min(m,n)/2*0.8;
phi(xi,yi) = sqrt((x-xc).^2+(y-yc).^2)-r; % 这里应该是做初始化
        
% g(I)
    
C = gaussian(B,2,5); % 做高斯滤波去噪
[gn,ge,gs,gw]=g(C,b);
    
imshow(uint8(A))
hold on
contour(phi(xi,yi),[0 0],'b')

while time < tfinal
    
    % du/dn=0 on the boundary
    
    phi(1,yi)= phi(2,yi);
    phi(end,yi)= phi(end-1,yi);
    phi(xi,1)= phi(xi,2);
    phi(xi,end)= phi(xi,end-1);
    
    Bn = phi(xi,yi+1);
    Be = phi(xi+1,yi);
    Bs = phi(xi,yi-1);
    Bw = phi(xi-1,yi);
    Bc = phi(xi,yi);
    
    [Gn,Ge,Gs,Gw] = NormGrad(phi,1e-5);    
    Gn = gn./Gn;
    Ge = ge./Ge;
    Gs = gs./Gs;
    Gw = gw./Gw;    
        
    % phi^k+1 = phi^k + dt * div[ g(I) grad(phi) / |grad(phi)| ] * |grad(phi)|
    
    temp = Gn.*Bn + Ge.*Be + Gs.*Bs + Gw.*Bw ...
        - (Gn+Ge+Gs+Gw).*Bc;
    
    temp = temp.* sqrt( (Bn-Bs).^2+(Be-Bw).^2 )/2;
    
    dt = 0.45;        
    dt = min(dt,tfinal-time);    
    time=time+dt;
    [time dt tfinal]
    
    phi(xi,yi)=phi(xi,yi)+dt*temp;

    % check for intermediate output
    
    if (floor(time/dtout) ~= floor((time-dt)/dtout) )
        figure
        imshow(uint8(A))
        hold on
        contour(phi(xi,yi),[0 0],'b')
    end
    
end

phi=phi(2:end-1,2:end-1);

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [gn,ge,gs,gw]=g(u,b)

% g(u)= 1/ (1+b |grad(u)|^2)

[m,n]=size(u);
xi=[2:m-1];
yi=[2:n-1];

uxp=u(xi+1,yi);
uxm=u(xi-1,yi);
uyp=u(xi,yi+1);
uym=u(xi,yi-1);
uc=u(xi,yi);

duxp=uxp-uc;
duxm=uc-uxm;
duyp=uyp-uc;
duym=uc-uym;

duxc=(duxp+duxm)/2;
duyc=(duyp+duym)/2;

gn=1./(1+b*(duxc.^2+duyp.^2));
ge=1./(1+b*(duxp.^2+duyc.^2));
gs=1./(1+b*(duxc.^2+duym.^2));
gw=1./(1+b*(duxm.^2+duyc.^2));

return

function [gn,ge,gs,gw]=NormGrad(u,eps)

% NormGrad(u)= |grad(u)|_eps

[m,n]=size(u);
xi=[2:m-1];
yi=[2:n-1];

uxp=u(xi+1,yi);
uxm=u(xi-1,yi);
uyp=u(xi,yi+1);
uym=u(xi,yi-1);
uc=u(xi,yi);

duxp=uxp-uc;
duxm=uc-uxm;
duyp=uyp-uc;
duym=uc-uym;

duxc=(duxp+duxm)/2;
duyc=(duyp+duym)/2;

gn=sqrt(duxc.^2+duyp.^2+eps^2);
ge=sqrt(duxp.^2+duyc.^2+eps^2);
gs=sqrt(duxc.^2+duym.^2+eps^2);
gw=sqrt(duxm.^2+duyc.^2+eps^2);

return

三、CV 模型

在图像中，对象与背景的区别可能会表现在平均灰度值的明显不同，因此 T.Chan 和 L.Vese 提出了以下能量泛函：
$$E(c_1.c_2.C)=\mu {\oint }_{C}ds+{\lambda }_1{\iint }_{{\Omega }_1}(I-c_1{)}^{2}dxdy+{\lambda }_2{\iint }_{{\Omega }_2}(I-c_2{)}^{2}dxdy$$
第二项、第三项分别计算实际图像与假定的“分片图像”之间的平均灰度值差异。下图显示了曲线 $C$ 在不同位置时第二项与第三项的取值，可以看到只有当曲线 $C$ 在对象边缘附近时，这两项同时达到最小值。

此图像来自知乎，也是 CV 模型进行说明时的经典图像。
CV 模型是基于区域的图像分割方法，也称为测地线活动区域（Geodesic activate region, GAR）模型.

CV 模型的梯度下降流

根据变分水平集方法，引入 Heaviside 函数
$$\begin{array}{rl}E(c_1,c_2,C)& =\mu {\iint }_{\Omega }\delta (u)|\nabla u|dxdy\\ & +{\lambda }_1{\iint }_{{\Omega }_1}(I-c_1{)}^{2}H(u)dxdy\\ & +{\lambda }_2{\iint }_{{\Omega }_2}(I-c_2{)}^2(1-H(u))dxdy\end{array}$$
根据 $c_1,c_2$ 最小化上式
$$c_i=\frac{{\iint }_{{\Omega }_i}Idxdy}{{\iint }_{{\Omega }_i}dxdy},i=1,2$$
每次迭代，$c_1,c_2$是已知的，因此 $c_1,c_2$固定时，上式的梯度下降流为
$$\begin{array}{rl}E(u+ϵv)& =\mu {\iint }_{\Omega }\delta (u+ϵv)|\nabla (u+ϵv)|dxdy\\ & +{\lambda }_1{\iint }_{{\Omega }_1}(I-c_1{)}^{2}H(u+ϵv)dxdy\\ & +{\lambda }_2{\iint }_{{\Omega }_2}(I-c_2{)}^2(1-H(u+ϵv))dxdy\end{array}$$
$$\begin{array}{rl}{E}_{ϵ}^{{}^{\prime }}& =\mu {\iint }_{\Omega }\delta (u)\frac{\nabla u+ϵ\nabla v}{|\nabla u+ϵ\nabla v|}\cdot \nabla vdxdy\\ & +{\lambda }_1{\iint }_{{\Omega }_1}(I-c_1{)}^2\delta (u+ϵv)vdxdy\\ & -{\lambda }_2{\iint }_{{\Omega }_2}(I-c_2{)}^2\delta (u+ϵv))vdxdy\end{array}$$
$$\begin{array}{rl}{E}_{ϵ=0}^{{}^{\prime }}& =\mu {\iint }_{\Omega }-\delta (u)div(\frac{\nabla u}{|\nabla u|})\cdot vdxdy\\ & +{\lambda }_1{\iint }_{{\Omega }_1}(I-c_1{)}^2\delta (u)vdxdy\\ & -{\lambda }_2{\iint }_{{\Omega }_2}(I-c_2{)}^2\delta (u))vdxdy\end{array}$$
对应的梯度下降流为
$$\frac{\partial u}{\partial t}={\delta }_{ϵ}(\mu div(\frac{\nabla u}{|\nabla u|})-{\lambda }_1(I-c_1{)}^2+{\lambda }_2(I-c_2{)}^2)$$
最后同 GAC 模型一样，数值离散后迭代求解上式即可。变分水平集介绍即使用参见水平集方法。

参考文献：《图像处理的偏微分方程方法》，王大凯，候榆青等。