图像处理偏微分方法实验报告.docx

图像处理偏微分方法实验报告

———————————————————————————————— 作者：
———————————————————————————————— 日期：

14

图像处理偏微分方法
实验报告

偏微分方程在图像处理领域的研究开始于六、七十年代，从最初的去噪角度和图像恢复的角度出发，相应地引入了偏微分方程。但直到九十年代才比较系统地将偏微分方程引入图像处理领域，结合其他一些数学工具如数学形态学和仿射几何等，形成了比较完整的理论体系。由s．Osher等提出的水平集方法具有良好的几何插值性，将灰度插值转化为曲线插值，在图像处理中有很大的影响，在图像修复与去噪、边缘检测、图像匹配、图像识别等方面都取得了相对较好的结果。

一、偏微分方程在图像处理中的应用
偏微分方程在图像处理领域的研究开始于六、七十年代，Gabor，Rudin，Osher等人分别从去噪的角度出发，把偏微分方程引入了图像处理领域。同时，Osher还从图像恢复的角度出发，使用最小化全局变差(Total Vadation)方法和变分法，也使用了相应的偏微分方程。但直到上世纪九十年代，Alvarez，Lions，Morel等才首次系统地将偏微分方程引入了图像处理领域，在理论上有比较大的突破，并结合其他数学工具，如数学形态学、变分法、逼近论、仿射几何等，形成了比较完整的理论体系。
由于理论体系本身的优点，从二维图像(静止图像，如照片)到三维图像(运动图像，如电影)应用的拓展相对容易一些。这也是偏微分方程在图像处理研究领域中的一个优势。当然，时间这一维数有它本身的特点，与图像平面有本质的区别，在应用中需要加以注意。概括地说偏微分方程的图像处理技术具有以下的特点： (1)基于PDE的方法具有良好的数学基础，可以提供深刻的理论结果，并且算法具备良好的稳定性； (2)一些经典的方法如高斯滤波、中值滤波、膨胀和腐蚀在PDE的统一框架下得到全新的解释； (3)这种视角产生了新的方法，它们比经典的方法包容更多的不变性，如保持结构的滤波、线形增强等；(4)PDE是连续的模型，与具体的离散网格无关，并且具有旋转不变性。
偏微分方程在图像处理中的研究主要集中在图像处理的两个基础部分：图像滤波和图像边缘提取。在实际应用方面，它在以下的图像处理的领域，如图像匹配、图像分割、图像重建、图像检索、图像边缘检测等，也取得了令人欣喜的成果，并表现出很大的发展潜力。以上这些方面及其分支覆盖了图像处理的绝大多数领域。例如，Cohen从Total Variation的角度出发结合相应的偏微分方程，在图像分割、医学图像处理中取得了良好的理论和试验结果，并将其研究领域拓展到了对图像非线性形变的分析。同时，对于动态图像的研究和分析也得到了广泛的应用，在动态目标跟踪、运动分析、形变分析方面都取得了令人满意的结果。由于发展迅速，近年来这方面的研究已经成为图像处理中的又一个热点。
二、二值图像及BMO算法
1．二值图像
在灰度图像中的一类特殊的图像，灰度值只由0，1（连续模型）或者0，255（离散模型）两个值组成，称为二值图像，也被称为像形状，因为二值图像的信息实际等同于图像支撑集W中的一个子集X∈W。
任何二值图像都具有1x这样的形式，其中1x表示定义在X上的特征函数，具体定义是

为去除二值图像中的多余信息，使经过处理后可以使形状变得规则，Koenderind和Van Dorn于1986年提出了动态边界（dynamic shape）方法，总体思想是先求出图像1x经过热传导方程滤波器作用后的结果，选择参数为t，得到的解为1x，t 。很显然1x，t不在是二值图像，再采用二值化的方法将其生成二值图像，仍然记为Xt，定义如下

图1 二值化例

图2 N=20,50,150,250,400处理的结果

Matlab源码：

I = imread('dynamic.bmp');
subplot(1,2,1),imshow(I);
I = double(I);
h = [0,1,0;1,-4,1;0,1,0]*0.1;
for N = 1:400
I = I + filter2(h,I);
End
I = im2bw(uint8(I),0.5);
subplot(1,2,2),imshow(I);

从结果可以看出，随着t的增加，图像的细节信息越来越少，但没有实现保留重要信息的目的，因为图像的内容的拓扑结构发生了变化，两个物体变成了一个物体。

2．BMO算法
为解决上面的问题，介绍Bence-Merrman-Osher算法（BMO），整个过程是一个迭代过程，描述如下：
（1）求1x，t，选择参数t足够小至于作用后不影响形状的拓扑结构；
（2）定义X1={ x | 1x,t(x)≥1/2 }，得到二值图像 1X1 ；
（3）令X = X1，再回到（1）。
随着迭代次数的增加，细节信息越来越少，同时还保留了形状的拓扑结构。

Matlab源码：

Iy=imread('bin_3.jpg');
figure(1);imshow(Iy);xlabel(原图');
[w,h]=size(Iy);
I=[255*ones(w,1),Iy,255*ones(w,1)];
I=[255*ones(1,h+2);I;255*ones(1,h+2)];
I=double(I);
r=0.1;
n=30;
N=[2,5,10,20,30];
for p=1:length(N)
for m=1:n
for k=1:N(p)
D=255*ones(w+2,h+2);
for i=2:w
for j=2:h
D(i,j)=-4*I(i,j)+I(i-1,j)+I(i,j+1)+I(i+1,j)+I(i,j-1);
end
end
I=I+r*D;
end
I=im2bw(uint8(I),0.5)*255;
end
figure(N(p));imshow(I);xlabel(['迭代',int2str(N(p)),'次的结果']);
end

三、水平集与Extrema killer算法
1．水平集方程
为了使演化曲线成为引进的曲面的零-水平集，定义曲面方程，
其中d为点(x．Y)到曲线C(s，t)的距离，正负表示点(x，y)在曲线的内部和外部，由此方程可以满足

表示在任何时刻的演化曲线都可看作运动曲面的零-水平集。我们感兴趣的正是零-水平集的演化情况，对上式关于t微分得

再关于S微分

于是

综合以上，可以得到
在已知时
上述方程被称为Hamilton．Jacobi方程。
图2可以帮助我们理解水平集的思想：假定圆曲线以依赖曲率的速度演化，初始时刻的曲线如图2(a)所示，嵌入在曲面的情形如图2(b)中；过了一个时刻，曲线演化为2(c)，此时嵌入在演化曲面的情形如2(d) ，从中可以看到曲面在不断演化膨胀。

图3 演化的水平集及其嵌入的曲面

2．形态学算子
形态学算子属于数学形态学研究的范畴，其历史可追溯到19世纪的Euler，Steiner，Crofton以及本世纪初Minlowski的论述中。数学形态学是一门新兴学科，1964年，法国的G．Mathem和J．Serra在积分几何的基础上首次创立了这门学科。1982年，J．Serra的专著《图像分析与数学形态学》问世后，数学形态学才在图像处理、模式识别和计算机视觉等领域引起广泛的重视和应用，这些应用反过来又促进它的进一步发展。数学形态学以图像的形态特征为研究对象，它的主要内容是设计一整套概念、变换和算法，用来描述图像的基本特征和基本结构，也就是描述图像中元素与元素、部分与部分间的关系。
最初，由Maheron和Serra提出的数学形态学研究以二值图像为对象，称为二值形态学。此后，Serra和Sternberg等借助于伞理论，把二值形态算子推广到灰度图像，因而使灰度形态学的理论和应用研究也得到很大的发展，已经成为数字图像信号处理和计算机视觉领域中的一种有效方法。而今，数学形态学已经成为图像处理理论的一个重要方面。广泛地应用到图像处理的增强、分割、边缘检测、结构分析、形态分析等众多领域之中。
把同时满足单调性和对比不变的图像变换被称为形态学算子。形态学算子是分析几何形状和结构的数学方法，是建立在集合代数基础上，用集合论方法定量描述几何结构的科学。
Extrema killer算子是一个基于水平集面积的图像光滑算法，这个算法的作用是去除图像中的“峰”，即孤立的水平集。图像中的峰通常是由脉冲噪声引进的，因此该算法在对椒盐噪声的去噪中表现出非常显著的效果。

图4 Vincent, Serra（1993）Extrema killer算法的结果

其算法如下：
（1）给定一个阈值t；
（2）记图像u在灰度λ的水平集为，它由若干个连通区域组成= ；
（3）去除面积比阈值t小的区域，将会改变图像的水平集，则新的水平集为 = ，其中meas是lebesgue测度；
（4）图像变换的结果由下式得到
。

下面是利用Extrema killer算子的一个平滑实例。

Matlab源码：

%Emtra killer
tic
Iy=imread('lena.bmp');
subplot(131);imshow(uint8(Iy));xlabel('原图');
I=imnoise(Iy,'salt & pepper',0.1);
[w,h]=size(I);
I=double(I);
Ij=I;
Ix=[];
n=10;
subplot(132);imshow(uint8(I));xlabel('加噪后的图像');
for lamda=1:256
for i=1:w
for j=1:h
if I(i,j)>=lamda-1
I(i,j)=1;
else I(i,j)=0;
end
end
end
for m=1:2
[L,num]=bwlabel(I,4);
for k=1:num
[r,c]=find(L==k);
rr=length(r);
if rr