snake算法

获取手动取点坐标%读取显示图像%I=imread('Coronary_MPR.jpg');I=imread('test1.png');%转化为双精度型%I=im2double(I);%若为彩色，转化为灰度%if(size(I,3)==3),I=rgb2gray(I);endfigure(1),imshow(I);%---------------------------%高斯滤波%---------------------------sigma=1;%创建特定形式的二维高斯滤波器HH=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]=getptswhilecN[xi,yi,button]=ginput(1);%获取坐标向量x=[xxi];y=[yyi];holdon%text(xi,yi,'o','FontSize',10,'Color','red');plot(xi,yi,'ro');%若为右击，则停止循环if(button==3),break;endc=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);xs=xys(1,:);ys=xys(2,:);%样条差值效果holdontemp=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;[rowcol]=size(Igs);%图像力-线函数Eline=Igs;%图像力-边函数[gx,gy]=gradient(Igs);Eedge=-1*sqrt((gx.*gx+gy.*gy));%图像力-终点函数%卷积是为了求解偏导数，而离散点的偏导即差分求解m1=[-11];m2=[-1;1];m3=[1-21];m4=[1;-2;1];m5=[1-1;-11];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');fori=1:rowforj=1:colEterm(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);endend%figure(3),imshow(Eterm);%figure(4),imshow(abs(Eedge));%外部力Eext=Eimage+EconEext=wl*Eline+we*Eedge+wt*Eterm;%计算梯度[fx,fy]=gradient(Eext);xs=xs';ys=ys';[mn]=size(xs);[mmnn]=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);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);%计算矩阵的逆[LU]=lu(A+gamma.*eye(m));Ainv=inv(U)*inv(L);%=========================================================================%画图部分%=========================================================================%text(10,10,'+','FontName','Time','FontSize',20,'Color','red');%迭代计算figure(3)fori=1:NIter;ssx=gamma*xs-kappa*interp2(fx,xs,ys);ssy=gamma*ys-kappa*interp2(fy,xs,ys);%计算snake的新位置xs=Ainv*ssx;ys=Ainv*ssy;%显示snake的新位置imshow(I);holdon;plot([xs;xs(1)],[ys;ys(1)],'r-');holdoff;pause(0.001)end