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

基于偏微分方程的图像分割(二)Snake模型 Matlab实现一、Snake模型[3]的数字原理来自[1]1987年Kass等人共同发表了题为“Snake:Activecontourmodels”的论文,首次引进了变分法,提出了运用活动轮廓模型进行图像分割的思想….

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

来自[1]

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

v2-1d1a82c4cddede8c9961590d45266bfb_b.jpg

v2-63a0e16d4d0047e99f9663b0cd6a305f_b.jpg

v2-064eb7329a3044f9b70f26303f36ed5d_b.jpg

v2-e01c160a0256f4f2fefdf5d232b1401a_b.jpg

二、原理解释

来自[2]

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

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

v2-70ccd116f58cadee150cb358d56f346b_b.jpg

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

v2-4a9494ba55b4c82c48a3a4d5399466a2_b.png

也称图像力。(当轮廓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] CSDN 博客 - 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  
      

四、演示效果

reference

1、doc88.com/p-67417855986

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

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

今天的文章基于偏微分方程的图像分割(二)Snake模型 Matlab实现分享到此就结束了,感谢您的阅读,如果确实帮到您,您可以动动手指转发给其他人。

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。
如需转载请保留出处:https://bianchenghao.cn/27595.html

(0)
编程小号编程小号

相关推荐

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注