偏微分方程与图像处理

(GAC的水平集方法)

实验二 GAC的水平集方法

一 实验目的

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

二 原理分析

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

ugcugugkut• （3.31）

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

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

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

第二种力来自于g的梯度(1,2)g，它是一种不论当前C的局部是在对象内部或外部，都能将曲线引向边界的“吸引力”。从而gu•总是使曲线向着更接近于边界线的方向运动，最终达到贴近对象边界的稳定状态。

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

采用单边迎风方案，根据（1.76）式的数值方案实现上式：

考虑到 0g,0c

可得：

(1)()(){nnijijijuutgc

()()()()max(1,0)min(1,0)max(2,0)min(2,0)xijxijyijyijDuDuDuDu (0)2(0)212[()()]}nijijxijyijgkDuDu (2.1)

其中

()2222[(max(,0))(min(,0))(max(,0))(min(,0))]xijxijyijyijDuDuDuDu

(2.3)

,1,1(0)2ijijxijuuDu 中心差分 (2.2) ,1,xijijijDuuu 向前单边差分 (2.3)

,,1xijijijDuuu 向后单边差分 (2.4)

三 编程过程

1 准备工作

1）读入图像I，将其转化为灰度图象，重新调整图象的大小为[100，100]。

2）进行预平滑，即对图象进行高斯卷积滤波。

3）计算图像梯度模值I，rI,代入函数exp()grk计算()gI，然后计算g的梯度(1,2)g。

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

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

2 迭代运算

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

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_central).*(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