图像复原
一、目的
1了解图像退化原因与复原技术分类化的数学模型;
2熟悉图像复原的经典与现代方法;
3热练掌握图像复原的应用;
4、通过本实验掌握利用MATLAB编程实现数字图像的图像复原。
二、设计原理:
图像复原处理是建立在图像退化的数学模型基础上的,这个退化数学模型能够反映图像退化的原因。图像的退化过程可以理解为施加于原图像上的运算和噪声两者联合作用的结果,图像退化模型如图1所示,可以表示为:
g ( x, y) H [ f ( x, y)] n( x, y) f ( x, y)h( x, y) n( x, y) (1)
![v2-1a4c7b6fe6624c5d5d3b5feac4eb36cd_b.jpg](https://img-blog.csdnimg.cn/img_convert/c7d9dbd634efd36958fa68962081f453.png)
图1 图像退化模型
- 在测试图像上产生高斯噪声lena图-需能指定均值和方差;并用滤波器(自选)恢复图像;
噪声是最常见的退化因素之一,也是图像恢复中重点研究的内容,图像中的噪声可定义为图像中不希望有的部分。噪声是一种随机过程,它的波形和瞬时振幅以及相位都随时间无规则变化,因此无法精确测量,所以不能当做具体的处理对象,而只能用概率统计的理论和方法进行分析和处理。本文中研究高斯噪声对图像的影响及其去噪过程。
①高斯噪声的产生:
所谓高斯噪声是指它的概率密度函数服从高斯分布(即正态分布)的一类噪声。一个高斯随机变量z的PDF可表示为:
P(z)=
![v2-a456d9a3e750e4da39ba5bad66370ee0_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-a456d9a3e750e4da39ba5bad66370ee0_b.jpg)
(2)
其中z代表灰度,u是z的均值,
![v2-43fca9fa39f35b19c71998bcfc974f52_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-43fca9fa39f35b19c71998bcfc974f52_b.jpg)
是z的标准差。高斯噪声的灰度值多集中在均值附近。
![v2-cb789e9e749fdf40d17767681dda7c06_b.jpg](https://img-blog.csdnimg.cn/img_convert/377b268adcddfe116478744c5942b061.png)
图2 高斯函数
可以通过不同的算法用matlab来产生高斯噪声。
②高斯噪声对信号的影响
噪声影响图像处理的输入、采集、处理的各个环节以及输出结果的全过程,在图像中加高斯噪声通常会使图像变得模糊并且会出现细小的斑点,使图像变得不清晰。
③去除高斯噪声的一些方法
去除高斯噪声的方法有直方图变换,低通滤波,高通滤波,逆滤波,维纳滤波,中值滤波等。本文应用高斯平滑滤波进行去噪处理。
维纳滤波对高斯白噪声的图像滤波效果较好,具有比较好的选择性,可以更好地保存图像的边缘和高频细节信息。所以,维纳滤波在大多数情况下都可以获得满意的结果,尤其对含有高斯噪声的图像。
三、实现步骤
1打开计算机,安装和启动MATLAB程序;程序组中“work”文件夹中应有待处理的图像文件;
2利用MatLab工具箱中的函数编制图像复原的函数;
3 在测试图像上产生高斯噪声lena图-需能指定均值和方差;并用滤波器(自选)恢复图像;
4 推导维纳滤波器并实现下边要求;
维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理,其目标是寻找一个滤波器,使得复原后图像 与原始图像 的均方误差最小:
因此维纳滤波器又称为最小均方误差滤波器。
在频率中用下式表达:
![v2-1e852c24bc05bf7e37d3ac33dbbb57c2_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-1e852c24bc05bf7e37d3ac33dbbb57c2_b.jpg)
其中,
![v2-d79fbb824eb3f7f682a431aa7a5317c7_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-d79fbb824eb3f7f682a431aa7a5317c7_b.jpg)
是退化图像的傅立叶变换,
![v2-f847c3b1d68ded335e17727a7c4542f7_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-f847c3b1d68ded335e17727a7c4542f7_b.jpg)
是退化函数。
![v2-7b400cad5a361a860d5e4997c72abf00_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-7b400cad5a361a860d5e4997c72abf00_b.jpg)
,
![v2-7d9c6f7c4d95ed7dfa2beab5b5179581_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-7d9c6f7c4d95ed7dfa2beab5b5179581_b.jpg)
是
![v2-05a22e78b0a2cfc10b2bb18f2d1714af_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-05a22e78b0a2cfc10b2bb18f2d1714af_b.jpg)
的复共轭。
![v2-2cbb8802be7a93c624c02e050ff60863_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-2cbb8802be7a93c624c02e050ff60863_b.jpg)
,为噪声的功率谱。
![v2-0f604318d834497e049822c3b66bde21_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-0f604318d834497e049822c3b66bde21_b.jpg)
,为未退化图像的功率谱。
维纳滤波器的推导:
合理假设要估计的信号f(x,y)为0均值平稳随机过程,噪声为0均值平稳随机过程。
g ( x, y) h( x, y) * f ( x, y)
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
( x, y)
![v2-88e8252e6365bb6d5ae46ff870c984de_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-88e8252e6365bb6d5ae46ff870c984de_b.jpg)
根据已退化图像g(x,y),利用线性估计器
![v2-6fef1070e048291e7b2b021b82472905_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6fef1070e048291e7b2b021b82472905_b.jpg)
来估计原始图像。
先证: g ( x, y) f ( x, y)
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
( x, y),
最小化均方误差:
![v2-64645e8da5206aef0645c18cdadf7bf5_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-64645e8da5206aef0645c18cdadf7bf5_b.jpg)
=
![v2-ca50a8cfe62ba2afa826b6a28228b89a_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-ca50a8cfe62ba2afa826b6a28228b89a_b.jpg)
,
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 {
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
w(k1, k2 )[g (m−k1, n−k2)g*(x,y)]}
= E {
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
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).[
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
h (k1, k2)f(x-k1,y-k2)+
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
(x, y)]*
=
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
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 [
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
h (k1, k2) f(m-k1,n-k2)+
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
(m,n)]
*[
![v2-be83ddcff400877757f5bab8438e397a_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-be83ddcff400877757f5bab8438e397a_b.jpg)
h *(
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,
![v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg)
2) f*(x-
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,y-
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
2)+
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
*(x, y)]
=
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
![v2-be83ddcff400877757f5bab8438e397a_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-be83ddcff400877757f5bab8438e397a_b.jpg)
h (k1, k2) h *(
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,
![v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg)
2) E[f(m-k1,n-k2)
f*(x-
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,y-
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
2)]+ E{
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
(m,n)
![v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-df6ce3893b09bb04e7896f09abc4869f_b.jpg)
*(x, y)}
得: Rg (m-x, n-y)
=
![v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-e05f135dc4c51b72dbe10fae9dc53071_b.jpg)
![v2-be83ddcff400877757f5bab8438e397a_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-be83ddcff400877757f5bab8438e397a_b.jpg)
h(k1,k2)h*(
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,
![v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg](https://ss.csdn.net/p?https://pic4.zhimg.com/v2-bff47a54059129c26cb5c0b7aef91fb7_b.jpg)
2) Rf (m-k1-x+
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
1,n-k2-y+
![v2-6642c095384106da636366ed218defbd_b.jpg](https://ss.csdn.net/p?https://pic2.zhimg.com/v2-6642c095384106da636366ed218defbd_b.jpg)
2) + R
![v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg)
(m-x,n-y)
= Rf (m-x,n-y)* h(m-x,n-y)* h* (-(m-x),-(n-y))+ R
![v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg)
(m-x,n-y)
所以,
Rg (m, n)= Rf (m,n)* h(m,n)* h*(-m,-n)+ R
![v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg)
(m-x,n-y)
上式作Flouier变换得:
Pg (w1,w2)= Pf (w1, w2)*H(w1, w2)*H*(-w1,-w2)+P
![v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg)
(w1, w2)
=
![v2-bc5034f2d5f7c76d29215edaca72a5fe_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-bc5034f2d5f7c76d29215edaca72a5fe_b.jpg)
Pf (w1, w2)+ P
![v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg](https://ss.csdn.net/p?https://pic1.zhimg.com/v2-c5f77b93b32569a52de7e2fdf3c6b9b8_b.jpg)
(w1, w2) (5)
综合以上式(3)(4)(5)可得:
![v2-929f10624e21d3b955ac6eed12ff0676_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-929f10624e21d3b955ac6eed12ff0676_b.jpg)
推导完毕。
- 实现模糊滤波器如方程Eq.(5.6-11).
在图像复原中,有3种主要的估计退化函数的方法:(1)图像观察估计法,(2) 试验估计法,(3)模型估计法。本文实现从基本原理开始推导数学模型的模型估计法。其SFR为:
![v2-dc0c33664ff8f9aed0ad42e3babb4c8e_b.jpg](https://ss.csdn.net/p?https://pic3.zhimg.com/v2-dc0c33664ff8f9aed0ad42e3babb4c8e_b.jpg)
根据上式设计滤波器,实现的程序代码。
I=imread('E:\test1.bmp');
figure(1);imshow(I,[]);
title('原图像');
%设置PSF参数
MF=imfilter(I,PSF,'circular');
noise=imnoise(zeros(size(I)),'gaussian',0,0.001);
MFN=imadd(MF,im2uint8(noise));
figure(2);imshow(MFN,[]);
title('运动模糊图像');
NSR=sum(noise(:).^2)/sum(MFN(:).^2);
figure(3);
imshow(deconvwnr(MFN,PSF,NSR),[]);
title('逆滤波复原');
figure(4);
imshow(deconvwnr(MFN,PSF,NSR),[]);
title('维纳滤波复原');
NP=0.002*prod(size(I));
[reg1 LAGRA]=deconvreg(MFN,PSF,NP/3.0);
figure(5);imshow(reg1);
title('最小二乘滤波复原');
五、实验示例
![v2-69ad753126b05e5d37b018628396e930_b.jpg](https://img-blog.csdnimg.cn/img_convert/1607f919a0a377e6d1379bcc84b825c9.png)
![v2-b6be3f322bfe75b5bd22654015bbe8b9_b.jpg](https://img-blog.csdnimg.cn/img_convert/65a4383e0b3a18abfb0e54bba1bf3fd8.png)
![v2-8789a81aba1d8bb7e60b15c7257806a3_b.jpg](https://img-blog.csdnimg.cn/img_convert/e0e02f755defe0f559725f635a429425.png)
![v2-f6ae49fd8adf68419712f88df469b7c0_b.jpg](https://img-blog.csdnimg.cn/img_convert/4b94efab363a678368e2eee46c1315fd.png)
![v2-e4ddaf690376e02037a6c9c26c79761c_b.jpg](https://img-blog.csdnimg.cn/img_convert/3282081f1ed0acabce4da68531dbb137.png)
六、结果分析
实验中发现,采用维纳滤波恢复可以取得比较好的效果,这个算法可以使估计的点扩散函数值PSF更加接近它的真实值。在我们知道模糊图像的点扩展函数的情况下,可以调用常规的图像复原算法;而现实里还会遇见不知道点扩展函数的情况,这个时候我们就可以利用盲卷积复原算法。它是利用原始图像模糊,同时进行清晰图像的恢复和点扩展函数计算的一种方法,因此,盲卷积复原算法的优点就是,对失真情况还未知的情形下,仍然能够操作恢复模糊图像。
利用约束最小二乘方法实现对受到噪声等因素所干扰的数字图像其恢复的效果和原始图像相比还有一定的差距。