1、空间滤波基础概念
1、空间滤波基础
空间滤波一词中滤波取自数字信号处理，指接受或拒绝一定的频率成分，但是空间滤波学习内容实际上和通过傅里叶变换实现的频域的滤波是等效的，故而也称为滤波。空间滤波主要直接基于领域（空间域）对图像中的像素执行计算，用滤波器（也成为空间掩膜、核、模板和窗口）直接作用于图像本身完成类似的平滑。
2、空间滤波机理
对空间域中的每一点(x,y)，重复如下操作：
- 对预先定义的以（x,y)为中心的领域内的像素进行预定义运算。
- 将（1）中运算的结果作为(x,y)点新的响应。
上述过程就称为邻域处理或空间域滤波。
3、空间滤波分类
根据预定义的操作，可以将滤波器分为：
- 线性滤波器
- 非线性滤波器
而根据滤波器最终对图像造成的影响，可以将滤波器分为：
- 平滑滤波器 ，通常用于模糊图像或者去除图像中的噪声
- 锐化滤波器，突出图像中的边缘细节部分
4、空间相关和卷积
在执行线性空间vlbo时，有两个相近的概念需要理解，相关和卷积。相关是滤波器模板移过图像并计算每个位置乘积之和的处理。卷积机理类似，但是滤波器首先需要旋转180°。关于卷积的详细说明，可以参考文章
卷积在图像处理中的应用（转自https://medium.com/@irhumshafkat/intuitively-understanding-convolutions-for-deep-learning-1f6f42faee1）
一些算子及不同类型的卷积
这里不再做重复说明。
2、为图片添加噪声
使用摄像头拍摄图片，噪声是难以避免的，这里重点提下两种噪声：
- 椒盐噪声：噪声的幅值基本上相同，但是噪声出现的位置是随机的；（中值滤波效果好）
- 高斯噪声：每一点都存在噪声，但噪声的幅值是随机分布的。
但是，我们网上找到的图片大多是经过处理的，或者说噪声不是明显的，这里为了研究，为图片添加噪声。
椒盐噪声表现在图片上是图片中离散分布的黑点和白点
// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
高斯噪声是一种加性噪声，为图像添加高斯噪声的代码如下：
// 添加Gussia噪声
// 使用指针访问
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}示例：
为图像添加两种噪声的程序如下：
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>
using namespace cv;
using namespace std;
// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}int main()
{Mat srcImage = imread("111.jpg");Mat src_salt, src_gaussian;src_salt = imread("111.jpg");src_gaussian = imread("111.jpg");if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);imshow("原图", srcImage);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);waitKey();return 0;
}添加完噪声后的图像如下：
3、线性平滑空间滤波
平滑滤波器主要用于模糊处理和降低噪音。
1、平滑线性滤波器(又称均值滤波）
首先需要说明的是图像平滑是一种减少和抑制图像噪声的实用数字图像处理技术，一般来说，图像具有局部连续性质，即相邻像素的数值相近，而噪声的存在使得在噪声点处产生灰度跳跃，但一般可以合理的假设偶尔出现的噪声并没有改变图像局部连续的性质。其中，线性滤波是现实将空域模板对应邻域内像素的灰度值加权之和作为邻域内中心像素的相应输出。线性平滑模板的权系数全为正值而且系数之和等于1，如下图所示：
因此这种方法不会增加图像中总体的灰度程度。也就是说在灰度一致的区域，线性平滑滤波的相应输出不变。此外，可以知道，线性平滑滤波等效于低通滤波。但是同样的我们可以看出，线性平滑模板在模糊和降噪的同时，图像中的边缘和细节的锐度也都丢失了，也就是说，在平滑的过程中，使得图像的一部分细节信息丢失。
在OpenCV中，函数blur表示使用该模板的均值滤波器，其声明如下：
void blur( InputArray src, OutputArray dst, Size ksize, Point anchor = Point(-1,-1), int borderType = BORDER_DEFAULT );src是输入图像，dst为输出图像；ksize是滤波器模板窗口的大小；后两个参数分别表示，待处理像素在模板窗口的位置，默认值是窗口的中心位置，所以窗口的大小一般为奇数，最后一个参数表示对编解类型的处理，使用默认值即可。其调用示例blur(src,dst,Size(5,5)，模板窗口的大小为5×55×5。
示例：
下面通过程序对两种含噪声的图片进行线性平滑滤波处理
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>using namespace cv;
using namespace std;// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}int main()
{Mat srcImage = imread("111.jpg");Mat src_salt, src_gaussian;src_salt = imread("111.jpg");src_gaussian = imread("111.jpg");if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);//平滑滤波Mat aft_mean_filter_salt, aft_mean_filter_gaussian;blur(src_salt, aft_mean_filter_salt, Size(5, 5), Point(-1, -1));blur(src_gaussian, aft_mean_filter_gaussian, Size(5, 5), Point(-1, -1));//增强图像对比度//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));imshow("原图", srcImage);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);imshow("经过平滑滤波后含椒盐噪声的的结果", aft_mean_filter_salt);imshow("经过平滑滤波后含高斯噪声的的结果", aft_mean_filter_gaussian);waitKey();return 0;
}
2、高斯平滑
均值平滑对于领域内的像素一视同仁，为了减少平滑处理中的模糊，得到更自然的平滑效果，我们可以很自然的想到适当加大模板中心点的权重，随着距离中心点的距离增加，权重迅速减小，从而可以确保中心点看起来更接近于与它距离更近的点，基于这样的考虑得到的模板即为高斯模板。常用的3X3高斯模板如下所示：
高斯模板名字的由来是高斯函数，即二维正态分布密度函数，一个二维的高斯函数如下,
$$
h(x, y)=e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}}
$$
公式中(x,y)是点坐标，在图像处理中可认为是整数，σ是标准差。
对于窗口模板的大小为 (2K+1)X(2K+1)的矩阵，其(i，j）位置的元素值可如下确定：
高斯模板的生成方法如下所示：
void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{static const double pi = 3.1415926;int center = ksize / 2; // 模板的中心位置，也就是坐标的原点double x2, y2;for (int i = 0; i < ksize; i++){x2 = pow(i - center, 2);for (int j = 0; j < ksize; j++){y2 = pow(j - center, 2);double g = exp(-(x2 + y2) / (2 * sigma * sigma));g /= 2 * pi * sigma;window[i][j] = g;}}double k = 1 / window[0][0]; // 将左上角的系数归一化为1for (int i = 0; i < ksize; i++){for (int j = 0; j < ksize; j++){window[i][j] *= k;}}
}
需要一个二维数组，存放生成的系数（这里假设模板的最大尺寸不会超过11）；第二个参数是模板的大小（不要超过11）；第三个参数就比较重要了，是高斯分布的标准差。
生成的过程，首先根据模板的大小，找到模板的中心位置ksize/2。 然后就是遍历，根据高斯分布的函数，计算模板中每个系数的值。
需要注意的是，最后归一化的过程，使用模板左上角的系数的倒数作为归一化的系数（左上角的系数值被归一化为1），模板中的每个系数都乘以该值（左上角系数的倒数），然后将得到的值取整，就得到了整数型的高斯滤波器模板。
至于小数形式的生成也比较简单，去掉归一化的过程，并且在求解过程后，模板的每个系数要除以所有系数的和。具体代码如下：
void generateGaussianTemplate(double window[][11], int ksize, double sigma)
{static const double pi = 3.1415926;int center = ksize / 2; // 模板的中心位置，也就是坐标的原点double x2, y2;double sum = 0;for (int i = 0; i < ksize; i++){x2 = pow(i - center, 2);for (int j = 0; j < ksize; j++){y2 = pow(j - center, 2);double g = exp(-(x2 + y2) / (2 * sigma * sigma));g /= 2 * pi * sigma;sum += g;window[i][j] = g;}}//double k = 1 / window[0][0]; // 将左上角的系数归一化为1for (int i = 0; i < ksize; i++){for (int j = 0; j < ksize; j++){window[i][j] /= sum;}}
}σ选取技巧
当标准差σ选取不同的值是，二维高斯函数形状会发生很大的变化，
- 如果σ过小，偏离中心的所有像素权杖重将会非常小，相当于加权和响应基本不考虑领域像素的作用，这样无法起到滤波平滑的效果。
- 如果σ过大，而领域相对较小，这样在领域内高斯模板则退化为平均模板。
在matlab中，σ的默认值为0.5，在实际应用中，通常3X3的模板取σ为0.8左右，更大的模板可以适当增加σ的值。
在OpenCV中使用GuassianBlur()函数来实现高斯平滑，函数的声明如下：
void GaussianBlur( InputArray src, OutputArray dst, Size ksize, double sigmaX, double sigmaY=0, int borderType=BORDER_DEFAULT );
这里由于高斯函数的可分离性，尺寸较大的高斯滤波器可以分成两步进行：首先将图像在水平（竖直）方向与一维高斯函数进行卷积；然后将卷积后的结果在竖直（水平）方向使用相同的一维高斯函数得到的模板进行卷积运算。
你也可以用如下源码进行实现：
void separateGaussianFilter(const Mat &src, Mat &dst, int ksize, double sigma)
{CV_Assert(src.channels()==1 || src.channels() == 3); // 只处理单通道或者三通道图像// 生成一维的高斯滤波模板double *matrix = new double[ksize];double sum = 0;int origin = ksize / 2;for (int i = 0; i < ksize; i++){// 高斯函数前的常数可以不用计算，会在归一化的过程中给消去double g = exp(-(i - origin) * (i - origin) / (2 * sigma * sigma));sum += g;matrix[i] = g;}// 归一化for (int i = 0; i < ksize; i++)matrix[i] /= sum;// 将模板应用到图像中int border = ksize / 2;copyMakeBorder(src, dst, border, border, border, border, BorderTypes::BORDER_REFLECT);int channels = dst.channels();int rows = dst.rows - border;int cols = dst.cols - border;// 水平方向for (int i = border; i < rows; i++){for (int j = border; j < cols; j++){double sum[3] = { 0 };for (int k = -border; k <= border; k++){if (channels == 1){sum[0] += matrix[border + k] * dst.at<uchar>(i, j + k); // 行不变，列变化；先做水平方向的卷积}else if (channels == 3){Vec3b rgb = dst.at<Vec3b>(i, j + k);sum[0] += matrix[border + k] * rgb[0];sum[1] += matrix[border + k] * rgb[1];sum[2] += matrix[border + k] * rgb[2];}}for (int k = 0; k < channels; k++){if (sum[k] < 0)sum[k] = 0;else if (sum[k] > 255)sum[k] = 255;}if (channels == 1)dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);else if (channels == 3){Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };dst.at<Vec3b>(i, j) = rgb;}}}// 竖直方向for (int i = border; i < rows; i++){for (int j = border; j < cols; j++){double sum[3] = { 0 };for (int k = -border; k <= border; k++){if (channels == 1){sum[0] += matrix[border + k] * dst.at<uchar>(i + k, j); // 列不变，行变化；竖直方向的卷积}else if (channels == 3){Vec3b rgb = dst.at<Vec3b>(i + k, j);sum[0] += matrix[border + k] * rgb[0];sum[1] += matrix[border + k] * rgb[1];sum[2] += matrix[border + k] * rgb[2];}}for (int k = 0; k < channels; k++){if (sum[k] < 0)sum[k] = 0;else if (sum[k] > 255)sum[k] = 255;}if (channels == 1)dst.at<uchar>(i, j) = static_cast<uchar>(sum[0]);else if (channels == 3){Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };dst.at<Vec3b>(i, j) = rgb;}}}delete[] matrix;
}示例：
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>using namespace cv;
using namespace std;// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}int main()
{Mat srcImage = imread("111.jpg");Mat src_salt, src_gaussian;src_salt = imread("111.jpg");src_gaussian = imread("111.jpg");if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);//平滑滤波Mat aft_mean_filter_salt, aft_mean_filter_gaussian;blur(src_salt, aft_mean_filter_salt, Size(5, 5), Point(-1, -1));blur(src_gaussian, aft_mean_filter_gaussian, Size(5, 5), Point(-1, -1));//增强图像对比度//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));imshow("原图", srcImage);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);imshow("经过平滑滤波后含椒盐噪声的的结果", aft_mean_filter_salt);imshow("经过平滑滤波后含高斯噪声的的结果", aft_mean_filter_gaussian);waitKey();return 0;
}
4、非线性平滑空间滤波
1、中值滤波
中值滤波本质上是一种统计排序滤波器，对于原图像中某点(i,j)，中值滤波以该点为中心的领域内的所有像素的统计排序中值作为(i,j)的响应，即它将每一像素点的灰度值设置为该点某邻域窗口内的所有像素点灰度值的中值，也就是将中心像素的值用所有像素值的中间值（不是平均值）替换。中值滤波不同于均值，是指排序队列中位于中间位置的元素的值。
中值滤波对于某些类型的随机噪声具有非常理想的降噪能力，对于线性平滑滤波而言，在处理的像素领域之内包含噪声点时，噪声的存在总会或多或少的影响该点的像素值的计算(对于高斯平滑，影响程度同噪声点到中点的距离成正比），但在中值滤波中，噪声点则常常是直接忽略掉的；而且在同线性平滑滤波器相比，中值滤波在降噪的同时引起的模糊效应较低。
OpenCV提供了中值滤波的API函数如下：
void medianBlur(InputArray src,OutputArray dst,int ksize)
参数解释：原图像，目标图像，模板尺寸。
示例：
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>using namespace cv;
using namespace std;// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}int main()
{Mat srcImage = imread("111.jpg");Mat src_salt, src_gaussian;src_salt = imread("111.jpg");src_gaussian = imread("111.jpg");if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);//平滑滤波Mat aft_median_filter_salt, aft_median_filter_gaussian;medianBlur(src_salt, aft_median_filter_salt, 3);medianBlur(src_gaussian, aft_median_filter_gaussian, 3);//增强图像对比度//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));imshow("原图", srcImage);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);imshow("经过平滑滤波后含椒盐噪声的的结果", aft_median_filter_salt);imshow("经过平滑滤波后含高斯噪声的的结果", aft_median_filter_gaussian);waitKey();return 0;
}
2、自适应中值滤波
中值滤波的效果依赖于滤波窗口的大小，太大会使图像边缘模糊，太小则去噪效果不好。因为噪声点和边缘点同样是灰度变化较为剧烈的像素，普通中值滤波在改变噪声点灰度值的时候，会一定程度上改变边缘像素的灰度值。但是噪声点几乎就是邻近像素的极值，而边缘往往不是，我们可以通过这个特性来限制中值滤波。
具体的改进方法如下：追行扫描图像， 当处理每一个像素时，判断该像素是否是滤波窗口所覆盖下领域像素的极大或者极小值，如果是，则采用正常的中值滤波处理该像素，如果不是，则不予处理。这种方法能够非常有效的去除突发噪声点，尤其是椒盐噪声，而且几乎不影响边缘。这种方法又称为自适应中值滤波。
算法实现如下所示：
uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize)
{vector<uchar> pixels;for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)for (int b = -kernelSize / 2; b <= kernelSize / 2; b++){pixels.push_back(im.at<uchar>(row + a, col + b));}sort(pixels.begin(), pixels.end());auto min = pixels[0];auto max = pixels[kernelSize * kernelSize - 1];auto med = pixels[kernelSize * kernelSize / 2];auto zxy = im.at<uchar>(row, col);if (med > min && med < max){// to Bif (zxy > min && zxy < max)return zxy;elsereturn med;}else{kernelSize += 2;if (kernelSize <= maxSize)return adaptiveProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，继续A过程。elsereturn med;}
}
示例：
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>using namespace cv;
using namespace std;
uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize);
// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}uchar adaptiveProcess(const Mat &im, int row, int col, int kernelSize, int maxSize)
{vector<uchar> pixels;for (int a = -kernelSize / 2; a <= kernelSize / 2; a++)for (int b = -kernelSize / 2; b <= kernelSize / 2; b++){pixels.push_back(im.at<uchar>(row + a, col + b));}sort(pixels.begin(), pixels.end());auto min = pixels[0];auto max = pixels[kernelSize * kernelSize - 1];auto med = pixels[kernelSize * kernelSize / 2];auto zxy = im.at<uchar>(row, col);if (med > min && med < max){// to Bif (zxy > min && zxy < max)return zxy;elsereturn med;}else{kernelSize += 2;if (kernelSize <= maxSize)return adaptiveProcess(im, row, col, kernelSize, maxSize); // 增大窗口尺寸，继续A过程。elsereturn med;}
}int main()
{Mat srcImage = imread("111.jpg",0);Mat src_salt, src_gaussian;src_salt = imread("111.jpg",0);src_gaussian = imread("111.jpg",0);if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);int minSize = 3; // 滤波器窗口的起始尺寸int maxSize = 7; // 滤波器窗口的最大尺寸//平滑滤波Mat aft_admedian_filter_salt, aft_admedian_filter_gaussian;// 扩展图像的边界copyMakeBorder(src_salt, aft_admedian_filter_salt, maxSize / 2, maxSize / 2, maxSize / 2, maxSize / 2, BorderTypes::BORDER_REFLECT);for (int j = maxSize / 2; j < aft_admedian_filter_salt.rows - maxSize / 2; j++){for (int i = maxSize / 2; i < aft_admedian_filter_salt.cols * aft_admedian_filter_salt.channels() - maxSize / 2; i++){aft_admedian_filter_salt.at<uchar>(j, i) = adaptiveProcess(aft_admedian_filter_salt, j, i, minSize, maxSize);}}copyMakeBorder(src_salt, aft_admedian_filter_gaussian, maxSize / 2, maxSize / 2, maxSize / 2, maxSize / 2, BorderTypes::BORDER_REFLECT);for (int j = maxSize / 2; j < aft_admedian_filter_gaussian.rows - maxSize / 2; j++){for (int i = maxSize / 2; i < aft_admedian_filter_gaussian.cols * aft_admedian_filter_gaussian.channels() - maxSize / 2; i++){aft_admedian_filter_gaussian.at<uchar>(j, i) = adaptiveProcess(aft_admedian_filter_gaussian, j, i, minSize, maxSize);}}medianBlur(src_gaussian, aft_admedian_filter_gaussian, 3);//增强图像对比度//Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);//filter2D(src_salt, src_salt, src_salt.depth(), kernel, Point(-1, -1));imshow("原图", srcImage);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);imshow("经过平滑滤波后含椒盐噪声的的结果", aft_admedian_filter_salt);imshow("经过平滑滤波后含高斯噪声的的结果", aft_admedian_filter_gaussian);waitKey();return 0;
}
3、双边滤波
双边滤波（Bilateral filter）是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的。双边滤波能够提供一种不会将边缘平滑掉的方法，但作为代价，需要更多的处理时间。与高斯滤波类似，双边滤波会依据每个像素及其领域构造一个加权平均值，加权计算包括两个部分，其中第一部分加权方式与高斯平滑中相同，第二部分也属于高斯加权，但不是基于中心像素点与其他像素点的空间距离之上的加权，而是基于其他像素与中心像素的亮度差值的加权。可以将双边滤波视为高斯平滑，对相似的像素赋予较高的权重，不相似的像素赋予较小的权重，也可用于图像分割之中。
双边滤波器之所以能够做到在平滑去噪的同时还能够很好的保存边缘（Edge Preserve），是由于其滤波器的核由两个函数生成：
- 一个函数由像素欧式距离决定滤波器模板的系数
- 另一个函数由像素的灰度差值决定滤波器的系数
其综合了高斯滤波器（Gaussian Filter）和α-截尾均值滤波器（Alpha-Trimmed mean Filter）的特点。高斯滤波器只考虑像素间的欧式距离，其使用的模板系数随着和窗口中心的距离增大而减小；α-截尾均值滤波器则只考虑了像素灰度值之间的差值，去掉α%的最小值和最大值后再计算均值。
双边滤波器使用二维高斯函数生成距离模板，使用一维高斯函数生成值域模板。
距离模板系数的生成公式如下：
$$
d(i, j, k, l)=\exp \left(-\frac{(i-k)^{2}+(j-l)^{2}}{2 \sigma_{d}^{2}}\right)
$$
其中,(k,l)为模板窗口的中心坐标;(i,j)为模板窗口的其他系数的坐标；σd 为高斯函数的标准差。 使用该公式生成的滤波器模板和高斯滤波器使用的模板是没有区别的。
值域模板系数的生成公式如下：
$$
r(i, j, k, l)=\exp \left(-\frac{\|f(i, j)-f(k, l)\|^{2}}{2 \sigma_{r}^{2}}\right)
$$
其中，函数f(x,y)表示要处理的图像，f(x,y)表示图像在点(x,y)处的像素值；(k,l)为模板窗口的中心坐标；(i,j)为模板窗口的其他系数的坐标；σr为高斯函数的标准差。
将上述两个模板相乘就得到了双边滤波器的模板
$$
w(i, j, k, l)=d(i, j, k, l) * r(i, j, k, l)=\exp \left(-\frac{(i-k)^{2}+(j-l)^{2}}{2 \sigma_{d}^{2}}-\frac{\|f(i, j)-f(k, l)\|^{2}}{2 \sigma_{r}^{2}}\right)
$$
其中定义域滤波和值域滤波可以用下图形象表示：
代码实现如下：
void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma) {int channels = src.channels();CV_Assert(channels == 1 || channels == 3);double space_coeff = -0.5 / (space_sigma * space_sigma);double color_coeff = -0.5 / (color_sigma * color_sigma);int radius = ksize / 2;Mat temp;copyMakeBorder(src, temp, radius, radius, radius, radius, BorderTypes::BORDER_REFLECT);vector<double> _color_weight(channels * 256); // 存放差值的平方vector<double> _space_weight(ksize * ksize); // 空间模板系数vector<int> _space_ofs(ksize * ksize); // 模板窗口的坐标double *color_weight = &_color_weight[0];double *space_weight = &_space_weight[0];int *space_ofs = &_space_ofs[0];for (int i = 0; i < channels * 256; i++)color_weight[i] = exp(i * i * color_coeff);// 生成空间模板int maxk = 0;for (int i = -radius; i <= radius; i++){for (int j = -radius; j <= radius; j++){double r = sqrt(i*i + j * j);if (r > radius)continue;space_weight[maxk] = exp(r * r * space_coeff); // 存放模板系数space_ofs[maxk++] = i * temp.step + j * channels; // 存放模板的位置，和模板系数相对应 }}// 滤波过程for (int i = 0; i < src.rows; i++){const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;uchar *dptr = dst.data + i * dst.step;if (channels == 1){for (int j = 0; j < src.cols; j++){double sum = 0, wsum = 0;int val0 = sptr[j]; // 模板中心位置的像素for (int k = 0; k < maxk; k++){int val = sptr[j + space_ofs[k]];double w = space_weight[k] * color_weight[abs(val - val0)]; // 模板系数 = 空间系数 * 灰度值系数sum += val * w;wsum += w;}dptr[j] = (uchar)cvRound(sum / wsum);}}else if (channels == 3){for (int j = 0; j < src.cols * 3; j+=3){double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;int b0 = sptr[j];int g0 = sptr[j + 1];int r0 = sptr[j + 2];for (int k = 0; k < maxk; k++){const uchar *sptr_k = sptr + j + space_ofs[k];int b = sptr_k[0];int g = sptr_k[1];int r = sptr_k[2];double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];sum_b += b * w;sum_g += g * w;sum_r += r * w;wsum += w;}wsum = 1.0f / wsum;b0 = cvRound(sum_b * wsum);g0 = cvRound(sum_g * wsum);r0 = cvRound(sum_r * wsum);dptr[j] = (uchar)b0;dptr[j + 1] = (uchar)g0;dptr[j + 2] = (uchar)r0;}}} }
OpenCV提供了对应双边滤波的API接口
void bilateralFilter(InputArray src,OutputArray dst,int d,double sigmaColor,double sigmaSpace,int borderType=BORDER_DEFAULT)第一个参数，InputArray类型的src，输入图像，即源图像，需要为8为或者浮点型单通道、三通道的图像；
第二个参数，OutputArray类型的dst，即目标图像，需要和源图像有一样的尺寸和类型；
第三个参数，int类型的d，表示在过滤过程中每个像素邻域的直径。如果这个参数被设置为负值，那么OpenCV会从第五个参数sigmaSpace来计算出它；
第四个参数，double类型的sigmaColor，颜色空间滤波器的sigma值，这个参数的值越大，就表明该像素邻域内有越宽广的颜色会被混合到一起，产生较大的半相等颜色区域；
第五个参数，double类型的 sigmaSpace，坐标空间中的sigma值，坐标空间的标注方差。它的数值越大，意味着越远的像素会相互影响，从而使更大区域中足够相似的颜色获取相同的颜色，当d>0时，d制定了邻域大小且与sigmaSpace无关否则，d正比于sigmaSpace；
第六个参数，int类型的borderType，用于推断图像外部像素的某种边界模式，有默认值BORDER_DEFAULT。
示例：
//实现直方图的反向投影
#include "stdafx.h"
#include <math.h>
#include <random>
#include <iostream>
#include <opencv2\core\core.hpp>
#include <opencv2\highgui\highgui.hpp>
#include <opencv2\imgproc\imgproc.hpp>using namespace cv;
using namespace std;
// 添加椒盐噪声
void addSaltNoise(Mat &m, int num)
{// 随机数产生器std::random_device rd; //种子std::mt19937 gen(rd()); // 随机数引擎auto cols = m.cols * m.channels();for (int i = 0; i < num; i++){auto row = static_cast<int>(gen() % m.rows);auto col = static_cast<int>(gen() % cols);auto p = m.ptr<uchar>(row);p[col++] = 255;p[col++] = 255;p[col] = 255;}
}
void addGaussianNoise(Mat &m, int mu, int sigma)
{// 产生高斯分布随机数发生器std::random_device rd;std::mt19937 gen(rd());std::normal_distribution<> d(mu, sigma);auto rows = m.rows; // 行数auto cols = m.cols * m.channels(); // 列数for (int i = 0; i < rows; i++){auto p = m.ptr<uchar>(i); // 取得行首指针for (int j = 0; j < cols; j++){auto tmp = p[j] + d(gen);tmp = tmp > 255 ? 255 : tmp;tmp = tmp < 0 ? 0 : tmp;p[j] = tmp;}}
}
void myBilateralFilter(const Mat &src, Mat &dst, int ksize, double space_sigma, double color_sigma)
{int channels = src.channels();CV_Assert(channels == 1 || channels == 3);double space_coeff = -0.5 / (space_sigma * space_sigma);double color_coeff = -0.5 / (color_sigma * color_sigma);int radius = ksize / 2;Mat temp;copyMakeBorder(src, temp, radius, radius, radius, radius, BorderTypes::BORDER_REFLECT);vector<double> _color_weight(channels * 256); // 存放差值的平方vector<double> _space_weight(ksize * ksize); // 空间模板系数vector<int> _space_ofs(ksize * ksize); // 模板窗口的坐标double *color_weight = &_color_weight[0];double *space_weight = &_space_weight[0];int *space_ofs = &_space_ofs[0];for (int i = 0; i < channels * 256; i++)color_weight[i] = exp(i * i * color_coeff);// 生成空间模板int maxk = 0;for (int i = -radius; i <= radius; i++){for (int j = -radius; j <= radius; j++){double r = sqrt(i*i + j * j);if (r > radius)continue;space_weight[maxk] = exp(r * r * space_coeff); // 存放模板系数space_ofs[maxk++] = i * temp.step + j * channels; // 存放模板的位置，和模板系数相对应}}// 滤波过程for (int i = 0; i < src.rows; i++){const uchar *sptr = temp.data + (i + radius) * temp.step + radius * channels;uchar *dptr = dst.data + i * dst.step;if (channels == 1){for (int j = 0; j < src.cols; j++){double sum = 0, wsum = 0;int val0 = sptr[j]; // 模板中心位置的像素for (int k = 0; k < maxk; k++){int val = sptr[j + space_ofs[k]];double w = space_weight[k] * color_weight[abs(val - val0)]; // 模板系数 = 空间系数 * 灰度值系数sum += val * w;wsum += w;}dptr[j] = (uchar)cvRound(sum / wsum);}}else if (channels == 3){for (int j = 0; j < src.cols * 3; j += 3){double sum_b = 0, sum_g = 0, sum_r = 0, wsum = 0;int b0 = sptr[j];int g0 = sptr[j + 1];int r0 = sptr[j + 2];for (int k = 0; k < maxk; k++){const uchar *sptr_k = sptr + j + space_ofs[k];int b = sptr_k[0];int g = sptr_k[1];int r = sptr_k[2];double w = space_weight[k] * color_weight[abs(b - b0) + abs(g - g0) + abs(r - r0)];sum_b += b * w;sum_g += g * w;sum_r += r * w;wsum += w;}wsum = 1.0f / wsum;b0 = cvRound(sum_b * wsum);g0 = cvRound(sum_g * wsum);r0 = cvRound(sum_r * wsum);dptr[j] = (uchar)b0;dptr[j + 1] = (uchar)g0;dptr[j + 2] = (uchar)r0;}}}
}
int main()
{Mat srcImage = imread("111.jpg");Mat src_salt, src_gaussian;src_salt = imread("111.jpg");src_gaussian = imread("111.jpg");if (!srcImage.data){cout << "图像打开失败！" << endl;return -1;}addSaltNoise(src_salt, 100);addGaussianNoise(src_gaussian,100,1);//平滑滤波Mat aft_bilateral_filter_salt, aft_bilateral_filter_gaussian, srcImage1;bilateralFilter(srcImage, srcImage1, 5, 100, 3);bilateralFilter(src_salt, aft_bilateral_filter_salt, 5, 100, 3);bilateralFilter(src_gaussian, aft_bilateral_filter_gaussian, 5, 100, 3);//增强图像对比度Mat kernel = (Mat_<int>(3, 3) << 0, -1, 0, -1, 5, -1, 0, -1, 0);filter2D(srcImage1, srcImage1, srcImage1.depth(), kernel, Point(-1, -1));imshow("原图", srcImage);imshow("双边滤波处理过后的原图", srcImage1);imshow("添加椒盐噪声后", src_salt);imshow("添加高斯噪声后", src_gaussian);imshow("经过平滑滤波后含椒盐噪声的的结果", aft_bilateral_filter_salt);imshow("经过平滑滤波后含高斯噪声的的结果", aft_bilateral_filter_gaussian);waitKey();return 0;
}
程序运行结果如下，双边滤波很好的对原图像进行了优化，但是对应椒盐噪声之类的高频噪声没有滤波效果。
参考资料：
图像处理基础(1)：噪声的添加和过滤
图像处理基础(2)：自适应中值滤波器(基于OpenCV实现)
图像处理基础(3)：均值滤波器及其变种
图像处理基础(4)：高斯滤波器详解
图像处理基础(5)：双边滤波器
图像处理基础(6)：锐化空间滤波器
OpenCV--Python 图像平滑之高斯平滑、均值平滑
【OpenCV图像处理】十五、图像空域滤波（上）
【OpenCV图像处理】十六、图像空域滤波（下）
【OpenCV】邻域滤波：方框、高斯、中值、双边滤波
【OpenCV学习笔记】之图像平滑（线性/非线性滤波器）
[Python图像处理] 四.图像平滑之均值滤波、方框滤波、高斯滤波及中值滤波
C++图像处理 -- 表面模糊
选择性模糊及其算法的实现。
