基于偏微分方程的图像分割 Snake模型 Matlab实现

一、Snake模型[3]的数字原理

来自[1]

1987年Kass等人共同发表了题为“Snake:Active contour models”的论文,首次引进了变分法,提出了运用活动轮廓模型进行图像分割的思想.

二、原理解释

来自[2]

Snake模型首先需要在感兴趣区域的附近给出一条初始曲线，接下来最小化能量泛函，让曲线在图像中发生变形并不断逼近目标轮廓。

Kass等提出的原始Snakes模型由一组控制点：v(s)=[x(s), y(s)] s∈[0, 1] 组成，这些点首尾以直线相连构成轮廓线。其中x(s)和y(s)分别表示每个控制点在图像中的坐标位置。 s 是以傅立叶变换形式描述边界的自变量。在Snakes的控制点上定义能量函数（反映能量与轮廓之间的关系）：

其中第1项称为弹性能量是v的一阶导数的模，第2项称为弯曲能量，是v的二阶导数的模，第3项是外部能量（外部力），在基本Snakes模型中一般只取控制点或连线所在位置的图像局部特征例如梯度：

也称图像力。（当轮廓C靠近目标图像边缘，那么C的灰度的梯度将会增大，那么上式的能量最小，由曲线演变公式知道该点的速度将变为0，也就是停止运动了。这样，C就停在图像的边缘位置了，也就完成了分割。那么这个的前提就是目标在图像中的边缘比较明显了，否则很容易就越过边缘了。）

弹性能量和弯曲能量合称内部能量（内部力），用于控制轮廓线的弹性形变，起到保持轮廓连续性和平滑性的作用。而第三项代表外部能量，也被称为图像能量，表示变形曲线与图像局部特征吻合的情况。内部能量仅仅跟snake的形状有关，而跟图像数据无关。而外部能量仅仅跟图像数据有关。在某一点的α和β的值决定曲线可以在这一点伸展和弯曲的程度。

最终对图像的分割转化为求解能量函数Etotal(v)极小化（最小化轮廓的能量）。在能量函数极小化过程中，弹性能量迅速把轮廓线压缩成一个光滑的圆，弯曲能量驱使轮廓线成为光滑曲线或直线，而图像力则使轮廓线向图像的高梯度位置靠拢。基本Snakes模型就是在这3个力的联合作用下工作的。

因为图像上的点都是离散的，所以我们用来优化能量函数的算法都必须在离散域里定义。所以求解能量函数Etotal(v)极小化是一个典型的变分问题（微分运算中，自变量一般是坐标等变量，因变量是函数；变分运算中，自变量是函数，因变量是函数的函数，即数学上所谓的泛函。对泛函求极值的问题，数学上称之为变分法）。

在离散化条件（数字图像）下，由欧拉方程可知最终问题的答案等价于求解一组差分方程：（欧拉方程是泛函极值条件的微分表达式，求解泛函的欧拉方程，即可得到使泛函取极值的驻函数，将变分问题转化为微分问题。）

记外部力 F = −∇ P， Kass等将上式离散化后，对x(s)和y(s)分别构造两个五对角阵的线性方程组，通过迭代计算进行求解。在实际应用中一般先在物体周围手动点出控制点作为Snakes模型的起始位置，然后对能量函数迭代求解。

三、Matlab实现

% ========================================================================= % Snakes：Active Contour Models % ========================================================================= % By gujinjin 2012/12/10-12 Sunny % 基于KASS等的论文思想 % 参考文献： % [1] KASS etc.Snakes：Active Contour Models % [2] 博客 - Author：乐不思蜀Tone % [3] Ritwik Kumar(Harvard University),D.Kroon（Twente University)的程序 % [4] 《数学建模与数学实验》 % ------ clc;clf;clear all; % ========================================================================= % 获取手动取点坐标 % ========================================================================= % 读取显示图像 %I = imread('Coronary_MPR.jpg'); I = imread('t2.bmp'); % 转化为双精度型 %I = im2double(I); % 若为彩色，转化为灰度 if(size(I,3)==3), I=rgb2gray(I); end figure(1),imshow(I); %--------------------------- % 高斯滤波 %--------------------------- sigma=1.0; % 创建特定形式的二维高斯滤波器H H = 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]=getpts while c<N [xi,yi,button]=ginput(1); % 获取坐标向量 x=[x xi]; y=[y yi]; hold on % text(xi,yi,'o','FontSize',10,'Color','red'); plot(xi,yi,'ro'); % 若为右击，则停止循环 if(button==3), break; end c=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,:); % 样条差值效果 hold on temp=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; [row col] = size(Igs); % 图像力-线函数 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'); for i = 1:row for j= 1:col Eterm(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); end end %figure(3),imshow(Eterm); %figure(4),imshow(abs(Eedge)); % 外部力 Eext = Eimage + Econ Eext = wl*Eline + we*Eedge + wt*Eterm; % 计算梯度 [fx,fy]=gradient(Eext); xs=xs'; ys=ys'; [m n] = size(xs); [mm nn] = 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); % 计算矩阵的逆 [L U] = lu(A + gamma.* eye(m)); Ainv = inv(U) * inv(L); % ========================================================================= % 画图部分 % ========================================================================= %text(10,10,'+','FontName','Time','FontSize',20,'Color','red'); % 迭代计算 figure(3) for i=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); hold on; plot([xs; xs(1)], [ys; ys(1)], 'r-'); hold off; pause(0.001) end

四、演示效果

五、不足

Kass 等提出的基本 Snakes 模型，在没有图像力平衡的条件下，内部力将把所有控制点收缩为一点或一条直线。也就是说，被分割物体必须完全包含在 Snakes 的初始位置之内，否则陷在内部的控制点将无法回到物体的边界（下图所示)

reference

1、http://www.doc88.com/p-6741785598618.html

2、《Matlab图像处理》part1 Snakes：Active Contour Models 主动轮廓模型

3、M. Kass, A. Witkin, and D. Terzopoulos. Snakes Active contour models. Int.