高斯滤波(GaussianFilter)原理及C++实现_高斯滤波 c++

写在前面

首先，搞清楚几个概念：滤波(高通、低通、带通、带阻) 、模糊、去噪、平滑，看下图：

                                              

滤波是对输入信号进行卷积处理的一个过程，写成一个表达式的形式是这样的：
滤波 = 卷积（ 输入信号 ，卷积模板 ）， 卷积模板/掩膜 的不同决定了不同的滤波方式，也因此产生了高通、低通、带通、带阻等基本的滤波方式。

针对低通滤波，就是保留将信号中的低频部分，抑制高频部分。要达到这个目的，可以利用均值掩膜、高斯掩膜等对输入信号进行处理。
采用均值掩膜对输入信号进行卷积的滤波方式叫均值滤波；
采用高斯掩膜对输入信号进行卷积的滤波方式叫高斯滤波；

这里高斯滤波，也叫高斯模糊。

应用：高斯滤波是一种线性平滑滤波器，对于服从正态分布的噪声有很好的抑制作用。在实际场景中，我们通常会假定图像包含的噪声为高斯白噪声，所以在许多实际应用的预处理部分，都会采用高斯滤波抑制噪声，如车牌识别等。

opencv函数声明:

void GaussianBlur(InputArray src, OutputArray dst, Size ksize, 
                  double sigmaX, double sigmaY=0,
                  int borderType=BORDER_DEFAULT )

Github: https://github.com/2209520576/Image-Processing-Algorithm/blob/master/Image%20Filtering/GaussianFilter.cpp

原理

高斯滤波和均值滤波一样，都是利用一个掩膜和图像进行卷积求解。不同之处在于：均值滤波器的模板系数都是相同的为1，而高斯滤波器的模板系数，则随着距离模板中心的增大而系数减小（服从二维高斯分布）。所以，高斯滤波器相比于均值滤波器对图像个模糊程度较小，更能够保持图像的整体细节。

二维高斯分布

高斯分布公式终于要出场了！

                                                   $f(x,y)=\frac{1}{(\sqrt{2\pi }\sigma )^{^{2}}}e^{-((x-ux)^{^{2}}+(y-uy)^{2})/2\sigma ^{^{2}}}$

其中不必纠结于系数$\frac{1}{(\sqrt{2\pi }\sigma )^{^{2}}}$,因为它只是一个常数！并不会影响互相之间的比例关系，并且最终都要进行归一化，所以在实际计算时我们是忽略它而只计算后半部分:

                                                  $f(x,y)=e^{-((x-ux)^{^{2}}+(y-uy)^{2})/2\sigma ^{^{2}}}$

其中(x,y)为掩膜内任一点的坐标，(ux,uy)为掩膜内中心点的坐标，在图像处理中可认为是整数；σ是标准差。

例如：要产生一个3×3的高斯滤波器模板，以模板的中心位置为坐标原点进行取样。模板在各个位置的坐标，如下所示（x轴水平向右，y轴竖直向下）  

这样，将各个位置的坐标带入到高斯函数中，得到的值就是模板的系数。
对于窗口模板的大小为 (2k+1)×(2k+1)，模板中各个元素值的计算公式如下：

                                                         
这样计算出来的模板有两种形式：小数和整数。

• 小数形式的模板，就是直接计算得到的值，没有经过任何的处理；
• 整数形式的，则需要进行归一化处理，将模板左上角的值归一化为1，具体介绍请看这篇博文。使用整数的模板时，需要在模板的前面加一个系数，系数为模板系数和的倒数。

生成高斯掩膜（小数形式）

知道了高斯分布原理，实现起来也就不困难了。

首先我们要确定我们生产掩模的尺寸wsize，然后设定高斯分布的标准差。生成的过程，我们首先根据模板的大小，找到模板的中心位置center。 然后就是遍历，根据高斯分布的函数，计算模板中每个系数的值。

最后模板的每个系数要除以所有系数的和。这样就得到了小数形式的模板。 

///
//x，y方向联合实现获取高斯模板
//
void generateGaussMask(cv::Mat& Mask,cv::Size wsize, double sigma){
	Mask.create(wsize,CV_64F);
	int h = wsize.height;
	int w = wsize.width;
	int center_h = (h - 1) / 2;
	int center_w = (w - 1) / 2;
	double sum = 0.0;
	double x, y;
	for (int i = 0; i < h; ++i){
		y = pow(i - center_h, 2);
		for (int j = 0; j < w; ++j){
			x = pow(j - center_w, 2);
			//因为最后都要归一化的，常数部分可以不计算，也减少了运算量
			double g = exp(-(x + y) / (2 * sigma*sigma));
			Mask.at<double>(i, j) = g;
			sum += g;
		}
	}
	Mask = Mask / sum;
}

3×3,σ=0.8的小数型模板：

 

σ的意义及选取

通过上述的实现过程，不难发现，高斯滤波器模板的生成最重要的参数就是高斯分布的标准差σ。标准差代表着数据的离散程度，如果σ较小，那么生成的模板的中心系数较大，而周围的系数较小，这样对图像的平滑效果就不是很明显；反之，σ较大，则生成的模板的各个系数相差就不是很大，比较类似均值模板，对图像的平滑效果比较明显。

来看下一维高斯分布的概率分布密度图：

于是我们有如下结论：σ越小分布越瘦高，σ越大分布越矮胖。

σ越大，分布越分散，各部分比重差别不大，于是生成的模板各元素值差别不大，类似于平均模板；
σ越小，分布越集中，中间部分所占比重远远高于其他部分，反映到高斯模板上就是中心元素值远远大于其他元素值，于是自然而然就相当于中间值得点运算。

 

基于OpenCV的c++实现

 
//按二维高斯函数实现高斯滤波
///
void GaussianFilter(cv::Mat& src, cv::Mat& dst, cv::Mat window){
	int hh = (window.rows - 1) / 2;
	int hw = (window.cols - 1) / 2;
	dst = cv::Mat::zeros(src.size(),src.type());
	//边界填充
	cv::Mat Newsrc;
	cv::copyMakeBorder(src, Newsrc, hh, hh, hw, hw, cv::BORDER_REPLICATE);//边界复制
	
	//高斯滤波
	for (int i = hh; i < src.rows + hh;++i){
		for (int j = hw; j < src.cols + hw; ++j){
			double sum[3] = { 0 };
 
			for (int r = -hh; r <= hh; ++r){
				for (int c = -hw; c <= hw; ++c){
					if (src.channels() == 1){
						sum[0] = sum[0] + Newsrc.at<uchar>(i + r, j + c) * window.at<double>(r + hh, c + hw);
					}
					else if (src.channels() == 3){
						cv::Vec3b rgb = Newsrc.at<cv::Vec3b>(i+r,j + c);
						sum[0] = sum[0] + rgb[0] * window.at<double>(r + hh, c + hw);//B
						sum[1] = sum[1] + rgb[1] * window.at<double>(r + hh, c + hw);//G
						sum[2] = sum[2] + rgb[2] * window.at<double>(r + hh, c + hw);//R
					}
				}
			}
 
			for (int k = 0; k < src.channels(); ++k){
				if (sum[k] < 0)
					sum[k] = 0;
				else if (sum[k]>255)
					sum[k] = 255;
			}
			if (src.channels() == 1)
			{
				dst.at<uchar>(i - hh, j - hw) = static_cast<uchar>(sum[0]);
			}
			else if (src.channels() == 3)
			{
				cv::Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
				dst.at<cv::Vec3b>(i-hh, j-hw) = rgb;
			}
 
		}
	}
 
}

只处理单通道或者三通道图像，模板生成后，其滤波（卷积过程）就比较简单了。不过，这样的高斯滤波过程，其循环运算次数为m×n×ksize×ksize，其中m，n为图像的尺寸；ksize为高斯滤波器的尺寸。这样其时间复杂度为O(ksize2)，随滤波器的模板的尺寸呈平方增长，当高斯滤波器的尺寸较大时，其运算效率是极低的。为了，提高滤波的运算速度，可以将二维的高斯滤波过程分解开来。

 

分离实现高斯滤波

由于高斯函数的可分离性，尺寸较大的高斯滤波器可以分成两步进行：首先将图像在水平（竖直）方向与一维高斯函数进行卷积；然后将卷积后的结果在竖直（水平）方向使用相同的一维高斯函数得到的模板进行卷积运算。具体实现代码如下：

 
//分离计算实现高斯滤波，更加高效
///
void separateGaussianFilter(cv::Mat& src, cv::Mat& dst, int wsize, double sigma){
	//获取一维高斯滤波模板
	cv::Mat window;
	window.create(1 , wsize, CV_64F);
	int center = (wsize - 1) / 2;
	double sum = 0.0;
	for (int i = 0; i < wsize; ++i){
		double g = exp(-(pow(i - center, 2)) / (2 * sigma*sigma));
		window.at<double>(0, i) = g;
		sum += g;
	}
	window = window / sum;
	//std::cout << window << std::endl;
 
	//边界填充
	int boder = (wsize - 1) / 2;
	dst = cv::Mat::zeros(src.size(), src.type());
	cv::Mat Newsrc;
	cv::copyMakeBorder(src, Newsrc, boder, boder, boder, boder, cv::BORDER_REPLICATE);//边界复制
 
	//高斯滤波--水平方向
	for (int i = boder; i < src.rows + boder; ++i){
		for (int j = boder; j < src.cols + boder; ++j){
			double sum[3] = { 0 };
 
			for (int r = -boder; r <= boder; ++r){
					if (src.channels() == 1){
						sum[0] = sum[0] + Newsrc.at<uchar>(i, j + r) * window.at<double>(0, r + boder); //行不变列变
					}
					else if (src.channels() == 3){
						cv::Vec3b rgb = Newsrc.at<cv::Vec3b>(i, j +r);
						sum[0] = sum[0] + rgb[0] * window.at<double>(0, r + boder);//B
						sum[1] = sum[1] + rgb[1] * window.at<double>(0, r + boder);//G
						sum[2] = sum[2] + rgb[2] * window.at<double>(0, r + boder);//R
				}
			}
			for (int k = 0; k < src.channels(); ++k){
				if (sum[k] < 0)
					sum[k] = 0;
				else if (sum[k]>255)
					sum[k] = 255;
			}
			if (src.channels() == 1){
				dst.at<uchar>(i - boder, j - boder) = static_cast<uchar>(sum[0]);
			}
			else if (src.channels() == 3){
				cv::Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
				dst.at<cv::Vec3b>(i - boder, j - boder) = rgb;
			}
		}
	}
 
	//高斯滤波--垂直方向
	//对水平方向处理后的dst边界填充
	cv::copyMakeBorder(dst, Newsrc, boder, boder, boder, boder, cv::BORDER_REPLICATE);//边界复制
	for (int i = boder; i < src.rows + boder; ++i){
		for (int j = boder; j < src.cols + boder; ++j){
			double sum[3] = { 0 };
 
			for (int r = -boder; r <= boder; ++r){
					if (src.channels() == 1){
						sum[0] = sum[0] + Newsrc.at<uchar>(i + r, j ) * window.at<double>(0, r + boder); //列不变行变
					}
					else if (src.channels() == 3){
						cv::Vec3b rgb = Newsrc.at<cv::Vec3b>(i + r, j);
						sum[0] = sum[0] + rgb[0] * window.at<double>(0, r + boder);//B
						sum[1] = sum[1] + rgb[1] * window.at<double>(0, r + boder);//G
						sum[2] = sum[2] + rgb[2] * window.at<double>(0, r + boder);//R
					}
				}
			for (int k = 0; k < src.channels(); ++k){
				if (sum[k] < 0)
					sum[k] = 0;
				else if (sum[k]>255)
					sum[k] = 255;
			}
			if (src.channels() == 1){
				dst.at<uchar>(i - boder, j - boder) = static_cast<uchar>(sum[0]);
			}
			else if (src.channels() == 3){
				cv::Vec3b rgb = { static_cast<uchar>(sum[0]), static_cast<uchar>(sum[1]), static_cast<uchar>(sum[2]) };
				dst.at<cv::Vec3b>(i - boder, j - boder) = rgb;
			}
	   	}	
	}
 }

首先得到一维高斯函数的模板，在卷积（滤波）的过程中，保持行不变，列变化，在水平方向上做卷积运算；接着在上述得到的结果上，保持列不边，行变化，在竖直方向上做卷积运算。 这样分解开来，算法的时间复杂度为O(ksize)，运算量和滤波器的模板尺寸呈线性增长。

 

测试：

int main(){
	cv::Mat src = cv::imread("I:\\Learning-and-Practice\\2019Change\\Image process algorithm\\Img\\4.bmp");
	if (src.empty()){
		return -1;
	}
	cv::Mat Mask;
	cv::Mat dst1;
	cv::Mat dst2;
 
	//暴力实现
	generateGaussMask(Mask, cv::Size(3, 3), 0.8);//获取二维高斯滤波模板
	std::cout << Mask << std::endl;
	GaussianFilter(src, dst1, Mask);
 
	//分离实现
	separateGaussianFilter( src, dst2, 7, 2);
	
	
	cv::namedWindow("src");
	cv::imshow("src", src);
	cv::namedWindow("暴力实现",CV_WINDOW_NORMAL);
	cv::imshow("暴力实现", dst1);
	cv::namedWindow("分离实现", CV_WINDOW_NORMAL);
	cv::imshow("分离实现", dst2);
	//cv::imwrite("I:\\Learning-and-Practice\\2019Change\\Image process algorithm\\Image Filtering\\GaussianFilter\\txt.jpg", dst2);
	cv::waitKey(0);
	return 0;
}

效果

参考：

https://www.zhihu.com/question/54918332（讲解很幽默，很通俗易懂）

https://www.cnblogs.com/wangguchangqing/p/6407717.html（讲解非常专业，很全面，有公式、图像、代码、效果展示）