由于受到光学系统的像差,成像设备与物体的相对运动等因素的影响,图像会出现一定的失真。要想得到高质量的图像,需要对已经退化后的图像进行复原。
图像退化模型
一般来说,图像的退化模型可以表示为
g
(
x
,
y
)
=
h
(
x
,
y
)
∗
f
(
x
,
y
)
+
n
(
x
,
y
)
g\left( x,y \right)=h\left( x,y \right)*f\left( x,y \right)+n\left( x,y \right)
g(x,y)=h(x,y)∗f(x,y)+n(x,y)
其中,
g
(
x
,
y
)
g\left( x,y \right)
g(x,y)表示退化后的图像,
h
(
x
,
y
)
h\left( x,y \right)
h(x,y)为点扩散函数,
f
(
x
,
y
)
f\left( x,y \right)
f(x,y)为原始图像,
n
(
x
,
y
)
n\left( x,y \right)
n(x,y)为引入的噪声。
在频域上面可以表示为
G
(
u
,
v
)
=
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
G\left( u,v \right)=H\left( u,v \right)F\left( u,v \right)+N\left( u,v \right)
G(u,v)=H(u,v)F(u,v)+N(u,v)
逆滤波
逆滤波是一种常见且直观的图像恢复方法,它的主要过程是将退化后的图像从空间域变换到频域,进行逆滤波后在变换到空间域上从而实现图像的复原。
从上面的介绍可以得到,图像的退化模型为
g
(
x
,
y
)
=
h
(
x
,
y
)
∗
f
(
x
,
y
)
+
n
(
x
,
y
)
=
∫
−
∞
+
∞
f
(
α
,
β
)
h
(
x
−
α
,
y
−
β
)
d
α
d
β
+
n
(
x
,
y
)
\begin{aligned} & g\left( x,y \right)=h\left( x,y \right)*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{aligned}
g(x,y)=h(x,y)∗f(x,y)+n(x,y)=∫−∞+∞f(α,β)h(x−α,y−β)dαdβ+n(x,y)
将其变换到频域可以得到
G
(
u
,
v
)
=
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
G\left( u,v \right)=H\left( u,v \right)F\left( u,v \right)+N\left( u,v \right)
G(u,v)=H(u,v)F(u,v)+N(u,v)
那么
F
′
(
u
,
v
)
=
G
(
u
,
v
)
H
(
u
,
v
)
−
N
(
u
,
v
)
H
(
u
,
v
)
=
F
(
u
,
v
)
−
N
(
u
,
v
)
H
(
u
,
v
)
\begin{aligned} & {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\left( u,v \right)-\frac{N\left( u,v \right)}{H\left( u,v \right)} \end{aligned}
F′(u,v)=H(u,v)G(u,v)−H(u,v)N(u,v)=F(u,v)−H(u,v)N(u,v)
F
(
u
,
v
)
=
G
(
u
,
v
)
H
(
u
,
v
)
+
N
(
u
,
v
)
H
(
u
,
v
)
=
F
′
(
u
,
v
)
+
N
(
u
,
v
)
H
(
u
,
v
)
\begin{aligned} & {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'\left( u,v \right)+\frac{N\left( u,v \right)}{H\left( u,v \right)} \end{aligned}
F(u,v)=H(u,v)G(u,v)+H(u,v)N(u,v)=F′(u,v)+H(u,v)N(u,v)
将其进行傅里叶反变换可以得到复原图像
f
⌢
(
x
,
y
)
=
F
−
1
[
F
(
u
,
v
)
]
\overset{\scriptscriptstyle\frown}{f}\left( x,y \right)={{F}^{-1}}\left[ {F}\left( u,v \right) \right]
f⌢(x,y)=F−1[F(u,v)]
但是仅仅知道
H
(
u
,
v
)
H\left( u,v \right)
H(u,v)时,逆滤波复原存在一定的难度,而且逆滤波复原对信噪比的要求比较高。
维纳滤波
维纳滤波是另外一种比较常见的图像复原方法,计算复杂度相对较小并且考虑了噪声的影响。通常情况下,假设图像和噪声是相互独立的且至少存在一个均值为零,恢复图像和模糊图像两者的灰度级呈线性关系时,维纳滤波可表示为
F
′
(
u
,
v
)
=
(
1
H
(
u
,
v
)
∣
H
(
u
,
v
)
∣
2
∣
H
(
u
,
v
)
∣
2
+
K
)
G
(
u
,
v
)
{F}'\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)
F′(u,v)=(H(u,v)1∣H(u,v)∣2+K∣H(u,v)∣2)G(u,v)
∣
H
(
u
,
v
)
∣
2
=
H
ˉ
(
u
,
v
)
H
(
u
,
v
)
|H\left( u,v \right){{|}^{2}}=\bar{H}\left( u,v \right)H\left( u,v \right)
∣H(u,v)∣2=Hˉ(u,v)H(u,v)
K
=
S
n
(
u
,
v
)
/
S
f
(
u
,
v
)
K={{S}_{n}}\left( u,v \right)/{{S}_{f}}\left( u,v \right)
K=Sn(u,v)/Sf(u,v)
其中,
S
n
(
u
,
v
)
{{S}_{n}}\left( u,v \right)
Sn(u,v)为噪声功率谱,
S
f
(
u
,
v
)
{{S}_{f}}\left( u,v \right)
Sf(u,v)为原图像功率谱。当信噪比较大时,
K
→
0
K\to 0
K→0,此时维纳滤波可转换为逆滤波;当信噪比较小时,
K
→
+
∞
K\to +\infty
K→+∞,图像存在难以复原的问题。
仿真结果
从图中可以看出,经过逆滤波后的效果不如维纳滤波效果好,同时也能看出由于只知道点扩散函数,不知道噪声对最终的结果有较大的影响,相比之下,维纳滤波表现出更好的性能。
代码如下
主函数
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