实 验 报 告
实验课程名称 数字图像处理
实验项目名称 图像复原
年 级 2010
专 业 光信息科学与技术
学生姓名 XXX
学 号 XXXXX
理 学 院
实验时间: 2013 年 5 月 2 日
学院:理学院 专业:光信息科学与技术 班级:光信XXX
数字图像处理作业
——图像恢复
摘要
数字图像恢复是数字图像处理的一个基本的和重要的课题,它是后期图像处理(分析和理解)的前提。图像在摄取、传输、储存的过程中不可避免地引起图像质量的下降(图像退化),图像恢复就是试图利用退化过程的先验知识使已退化的图像恢复本来面貌,即根据退化的原因,分析引起退化的环境因素,建立相应的数学模型,并沿着使图像降质的逆过程恢复图像。本文首先对测试图像进行模糊及加噪处理,然后用不同的图像恢复方法,如维纳滤波恢复、约束最小二乘滤波进行图像恢复,并比较它们的处理效果。发现维纳滤波较约束最小二乘法滤波效果要好,这是因为前者利用了原图像的统计信息,采用了真实的PSF函数来恢复。无论何种算法,它们都要依据获取的相关信息才能有效地实施,算法利用的信息越多,信息的准确性越高,复原图像的质量也就越高。
实验原理:
图像复原处理是建立在图像退化的数学模型基础上的,这个退化数学模型能够反映图像退化的原因。图像的退化过程可以理解为施加于原图像上的运算和噪声两者联合作用的结果,图像退化模型如图1所示,可以表示为:
g ( x, y) = H [ f ( x, y)] + n( x, y) = f ( x, y) *h( x, y) + n( x, y) (1)
图1 图像退化模型
(1)在测试图像上产生高斯噪声lena图-需能指定均值和方差;并用滤波器(自选)恢复图像;
实验原理:
噪声是最常见的退化因素之一,也是图像恢复中重点研究的内容,图像中的噪声可定义为图像中不希望有的部分。噪声是一种随机过程,它的波形和瞬时振幅以及相位都随时间无规则变化,因此无法精确测量,所以不能当做具体的处理对象,而只能用概率统计的理论和方法进行分析和处理。本文中研究高斯噪声对图像的影响及其去噪过程。
①高斯噪声的产生:
所谓高斯噪声是指它的概率密度函数服从高斯分布(即正态分布)的一类噪声。一个高斯随机变量z的PDF可表示为:
P(z)= (2)
其中z代表灰度,u是z的均值,是z的标准差。高斯噪声的灰度值多集中在均值附近。
图2 高斯函数
可以通过不同的算法用matlab来产生高斯噪声。
②高斯噪声对信号的影响
噪声影响图像处理的输入、采集、处理的各个环节以及输出结果的全过程,在图像中加高斯噪声通常会使图像变得模糊并且会出现细小的斑点,使图像变得不清晰。
③去除高斯噪声的一些方法
去除高斯噪声的方法有直方图变换,低通滤波,高通滤波,逆滤波,维纳滤波,中值滤波等。本文应用高斯平滑滤波进行去噪处理。
处理结果如下图:
维纳滤波对高斯白噪声的图像滤波效果较好,具有比较好的选择性,可以更好地保存图像的边缘和高频细节信息。所以,维纳滤波在大多数情况下都可以获得满意的结果,尤其对含有高斯噪声的图像。
(2)推导维纳滤波器并实现下边要求;
实验原理:
维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理,其目标是寻找一个滤波器,使得复原后图像 与原始图像 的均方误差最小:
因此维纳滤波器又称为最小均方误差滤波器。
在频率中用下式表达:
其中,是退化图像的傅立叶变换,是退化函数。,是的复共轭。,为噪声的功率谱。,为未退化图像的功率谱。
维纳滤波器的推导:
合理假设要估计的信号f(x,y)为0均值平稳随机过程,噪声为0均值平稳随机过程。
g ( x, y) = h( x, y) * f ( x, y) + ( x, y)
根据已退化图像g(x,y),利用线性估计器来估计原始图像。
先证: g ( x, y) = f ( x, y) + ( x, y),
最小化均方误差:
=,
e(x,y)=f(x,y)? f?(x,y), f?(x,y)=w(x,y)*g(x,y)
最小化均方误差可以由正交原则求解,即最优求解的误差与观测的信号值不相关。
得: E[e(m, n ) g * (x, y)]= E[f(m,n)? f?(m,n)] g*(x, y)]=0
进一步分解得:
E [f(m,n)g * (x, y)]=E[f?(m,n)g * (x, y)]
? = E {w(k1, k2 )[g (m?k1, n?k2)g*(x,y)]}
= E {w(k1, k2 )(Rg (m?k1-x, n?k2-y)}
=w(m-x,n-y)* Rg (m-x, n-y)
记: Rfg ( m-x, n-y)= E[f(m,n)g * (x, y)]
因此, Rfg ( m, n)= w(m,n)* Rg (m, n)
上式作Flouier变换得:
Pfg ( w1, w2)= W(w1, w2)* Pg (w1, w2)
即: W(w1, w2)= Pfg ( w1, w2)/ Pg (w1, w2). (3)
又有: Rfg ( m-x, n-y)= E[f(m,n)g * (x, y)]
= E [f(m,n).[h (k1, k2)f(x-k1,y-k2)+(x, y)]*
= h *(k1, k2) Rf (m-x+ k1, n-y +k2)
= h *(-(m-x),-(n-y))* Rf (m-x, n-y)
因此, Rfg ( m, n)= h *(-m,-n)* Rf (m, n)
对上式作Flouier变换得:
Pfg ( w1, w2)= H*( w1, w2)Pf (w1, w2) (4)
同时, Rg (m-x, n-y)= E[g(m,n)g * (x, y)]
= E [h (k1, k2) f(m-k1,n-k2)+ (m,n)]
*[ h *( 1,2) f*(x-1,y-2)+ *(x, y)]
=h (k1, k2) h *( 1,2) E[f(m-k1,n-k2)
f*(x-1,y-2)]+ E{(m,n) *(x, y)}
得: Rg (m-x, n-y)
= h(k1,k2)h*(1,2) Rf (m-k1-x+1,n-k2-y+2) + R (m-x,n-y)
= Rf (m-x,n-y)* h(m-x,n-y)* h* (-(m-x),-(n-y))+ R (m-x,n-y)
所以,
Rg (m, n)= Rf (m,n)* h(m,n)* h*(-m,-n)+ R (m-x,n-y)
上式作Flouier变换得:
Pg (w1,w2)= Pf (w1, w2)*H(w1, w2)*H*(-w1,-w2)+P(w1, w2)
= Pf (w1, w2)+ P(w1, w2) (5)
综合以上式(3)(4)(5)可得:
推导完毕。
(a) 实现模糊滤波器如方程Eq.(5.6-11).
在图像复原中,有3种主要的估计退化函数的方法:(1)图像观察估计法,(2) 试验估计法,(3)模型估计法。本文实现从基本原理开始推导数学模型的模型估计法。其SFR为:
根据上式设计滤波器,实现的程序代码见附录。
(b) 模糊lena图像:45度方向,T=1;
根据(a)所设计的滤波器处理图像,效果如下图示:
(c) 在模糊的lena图像中增加高斯噪声,均值=0,方差=10 pixels 以产生模糊图像;
处理结果如下图示(程序代码见附录):
(d) 分别利用方程Eq. (5.8-2)和(5.9-4),恢复图像;
考虑到设计的复杂性以及时间的有限性,本次实验的图像复原利用MATLAB数字图像处理工具箱来完成。
①维纳滤波复原
程序中仿真了一个运动模糊PSF,对原图像进行模糊操作,并指定运动位移为50个像素,运动角度为45度。利用函数deconvwnr进行维纳滤波复原图像,调用格式为:
Wnr=deconvwnr(blurredNoisy,PSF);
MATLAB代码见附录,其处理结果如下图示:
可见,直接用无信噪比NSR参数的函数deconvwnr进行复原的滤波效果很差,因此引人信噪比NSR作为噪声参数进行图像复原,即:
noise=imnoise(zeros(size(I)),’gaussian’,0,0.01);
NSR=sum(noise(:).^2)/sum(im2double(I(:)).^2);
Wnr=deconvwnr(blurredNoisy,PSF,NSR);
处理结果如图所示:
②最小二乘滤波复原
维纳滤波建立在最小化统计准则的基础上,在平均意义上它是最优的。本节所提供的算法具有显著的特点对于处理的每一幅图像它都能产生最优结果。根据卷积定义,有
H对噪声敏感,为了减少噪声敏感性问题,以平滑措施的最佳复原为基础,如一幅图像的二阶导数(“拉普拉斯变换”算子)。但是,复原必须被所用到的参量约束。因此,定义一个最小的准则函数C如下:
(6)
约束为:
(7)
这里,是欧几里得矢量范数,是未退化图像的估计值。▽2 为拉普拉斯算子。
这个最佳化问题的频域解决方法由下面的表达式给出:
(8)
这里,γ是一个参数,它必须被调整以使式(7)满足条件,P(u,v)是函数p(x,y)的傅里叶变换。
(9)
p(x,y)和所有其他相关的空间域函数,在用式(8)计算它们的傅里叶变换之前要用零进行适当的延拓,注意,当γ为0时,式(8)变为逆滤波。
根据式(8)设计约束最小二乘滤波器,处理结果如下图:
结果分析:
实验中发现,采用维纳滤波恢复可以取得比较好的效果,这个算法可以使估计的点扩散函数值PSF更加接近它的真实值。在我们知道模糊图像的点扩展函数的情况下,可以调用常规的图像复原算法;而现实里还会遇见不知道点扩展函数的情况,这个时候我们就可以利用盲卷积复原算法。它是利用原始图像模糊,同时进行清晰图像的恢复和点扩展函数计算的一种方法,因此,盲卷积复原算法的优点就是,对失真情况还未知的情形下,仍然能够操作恢复模糊图像。
利用约束最小二乘方法实现对受到噪声等因素所干扰的数字图像其恢复的效果和原始图像相比还有一定的差距。建立在该方法的基础上,已经有不少新的恢复算法不断地被提出,而且使得对数字图像的恢复有了越来越好的效果。恢复的图像存在一定的“环”,这些环是由图像灰度变换较大的部分或图像边界产生的。在图像恢复处理中使用的方法很多,但具体使用哪一种,要按照图像的情况做具体分析,然后再决定采用哪种方法进行图像恢复。
附录
一、参考文献
[1] 冈萨雷斯著.数字图像处理(第三版).北京:电子工业出版社,2010
[2] 朱冠南著.基于MATLAB的图像复原设计. 技术交流,2009
[3] 孟永定 马佳著.基于MATLAB实现数字图像恢复. 北京:电子工业出版社,2009
[4] 刘红岩 徐志鹏著.基于MATLAB的数字图像恢复.高校理科研究,2008
[5] 李国立 段汕著.基于约束最小二乘数字图像恢复.人工智能及识别技术,2008
[6] 百度文库. http://wenku.baidu.com/view/c6f69cdb6f1aff00bed51ee0.html,20##-6
二、源代码:
1、高斯噪声的添加以及滤波处理
I= imread('E:\大三下\图像处理英文课件\作业\第六次作业\lena.bmp','bmp');
J=imnoise(I,'gaussian',0,0.01);
figure;
subplot(1,2,1);
imshow(I);
title('源图像lena.bmp');
subplot(1,2,2);
imshow(J);
title('加入gaussian噪声后的lena.bmp');
n1=7;sigma1=1.5;n2=3;sigma2=1.5;theta=0;
r=[cos(theta) -sin(theta); sin(theta) cos(theta)];
for i = 1 : n2
for j = 1 : n1
u = r*[j-(n1+1)/2 i-(n2+1)/2]';
h(i,j)=exp(-u(1)^2/(2*sigma1^2))/(sigma1*sqrt(2*pi))*exp(-u(2)^2/
(2*sigma2^2))/(sigma2*sqrt(2*pi));
end
end
h = h / sqrt(sum(sum(h.*h)));
f1=conv2(J,h,'same');
subplot(1,2,2);
figure;
imagesc(f1);
title('高斯平滑后的lena.bmp(7x7)');
colormap(gray);
2、图像模糊及添加噪声
①图像的运动模糊
I= imread('E:\大三下\图像处理英文课件\作业\第六次作业\lena.bmp','bmp');
figure;
subplot(1,2,1);
imshow(I);
title('源图像lena.bmp');
f=double(I); % 数据类型转换,MATLAB不支持图像的无符号整型的计算
g=fft2(f); % 傅立叶变换
g=fftshift(g); % 转换数据矩阵
[M,N]=size(g);
a=0.5;b=0.5;T=0.1;
m=fix(M/2);
n=fix(N/2);
j=sqrt(-1);
for i=1:M
for k=1:N
h=(T/(pi*(i*a+k*b)))*sin(pi*(i*a+k*b))*exp(-j*pi*(i*a+k*b));
end
result(i,k)=h*g(i,k);
end
result=ifftshift(result);
J1=ifft(result);
J2=uint8(real(J1));
subplot(1,2,2);
imshow(J2);
title('模糊化lena.bmp');
②图像加噪
figure;subplot(1,2,1);
imshow(J2);
title('运动模糊后的lena.bmp(角度为45)');
J3=imnoise(J2,'gaussian',0,0.01);
subplot(1,2,2);
imshow(J3);
title('加噪并模糊的lena.bmp');
3、图像恢复
①维纳滤波恢复图像
I= imread('E:\大三下\图像处理英文课件\作业\第六次作业\lena.bmp','bmp');
H=fspecial('motion',50,45);
J=imfilter(I,H,'circular','conv');
figure;
subplot(1,2,1);
imshow(J);
title('运动模糊后的lena.bmp(角度为45)');
J1=imnoise(J,'gaussian',0,0.01);
subplot(1,2,2);
imshow(J1);
title('加噪并模糊的lena.bmp');
%figure;
J2=deconvwnr(J1,H),[]);
imshow(J2);
title('模糊噪声图像的维纳滤波复原');
figure;
noise=imnoise(zeros(size(I)),'gaussian',0,0.01);
NSR=sum(noise(:).^2)/sum(im2double(I(:)).^2);
J3=deconvwnr(J1,H,NSR);
imshow(J3);
title('引入SNR的维纳滤波复原');
②约束最小二乘方法恢复图像
I= imread('E:\大三下\图像处理英文课件\作业\第六次作业\lena.bmp','bmp');
I1=checkerboard(8);
PSF=fspecial('motion',50,45);
V=0.0001;
J=imfilter(I,PSF,'circular','conv');
J1=imnoise(J,'gaussian',0,0.01);
figure;subplot(1,2,1);
imshow(J1);
title('模糊加噪图像');
NoisePower=V*prod(size(I));
[G,LAGRA]=deconvreg(J,PSF,NoisePower);
subplot(1,2,2);
imshow(G);
title('约束最小二乘滤波复原');