偏微分方程与图像处理.

偏微分方程与图像处理(曲线的演化)实验名称： 平面曲线的演化实验内容：1.用水平集方法对曲线进行演化；2.用离散中值滤波方法进行演化。

理论分析：我们已知道：曲线演化方程式（平均曲率运动方程MCM )ck N t∂=∂； 1. 曲线演化水平集方法平面封闭曲线可以表达为一个二维函数u(x,y)的水平（线）集(,,){(,,):(,,)}c L x y t x y t u x y t c ==这样就可将曲线演化问题嵌入到函u(x,y,t)的演化问题。

即转化为水平集演化问题 曲线演化水平集方法的基本方程式如下：||uk u t∂=∇∂其中，||u ∇=()223/2222xx y x y xy yy x xyu u u u u u u k uu -+=+进而推得：22222xx y x y xy yy xx yu u u u u u u u t u u -+∂=∂+；其中x u ,xy u ,xx u 可采用中心差分近似 ()()1,1,1,,1,21,11,11,11,12(,)22(,)(,)4i j i jx i j i j i jxx i j i j i j i j xy u u u i j xu u u u i j x u u u u u i j x +-+-++--+--+-=∆-+=∆+--=∆对于y u ,yy u 有类似的表达式。

x ∆表示相邻几个点。

从而完整的演化公式为：221,,222xx y x y xy yy x n ni ji j x yu u u u u u u u u tu u +-+=+∆+ （1）其中，t ∆为演化步长，在本程序中取为1。

这样就涉及到两个问题： （1）.嵌入函数的选用嵌入函数为—令u(x,y)表示平面上(x,y)点到曲线C 的带有符号的距离(见课本)。

因此研究的曲线总对应于零水平集，这样只要检测过零点条件,1,.0i j i j u u +< 或 ,,1.0i j i j u u +<就可决定曲线C 目前所处的位置，事实上，我们在程序中也是这样做的。

（2）.初始化嵌入函数：计算平面上每一网格点到初始曲线0C 的距离，即到曲线各点距离的最小值，再根据此网格点在0C 的内部或外部赋以正号或符号。

只要前两步准备好，再按式（1）对某一闭合曲线进行演化，该曲线就会逐步光滑，转成椭圆，圆，直至点，然后消失。

（见后续的实验结果）； 2. 离散中值滤波方法根据课本分析，我们知道，曲线演化方程式对应于离散中值滤波。

因此，利用该方法也可以看到曲线随着时间的演化。

具体利用下式实现：{}{}1{,,#1#,,}2xy M x y x b B x X b B =⊗=≥∈∈ 其中，⊗为“逻辑与”运算，xy B 是将结构元素的中心平移至(),x y 的集合。

实验主程序：1. 嵌入函数初始化程序该函数主要完成嵌入函数(,,)u x y t 的初始化，以便进行后续演化工作。

该函数的输出为嵌入函数（距离）和闭合曲线矩阵。

function [dis_u, fish_line]=fish_initialize()clear;[DataRange,datalen]=data_get(); %获得文本文档kk99.c 中的闭合曲线相应坐标信息。

max_row=max(DataRange(:,1)); %max_row=172 max_column=max(DataRange(:,2));%max_column=305%% -------------------------------------------完成闭合曲线的显示 m=max_row+10; n=max_column+10; fish_line=ones(m,n); for i=6:m-5 for j=6:n-5for k=1:datalenif (i==DataRange(k,1)+5)&&(j==DataRange(k,2)+5)==1 fish_line(i,j)=0; end end end endfigure;imshow(fish_line);f=full_fill(fish_line); figure;imshow(f); %填充闭合曲线，以便将曲线内的距离赋为负值；dis_u=distance(f,DataRange,datalen,5); %计算网格上的每一点到闭合曲线的距离% lever_set_display(dis_u); %4，0，-4 水平集的显示2.水平集演化程序该函数对水平集演化3000次，每300次输出一幅0水平集function yan_hua()clear;[dis_u, fish_line]=fish_initialize();u0=dis_u; det=1;[m,n]=size(u0); u1=ones(m,n);newpic=fish_line; figure;imshow(newpic);for num=1:3000for i=3:m-3for j=3:n-3u_y=(u0(i,j+2)-u0(i,j-2))/4;u_x=(u0(i+2,j)-u0(i-2,j))/4;u_yy=(u0(i,j+2)-2*u0(i,j)+u0(i,j-2))/4;u_xx=(u0(i+2,j)-2*u0(i,j)+u0(i-2,j))/4;u_xy=(u0(i+2,j+2)-u0(i-2,j+2)-u0(i+2,j-2)+u0(i-2,j-2))/16;g=(u_xx*u_y*u_y-2*u_x*u_y*u_xy+u_yy*u_x*u_x)/(u_x*u_x+u_y*u_y+eps);u1(i,j)=u0(i,j)+det*g;endendu0=u1;if mod(num,300)==0 %每300次输出一条0水平集for i=2:m-1for j=2:n-1if (u0(i,j)*u0(i+1,j)<0)|(u0(i,j)*u0(i,j+1)<0)==1newpic(i,j)=0;endendendfigure; imshow(newpic);endend3.离散中值滤波该函数对经过填充的闭合曲线（图像）应用离散中值算子进行滤波。

滤波240次，每30次输出一次图片function mid_smooth()clear;c=imread('w_fish.jpg'); %读入填充好的闭合曲线图像b=rgb2gray(c); %将读入图像转换为灰度图像并二值化，1（white）,0(black)[m,n]=size(b);a=im2bw(b);%figure;imshow(a);newpic=ones(m,n); %提取二值化后的图像边缘线，以便嵌套滤波结果边缘线，方便对照for i=1:m-1for j=1:n-1if (a(i,j)+a(i+1,j)==1)|(a(i,j)+a(i,j+1)==1)==1newpic(i,j)=0;endendend% figure; imshow(newpic);for k=1:180 %采用9*9窗口对图像滤波180次，每30次输出一条边缘线c=ones(m,n);for i=5:m-4for j=5:n-4num=0;%----------------------------------------计算9*9窗口中白点的个数ass1=a(i-4,j-4)+a(i-4,j-3)+a(i-4,j-2)+a(i-4,j-1)+a(i-4,j)+a(i-4,j+1)+a(i-4,j+2)+a(i-4,j+3)+a(i-4,j+4); ass2=a(i-3,j-4)+a(i-3,j-3)+a(i-3,j-2)+a(i-3,j-1)+a(i-3,j)+a(i-3,j+1)+a(i-3,j+2)+a(i-3,j+3)+a(i-3,j+4); ass3=a(i-2,j-4)+a(i-2,j-3)+a(i-2,j-2)+a(i-2,j-1)+a(i-2,j)+a(i-2,j+1)+a(i-2,j+2)+a(i-2,j+3)+a(i-2,j+4);ass4=a(i-1,j-4)+a(i-1,j-3)+a(i-1,j-2)+a(i-1,j-1)+a(i-1,j)+a(i-1,j+1)+a(i-1,j+2)+a(i-1,j+3)+a(i-1,j+4); ass5=a(i,j-4)+a(i,j-3)+a(i,j-2)+a(i,j-1)+a(i,j)+a(i,j+1)+a(i,j+2)+a(i,j+3)+a(i,j+4);ass6=a(i+1,j-4)+a(i+1,j-3)+a(i+1,j-2)+a(i+1,j-1)+a(i+1,j)+a(i+1,j+1)+a(i+1,j+2)+a(i+1,j+3)+a(i+1, j+4);ass7=a(i+2,j-4)+a(i+2,j-3)+a(i+2,j-2)+a(i+2,j-1)+a(i+2,j)+a(i+2,j+1)+a(i+2,j+2)+a(i+2,j+3)+a(i+2, j+4);ass8=a(i+3,j-4)+a(i+3,j-3)+a(i+3,j-2)+a(i+3,j-1)+a(i+3,j)+a(i+3,j+1)+a(i+3,j+2)+a(i+3,j+3)+a(i+3, j+4);ass9=a(i+4,j-4)+a(i+4,j-3)+a(i+4,j-2)+a(i+4,j-1)+a(i+4,j)+a(i+4,j+1)+a(i+4,j+2)+a(i+4,j+3)+a(i+4, j+4);num=ass1+ass2+ass3+ass4+ass5+ass6+ass7+ass8+ass9;%----------------------------------------计算9*9窗口中白点的个数% % 滤波过程 if num>=41c(i,j)=1; %曲线回缩 elsec(i,j)=0; %曲线外放 end end end ka=c;%%----------------------------------------每30次提取一次边缘线，并将其套在前面的边缘线中if mod(k,30)==0for i=1:m-1 for j=1:n-1if a(i,j)+a(i+1,j)==1|a(i,j)+a(i,j+1)==1 newpic(i,j)=0; end end end end endfigure; imshow(newpic); %在一幅图newpic 中显示6条滤波曲线实验结果及分析：图1 给定的闭合曲线及填充结果图2 4，0，-4的水平集显示1.用水平集方法演化结果图3 演化3000次，每300次输出一条曲线2.离散中值滤波方法结果图4 滤波180次，每30次输出一次边缘线，第一次为原闭合曲线结果分析：由图3和图4可以看到（1）水平集方法和离散中值滤波方法都可以实现曲线的演化（对后者来讲是图像的轮廓）（2）然而，两者有以下不同：①前者比后者演化推进速度慢得多。