图像复原之直接逆滤波

由退化函数H退化的图像复原的最简单的方法是直接做逆滤波，设图像退化前的傅里叶变换为F(u,v)，退化后的傅里叶变换为G(u,v)，系统函数即退化函数的傅里叶变换为H(u,v)。

所谓直接逆滤波，就是用退化函数除退化图像的傅里叶变换，得到退化前图像的傅里叶变换的估计，

为F(u,v)的估计，则，同时，其中N(u,v)为噪声的傅里叶变换。因此，可得，由该式可知，即使知道退化函数，也不能准确的复原图像，因为N(u,v)未知，甚至有更糟的情况是如果退化函数是零或是非常小的值，则N(u,v)/H(u,v)的值比较大，很容易支配F(u,v)的估计值。解决这个问题的一种方法是限制滤波的频率，从频谱图可知，高频分量的值接近0，而H(0,0)在频率域中通常是H(u,v)的最高值。因此可缩短滤波半径，使通过的频率接近原点，减少遇到零值的概率。

直接逆滤波的MATLAB代码如下

image_o=imread('C:\Program Files\MATLAB\R2013a\bin\work\图像复原\lena.bmp');
subplot(1,3,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);
subplot(1,3,2);
imshow(x);
%
title('退化图像');

image_d=imread('C:\Program Files\MATLAB\R2013a\bin\work\图像复原\lena_deterioration.bmp');
%直接逆滤波图像复原

ff=im2double(image_d);%将图像灰度值归一化到0-1之间

% 傅里叶变换
f_Id=fft2(ff);
f_Id=fftshift(f_Id);
fH_Id=f_Id;
[M,N]=size(fH_Id);
% 逆滤波
threshold=78;
if threshold>M/2
        %全滤波
        fH_Id=fH_Id./(H+eps);
else
        %对一定半径范围内进行滤波
        for i=1:M
            for j=1:N
                if sqrt((i-M/2).^2+(j-N/2).^2)<threshold
                    fH_Id(i,j)=fH_Id(i,j)./(H(i,j)+eps);
                end
            end
        end
end

% 执行傅立叶逆变换
fH_Id1=ifftshift(fH_Id);
f_new=ifft2(fH_Id1);
f_new=uint8(abs(f_new)*255);
subplot(1,3,3);
imshow(f_new);
title('滤波半径=78的逆滤波复原图像');

效果如下

从实验结果可知，滤波为78比较适宜。