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

1、主动轮廓模型Snake分割算法matlab源码学习图像分割算法，在网上找到的关于主动轮廓模型的实现代码，自己简化总结了一下，在这里和大家分享，欢迎提问进入正题：snake是一种能量最小的曲线，表示为v(s)=(x(s),y(s),s为归一化的曲线长度，swO,1。能量函数由曲线内部能量和外部约束力(图像力)组成，表示为Esnake=JEsnake(v(s)ds=J(Eint(v(s)+Eimage(v(s)ds曲线内部能量图像力内部能量分为弹性势能和弯曲势能两部分：Eint=(a(s)|vs(s)|人2+B(s)|Vss(s)|人2)/2弹性势能弯曲势能vs(s)是曲线的一阶导数;Vss(s)

2、是曲线的二阶导数。假设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=艺a|vi+1-vi|+B|vi+1-2vi+vi-1|A2;图像力分为三部分，分别驱使snake趋向于lines(线),edges(边),termination(终端)Eimage=wlineEline+wedgeEedge+wtermEterm般设定Eline为图像强度，Eedge为亮度的梯度变化；Eline=I(x,y);Eedge=-l(x,y)|A2;C(x,y)为高斯滤波后的图像，8是(

3、x,y)处的梯度角度；C(x,y)=Go(x,y)*I(x,y);tan8=Cy/Cx;规定n=(cos8,sin8),门丄=(-sin8,cos8)则，终端的能量函数定义为：E=im弘丄_护cmSC/dn综上，目标轮廓的确定就转化为极小化如下的能量泛函的问题Esnake=J(a(s)|vs(s)|A2+B(s)|Vss(s)|A2)/2+Eimage)ds求解能量的极小化是一个典型的变分问题，依据变分法的原理将其转化为欧拉公式，将变分问题转化为微分问题，进而求得极小值%读入图像I=imread(sample.tif);I=imread(sample.tif);Igs=im2double(I)

4、;figure,imshow(Igs)%手动获取snake轮廓点x=;y=;c=1;N=100;whilecNxi,yi,button=ginput(1);%精确获取轮廓点x=x,xi;%将获取的点存入x,y集合y=y,yi;holdon;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=x

5、ys(2,:);%初始取点纵坐标%查看插值效果holdontemp=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);%梯度图像%图像力终点函数ml=-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)

6、;cyy=conv2(Igs,m4,same);cxy=conv2(Igs,m5,same);row,col=size(Igs);fori=1:rowforj=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)U5);endendwl=0;we=0.4;wt=0;%计算外部力Eext=wl*Eline+we*Eedge+wt*Eterm;%计算梯度fx,fy=gradient(Eext);%

7、计算五对角状矩阵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),O);A=A+b*circshift(eye(m),-1);A=A+b(5)*circshift(eye(m),-2);%计算矩阵的逆LU=lu(A+gama.*eye(m);Ainv=inv(U)*inv(L);%画图部分NIter=1000;figurefori=1:NIter;ssx=gama*xs-kappa*interp2(fx,xs,ys);ssy=gama*ys-kappa