图像降噪算法——维纳滤波
维纳滤波是在频域中处理图像的一种算法,是一种非常经典的图像增强算法,不仅可以进行图像降噪,还可以消除由于运动等原因带来的图像模糊。
1. 基本原理
在图像拍摄过程中由于各种原因会造成图像退化,图像退化模型如下:
g
(
x
,
y
)
=
h
(
x
,
y
)
⋆
f
(
x
,
y
)
+
η
(
x
,
y
)
g(x, y)=h(x, y) \star f(x, y)+\eta(x, y)
g(x,y)=h(x,y)⋆f(x,y)+η(x,y)其中,
⋆
\star
⋆为卷积符号,
f
(
x
,
y
)
f(x,y)
f(x,y)为输入图像,
g
(
x
,
y
)
g(x,y)
g(x,y)为退化图像,
h
(
x
,
y
)
h(x,y)
h(x,y)为退化函数,
η
(
x
,
y
)
\eta(x,y)
η(x,y)为加性噪声,将上式进行傅里叶变换有:
G
(
u
,
v
)
=
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
G(u, v)=H(u, v) F(u, v)+N(u, v)
G(u,v)=H(u,v)F(u,v)+N(u,v)根据傅里叶变换的特性,空间域中的卷积相当于频率域中的乘积。
(1) 如果不考虑退化函数,图像退化模型就简化为图像噪声模型:
g
(
x
,
y
)
=
f
(
x
,
y
)
+
η
(
x
,
y
)
g(x, y)=f(x, y)+\eta(x, y)
g(x,y)=f(x,y)+η(x,y)图像增强问题成为单纯的图像去噪问题,可以通过空间域滤波等众多方法解决。
(2) 如果不考虑加性噪声,图像退化模型就简化为:
g
(
x
,
y
)
=
h
(
x
,
y
)
⋆
f
(
x
,
y
)
g(x, y)=h(x, y) \star f(x, y)
g(x,y)=h(x,y)⋆f(x,y)这种问题可以通过逆滤波解决,即通过傅里叶变化以及阵列除法即可获得恢复后的图像频谱
F
^
(
u
,
v
)
\hat{F}(u, v)
F^(u,v):
F
^
(
u
,
v
)
=
G
(
u
,
v
)
H
(
u
,
v
)
\hat{F}(u, v)=\frac{G(u, v)}{H(u, v)}
F^(u,v)=H(u,v)G(u,v)那么
H
(
u
,
v
)
H(u,v)
H(u,v)怎么获得呢?《数字图像处理》中的方法有观察估计、实验估计和建模估计,例如建模估计中可以通过运动数学模型将退化函数构造为:
H
(
u
,
v
)
=
T
π
(
u
a
+
v
b
)
sin
[
π
(
u
a
+
v
b
)
]
e
−
j
π
(
u
a
+
v
b
)
H(u, v)=\frac{T}{\pi(u a+v b)} \sin [\pi(u a+v b)] \mathrm{e}^{-\mathrm{j} \pi(u a+v b)}
H(u,v)=π(ua+vb)Tsin[π(ua+vb)]e−jπ(ua+vb)
(3) 如果退化函数和加性噪声都考虑,空域滤波器无法解决图像退化问题,逆滤波效果因为噪声的存在会变得非常差,这个时候就需要用到维纳滤波,(维纳滤波的推导写在结论中)维纳滤波公式如下:
F
^
(
u
,
v
)
=
[
1
H
(
u
,
v
)
∣
H
(
u
,
v
)
∣
2
∣
H
(
u
,
v
)
∣
2
+
S
η
(
u
,
v
)
/
S
f
(
u
,
v
)
]
G
(
u
,
v
)
\hat{F}(u, v)=\left[\frac{1}{H(u, v)} \frac{|H(u, v)|^{2}}{|H(u, v)|^{2}+S_{\eta}(u, v) / S_{f}(u, v)}\right] G(u, v)
F^(u,v)=[H(u,v)1∣H(u,v)∣2+Sη(u,v)/Sf(u,v)∣H(u,v)∣2]G(u,v)其中,
S
η
(
u
,
v
)
=
∣
N
(
u
,
v
)
∣
2
S_{\eta}(u, v)=|N(u, v)|^{2}
Sη(u,v)=∣N(u,v)∣2为噪声的功率谱,这个我们可以通过用户输入的方差构造一个噪声图像
N
(
u
,
v
)
N(u, v)
N(u,v)并计算功率谱;
S
f
(
u
,
v
)
=
∣
F
(
u
,
v
)
∣
2
S_{f}(u, v)=|F(u, v)|^{2}
Sf(u,v)=∣F(u,v)∣2为输入图像的功率谱,这里乍一看会觉得有点问题,我们如果知道输入图像还需要滤波干嘛?我们当然不知道输入图像,因为真实图像的功率谱都是类似的,因此我们使用一个参考图像计算功率谱即可,在下面的例子中就是使用的lena的灰度图作为参考图像;
下面的例子中我们还对退化函数进行了简化,将退化函数置为1,因此维纳滤波公式简化为:
F
^
(
u
,
v
)
=
[
S
f
(
u
,
v
)
S
f
(
u
,
v
)
+
S
η
(
u
,
v
)
]
G
(
u
,
v
)
\hat{F}(u, v)=\left[\frac{S_{f}(u, v)}{S_{f}(u, v)+S_{\eta}(u, v)}\right] G(u, v)
F^(u,v)=[Sf(u,v)+Sη(u,v)Sf(u,v)]G(u,v)
2. C++代码实现
下面基于OpenCV实现维纳滤波
Mat Denoise::WienerFilter(const Mat &src, const Mat &ref, int stddev)
{
Mat pad, cpx, dst;
int m = getOptimalDFTSize(src.rows);
int n = getOptimalDFTSize(src.cols);
copyMakeBorder(src, pad, 0, m-src.rows, 0, n-src.cols, BORDER_CONSTANT, Scalar::all(0));
Mat tmpR(pad.rows, pad.cols, CV_8U);
resize(ref, tmpR, tmpR.size());
Mat refSpectrum = GetSpectrum(tmpR);
Mat tmpN(pad.rows, pad.cols, CV_32F);
randn(tmpN, Scalar::all(0), Scalar::all(stddev));
Mat noiseSpectrum = GetSpectrum(tmpN);
Mat planes[] = {Mat_<float>(pad), Mat::zeros(pad.size(), CV_32F)};
merge(planes, 2, cpx);
dft(cpx, cpx);
split(cpx, planes);
Mat factor = refSpectrum / (refSpectrum + noiseSpectrum);
multiply(planes[0], factor, planes[0]);
multiply(planes[1], factor, planes[1]);
merge(planes, 2, cpx);
idft(cpx, dst, DFT_SCALE | DFT_REAL_OUTPUT);
dst.convertTo(dst, CV_8UC1);
return dst;
}
Mat Denoise::GetSpectrum(const Mat &src)
{
Mat dst, cpx;
Mat planes[] = {Mat_<float>(src), Mat::zeros(src.size(), CV_32F)};
merge(planes, 2, cpx);
dft(cpx, cpx);
split(cpx, planes);
magnitude(planes[0], planes[1], dst);
multiply(dst, dst, dst);
return dst;
}
下面是运行结果:
首先是原图:
添加上高斯噪声后:
维纳滤波的效果:
这个效果看上去有点神奇…
3. 结论
-
从滤波结果看,和频域的高斯低通滤波是不同的,从频域图上看其实就是筛掉了不同的高频信号
-
咱学知识就要学明白,这里我们把维纳滤波的推导再过一遍:
维纳滤波是采用优化的方法推导的的,优化的目标是是的估计图像
F
^
(
u
,
v
)
\hat{F}(u,v)
F^(u,v)和参考图像
F
(
u
,
v
)
F(u,v)
F(u,v)(或者说是输入图像)频谱的均方差最小,即:
e
=
E
{
∣
F
(
u
,
v
)
−
F
^
(
u
,
v
)
∣
2
}
e=E\{|F(u,v)-\hat{F}(u,v)|^{2}\}
e=E{∣F(u,v)−F^(u,v)∣2}将估计图像频谱进行替换:
e
=
E
{
∣
F
(
u
,
v
)
−
X
(
u
,
v
)
G
(
u
,
v
)
∣
2
}
=
E
{
∣
F
(
u
,
v
)
−
X
(
u
,
v
)
[
H
(
u
,
v
)
F
(
u
,
v
)
+
N
(
u
,
v
)
]
∣
2
}
=
E
{
∣
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
F
(
u
,
v
)
−
X
(
u
,
v
)
N
(
u
,
v
)
∣
2
}
\begin{aligned} e &=E\{|F(u,v)-X(u,v) G(u,v)|^{2}\} \\ &=E\{|F(u,v)-X(u,v) [H(u, v) F(u, v)+N(u, v)]|^{2}\} \\ &=E\{|[1-X(u,v)H(u, v) ] F(u, v)-X(u,v)N(u, v)|^{2}\} \end{aligned}
e=E{∣F(u,v)−X(u,v)G(u,v)∣2}=E{∣F(u,v)−X(u,v)[H(u,v)F(u,v)+N(u,v)]∣2}=E{∣[1−X(u,v)H(u,v)]F(u,v)−X(u,v)N(u,v)∣2}其中,
X
(
u
,
v
)
X(u,v)
X(u,v)就是我们需要估计的维纳滤波系数,然后对平方进行展开
e
=
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
∗
×
E
{
∣
F
(
u
,
v
)
∣
2
}
−
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
X
∗
(
u
,
v
)
×
E
{
F
(
u
,
v
)
N
∗
(
u
,
v
)
}
−
X
(
u
,
v
)
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
∗
×
E
{
F
(
u
,
v
)
∗
N
(
u
,
v
)
}
+
X
(
u
,
v
)
X
∗
(
u
,
v
)
×
E
{
∣
N
(
u
,
v
)
∣
2
}
\begin{aligned} e &=[1-X(u,v)H(u, v) ][1-X(u,v)H(u, v) ]^*×E\{|F(u,v)|^{2}\} \\ &-[1-X(u,v)H(u, v) ]X^*(u,v)× E\left\{F(u,v)N^*(u, v)\right\} \\ &-X(u,v)[1-X(u,v)H(u, v)]^*×E\left\{F(u,v)^*N(u, v)\right\} \\ &+X(u,v)X^*(u,v)×E\{|N(u, v)|^{2}\} \end{aligned}
e=[1−X(u,v)H(u,v)][1−X(u,v)H(u,v)]∗×E{∣F(u,v)∣2}−[1−X(u,v)H(u,v)]X∗(u,v)×E{F(u,v)N∗(u,v)}−X(u,v)[1−X(u,v)H(u,v)]∗×E{F(u,v)∗N(u,v)}+X(u,v)X∗(u,v)×E{∣N(u,v)∣2}其中,由于噪声和信号是独立无关的,因此
E
{
F
(
u
,
v
)
N
(
u
,
v
)
}
=
E
{
F
(
u
,
v
)
N
(
u
,
v
)
}
=
0
E\left\{F(u,v)N(u, v)\right\}=E\left\{F(u,v) N(u, v)\right\}=0
E{F(u,v)N(u,v)}=E{F(u,v)N(u,v)}=0定义如下功率谱有:
S
f
(
u
,
v
)
=
E
{
∣
F
(
u
,
v
)
∣
2
}
S
η
(
u
,
v
)
=
E
{
∣
N
(
u
,
v
)
∣
2
}
\begin{array}{l} S_f(u,v)=E\{|F(u,v)|^{2}\} \\ S_\eta(u,v)=E\{|N(u, v)|^{2}\} \end{array}
Sf(u,v)=E{∣F(u,v)∣2}Sη(u,v)=E{∣N(u,v)∣2}于是有:
e
=
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
]
∗
S
f
(
u
,
v
)
+
X
(
u
,
v
)
X
∗
(
u
,
v
)
(
f
)
S
η
(
u
,
v
)
e=[1-X(u,v)H(u, v) ][1-X(u,v)H(u, v) ]^*S_f(u,v)+X(u,v)X^*(u,v)(f) S_\eta(u,v)
e=[1−X(u,v)H(u,v)][1−X(u,v)H(u,v)]∗Sf(u,v)+X(u,v)X∗(u,v)(f)Sη(u,v)然后对
X
(
u
,
v
)
X(u,v)
X(u,v)求导
d
(
e
)
d
X
(
u
,
v
)
=
X
∗
(
u
,
v
)
S
η
(
u
,
v
)
−
X
(
u
,
v
)
[
1
−
X
(
u
,
v
)
H
(
u
,
v
)
)
]
∗
S
f
(
u
,
v
)
=
0
\frac{\mathrm{d}(e)}{\mathrm{d} X(u,v)}=X^{*}(u,v)S_\eta(u,v)-X(u,v)[1-X(u,v)H(u, v))]^{*}S_f(u,v)=0
dX(u,v)d(e)=X∗(u,v)Sη(u,v)−X(u,v)[1−X(u,v)H(u,v))]∗Sf(u,v)=0求解可获得维纳滤波公式。
有问题欢迎交流~
此外,这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结,对图像降噪算法感兴趣的同学欢迎参考
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)