图像复原-------逆滤波和维纳滤波介绍及MATLAB实现

由于受到光学系统的像差，成像设备与物体的相对运动等因素的影响，图像会出现一定的失真。要想得到高质量的图像，需要对已经退化后的图像进行复原。

图像退化模型

一般来说，图像的退化模型可以表示为
$$g\left(x,y\right)=h\left(x,y\right)\ast f\left(x,y\right)+n\left(x,y\right)$$
其中，$g\left(x,y\right)$表示退化后的图像，$h\left(x,y\right)$为点扩散函数，$f\left(x,y\right)$为原始图像，$n\left(x,y\right)$为引入的噪声。
在频域上面可以表示为
$$G\left(u,v\right)=H\left(u,v\right)F\left(u,v\right)+N\left(u,v\right)$$

逆滤波

逆滤波是一种常见且直观的图像恢复方法，它的主要过程是将退化后的图像从空间域变换到频域，进行逆滤波后在变换到空间域上从而实现图像的复原。
从上面的介绍可以得到，图像的退化模型为
$$\begin{array}{rl}& g\left(x,y\right)=h\left(x,y\right)\ast f\left(x,y\right)+n\left(x,y\right)\\ & \text{=}{\int }_{-\infty }^{+\infty }f\left(\alpha ,\beta \right)h\left(x-\alpha ,y-\beta \right)d\alpha d\beta +n\left(x,y\right)\end{array}$$
将其变换到频域可以得到
$$G\left(u,v\right)=H\left(u,v\right)F\left(u,v\right)+N\left(u,v\right)$$
那么
$$\begin{array}{rl}& F^{\prime }\left(u,v\right)\text{=}\frac{G\left(u,v\right)}{H\left(u,v\right)}-\frac{N\left(u,v\right)}{H\left(u,v\right)}\\ & =F\left(u,v\right)-\frac{N\left(u,v\right)}{H\left(u,v\right)}\end{array}$$
$$\begin{array}{rl}& F\left(u,v\right)\text{=}\frac{G\left(u,v\right)}{H\left(u,v\right)}+\frac{N\left(u,v\right)}{H\left(u,v\right)}\\ & =F^{\prime }\left(u,v\right)+\frac{N\left(u,v\right)}{H\left(u,v\right)}\end{array}$$
将其进行傅里叶反变换可以得到复原图像$\stackrel{⌢}{f}\left(x,y\right)=F^{-1}\left[F\left(u,v\right)\right]$
但是仅仅知道$H\left(u,v\right)$时，逆滤波复原存在一定的难度，而且逆滤波复原对信噪比的要求比较高。

维纳滤波

维纳滤波是另外一种比较常见的图像复原方法，计算复杂度相对较小并且考虑了噪声的影响。通常情况下，假设图像和噪声是相互独立的且至少存在一个均值为零，恢复图像和模糊图像两者的灰度级呈线性关系时，维纳滤波可表示为
$$F^{\prime }\left(u,v\right)=\left(\frac{1}{H\left(u,v\right)}\frac{|H\left(u,v\right){|}^2}{|H\left(u,v\right){|}^2+K}\right)G\left(u,v\right)$$
$$|H\left(u,v\right){|}^2=\bar{H}\left(u,v\right)H\left(u,v\right)$$
$$K=S_n\left(u,v\right)/S_f\left(u,v\right)$$
其中，$S_n\left(u,v\right)$为噪声功率谱，$S_f\left(u,v\right)$为原图像功率谱。当信噪比较大时，$K\to 0$，此时维纳滤波可转换为逆滤波；当信噪比较小时，$K\to +\infty$，图像存在难以复原的问题。

仿真结果

从图中可以看出，经过逆滤波后的效果不如维纳滤波效果好，同时也能看出由于只知道点扩散函数，不知道噪声对最终的结果有较大的影响，相比之下，维纳滤波表现出更好的性能。
代码如下
主函数

clear;
close all;
clc;
N = 256;
x = double(imread('lenna.jpg'));
figure(1)
subplot(221)
imshow(x,gray(256));title('原始图像');
h = ones(4,4)/16;  %点扩散函数
sigma = 15;         %噪声方差

Xf = fft2(x);
Hf = fft2(h,N,N);
y = real(ifft2(Hf.*Xf))+sigma*randn(N,N); %退化图像
%y = filter2(h,x)+sigma*randn(N,N);	  %退化图像
subplot(222)
imshow(y,gray(256));title('退化图像');
gamma = 2;
eix = inverseFilter(y,h,gamma);
subplot(223)
imshow(eix,gray(256));title('逆滤波复原');
%% 广义维纳滤波
gamma = 1;
alpha = 1;
ewx = wienerFilter(y,h,sigma,gamma,alpha);
subplot(224)
imshow(ewx,gray(256));title('维纳滤波复原');

inverseFilter.m

function restored_x = inverseFilter(y,h,gamma)
%% 广义逆滤波
N = size(y,1);  %大小
Yf = fft2(y);   %频域变换
Hf = fft2(h,N,N); %点扩散函数
sHf = Hf.*(abs(Hf)>0)+1/gamma*(abs(Hf)==0);  %处理为零的情况
iHf = 1./sHf;     %求逆
iHf = iHf.*(abs(Hf)*gamma>1)+gamma*abs(sHf).*iHf.*(abs(sHf)*gamma<=1);
restored_x = real(ifft2(iHf.*Yf));   %复原图像
end

wienerFilter.m

function ex = wienerFilter(y,h,sigma,gamma,alpha)
%广义维纳滤波，当alpha=1时退化为维纳滤波

N = size(y,1);
Yf = fft2(y);
Hf = fft2(h,N,N);
Pyf = abs(Yf).^2/N^2;
sHf = Hf.*(abs(Hf)>0)+1/gamma*(abs(Hf)==0);
iHf = 1./sHf;
iHf = iHf.*(abs(Hf)*gamma>1)+gamma*abs(sHf).*iHf.*(abs(sHf)*gamma<=1);

Pyf = Pyf.*(Pyf>sigma^2)+sigma^2*(Pyf<=sigma^2);
Gf = iHf.*(Pyf-sigma^2)./(Pyf-(1-alpha)*sigma^2);
eXf = Gf.*Yf;
ex = real(ifft2(eXf));

end