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

偏微分方程与图像处理

(GAC的水平集方法)

实验二 GAC的水平集方法

一 实验目的

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

二 原理分析

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

?u?t?????gc?u??g??u?gk?u （3.31）

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

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

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

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

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

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

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

uij?min(?1,0)Dx(0)212(?)?max(?1,0)Dxn(0)(?)uij?max(?2,0)Dy(?)uij?min(?2,0)Dy(?)uij ?gijkij[(Dxuij)?(Dyuij)]} (2.1) 其中 ?(?)2?[(max(Dxuij,0))?(min(Dxuij,0))?(max(Dyuij,0))?(min(Dyuij,0))]?2?2?2?2

(2.3) D(0)xuij?ui,j?1?ui,j?12 中心差分 (2.2)

Dxuij?ui,j?1?ui,j 向前单边差分 (2.3) Dxuij?ui,j?ui,j?1 向后单边差分 (2.4)

??三 编程过程

1 准备工作

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

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

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

5）根据C0初始化水平集u0。

????2 迭代运算

1）根据（2.1）式进行迭代运算，由uij(n)计算uij(n?1)。

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

[Ix, Iy]=gradient(b); deltI=abs(sqrt(Ix.^2+Iy.^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( 'd', ( 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． g?exp(?rk)中k值选取过大，导致函数g 的下降速度过于平缓，不能在对象边缘有效的停止。

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，能够有效的进行分割。

本文来源：https://www.bwwdw.com/article/vb5t.html