主动轮廓模型——Snake分割算法 matlab源码

学习图像分割算法，在网上找到的关于主动轮廓模型的实现代码，自己简化总结了一下，在这里和大家分享，欢迎提问

进入正题：

snake是一种能量最小的曲线，表示为v(s) = (x(s), y(s)), s为归一化的曲线长度，s∈[0, 1]。

能量函数由曲线内部能量和外部约束力（图像力）组成， 表示为 Esnake =  ∫Esnake(v(s))ds

                                                                                = ∫( Eint(v(s))  +  Eimage(v(s)) )ds

                                                                                      曲线内部能量        图像力         

 内部能量分为弹性势能和弯曲势能两部分：

                                                                Eint = (α(s)|vs(s)|^2 + β(s)|Vss(s)|^2）/2

                                                                              弹性势能             弯曲势能

vs(s)是曲线的一阶导数；Vss(s)是曲线的二阶导数。假设Vi = (xi, yi),  i=0,1,…，n-1;

vs(s) = (vi+1 - vi-1)/2;

vss(s) ≈ ((vi+1 -vi) - (vi - vi-1))  = vi+1 - 2vi +vi-1;

所以 Eint = ∑α|vi+1 - vi| + β|vi+1 - 2vi +vi-1|^2 ;

图像力分为三部分，分别驱使snake趋向于lines(线),edges(边),termination(终端)

        Eimage = wlineEline + wedgeEedge + wtermEterm

           一般设定Eline为图像强度，Eedge为亮度的梯度变化；

        Eline = I(x,y);

        Eedge = -|▽I(x,y)|^2;

        C(x,y)为高斯滤波后的图像，θ是(x,y)处的梯度角度；

                    C(x,y) = Gσ(x,y) * I(x,y);

                    tanθ = Cy / Cx;

        规定 n = (cosθ,  sinθ ) , n⊥ = (-sinθ , cosθ)

        则， 终端的能量函数定义为：

综上，

目标轮廓的确定就转化为极小化如下的能量泛函的问题

Esnake = ∫((α(s)|vs(s)|^2 + β(s)|Vss(s)|^2）/2  + Eimage) ds

求解能量的极小化是一个典型的变分问题，依据变分法的原理将其转化为欧拉公式，将变分问题转化为微分问题，进而求得极小值

%读入图像

I = imread('sample.tif');
Igs = im2double(I);
figure,imshow(Igs)

%手动获取snake轮廓点
x=[];y=[];c=1;N=100;
while c<N
    [xi,yi,button] = ginput(1); %%精确获取轮廓点
    x = [x, xi];  %将获取的点存入x,y集合
    y = [y, yi];
    hold on;
    plot(xi,yi,'ro');
    if(button == 3), %当点击鼠标右键时，取点停止
        break; 
    end
    c = c+1;
end


%将第一个点复制到最后，构成完整的轮廓结构
xy = [x;y];
c = c+1;
xy(:,c) = xy(:,1);
%对轮廓线进行插值
t = 1:c;
ts = 1:0.1:c;
xys = spline(t,xy,ts);
xs = xys(1,:); %初始取点横坐标
ys= xys(2,:); %初始取点纵坐标
%查看插值效果
hold on
temp = plot(x(1),y(1),'ro',xs,ys,'b.');
legend(temp, '原点', '插值点');


%%snake算法主体部分
%图像力——线函数
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');
[row, col] = size(Igs);
for i = 1:row
    for j = 1:col
        Eterm(i,j) =(cyy(i,j)*cx(i,j)*cx(i,j) + cxx(i,j)*cy(i,j)*cy(i,j) -2*cxy(i,j)*cx(i,j)*cy(i,j))/(1+cx(i,j)*cx(i,j)+cy(i,j)*cy(i,j)^1.5);
    end
end


wl=0; we=0.4; wt=0;
%计算外部力
Eext = wl*Eline + we*Eedge + wt*Eterm;
%计算梯度
[fx, fy] = gradient(Eext);


%计算五对角状矩阵
xs = xs'; %初始取点横坐标集合转换为列向量
ys = ys';
[m,n] = size(xs);
[mm,nn] = size(fx);


alpha=0.2; beta=0.2; gama=1; kappa=0.1;
b(1)=beta;
b(2)=-(alpha + 4*beta);
b(3)=(2*alpha + 6*beta);%%b(i) 表示v(i)系数，从(i-2)到(i+2)
b(4)=b(2);
b(5)=b(1);


A = b(1)*circshift(eye(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);


%计算矩阵的逆
[L U] = lu(A+gama.*eye(m));
Ainv = inv(U) * inv(L);


% 画图部分
NIter = 1000;
figure
for i = 1:NIter;
    ssx = gama*xs - kappa*interp2(fx,xs,ys);
    ssy = gama*ys - kappa*interp2(fy,xs,ys);
    %计算新的轮廓点位置
    xs = Ainv*ssx;
    ys = Ainv*ssy;
    
    imshow(I)
    hold on;
    plot([xs; xs(1)],[ys; ys(1)], 'r-');
    hold off;
    pause(0.001)
end