写在前面

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

preview

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

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

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

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

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(讲解非常专业,很全面,有公式、图像、代码、效果展示)