图像退化与复原

一． 实验名称：图像退化与复原 二． 实验目的

1． 了解光电图像的退化原因；

2． 掌握和理解基本的噪声模型，并能对图像进行加噪处理；

3． 了解点扩展函数(PSF)与光学传递函数(OTF)的关系，熟悉几种经典的退化模 型的模拟试验和 OTF 估计方法；

4． 熟悉和掌握几种经典的图像复原方法及其基本原理；

5． 能熟练利用 MATLAB 或 C/C++工具进行图像的各种退化处理，并能编程实现 退化图像的复原。

三． 实验原理

光电成像系统出现图像退化的过程是复杂多变的，为了研究的需要，通常情况下都把退化简化为化为一个线性移不变过程，见下图 1 所示。

因此， 在空域中退化过程可以表示如下：

(x,y)(x,y)(x,y)(x,y)g f h h =*+ (1)

只有加性噪声不存在情况下，退化过程可以模型化如下表达式：

(x,y)(x,y)(x,y)g f h =+ (2)

其频域表达式为：

=(,)+(),)G ,(F u v N u v v u

(3)

图1光电图像退化与复原原理图

针对这种退化图像的复原，除了周期噪声以外，通常都可以采用空间域滤波的方法进行图像复原，此时图像复原与图像增强几乎是没有区别的。常见的空间域滤波方法有均值滤波器和统计排序滤波器。

当退化图像存在线性移不变退化时，图像的复原不能采用简单空间域滤波器来实现，要实现线性移不变退化图像的复原，必须知道退化系统的退化函数，即点扩展函数(x,y)h 。 在点扩展函数已知的情况下，常见图像复原方法有逆滤波和维纳滤波两种。

在考虑噪声的情况下，逆滤波的原理可以表示如下：

()()

()()()()

G u,v N u,v F u,v F u,v H u,v H u,v ù

=

=+

(4)

通常情况下，()N u,v 是未知的，因此即使知道退化模型也不能复原图像。 此外，当(),H u v 的任何元素为零或者值很小时，()(),/,N u v H u v 的比值决定着复原的结果，从而导致图像复原结果出现畸变。对于这种情况，通常采用限制滤波频率使其难以接近原点值，从而减少遇到零值的可能性。

维纳滤波则克服了逆滤波的缺点，其数学模型表示如下：

2*

2()1

()()()()(,)/(,)f H u,v F u,v G u,v H u,v H u,v S u v S u v h ù

轾犏=犏犏+犏臌

(5)

然而，为退化图像的功率谱很少是已知的，因此常常用下面表达式近似：

2*

2()1

()()()()H u,v F u,v G u,v H u,v H u,v k ù

轾犏=犏犏+犏臌

(6)

因此，本实验的内容就是利用上述经典图像复原的原理，对降质退化图像进行复原。

四． 实验步骤

本次实验主要包括光电图像的退化模型和复原方法实现两大部分内容。 (一) 图像的退化图像 1、 大气湍流的建模

1) 湍流引起图像退化的光学传递函数（OTF ）生成。已知湍流退化模型的 OTF 表达式如下：

225/6(,)exp[-()]H u v k u v =+

(7)

其中，k 为一个常数，反映了大气湍流的严重程度。(,)u v 分别代表了（x,y ） 方向的频率坐标。为了生成中心化的OTF ，可以考虑将式(7)改写为：

()

5/6

2

2(,)exp[-(/2)(/2)

]H u v k u M v N =-+- (8)

其中，M ，N 为图像的长和宽。

2) 读入一幅灰度图像，设定式(8)中0.0025k =，进行退化试验。分别显示 原始图像、退化模型和退化图像。

3) 设定0.0010.00025k =、重复上一步的试验。

对原图形进行灰度处理将上述结果进行fft 处理

得到FP

读入原始图像

设计退化湍流模型为H 结束

显示原图像

显示传递函数

由FP 与H 进行相关处理，

得到退化图像

显示退化图像

图2大气湍流的退化过程

2、 运动模糊的图像退化试验

1) 匀速直线运动引起图像退化的光学传递函数（ OTF ）生成。已知相机匀 速直线运动的 OTF 表达式如下：

()

(,)sin[()]()

j ua vb T

H u v ua vb e ua vb p p p -+=

++ (9)

其中，T 为相机曝光时间，a ,b 分别表示(,)x y 方向的速度；(,)u v 分别对 应(,)x y 方向的频率坐标。

2) 读入一幅灰度图像，设定式（ 9）中 T = 1.0， a=b=0.1，编写 MATLAB 代 码进行模糊退化试验。要求分别显示原始图像、退化模型和退化图像。

3) 设定不同的值，a ，b 值，重复上一步的试验。

4) 利用数字显微镜或其它图像采集设备，移动物体过程中，采集图像。

对原图形进行灰度处理将上述结果进行fft 处理

得到FP

读入原始图像

设计运动模糊模型为H 结束

显示原图像

显示传递函数

由FP 与H 进行相关处理，

得到退化图像

显示退化图像

图3运动模糊的图像退化

(二) 图像复原试验 1、 逆滤波

1) 根据试验（一） 设计一幅退化图像（包括噪声污染+模糊退化两部分），其中模糊退化可选高斯模糊、大气湍流模糊或运动模糊（ 方向可任意指定，如10 度、20度、45度等），噪声模型可自行设定。

2) 利用 MATLAB 编程实现利用全逆滤波方法对退化图像的复原。要求在同一个窗口下显示原始退化图像、复原结果及复原结果与理想图像的差值图共 3 个图，并对复原 结果进行必要的分析。 逆滤波复原公式如下：

()

()()

G u,v F u,v H u,v ù

=

(10)

其中，()G u,v 为退化图像的傅立叶变换，()H u,v 为退化系统的光学传递函数（OTF ）。

3) 伪逆滤波：为了防止逆滤波中(,)H u v 过小，使得复原后的图像数据过

大

和放大噪声，可采用频谱半径（阈值）限制下的逆滤波方法，即

1,(,)(,)0,R H u v P u v R

ì????=í

?????

(11)

其中，R 为中心化频谱(,)H u v 中某点到原点（零频）的距离或半径。另一种替代方法是直接限制(,)H u v 的值，即

1,

(,)(,)(,)0,(,)H u v H u v P u v H u v s s

ì??￡?=í

??>???

(12)

其中，σ为一个阈值，用于限制频谱的幅度值。这种方法被称为伪逆滤波。 实验要求利用式(11)方式的伪逆滤波重复实验步骤内容2）所涉及的图像。 2、 Wiener 滤波

1) 针对以上逆滤波设计的退化图，编程实现利用Wiener 滤波对其进行复原。 滤波原理如下：

2*

2()1

()()()()H u,v F u,v G u,v H u,v H u,v k ù

轾犏=犏犏+犏臌

(13)

其中，()G u,v 为退化图像的傅立叶变换，()H u,v 为退化系统的光学传递函数（OTF ），

k 为一个与信噪比有关的调节因子。要求在同一个窗口下显示理想图像（退化前）、

退化图像、复原结果等共3个图，并对复原结果进行必要的分析。

2) 改变k 值，重复试验内容 1）。

以上应根据原理自行编写代码，不允许直接调用MATLAB 自带的deconvwnr()函数。

对原图形进行灰度处理，

并加入椒盐噪声得到f1

由实验一对f1进行运动退化处理得f2（a=b=0.1

）

读入原始图像

对f2进行fft ，fftshift 处理得G1

结束

显示原始退化图像

分别利用逆滤波，伪逆滤波，wiener 滤波对G1进行复原处理

得到F1分别显示三种方式下的复原

结果F1

将复原结果F1与理想图像做减

法，得到差值F2

分别显示三种方式下的差值

F2

图4全逆，伪逆，wiener 滤波复原过程

五． 实验结果及分析

1、 大气湍流的建模

图5

分析：由上述结果可知，大气湍流会使图像变得模糊，而k 值越大，其模糊效果越明显。

2、 运动模糊的图像退化试验

图7

图6

图8

图

9

分析：由上述结果可知，随着a ，b 的值变大，图像模糊变得明显，人眼看起来好像是由于运动速度过快造成的模糊。 3、 图像复原试验

分析：图10，图11，图12 分别为全逆，伪逆，wiener 对运动模糊（a=b=0.1）滤波的结果,从中可以看出wiener 是三者中对运动模糊复原效果最好的滤波方式，且wiener 中k 值越小复原效果越好。

图10

图11

图

12

六．实验心得体会和建议

●心得体会：通过这次实验使我了解了图像退化的原因，以及相关的退化模型，并学会以matlab为平台利用退化模型对图像进行退化处理以及退化图像的复原处理。

●建议：可以要求利用C或C++进行图像的退化与复原处理。

七．程序源代码

% title：atmosphere

% explain：本程序利用大气湍流模型对理想图像进行退化

f=imread('3.jpg');

figure(1)

subplot(131),imshow(f),title('原始图像')

f=rgb2gray(f);

Fp=fft2(f);

[m,n]=size(f);%绘制网格点

[v,u]=meshgrid(1:n,1:m);

u=u-floor(m/2);

v=v-floor(n/2);

k=0.00025;

Duv=u.^2+v.^2;

H=exp(-k.*Duv.^(5/6));

G=H.*fftshift(Fp);

f1=abs(ifft2(G));

nchar = num2str(k);

ltext = strcat('k=', nchar);%标题注释

subplot(132),imshow(H),title(['传递函数',ltext]);

subplot(133),imshow(f1,[]),title('退化图像');

% title：move

% explain：本程序利用运动模糊模型对理想图像进行退化

f=imread('3.jpg');

figure(1)

subplot(131),imshow(f),title('原始图像')

f=rgb2gray(f);

[m,n]=size(f);

[v,u]=meshgrid(1:n,1:m);

u=u-floor(m/2);

v=v-floor(n/2);

T=1.0;

a=0.3;b=0.3;

% a=0.1,b=0.1;

% a=0.01,b=0.01;

z=pi*(u*a+v*b)+eps;

H=T./z.*sin(z).*exp(-1j*z);

Fp=fft2(f);

G=H.*fftshift(Fp);

f1=abs(ifft2(G));

nchar = num2str(a);

ltext = strcat('a=b=', nchar);

subplot(132),imshow(H),title(['传递函数',ltext]);

subplot(133),imshow(f1,[]),title('退化图像');

% title：recovery

% explain：本程序利用运动模糊对加入椒盐噪声的理想图像进行退化，之后分别用全逆滤波

%，伪逆滤波，wiener滤波对设计的退化图进行处理，观察三种滤波的复原效果。

clc,clearall,close all;

f=imread('4.jpg');

f=rgb2gray(f);

figure,imshow(f),title('原始图像')

f1=imnoise(f,'salt& pepper',0.02);%加入椒盐噪声后的图像

[m,n]=size(f);

[v,u]=meshgrid(1:n,1:m);%画网格点

u=u-floor(m/2);

v=v-floor(n/2);

%==================================================================== %运动模糊

%==================================================================== % T=1.0,a=0.1,b=0.1;

% T=1.0,a=0.01,b=0.01;

T=1.0；a=0.1；b=0.1;

z=pi*(u*a+v*b)+eps;

H=T./z.*sin(z).*exp(-1j*z);

G=H.*fftshift(fft2(f1)); %退化图像的频域

f2=real(ifft2(ifftshift(G)));%加了椒盐噪声的退化图像

G1=fftshift(fft2(f2));%求原始退化图像的频域表示

%==================================================================== %逆滤波

%==================================================================== F=G1./H;

F1=real(ifft2(ifftshift(F)));

%反中心化，反傅立叶变换取实部得到复原结果图（包含噪声）

F2=F1-double(f);%复原结果与理想图像的差值

figure,

subplot(131),imshow(f2,[]),title('原始退化图像');

subplot(132),imshow(F1,[]);title('复原结果');

subplot(133),imshow(F2,[]);title('差值图');

%==================================================================== %伪逆滤波

%==================================================================== Duv=u.^2+v.^2;

if Duv.^(1/2)<=1000

P=1./H;

else

P=0;

end

F=G1.*P;

F1=real(ifft2(ifftshift(F)));

F2=F1-double(f);

figure,

subplot(131),imshow(f2,[]),title('原始退化图像');

subplot(132),imshow(F1,[]);title('复原结果');

subplot(133),imshow(F2,[]);title('差值图');

%==================================================================== %wiener滤波

%==================================================================== k1=0.000025;

F=((1./H).*(abs(H).^2)./((abs(H).^2)+k1)).*G1;

F1=real(ifft2(ifftshift(F)));%反中心化，反傅立叶变换取实部得到复原结果图（包含噪声）

F2=F1-double(f);

figure,

subplot(131),imshow(f2,[]),title('原始退化图像'); subplot(132),imshow(F1,[]);title('复原结果'); subplot(133),imshow(F2,[]);title('差值图');

八． 思考题

1. 简要叙述图像退化的原因。 答：

● 原因：在图像的获取和传输过程中，由于成像系统各种因素的影响，使 得最终得到的图像相比于理想图像质量降低。 2. 图像复原与图像增强有什么异同？ 答：

表格1图像复原与图像增强的异同

3. 简述同态滤波的设计思想和处理流程，说明该方法适合于何种场合下的图 像恢复。

答：

● 设计思想：通过把频率过滤和灰度变换相结合进行图像处理，依靠图像 的照度/ 反射率模型作为频域处理的基础，利用压缩亮度范围和增强对比度来改善图像的质量。

● 处理流程：

1) 对原图取对数得到如下两个加性分量：

ln (,)ln (,)ln (,)i r f x y f x y f x y =+

(14)

2) 频域表示为：

()()()i r FFT ln f x,y =FFT ln f x,y +FFT ln f x,y ???????????? (15)

3) 设计一个滤波器()H u,v 进行滤波处理。

4)傅立叶反变换，即IFFT。

5)求指数结果，得到复原结果。

适用场合：作为非线性复原方法的同态滤波，它可以用来恢复在非均匀光照射下的图像。

报告评分：指导教师签字：