學習圖像分割算法，在網上找到的關于主動輪廓模型的實現代碼，自己簡化總結了一下，在這里和大家分享，歡迎提問

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

?

?

總結

以上是生活随笔為你收集整理的主动轮廓模型——Snake分割算法（MATLAB）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。