图像降噪算法——维纳滤波_维纳滤波算法

图像降噪算法——维纳滤波

维纳滤波是在频域中处理图像的一种算法，是一种非常经典的图像增强算法，不仅可以进行图像降噪，还可以消除由于运动等原因带来的图像模糊。

1. 基本原理

在图像拍摄过程中由于各种原因会造成图像退化，图像退化模型如下：$$g(x,y)=h(x,y)\star f(x,y)+\eta (x,y)$$其中，$\star$为卷积符号，$f(x,y)$为输入图像，$g(x,y)$为退化图像，$h(x,y)$为退化函数，$\eta (x,y)$为加性噪声，将上式进行傅里叶变换有：$$G(u,v)=H(u,v)F(u,v)+N(u,v)$$根据傅里叶变换的特性，空间域中的卷积相当于频率域中的乘积。

(1) 如果不考虑退化函数，图像退化模型就简化为图像噪声模型：$$g(x,y)=f(x,y)+\eta (x,y)$$图像增强问题成为单纯的图像去噪问题，可以通过空间域滤波等众多方法解决。

(2) 如果不考虑加性噪声，图像退化模型就简化为：$$g(x,y)=h(x,y)\star f(x,y)$$这种问题可以通过逆滤波解决，即通过傅里叶变化以及阵列除法即可获得恢复后的图像频谱$\hat{F}(u,v)$：$$\hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$$那么$H(u,v)$怎么获得呢？《数字图像处理》中的方法有观察估计、实验估计和建模估计，例如建模估计中可以通过运动数学模型将退化函数构造为：$$H(u,v)=\frac{T}{\pi (ua+vb)}\sin [\pi (ua+vb)]{\mathrm{e}}^{-\mathrm{j}\pi (ua+vb)}$$

(3) 如果退化函数和加性噪声都考虑，空域滤波器无法解决图像退化问题，逆滤波效果因为噪声的存在会变得非常差，这个时候就需要用到维纳滤波，（维纳滤波的推导写在结论中）维纳滤波公式如下：$$\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)$$其中，
$S_{\eta }(u,v)=|N(u,v){|}^2$为噪声的功率谱，这个我们可以通过用户输入的方差构造一个噪声图像$N(u,v)$并计算功率谱；
$S_f(u,v)=|F(u,v){|}^2$为输入图像的功率谱，这里乍一看会觉得有点问题，我们如果知道输入图像还需要滤波干嘛？我们当然不知道输入图像，因为真实图像的功率谱都是类似的，因此我们使用一个参考图像计算功率谱即可，在下面的例子中就是使用的lena的灰度图作为参考图像；
下面的例子中我们还对退化函数进行了简化，将退化函数置为1，因此维纳滤波公式简化为：$$\hat{F}(u,v)=\left[\frac{S_f(u,v)}{S_f(u,v)+S_{\eta }(u,v)}\right]G(u,v)$$

2. C++代码实现

下面基于OpenCV实现维纳滤波

Mat Denoise::WienerFilter(const Mat &src, const Mat &ref, int stddev)
{
    //这些图片是过程中会用到的，pad是原图像0填充后的图像，cpx是双通道频域图，mag是频域幅值图，dst是滤波后的图像
    Mat pad, cpx, dst;

    //获取傅里叶变化最佳图片尺寸，为2的指数
    int m = getOptimalDFTSize(src.rows);
    int n = getOptimalDFTSize(src.cols);

    //对原始图片用0进行填充获得最佳尺寸图片
    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);

    //对src进行傅里叶变换
    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]);

    //重新合并实部planes[0]和虚部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;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

下面是运行结果：
首先是原图：

添加上高斯噪声后：

维纳滤波的效果：

这个效果看上去有点神奇…

3. 结论

1. 从滤波结果看，和频域的高斯低通滤波是不同的，从频域图上看其实就是筛掉了不同的高频信号

2. 咱学知识就要学明白，这里我们把维纳滤波的推导再过一遍：
   维纳滤波是采用优化的方法推导的的，优化的目标是是的估计图像$\hat{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)$就是我们需要估计的维纳滤波系数，然后对平方进行展开 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){]}^{\ast }{S}_f(u,v)+X(u,v)X^{\ast }(u,v)(f)S_{\eta }(u,v)$$然后对$X(u,v)$求导$$\frac{\mathrm{d}(e)}{\mathrm{d}X(u,v)}=X^{\ast }(u,v)S_{\eta }(u,v)-X(u,v)[1-X(u,v)H(u,v)){]}^{\ast }{S}_f(u,v)=0$$求解可获得维纳滤波公式。

有问题欢迎交流～

此外，这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结，对图像降噪算法感兴趣的同学欢迎参考