Parametric Wiener Filter

摘   要：本实验主要通过维纳滤波法（也称最小均方误差滤波法）对图像进行恢复与重建，维纳滤波器没有逆滤波中的退化函数为零的问题，更容易把噪声与有用信息的图像分离出来。首先给图像中添加高斯噪声和运动污损滤波，然后从噪声中提取恢复图像，并与原图像进行对比。图像复原的大部分过程是一个客观的过程，是利用退化现象的某种先验知识来重建或复原被退化的图像。

 

 

一、实验目的

(a) 编写一个给图像中添加高斯噪声的程序，输入参数为噪声的均值与方差。

(b) 编写程序实现公式(5.6-11)所示的污损滤波；

(c) 如图 5.26(b)所示，对图像 5.26(a) 进行+45º方向，T = 1 的污损滤波；

(d) 对污损后的图像加入均值为０，方差为 10 的高斯噪声；

(e) 编写程序使用公式(5.8-6)所示的参数维纳滤波对图像进行恢复。

二、技术论述

1 高斯噪声

在空间域和频域中高斯噪声在数学上非常容易处理，所以在实践中常使用高斯噪声（又称正态噪声）模型，在多数场合高斯模型往往可以在一定程度上导致最好的处理结果。高斯随机变量z的概率密度函数是：

其中z表示灰度值， 表示z的均值， 为z的标准差， 为z的方差。

2 图像的退化和复原

图像复原试图利用退化现象的某种先验知识来复原被退化的图像，以预先确定的目标来改善图像。因此复原技术是面向退化模型的，并且采用相反的过程进行处理，以便恢复出原图像。如图2所示是图像退化/复原过程的模型，图像复原处理退化过程中可以被模型化为一个退化函数和一个加性噪声项，处理一幅输入图像f(x,y) 产生一幅退化图像。给定g(x,y) 和关于退化函数H(u,v) 的一些知识以及外加噪声项n(x,y), 图像复原的目的是获得关于原始图像的近似估计  (x,y) 。通常我们希望这一估计尽可能接近原始输入图像，并且H(u,v)和f(x,y) 的信息知道得越多，所得到的 (x,y)越与f(x,y)接近。

图2 图像退化/复原过程的模型

一个线性、位置不变性的系统函数H在空间域的退化函数为：

又因为空间域上的卷积等同于频域上的乘积，所以把退化图像的模型函数写成等价的频域得形式： 。

退化函数 为：

3 维纳（最小均方误差）滤波法

    维纳滤波法建立在图像和噪声函数都是随机变量的基础上，综合了退化函数和噪声统计特性两个方面进行复原处理，目标是找到未污染图像f的一个估计 ，让f和 之间的随机方差最小，误差量度为：

其中E{.}表示参数的期望值。设噪声和图像不相关，其中一个或另一个有零均值，且估计的灰度级是退化图像灰度级的线性函数。在这些条件下 中的误差函数的最小值在频域中的表达式为：

上式的推导过程中利用了一个复数量与其共轭的乘积等于该复数量幅度的平方。其中 是退化函数、 、 为噪声的功率谱、 为退化函数的功率谱

三、实验结果讨论

通过素材Fig5.07(a).jpg进行+45度方向，T=1的污损滤波，并对污损图像进行高斯噪声加噪处理，再验证维纳滤波法的有效性。通过编写的高斯噪声函数gaussian_noise.m为图像Fig5.07(a).jpg加入高斯噪声。对比观察加入的高斯噪声的均值为0，方差分别为10、500、2000的加噪后的图像效果,如图3.1所示，+45度方向，T=1的污损滤波处理后的图像如图3.2所示, 经过运动模糊、加噪处理和维纳滤波后的图像如图3.3(a)、(b)所示（程序见附录）

图3.1 方差分别为10、500、2000的加噪后的图像

图3.2 +45度方向，T=1的污损滤波处理后的图像

图3.3(a) 经过运动模糊、加噪处理（0均值，500方差）和维纳滤波后的图像

图3.3(b) 经过运动模糊、加噪处理（0均值，1000方差）和维纳滤波后的图像

四、实验结论

经过维纳滤波法（也称最小均方误差滤波法）对图像进行退化与重建的图像进行对比分析可知：

维纳滤波法对图像的恢复效果较好，在实验要求不高的情况下可以把有噪声且运动模糊的图像恢复出不错的效果，非常接近原图想的外观。但对于添加噪声方差较大的图像的滤波结果并不理想，恢复后会损失掉图像的一部分细节，要辨认图像中的文字还是有些困难的。基于维纳滤波法的这一缺点理论上有另一种方法——约束最小乘方滤波法，它的显著特点就是对于应用到每一幅图像都能产生最优的结果。

 

附录：程序清单

1 % 编写给图像加入高斯噪声的函数

% ima为输入图像，a和b分别为高斯函数的均值和方差

% 函数randn()产生正态分布的随机数或矩阵的函数

function  imt=gaussian_noise(ima,a,b)

%ima=imread('Fig5.26(a).jpg');

% a = 0;

% b = 10;

[m,n] =size(ima);

x = uint8(a +sqrt(b) * randn(m,n)); % 产生一个m×n的高斯噪声矩阵

imt = ima +x;                       % 给原图叠加高斯噪声

% imshow(imt); title(‘方差为10的加噪处理后图像’）

return

 

2 % 编写函数给图像加入模糊污损噪声

% ima为输入图像，a和b分别为x轴和y轴的运动速度参数，T为图像位移时间

functionimag=motion(ima,a,b,T)

%ima 为输入图像，a、b为方向，T为时间

%ima=imread('Fig5.26(a).jpg');

%a=0.1; b=0.1; T=1;

[m,n]=size(ima);

ima=fft2(ima);%原图像傅里叶变换

for u=1:m

    for v=1:n

       Y(u,v)=u*a+v*b;

       H(u,v)=(T*sin(pi*Y(u,v))*exp(-1j*pi*Y(u,v)))/(pi*Y(u,v));%退化传递函数

       G(u,v)=ima(u,v)*H(u,v);%频域相乘，得到频域的退化图像

    end

end

imag=ifft2(G);%傅里叶反变化到时域

imag=uint8(abs(imag));

% imshow(imag); title('+45度方向且T=1 的污损滤波图像'）

return

 

3 %维纳滤波函数

function F1 =Wiener(ima,imag)

%weina 维纳滤波函数

%ima是原始图像，imag是污损加噪声图像

F0=fft2(ima);   % 对原始图像作傅立叶变换

G0=fft2(imag);   % 对运动模糊加高斯噪声图像作傅立叶变换

[m,n]=size(ima);

K=0.0001;        

for u=1:m

  for v=1:n

     H(u,v) = G0(u,v)/F0(u,v);  % 加噪图像作傅立叶变换除以原始图像作傅立叶变换，得到退化函数

     % 计算维纳滤波函数

      H0(u,v)=(abs(H(u,v)))^2; 

      H1(u,v)=H0(u,v)/(H(u,v)*(H0(u,v)+K));

      F(u,v)=H1(u,v)*G0(u,v);  % 恢复、计算原始图像的傅立叶变换估计

  end

end

F1=ifft2(F);

F1=uint8(abs(F1));

end

 

4 %参数维纳滤波对图像进行恢复

ima=imread('Fig5.26(a).jpg');

[m,n]=size(ima);

figure(1);

subplot(2,2,1);imshow(ima);title('原始图像');

 

imag=motion(ima,0.1,0.1,1); %调用motion子函数模糊图像

subplot(2,2,2);imshow(imag);title('+45度运动模糊图像');

 

imt=gaussian_noise(imag,0,1000); %调用gussion_noise子函数添加高斯噪声

subplot(2,2,3);imshow(imt);title('加噪模糊图像');

 

F1=Wiener(ima,imt); %调用weina子函数进行维纳滤波，复原图像

subplot(2,2,4),imshow(F1),title('维纳滤波复原图像');