基于水平集的gac模型的图像分割报告

偏微分方程与图像处理(GAC的水平集方法)

实验二 GAC 的水平集方法

一 实验目的

采用GAC 模型的水平集方法检测图像中对象的轮廓，以便有效地进行分割。

二 原理分析

推广GAC 模型的水平集方法对应的PDE 为：

u gc u g u gk u t

?=?+???+??

（3.31）

按照上式，曲线运动将受两种“力”的支配，第一种力来自于曲率几何形变—曲率运动（gc u gk u ?+?），不过它的强弱还要受到因子()g I ?的影响。

I ?为图象I(x,y)的梯度模值，函数g (r) 是可以是任何具有单调减性的函数。

因为图象梯度模值I ?在图象的边缘附近有较大值，从而使g(I ?)取极小的值，故在图象边缘附近，该作用力将会变的很小，因此有时将边缘函数()g I ?称之为边缘停止函数。常数c 的作用是加速曲线向内部收缩。

第二种力来自于g 的梯度(1,2)g αα?=

，它是一种不论当前C 的局部是在对象内部或外部，都能将曲线引向边界的“吸引力”。从而g u ??? 总是使曲线向着更接近于边界线的方向运

动，最终达到贴近对象边界的稳定状态。

由于这两种作用使曲线演化可最终达到紧靠轮廓这一稳定状态而不再继续演化。

采用单边迎风方案，根据（1.76）式的数值方案实现上式： 考虑到 0g >,0c > 可得： (1)

()

()

{n n ij

ij

ij u u t g c +-=+??

()

()

()

()

max(1,0)min(1,0)max(2,0)min(2,0)x

ij x

ij y

ij y

ij

D u D u D u D u αααα-+-+++++ (0)

2

(0)

212

[()()]}n

ij ij x ij y ij g k D u D u ++ (2.1) 其中 ()

2222

[(max(,0))(min(,0))(max(,0))(min(,0))]

x ij x ij y ij y ij D u D u D u D u -+-+-?

=+++

(2.3) ,1,1

(0)

2

i j i j x

ij u u D

u +--=

中心差分 (2.2)

,1,x ij i j i j D u u u +

+=- 向前单边差分 (2.3) ,,1x ij i j i j D u u u -

-=- 向后单边差分 (2.4)

三 编程过程

1 准备工作

1）读入图像I ，将其转化为灰度图象，重新调整图象的大小为[100，100]。 2）进行预平滑，即对图象进行高斯卷积滤波。

3）计算图像梯度模值I ?，r I =?,代入函数exp()g r k =-计算()g I ?，然后计算g 的梯度(1,2)g αα?=

。

4）选取初始封闭曲线0C ，使其从外部全部包住对象，简单起见这里将0C 选取为以图像中心为圆心的封闭圆。

5）根据0C 初始化水平集0u 。

2 迭代运算

1）根据（2.1）式进行迭代运算，由()

n ij

u 计算(1)

n ij

u +。

2）每迭代5次，进行一次重新初始化，以免累计误差。具体的方法是根据当前演化得到的u 检测零水平集则为当前C ，根据当前C 重新初始化水平集u 。 3）每10次观察一次零水平集，当演化曲线迭代400次，则停止迭代。

3 参数取值：c=3~4, t ?=0.05~0.1, N=3~5, k=1~2.

四 程序

1、主程序：

% 用GAC 水平集演化方法检测图像轮廓 close all; clear all;

a=imread('3.bmp'); % figure;imshow(a); a=rgb2gray(a);

im=imresize(a,[100 100]); % figure;imshow(im); imsize=size(im); im_1=double(im);

% 对图像进行高斯滤波 sigma=2;

gauss_filter =fspecial('gaussian'，3，sigma); %默认值3*3,SIGMA=0.5

b=imfilter(im_1, gauss_filter ,'conv');

% 计算图像梯度[Ix,Iy]和梯度模值deltI

[I x, I y]=gradient(b);

deltI=abs(sqrt(I x.^2+I y.^2));

k=2;

g=exp(-deltI./k);

[gx,gy]=gradient(g);

% 初始化圆，定义中心和半径

center=[floor(size(im)/2)];

radius = min(center)-8;

u = init_u( imsize, center, radius);

% 调用迭代函数

filename = '3.bmp';

m_name = filename( 1 : strfind( filename, '.' ) - 1 );

num=400;

u_new=die_dai(im,u,g,num,m_name);

2、主要子程序

（1）初始化水平集函数：

function u = init_u( imsize, center, radius )

% 初始化水平集

m = imsize( 1 ); n = imsize( 2 );

u = zeros( imsize );

for i = 1 : m;

for j = 1 : n;

distance = sqrt( sum( ( center - [ i, j ] ).^2 ) );

u( i, j ) = distance - radius;

end

end

（2）迭代计算函数：

function u=die_dai(im,u,g,num,m_name)

[m,n]=size(u);

[gx,gy]=gradient(g);

u1=u;

newpic=im;

for i=2:m-1

for j=2:n-1

if (u(i,j)*u(i+1,j)<0)|(u(i,j)*u(i,j+1)<0)==1 newpic(i,j)=0;

end

end

end

figure;imshow(newpic);

[gx,gy]=gradient(g);

det=0.05;c=3;dx=1;dy=1;display_it=10;

for k=1:num

if mod(k,5)==0

u=re_init_u( u );

end

%x和y方向的前向差分和后向差分

diff_x_backward=( u - circshift( u, [ 0, 1 ] ) );

diff_x_forward=( circshift( u, [ 0, -1 ] ) - u );

diff_y_backward=( u - circshift( u, [ 1, 0 ] ) );

diff_y_forward=( circshift( u, [ -1, 0 ] ) - u );

du_1=g.*c.* ( (max( diff_x_forward,0 )).^2 + (min( diff_x_backward,0 )).^2

+(max( diff_y_forward,0 )).^2 + (min( diff_y_backward,0 )).^2);

%计算更新u的第二部分

du_2=max(gx,0) .* diff_x_backward + min(gx,0) .* diff_x_forward + max(gy,0) .*

diff_y_backward + min(gy,0) .* diff_y_forward ;

%--------

%计算更新u的第三部分

%中心差分

diff_y_central=( circshift( u, [ 0, -1 ] ) - circshift( u, [ 0, 1 ] ) ) / 2;

diff_x_central=( circshift( u, [ -1, 0 ] ) - circshift( u, [ 1, 0 ] ) ) / 2;

diff_yy_central=( circshift( u, [ 0, -1 ] ) - 2*u + circshift( u, [ 0, 1 ] ) )/1; diff_xx_central=( circshift( u, [ -1, 0 ] ) - 2*u + circshift( u, [ 1, 0 ] ) )/1; diff_xy_central=(circshift( u, [ -1, -1 ] ) + circshift( u, [ 1, 1] ) -circshift( u, [ -1, 1 ] )-circshift( u, [ 1, -1 ] ) )/4;

du_3=g.*( (diff_xx_central).*((diff_y_central).^2)-(2*(diff_x_central).*(diff_y_cen tral).*(diff_xy_central))+((diff_yy_central).*((diff_x_central).^2)))./((diff_ x_central).^2 +(diff_y_central).^2 + eps);

%计算更新的u1

u1=u1 + det.*( du_1 + du_2 + du_3 );

u=u1;

%每隔10次，显示一幅演化图片,并编号，存入当前路径

if mod(k,display_it)==0

fprintf( 1, '%d\n', k );

disp( 'Displaying segmented image' );

newpic=im;

for i=1:m-1

for j=1:n-1

if (u(i,j)*u(i+1,j)<0)|(u(i,j)*u(i,j+1)<0)==1

newpic(i,j)=0;

end

end

end

imshow( newpic ); drawnow;

filename = strcat( m_name, sprintf( '%06d', ( k / d_it ) + 1 ), '.jpg' );

imwrite( newpic, filename );

end

end

（3）重新初始化水平集函数：

function u = re_init_u( u )

[M,N] = size (u);

[x,y] = find( isfront( u ) ); %找出边缘位置,放在x,y矩阵内；

L = length ( x ); %x,y是对应的坐标,两者成对出现；所以只需求x长度；for i = 1 : M

for j = 1 : N

for k=1:L

dist(1,k)=sqrt((i-x(k)).^2+(j-y(k)).^2);

end

if u(i,j)>0

u(i,j)=min(dist);

else

u(i,j)=-min(dist);

end

end

end

五实验结果分析

以下为实验结果：

Fig.1 Fig.2 Fig.3 Fig.4

Fig.5 Fig.6 Fig.7 Fig.8

Fig.9 Fig.10 Fig.11

上面的图是每迭代40次选取一副图象，直至演化至稳态的结果（Fig.1为调整大小并滤波后的图）。

注意：

1．曲线最终演化基本稳定在对象边界。

2．由于c 选择过小，使曲线在对象凹陷部分（k<0）继续运动的作用力（gc u gk u

?+?）不足，但c过大又会使曲线穿透边界，因此c的选择要合适。

3．重新初始化间隔不要过大，会积累误差。

4．exp()

=-中k值选取过大，导致函数g 的下降速度过于平缓，不能在对象边缘

g r k

有效的停止。

5．时间步长delt_t选取过大，就会影响效果。

对不同的参数进行验证，实验结果表明各参数较为合理的值的范围：

c： 5 ~ 10

delt_t：0.05 ~ 0.1

re_init： 5 ~ 20

k：0.001 ~ 0.005

六实验小结：

通过本次实验，可以看到GAC（测地线活动轮廓模型）在图象分割中的理论运用。其数学思想是通过最小化一个能量泛函来实现的。这与前面所说的演化有些类似。只是在选取时采用什么样的函数来做不同。对于图象的分割，就是要找到边缘部分（即梯度变化在局部有最大值的地方）。所以用g(r)函数来做，使得在图像边缘处曲线停止演化，得到相应的水平集方法的PDE，能够有效的进行分割。