高斯滤波原理及实现

https://blog.csdn.net/linqianbi/article/details/78635941

高斯滤波器是一种线性滤波器，能够有效的抑制噪声，平滑图像。其作用原理和均值滤波器类似，都是取滤波器窗口内的像素的均值作为输出。其窗口模板的系数和均值滤波器不同，均值滤波器的模板系数都是相同的为1；而高斯滤波器的模板系数，则随着距离模板中心的增大而系数减小。所以，高斯滤波器相比于均值滤波器对图像个模糊程度较小。

什么是高斯滤波器

既然名称为高斯滤波器，那么其和高斯分布（正态分布）是有一定的关系的。一个二维的高斯函数如下：

其中为点坐标，在图像处理中可认为是整数；是标准差。要想得到一个高斯滤波器的模板，可以对高斯函数进行离散化，得到的高斯函数值作为模板的系数。例如：要产生一个的高斯滤波器模板，以模板的中心位置为坐标原点进行取样。模板在各个位置的坐标，如下所示（x轴水平向右，y轴竖直向下）

这样，将各个位置的坐标带入到高斯函数中，得到的值就是模板的系数。

对于窗口模板的大小为 ，模板中各个元素值的计算公式如下：

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

小数形式的模板，就是直接计算得到的值，没有经过任何的处理；

整数形式的，则需要进行归一化处理，将模板左上角的值归一化为1，下面会具体介绍。使用整数的模板时，需要在模板的前面加一个系数，系数为也就是模板系数和的倒数。

二维高斯模板的生成：

double **getTwoGuassionArray(int size, double sigma) {
	double sum = 0.0;
	
	int kerR = size / 2;
	
	
	double **arr = new double*[size];
	for (int i = 0; i < size; i++)
	{
		arr[i] = new double[size];
	}
	
	
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			
			arr[i][j] = exp(-((i-kerR)*(i-kerR) + (j-kerR)*(j-kerR)) / (2 * sigma*sigma));
			sum += arr[i][j];
		}
	}
	
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			arr[i][j] /= sum;
			cout << arr[i][j] << endl;
		}
	}
	
	return arr;
}

           

一维高斯模板的生成：

double *getOneGuassionArray(int size, double sigma)
{
	double sum = 0.0;
	
	int kerR = size / 2;
 
	
	double *arr = new double[size];
	for (int i = 0; i < size; i++)
	{
 
		
		arr[i] = exp(-((i - kerR)*(i - kerR)) / (2 * sigma*sigma));
		sum += arr[i];
 
	}
	
	for (int i = 0; i < size; i++)
	{
		arr[i] /= sum;
		cout << arr[i] << endl;
	}
	return arr;
}
           

值的意义及选取

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

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

横轴表示可能得取值x，竖轴表示概率分布密度F(x)，那么不难理解这样一个曲线与x轴围成的图形面积为1。（标准差）决定了这个图形的宽度，可以得出这样的结论：越大，则图形越宽，尖峰越小，图形较为平缓；越小，则图形越窄，越集中，中间部分也就越尖，图形变化比较剧烈。这其实很好理解，如果sigma也就是标准差越大，则表示该密度分布一定比较分散，由于面积为1，于是尖峰部分减小，宽度越宽（分布越分散）；同理，当越小时，说明密度分布较为集中，于是尖峰越尖，宽度越窄！

于是可以得到如下结论：

越大，分布越分散，各部分比重差别不大，于是生成的模板各元素值差别不大，类似于平均模板；

越小，分布越集中，中间部分所占比重远远高于其他部分，反映到高斯模板上就是中心元素值远远大于其他元素值，于是自然而然就相当于中间值得点运算。

#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include <iostream>  
#include <cmath>
using namespace cv;
using namespace std;
 
 
double **getTwoGuassionArray(int size, double sigma) {
	double sum = 0.0;
	
	int kerR = size / 2;
	
	
	double **arr = new double*[size];
	for (int i = 0; i < size; i++)
	{
		arr[i] = new double[size];
	}
	
	
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			
			arr[i][j] = exp(-((i-kerR)*(i-kerR) + (j-kerR)*(j-kerR)) / (2 * sigma*sigma));
			sum += arr[i][j];
		}
	}
	
	for (int i = 0; i < size; i++)
	{
		for (int j = 0; j < size; j++)
		{
			arr[i][j] /= sum;
			cout << arr[i][j] << endl;
		}
	}
	
	return arr;
}
 
 
void MyGaussianBlur(Mat &srcImage, Mat &dst,int size)
{
	CV_Assert(srcImage.channels() || srcImage.channels() == 3); 
	int kerR = size / 2;
	dst = srcImage.clone();
	int channels = dst.channels();
	double** arr;
	arr = getTwoGuassionArray(size, 1);
	
	
	for (int i = kerR; i < dst.rows - kerR; i++)
	{
		for (int j = kerR; j < dst.cols - kerR; j++)
		{
			double GuassionSum[3] = { 0 };
			
			for (int k = -kerR; k <= kerR; k++)
			{
				for (int n = -kerR; n <= kerR; n++)
				{
					if (channels == 1)
					{
						GuassionSum[0] += arr[kerR + k][kerR + n] * dst.at<uchar>(i + k, j + n);
					}
					else if(channels == 3)
					{
						Vec3b bgr = dst.at<Vec3b>(i + k, j + n);
						auto a = arr[kerR + k][kerR + n];
						GuassionSum[0] += a*bgr[0];
						GuassionSum[1] += a*bgr[1];
						GuassionSum[2] += a*bgr[2];
					}
					
				}
			}
			for (int k = 0; k < channels; k++)
			{
				if (GuassionSum[k] < 0)
					GuassionSum[k] = 0;
				else if (GuassionSum[k] > 255)
					GuassionSum[k] = 255;
			}
			if (channels == 1)
				dst.at<uchar>(i, j) = static_cast<uchar>(GuassionSum[0]);
			else if (channels == 3)
			{
				Vec3b bgr = { static_cast<uchar>(GuassionSum[0]), static_cast<uchar>(GuassionSum[1]), static_cast<uchar>(GuassionSum[2]) };
				dst.at<Vec3b>(i, j) = bgr;
			}		
			
		}
	}
	for (int i = 0; i < size; i++)
		delete[] arr[i];
	delete[] arr;
}
 
 
int main()
{
	Mat srcImage = imread("flower3.jpg");
	if (!srcImage.data)
	{
		printf("could not load image...\n");
		return -1;
	}
	//getGuassionArray(3, 1);
	
	Mat resMat;
	MyGaussianBlur(srcImage, resMat, 3);
	
	imshow("src", srcImage);
	imshow("resMat", resMat);
	waitKey(0);
	
	return 0;
}
           

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

分离实现高斯滤波

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

#include "opencv2/imgproc/imgproc.hpp"  
#include "opencv2/highgui/highgui.hpp"  
#include <iostream>  
#include <cmath>
using namespace cv;
using namespace std;
 
 
double *getOneGuassionArray(int size, double sigma)
{
	double sum = 0.0;
	
	int kerR = size / 2;
 
	double *arr = new double[size];
	for (int i = 0; i < size; i++)
	{
 
		
		arr[i] = exp(-((i - kerR)*(i - kerR)) / (2 * sigma*sigma));
		sum += arr[i];
 
	}
	
	for (int i = 0; i < size; i++)
	{
		arr[i] /= sum;
		cout << arr[i] << endl;
	}
	return arr;
}
 
 
 
 
void MyGaussianBlur(Mat &srcImage, Mat &dst, int size)
{
	CV_Assert(srcImage.channels() || srcImage.channels() == 3); 
	int kerR = size / 2;
	dst = srcImage.clone();
	int channels = dst.channels();
	double* arr;
	arr = getOneGuassionArray(size, 1);
	
	for (int i = kerR; i < dst.rows - kerR; i++)
	{
		for (int j = kerR; j < dst.cols - kerR; j++)
		{
			double GuassionSum[3] = { 0 };
			
			for (int k = -kerR; k <= kerR; k++)
			{
				
					if (channels == 1)
					{
						GuassionSum[0] += arr[kerR + k] * dst.at<uchar>(i , j + k);
					}
					else if (channels == 3)
					{
						Vec3b bgr = dst.at<Vec3b>(i, j + k);
						auto a = arr[kerR + k];
						GuassionSum[0] += a*bgr[0];
						GuassionSum[1] += a*bgr[1];
						GuassionSum[2] += a*bgr[2];
					}
			}
			for (int k = 0; k < channels; k++)
			{
				if (GuassionSum[k] < 0)
					GuassionSum[k] = 0;
				else if (GuassionSum[k] > 255)
					GuassionSum[k] = 255;
			}
			if (channels == 1)
				dst.at<uchar>(i, j) = static_cast<uchar>(GuassionSum[0]);
			else if (channels == 3)
			{
				Vec3b bgr = { static_cast<uchar>(GuassionSum[0]), static_cast<uchar>(GuassionSum[1]), static_cast<uchar>(GuassionSum[2]) };
				dst.at<Vec3b>(i, j) = bgr;
			}
 
		}
	}
 
	
	for (int i = kerR; i < dst.rows - kerR; i++)
	{
		for (int j = kerR; j < dst.cols - kerR; j++)
		{
			double GuassionSum[3] = { 0 };
			
			for (int k = -kerR; k <= kerR; k++)
			{
 
				if (channels == 1)
				{
					GuassionSum[0] += arr[kerR + k] * dst.at<uchar>(i + k, j);
				}
				else if (channels == 3)
				{
					Vec3b bgr = dst.at<Vec3b>(i + k, j);
					auto a = arr[kerR + k];
					GuassionSum[0] += a*bgr[0];
					GuassionSum[1] += a*bgr[1];
					GuassionSum[2] += a*bgr[2];
				}
			}
			for (int k = 0; k < channels; k++)
			{
				if (GuassionSum[k] < 0)
					GuassionSum[k] = 0;
				else if (GuassionSum[k] > 255)
					GuassionSum[k] = 255;
			}
			if (channels == 1)
				dst.at<uchar>(i, j) = static_cast<uchar>(GuassionSum[0]);
			else if (channels == 3)
			{
				Vec3b bgr = { static_cast<uchar>(GuassionSum[0]), static_cast<uchar>(GuassionSum[1]), static_cast<uchar>(GuassionSum[2]) };
				dst.at<Vec3b>(i, j) = bgr;
			}
 
		}
	}
	delete[] arr;
}
 
 
int main()
{
	Mat srcImage = imread("flower3.jpg");
	if (!srcImage.data)
	{
		printf("could not load image...\n");
		return -1;
	}
	//getGuassionArray(3, 1);
 
	Mat resMat;
	MyGaussianBlur(srcImage, resMat, 3);
 
	imshow("src", srcImage);
	imshow("resMat", resMat);
	waitKey(0);
 
	return 0;
}
           

代码没有重构较长，不过其实现原理是比较简单的。首先得到一维高斯函数的模板，在卷积（滤波）的过程中，保持行不变，列变化，在水平方向上做卷积运算；接着在上述得到的结果上，保持列不边，行变化，在竖直方向上做卷积运算。 这样分解开来，算法的时间复杂度为，运算量和滤波器的模板尺寸呈线性增长。