主動輪廓模型-Snake分割算法matlab源碼

1、主動輪廓模型Snake分割算法matlab源碼學習圖像分割算法，在網(wǎng)上找到的關于主動輪廓模型的實現(xiàn)代碼，自己簡化總結了一下，在這里和大家分享，歡迎提問進入正題：snake是一種能量最小的曲線，表示為v(s)=(x(s),y(s),s為歸一化的曲線長度，swO,1。能量函數(shù)由曲線內部能量和外部約束力(圖像力)組成，表示為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)是曲線的一階導數(shù);Vss(s)

2、是曲線的二階導數(shù)。假設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;規(guī)定n=(cos8,sin8),門丄=(-sin8,cos8)則，終端的能量函數(shù)定義為：E=im弘丄_護cmSC/dn綜上，目標輪廓的確定就轉化為極小化如下的能量泛函的問題Esnake=J(a(s)|vs(s)|A2+B(s)|Vss(s)|A2)/2+Eimage)ds求解能量的極小化是一個典型的變分問題，依據(jù)變分法的原理將其轉化為歐拉公式，將變分問題轉化為微分問題，進而求得極小值%讀入圖像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算法主體部分%圖像力線函數(shù)Eline=Igs;%原圖像%圖像力邊函數(shù)gx,gy=gradient(Igs);Eedge=-1*sqrt(gx.*gx+gy.*gy);%梯度圖像%圖像力終點函數(shù)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)系數(shù)，從(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