图像复原之维纳滤波

       维纳滤波也称最小均方误差滤波，它能处理被退化函数退化和噪声污染的图像。该滤波方法建立在图像和噪声都是随机变量的基础之上，目标是找到未污染图像f(x,y)的一个估计，使它们之间的均方误差最小，即，其中E{.}是参数期望值。在假设噪声和图像不相关，其中一个或另一个有零均值，且估计中的灰度级是退化图像中灰度级的线性函数的条件下，均方误差函数的最小值在频率域由如下表达式给出：

                                                 

其中，H(u,v)为退化函数，H*(u,v)为H(u,v)的复共轭，|H(u,v)|^2=H(u,v)H*(u,v)，Sη(u,v)=|N(u,v)|^2=噪声的功率谱，Sf(u,v)=|F(u,v)|^2=未退化图像的功率谱，H(u,v)是退化函数的傅里叶变换，G(u,v)是退化后图像的傅里叶变换。

        从上面的公式可以发现，如没有噪声，即Sη(u,v)=0，此时维纳滤波变为直接逆滤波，如有噪声，那么Sη(u,v)如何估计将成问题，同时Sf(u,v)的估计也成问题。在实际应用中假设退化函数已知，如果噪声为高斯白噪声，则Sη(u,v)为常数，但Sf(u,v)通常难以估计。一种近似的解决办法是用一个系数K代替Sη(u,v)/Sf(u,v)，因此上面的公式变为如下式所示：

                                                                                                   

在实际应用中，根据处理的效果选取合适的K值。

这里对256×256大小的灰度图像进行退化处理，然后添加均值为0，方差为0.005的高斯噪声，对得到的退化图像使用直接逆滤波和维纳滤波复原，并比较效果，MATLAB代码如下：

image_o=imread('C:\Program Files\MATLAB\R2013a\bin\work\图像复原\lena.bmp');
subplot(2,2,1);
imshow(image_o);
title('原图像');
%频率域退化图像，退化函数H(u,v)=exp(-0.0025*( (u-M/2).^2+(v-N/2).^2).^(5/6) )
%傅里叶变换
f=im2double(image_o);
F=fft2(f);
F=fftshift(F);
%执行退化
[M,N]=size(F);
[u,v]=meshgrid(1:M,1:N);%生成二维坐标系
H=exp(-0.0025* ( (u-M/2).^2+(v-N/2).^2).^(5/6) );
F=F.*H;
%添加高斯噪声
%傅里叶反变换
X=ifftshift(F);
x=ifft2(X);
x=uint8(abs(x)*256);
x1=imnoise(x,'gaussian',0,0.005);%添加高斯噪声
subplot(2,2,2);
imshow(x1);
%
title('退化图像');

%
Id=im2double(x1);
% 傅里叶变换
f_Id=fft2(Id);
f_Id=fftshift(f_Id);
fH_Id=f_Id;
D=abs(H);
D=D.^2;

[M1,N1]=size(fH_Id);
% 逆滤波
threshold=48;
K=0.08;
if threshold>M1/2
    % 全滤波
    fH_Id=fH_Id./(H+eps);
else
    %对一定半径范围内进行滤波
    for i=1:M1
        for j=1:N1
            if sqrt((i-M1/2).^2+(j-N1/2).^2)<threshold
                % 维纳滤波公式
                fH_Id(i,j)=fH_Id(i,j)./(H(i,j)) .* (D(i,j)./(D(i,j)+K));
            end
        end
    end
end
% 执行傅立叶逆变换
fH_Id1=ifftshift(fH_Id);
f_new=ifft2(fH_Id1);
f_new=uint8(abs(f_new)*255);
subplot(2,2,3);
imshow(f_new);
title('维纳滤波复原图像');

fH_Id2=f_Id;
[M2,N2]=size(fH_Id2);
if threshold>M2/2
        %全滤波
        fH_Id2=fH_Id2./(H+eps);
else
        %对一定半径范围内进行滤波
        for i=1:M2
            for j=1:N2
                if sqrt((i-M2/2).^2+(j-N2/2).^2)<threshold
                    fH_Id2(i,j)=fH_Id2(i,j)./(H(i,j)+eps);
                end
            end
        end
end

fH_Id3=ifftshift(fH_Id2);
f_new3=ifft2(fH_Id3);
f_new3=uint8(abs(f_new3)*255);
subplot(2,2,4);
imshow(f_new3);
title('直接逆滤波复原图像');

处理效果如下图所示

从运行结果来看，维纳滤波性能优于直接逆滤波