主動(dòng)輪廓模型-Snake分割算法matlab源碼

1、主動(dòng)輪廓模型Snake分割算法matlab源碼學(xué)習(xí)圖像分割算法，在網(wǎng)上找到的關(guān)于主動(dòng)輪廓模型的實(shí)現(xiàn)代碼，自己簡(jiǎn)化總結(jié)了一下，在這里和大家分享，歡迎提問(wèn)進(jìn)入正題：snake是一種能量最小的曲線，表示為v(s)=(x(s),y(s),s為歸一化的曲線長(zhǎng)度，swO,1。能量函數(shù)由曲線內(nèi)部能量和外部約束力(圖像力)組成，表示為Esnake=JEsnake(v(s)ds=J(Eint(v(s)+Eimage(v(s)ds曲線內(nèi)部能量圖像力內(nèi)部能量分為彈性勢(shì)能和彎曲勢(shì)能兩部分：Eint=(a(s)|vs(s)|人2+B(s)|Vss(s)|人2)/2彈性勢(shì)能彎曲勢(shì)能vs(s)是曲線的一階導(dǎo)數(shù);Vss(s)

2、是曲線的二階導(dǎo)數(shù)。假設(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;圖像力分為三部分，分別驅(qū)使snake趨向于lines(線),edges(邊),termination(終端)Eimage=wlineEline+wedgeEedge+wtermEterm般設(shè)定Eline為圖像強(qiáng)度，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弘丄_護(hù)cmSC/dn綜上，目標(biāo)輪廓的確定就轉(zhuǎn)化為極小化如下的能量泛函的問(wèn)題Esnake=J(a(s)|vs(s)|A2+B(s)|Vss(s)|A2)/2+Eimage)ds求解能量的極小化是一個(gè)典型的變分問(wèn)題，依據(jù)變分法的原理將其轉(zhuǎn)化為歐拉公式，將變分問(wèn)題轉(zhuǎn)化為微分問(wèn)題，進(jìn)而求得極小值%讀入圖像I=imread(sample.tif);I=imread(sample.tif);Igs=im2double(I)

4、;figure,imshow(Igs)%手動(dòng)獲取snake輪廓點(diǎn)x=;y=;c=1;N=100;whilecNxi,yi,button=ginput(1);%精確獲取輪廓點(diǎn)x=x,xi;%將獲取的點(diǎn)存入x,y集合y=y,yi;holdon;plot(xi,yi,ro);if(button=3),%當(dāng)點(diǎn)擊鼠標(biāo)右鍵時(shí)，取點(diǎn)停止break;endc=c+1;end%將第一個(gè)點(diǎn)復(fù)制到最后，構(gòu)成完整的輪廓結(jié)構(gòu)xy=x;y;c=c+1;xy(:,c)=xy(:,1);%對(duì)輪廓線進(jìn)行插值t=1:c;ts=1:0.1:c;xys=spline(t,xy,ts);xs=xys(1,:);%初始取點(diǎn)橫坐標(biāo)ys=x

5、ys(2,:);%初始取點(diǎn)縱坐標(biāo)%查看插值效果holdontemp=plot(x(1),y(1),ro,xs,ys,b.);legend(temp,原點(diǎn)插值點(diǎn));%snake算法主體部分%圖像力線函數(shù)Eline=Igs;%原圖像%圖像力邊函數(shù)gx,gy=gradient(Igs);Eedge=-1*sqrt(gx.*gx+gy.*gy);%梯度圖像%圖像力終點(diǎn)函數(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;%計(jì)算外部力Eext=wl*Eline+we*Eedge+wt*Eterm;%計(jì)算梯度f(wàn)x,fy=gradient(Eext);%

7、計(jì)算五對(duì)角狀矩陣xs=xs;%初始取點(diǎn)橫坐標(biāo)集合轉(zhuǎn)換為列向量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);%計(jì)算矩陣的逆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