学习图像分割算法，在网上找到的关于主动轮廓模型的实现代码，自己简化总结了一下，在这里和大家分享，欢迎提问
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; endc = 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:rowfor j = 1:colEterm(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参考博客：https://blog.csdn.net/special00/article/details/108768457?utm_medium=distribute.pc_relevant_bbs_down.none-task-blog-baidujs-1.nonecase&depth_1-utm_source=distribute.pc_relevant_bbs_down.none-task-blog-baidujs-1.nonecase
