偏微分方程與圖像處理

1、偏微分方程與圖像處理(曲線的演化)實驗名稱： 平面曲線的演化實驗內(nèi)容： 1.用水平集方法對曲線進行演化； 2.用離散中值濾波方法進行演化。理論分析：我們已知道：曲線演化方程式（平均曲率運動方程MCM)；1. 曲線演化水平集方法平面封閉曲線可以表達為一個二維函數(shù)u(x,y)的水平（線）集 這樣就可將曲線演化問題嵌入到函u(x,y,t)的演化問題。即轉(zhuǎn)化為水平集演化問題曲線演化水平集方法的基本方程式如下： 其中，; 進而推得：；其中,可采用中心差分近似對于,有類似的表達式。表示相鄰幾個點。從而完整的演化公式為： （1） 其中，為演化步長，在本程序中取為1。這樣就涉及到兩個問題：（1）.嵌入函數(shù)的選

2、用嵌入函數(shù)為令u(x,y)表示平面上(x,y)點到曲線C的帶有符號的距離(見課本)。因此研究的曲線總對應于零水平集，這樣只要檢測過零點條件 或 就可決定曲線C目前所處的位置，事實上，我們在程序中也是這樣做的。（2）.初始化嵌入函數(shù)：計算平面上每一網(wǎng)格點到初始曲線的距離，即到曲線各點距離的最小值，再根據(jù)此網(wǎng)格點在的內(nèi)部或外部賦以正號或符號。只要前兩步準備好，再按式（1）對某一閉合曲線進行演化，該曲線就會逐步光滑，轉(zhuǎn)成橢圓，圓，直至點，然后消失。（見后續(xù)的實驗結(jié)果）；2. 離散中值濾波方法根據(jù)課本分析，我們知道，曲線演化方程式對應于離散中值濾波。因此，利用該方法也可以看到曲線隨著時間的演化。具體利

3、用下式實現(xiàn)： 其中，為“邏輯與”運算，是將結(jié)構(gòu)元素的中心平移至的集合。實驗主程序：1. 嵌入函數(shù)初始化程序該函數(shù)主要完成嵌入函數(shù)的初始化，以便進行后續(xù)演化工作。該函數(shù)的輸出為嵌入函數(shù)（距離）和閉合曲線矩陣。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); %填充閉合曲線，以便將曲線內(nèi)的距離賦為負值；dis_u=distance(f,Da

5、taRange,datalen,5); %計算網(wǎng)格上的每一點到閉合曲線的距離% lever_set_display(dis_u); %4，0，-4 水平集的顯示 2. 水平集演化程序該函數(shù)對水平集演化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. 離散中值濾波 該函數(shù)對經(jīng)過填充的閉合曲線（圖像）應用離散中值算子進行濾波。濾波240次，每30次輸出一次圖片function mid_smooth()clear;c=imread('w_fish.jpg'); %讀入填充好的閉合曲線圖像b=

8、rgb2gray(c); %將讀入圖像轉(zhuǎn)換為灰度圖像并二值化，1（white）,0(black)m,n=size(b);a=im2bw(b);%figure;imshow(a); newpic=ones(m,n); %提取二值化后的圖像邊緣線，以便嵌套濾波結(jié)果邊緣線，方便對照 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窗口中白點的個數(shù) 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窗口中白點的個數(shù)% % 濾波過程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條濾波曲線實驗結(jié)果及分析：圖1 給定的閉合曲線及填充結(jié)果圖2 4，0，-4的水平集顯示1. 用水平集方法演化結(jié)果 圖3 演化3000次，每300次輸出一條曲線2. 離散中值濾波方法結(jié)果 圖4 濾波180次，每30次輸出一次邊緣線，第一次為原閉合曲線結(jié)果分析：由圖3和圖4可以看到（1） 水平集方法和離散中值濾波方法都可以實現(xiàn)曲線的演化（對后者來講是圖像的輪廓）（2） 然而，兩者有以下不同： 前者比后者演化推進速度慢得多。大

14、概前者75次迭代的結(jié)果才約等于后者4次的濾波結(jié)果（3000/120）。這是因為前者演化時用到的迭代未必能使該點之發(fā)生變化；但后者至少會對處于邊緣的點平滑一個像素點； 前者程序執(zhí)行時間遠遠大于后者，（前者：3000次,約50分鐘；后者：120次，約15s）這是由于前者的演化是在迭代的基礎(chǔ)上進行的，需要大量的運算，因而費時；而后者只是一些簡單的判斷語句，花費時間少。子程序附后：function DataRange,datalen=data_get()%獲得文本文檔kk99.c中的閉合曲線相應坐標信息，返回坐標數(shù)組及數(shù)據(jù)長度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) %填充函數(shù) %由于設(shè)計的部分填充函數(shù)效果不好，因此將原圖像分別旋轉(zhuǎn)0，90，180度，利用部分填%充函數(shù)填充，然后取三者的交集即為最終填充效果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) %部分填充函數(shù) 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