图像降噪算法——高斯低通滤波

1. 基本原理

通过离散傅里叶变换对图像进行滤波流程作非常简单，就如下图所示：

(1) 二维离散傅里叶变换的公式为：F(u,v)=∑x=0M−1∑y=0N−1f(x,y)e−j2π(ux/M+vy/N)F(u, v)=\sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x, y) \mathrm{e}^{-\mathrm{j} 2 \pi(u x / M+v y / N)}F(u,v)=x=0∑M−1y=0∑N−1f(x,y)e−j2π(ux/M+vy/N)其中，xxx和yyy分别为原始图像上的横纵坐标，MMM和NNN为为原始图像的高和宽，f(x,y)f(x,y)f(x,y)为原始图像上坐标(x,y)(x,y)(x,y)处的灰度值，其中，根据欧拉公式：e−j2π(ux/M+vy/N)=cos(2π(uxM+vyN))−jsin(2π(uxM+vyN))\mathrm{e}^{-\mathrm{j} 2 \pi(u x / M+v y / N)}=cos(2 \pi(\frac{u x} {M}+\frac{v y}{N}))-jsin(2 \pi(\frac{u x} {M}+\frac{v y}{N}))e−j2π(ux/M+vy/N)=cos(2π(Mux+Nvy))−jsin(2π(Mux+Nvy))因此原始图像通过二维离散傅里叶变换获得的二维频谱上每一个像素都是复数，在OpenCV中用Mat数据格式处理时，每个像素值都有两个通道。

(2) 二维离散傅里叶反变换的公式为：f(x,y)=1MN∑u=0M−1∑v=0N−1F(u,v)ej2π(ux/M+vy/N)f(x, y)=\frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) \mathrm{e}^{\mathrm{j} 2 \pi(u x / M+v y / N)}f(x,y)=MN1u=0∑M−1v=0∑N−1F(u,v)ej2π(ux/M+vy/N)其中，符号表示和上述二维离散傅里叶变换相同。

(3) 对二维频谱的各种操作的具体方式：将二维频谱与设计好的滤波模板进行乘积，为什么是乘积呢？

因为根据傅里叶变换的性质，频率域中的乘积相当于空间域中的卷积，即f(x,y)⋆h(x,y)⇔F(u,v)H(u,v)f(x, y) \star h(x, y) \Leftrightarrow F(u, v) H(u, v)f(x,y)⋆h(x,y)⇔F(u,v)H(u,v)，其中二维循环卷积公式如下f(x,y)⋆h(x,y)=∑m=0M−1∑n=0N−1f(m,n)h(x−m,y−n)f(x, y) \star h(x, y)=\sum_{m=0}^{M-1} \sum_{n=0}^{N-1} f(m, n) h(x-m, y-n)f(x,y)⋆h(x,y)=m=0∑M−1n=0∑N−1f(m,n)h(x−m,y−n)以高斯低通滤波为例，即H(u,v)=e−D2(u,v)/2σ2H(u, v)=\mathrm{e}^{-D^{2}(u, v) / 2 \sigma^{2}}H(u,v)=e−D2(u,v)/2σ2其中D(u,v)D(u,v)D(u,v)为二维频谱上像素(u,v)(u,v)(u,v)距离中心的距离，其通过傅里叶反变换得到：h(x,y)=2πσe−2π2σ2D2(x,y)h(x,y)=\sqrt{2 \pi} \sigma \mathrm{e}^{-2 \pi^{2} \sigma^{2} D^{2}(x,y)}h(x,y)=2πσe−2π2σ2D2(x,y)其中D(x,y)D(x,y)D(x,y)为空间域图像上像素(x,y)(x,y)(x,y)距离中心的距离。可以发现空间域中同样是一个高斯滤波器，两者之间的关系是频率中高斯低通滤波器越窄，说明其衰减的低频越多，引起的模糊就越大，在空间域中就相当于使用了更大的模板来增加模糊。

2. C++代码实现

下面是基于OpenCV实现的频域的高斯低通滤波器

//高斯低通滤波器（频域）
Mat Denoise::GaussianLowPassFilter(const Mat &src, double sigma)
{//这些图片是过程中会用到的，pad是原图像0填充后的图像，cpx是双通道频域图，mag是频域幅值图，dst是滤波后的图像Mat pad, cpx, mag, 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 gaussian(pad.size(),CV_32FC2);for(int i = 0; i<m; i++){float* p = gaussian.ptr<float>(i);for(int j = 0; j<n; j++){double d = pow(i-m/2, 2) + pow(j-n/2, 2);p[2*j] = expf(-d/sigma/sigma/2.0);p[2*j+1] = expf(-d/sigma/sigma/2.0);}}//建立双通道图片，其中planes[0]填充原始图片Mat planes[] = {Mat_<float>(pad), Mat::zeros(pad.size(), CV_32F)};merge(planes, 2, cpx);//进行傅里叶变换dft(cpx, cpx);//分离通道并进行象限变幻split(cpx, planes);planes[0] = ShiftQuadrant(planes[0]);planes[1] = ShiftQuadrant(planes[1]);
//    magnitude(planes[0], planes[1], mag);
//    mag += Scalar::all(1);
//    log(mag, mag);
//    normalize(mag, mag, 0, 1, CV_MINMAX);
//    imshow("mag1", mag);//进行滤波merge(planes, 2, cpx);multiply(cpx, gaussian, cpx);//分离通道并进行象限变幻split(cpx, planes);//计算幅值，并讲幅值存储再planes[0]中magnitude(planes[0], planes[1], mag);
//    mag += Scalar::all(1);
//    log(mag, mag);
//    normalize(mag, mag, 0, 1, CV_MINMAX);
//    imshow("mag2", mag);planes[0] = ShiftQuadrant(planes[0]);planes[1] = ShiftQuadrant(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::ShiftQuadrant(const Mat &src)
{// 交换前// ×××××××××××××××××××××××// ×    q1    ×    q2    ×// ×          ×          ×// ×××××××××××××××××××××××// ×    q3    ×    q4    ×// ×          ×          ×// ×××××××××××××××××××××××// 交换后// ×××××××××××××××××××××××// ×    q4    ×    q3    ×// ×          ×          ×// ×××××××××××××××××××××××// ×    q2    ×    q1    ×// ×          ×          ×// ×××××××××××××××××××××××Mat dst = src.clone();int m = src.rows, n = src.cols;Mat q1(dst, Rect(0,0,n/2,m/2));Mat q2(dst, Rect(n/2,0,n/2,m/2));Mat q3(dst, Rect(0,m/2,n/2,m/2));Mat q4(dst, Rect(n/2,m/2,n/2,m/2));//交换象限Mat temp;q1.copyTo(temp);q4.copyTo(q1);temp.copyTo(q4);q2.copyTo(temp);q3.copyTo(q2);temp.copyTo(q3);return dst;
}

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

添加上高斯噪声后：
进行二维离散傅里叶变换获得二维频谱（下面显示的幅值图mag，仅仅是为了演示频域滤波的过程，实际上用来操作的是两通道的频谱图cpx）：

对二维频谱进行操作，乘以高斯低通滤波器：

最后通过二维离散傅里叶反变换获得降噪图像：

3. 结论

本文是直接上手二维离散傅里叶变换进行滤波，啥是离散傅里叶变换？怎么就变换了？不太熟悉的同学建议复习下《数字信号处理》或者到网上看一些大佬的博客，这里推荐几篇：
1. 从菜鸟到完全学会二维傅立叶在图像处理算法中的应用【老司机教你学傅立叶】
2. 为什么用图像二维傅里叶变换的相位谱进行反变换，能够大致得到原图的形状，而幅度谱则不行呢？
3. 十分简明易懂的FFT（快速傅里叶变换）
4. 如果看了这篇文章你还不懂傅里叶变换，那就过来掐死我吧
本文介绍的二维离散傅里叶变换是最基本的形式，而实际应用中采用的肯定是快速傅里叶变换。
采用离散傅里叶进行滤波会有更直观，但是由于空间域可以进行并行计算，因此实际滤波是还是更加倾向于空间域滤波。

暂时分享这么多，有问题欢迎交流～
此外，这里我写一个各种算法的总结目录图像降噪算法——图像降噪算法总结，对图像降噪算法感兴趣的同学欢迎参考