学习图像分割算法,在网上找到的关于主动轮廓模型的实现代码,自己简化总结了一下,在这里和大家分享,欢迎提问
进入正题:
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; end c = 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:row for j = 1:col Eterm(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