偏微分方程与图像处理(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)2i 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.5b=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;endend（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-1for j=2:n-1if (u(i,j)*u(i+1,j)<0)|(u(i,j)*u(i,j+1)<0)==1 newpic(i,j)=0;endendendfigure;imshow(newpic);[gx,gy]=gradient(g);det=0.05;c=3;dx=1;dy=1;display_it=10;for k=1:numif mod(k,5)==0u=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);%计算更新的u1u1=u1 + det.*( du_1 + du_2 + du_3 );u=u1;%每隔10次，显示一幅演化图片,并编号，存入当前路径if mod(k,display_it)==0fprintf( 1, '%d\n', k );disp( 'Displaying segmented image' );newpic=im;for i=1:m-1for j=1:n-1if (u(i,j)*u(i+1,j)<0)|(u(i,j)*u(i,j+1)<0)==1newpic(i,j)=0;endendendimshow( newpic ); drawnow;filename = strcat( m_name, sprintf( '%06d', ( k / d_it ) + 1 ), '.jpg' );imwrite( newpic, filename );endend（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 : Mfor j = 1 : Nfor k=1:Ldist(1,k)=sqrt((i-x(k)).^2+(j-y(k)).^2);endif u(i,j)>0u(i,j)=min(dist);elseu(i,j)=-min(dist);endendend五实验结果分析以下为实验结果：Fig.1 Fig.2 Fig.3 Fig.4Fig.5 Fig.6 Fig.7 Fig.8Fig.9 Fig.10 Fig.11上面的图是每迭代40次选取一副图象，直至演化至稳态的结果（Fig.1为调整大小并滤波后的图）。
