偏微分方程与图像处理

1、偏微分方程与图像处理(曲线的演化)实验名称： 平面曲线的演化实验内容： 1.用水平集方法对曲线进行演化； 2.用离散中值滤波方法进行演化。理论分析：我们已知道：曲线演化方程式（平均曲率运动方程MCM)；1. 曲线演化水平集方法平面封闭曲线可以表达为一个二维函数u(x,y)的水平（线）集 这样就可将曲线演化问题嵌入到函u(x,y,t)的演化问题。即转化为水平集演化问题曲线演化水平集方法的基本方程式如下： 其中，; 进而推得：；其中,可采用中心差分近似对于,有类似的表达式。表示相邻几个点。从而完整的演化公式为： （1） 其中，为演化步长，在本程序中取为1。这样就涉及到两个问题：（1）.嵌入函数的选

2、用嵌入函数为令u(x,y)表示平面上(x,y)点到曲线C的带有符号的距离(见课本)。因此研究的曲线总对应于零水平集，这样只要检测过零点条件 或 就可决定曲线C目前所处的位置，事实上，我们在程序中也是这样做的。（2）.初始化嵌入函数：计算平面上每一网格点到初始曲线的距离，即到曲线各点距离的最小值，再根据此网格点在的内部或外部赋以正号或符号。只要前两步准备好，再按式（1）对某一闭合曲线进行演化，该曲线就会逐步光滑，转成椭圆，圆，直至点，然后消失。（见后续的实验结果）；2. 离散中值滤波方法根据课本分析，我们知道，曲线演化方程式对应于离散中值滤波。因此，利用该方法也可以看到曲线随着时间的演化。具体利

3、用下式实现： 其中，为“逻辑与”运算，是将结构元素的中心平移至的集合。实验主程序：1. 嵌入函数初始化程序该函数主要完成嵌入函数的初始化，以便进行后续演化工作。该函数的输出为嵌入函数（距离）和闭合曲线矩阵。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 % -

4、完成闭合曲线的显示 m=max_row+10; n=max_column+10; fish_line=ones(m,n); for i=6:m-5 for j=6:n-5 for k=1:datalen if (i=DataRange(k,1)+5)&&(j=DataRange(k,2)+5)=1 fish_line(i,j)=0; end end end end figure;imshow(fish_line);f=full_fill(fish_line); figure;imshow(f); %填充闭合曲线，以便将曲线内的距离赋为负值；dis_u=distance(f,Da

5、taRange,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:3000 for i=3:m-3 for j=3:n-

6、3 u_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; end end u0=u1;

7、if mod(num,300)=0 %每300次输出一条0水平集 for i=2:m-1 for j=2:n-1 if (u0(i,j)*u0(i+1,j)<0)|(u0(i,j)*u0(i,j+1)<0)=1 newpic(i,j)=0; end end end figure; imshow(newpic); endend3. 离散中值滤波 该函数对经过填充的闭合曲线（图像）应用离散中值算子进行滤波。滤波240次，每30次输出一次图片function mid_smooth()clear;c=imread('w_fish.jpg'); %读入填充好的闭合曲线图像b=

8、rgb2gray(c); %将读入图像转换为灰度图像并二值化，1（white）,0(black)m,n=size(b);a=im2bw(b);%figure;imshow(a); newpic=ones(m,n); %提取二值化后的图像边缘线，以便嵌套滤波结果边缘线，方便对照 for i=1:m-1 for j=1:n-1 if (a(i,j)+a(i+1,j)=1)|(a(i,j)+a(i,j+1)=1)=1 newpic(i,j)=0; end end end% figure; imshow(newpic); for k=1:180 %采用9*9窗口对图像滤波180次，每30次输出一条边缘

9、线 c=ones(m,n); for i=5:m-4 for j=5:n-4 num=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,

10、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

11、,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

12、-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>=41 c(i,j)=1; %曲线回缩 else c(i,j)=0; %曲线外放 end end end k a=c; %-每30次提取一次边缘线，并将其套在前面的边缘线中if mod(k,30)=0 for i=1:m-1 for j=1:n-1 if a(i,j)+a(i+1,j)=1|a(i

13、,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） 然而，两者有以下不同： 前者比后者演化推进速度慢得多。大

14、概前者75次迭代的结果才约等于后者4次的滤波结果（3000/120）。这是因为前者演化时用到的迭代未必能使该点之发生变化；但后者至少会对处于边缘的点平滑一个像素点； 前者程序执行时间远远大于后者，（前者：3000次,约50分钟；后者：120次，约15s）这是由于前者的演化是在迭代的基础上进行的，需要大量的运算，因而费时；而后者只是一些简单的判断语句，花费时间少。子程序附后：function DataRange,datalen=data_get()%获得文本文档kk99.c中的闭合曲线相应坐标信息，返回坐标数组及数据长度filename,pathname = uigetfile( . '

15、*.c', 'All c-Files (*.c)' . '*.dat', 'All tif-Files (*.dat)' . '*.*','All Files (*.*)', . 'Select image'); if isequal(filename,pathname,0,0) % If "Cancel" is selected then return return end File = fullfile(pathname,filename); strImage=str

16、cat(pathname,filename); % pathname+filenamefid = fopen(strImage,'rt'); tline = fgetl(fid); datalen=str2num(tline(3:length(tline); %length(tline)=5 DataRange=zeros(datalen,2); i=0; while 1 i=i+1; tline = fgetl(fid); if ischar(tline), break, end divpos=0; while 1 divpos=divpos+1; if tline(divp

17、os)=' ', break, end end DataRange(i,1)=str2num(tline(1:divpos); %datarange(i,1)=16 DataRange(i,2)=str2num(tline(divpos:length(tline); %datarange(i,2)=155 endfclose(fid);function full=full_fill(fish_line) %填充函数 %由于设计的部分填充函数效果不好，因此将原图像分别旋转0，90，180度，利用部分填%充函数填充，然后取三者的交集即为最终填充效果e=fish_line;m,n=s

18、ize(fish_line);part_1=part_fill(e,m,n); g=rot90(e); ass_2=part_fill(g,n,m);f1=rot90(ass_2); f2=rot90(f1); part_2=rot90(f2); h=rot90(g); ass_3=part_fill(h,m,n); g=rot90(ass_3); part_3=rot90(g); full=ones(m,n);for i=1:m for j=1:n if (part_1(i,j)=0)&&(part_2(i,j)=0)&&(part_3(i,j)=0)=1 full(i,j)=0; end endendfunction f=part_fill(e,m,n) %部分填充函数 f=ones(m,n); for i=2:m for j=2:n if e(i,j)=e(i,j-1) f(i,j)=f(i,j-1); elseif e(i,j)=e(i-1,j) f(i,j)=f(i-1,j); else f(i,j)=0; end end