matlab 均值滤波_数字图像处理基础 — 高斯滤波

高斯滤波，本文主要讲其如何通过 C语言实现。不太擅长写理论性质的文章，这里仅仅阐述自己怎么实现以及简单的优化过程。

通常我们对获取的图像进行进一步处理时，往往需要先进行一次降噪，而通常我们选择的是高斯滤波也称高斯模糊，其因为较好的降噪效果而被广泛使用，具体的原因大家可以做实验以及看相关的论文去理解。

如果你经常使用MATLAB或者Python，可能只需要调用一下库函数即可实现。以matlab为例，其高斯滤波的常见使用，如下：

% Generate Kernel[5, 5], sigma = 1.0.
gaussKernel = fspecial('gaussian', [5, 5], 1.0);
% Compute gaussian filter result.
gaussResult = imfilter(Image, gaussKernel, 'replicate');

真的是方便至极，但自己不实现一下，好像没什么意思。并且嵌入式设备没有相应的库函数让你使用，所有还是有必要对其进行理解和学习。以下为正文分割线。

理论公式与模板生成

高斯滤波作为一种线性滤波器，其本质是带权值的加权均值滤波，相对于均值滤波，高斯滤波除了中心位置为1，其他位置根据距离的远近，越远滤波系数权重越低。

对高斯滤波的公式进行理解学习，参考如下公式：

其中，

为模板中心坐标，即坐标系原点。有了该公式，我们就可以实现高斯模板的生成了，参考如下

C代码：

void kernelGenerate(float *gaussKernel, int k_size, float k_sigma) {
	
	float flag = 1 / (2 * PI * k_sigma * k_sigma);
	float coef;
	float sum  = 0;
	int center = (k_size - 1) / 2;

	for (int i = 0; i < k_size; i++) {
		float x2 = (float)pow((i - center), 2);
		for (int j = 0; j < k_size; j++) {
			float y2 = (float)pow((j - center), 2);
			coef = (-x2 - y2) * flag;
			gaussKernel[i * k_size + j] = flag * exp(coef);
			sum += gaussKernel[i * k_size + j];
		}
	}

	/* Normalize */
	if (sum) {
		for (int i = 0; i < k_size * k_size; i++) {
			gaussKernel[i] /= sum;
		}
	}
}

这个模板是生成的小数模板，要是生成整数模板呢？如下，将左上角置1。然后再进行取整操作。

void kernelGenerate(float *gaussKernel, int k_size, float k_sigma) {
	
	float flag = 1 / (2 * PI * k_sigma * k_sigma);
	float coef;
	int center = (k_size - 1) / 2;

	for (int i = 0; i < k_size; i++) {
		float x2 = (float)pow((i - center), 2);
		for (int j = 0; j < k_size; j++) {
			float y2 = (float)pow((j - center), 2);
			coef = (-x2 - y2) * flag;
			gaussKernel[i * k_size + j] = flag * exp(coef);
		}
	}

	/* Intergel Normalize */
	if (gaussKernel[0]) {
                float normal = 1.0 / gaussKernel[0];
		for (int i = 0; i < k_size * k_size; i++) {
			gaussKernel[i] *= normal;
		}
	}
}

有了模板之后，就需要进行滤波了。通过公式，可以发现滤波最关键的就是滤波的核的标准差，因此需要根据实际的需要对其进行相应的设置。

滤波处理

如果直接使用公式，进行滤波则可以写出如下代码：

void gaussian_filter(unsigned char *src, unsigned char *dst, int src_w, int src_h, int k_size, float k_sigma, int padd_type, int out_type) {

	int img_w, img_h;

	if (k_size == 1)
		dst = src;

	if (k_size % 2 == 0) {
		printf("The kernel size not odd and ksize > 1! Convert it to k_size + 1");
		k_size += 1;
	}
	int padd_size = (k_size - 1) / 2;

	img_w = src_w + padd_size * 2;
	img_h = src_h + padd_size * 2;

	unsigned char *pBlur = (unsigned char *)malloc(img_w * img_h);
	unsigned char *pDst  = dst;
	float *gaussKernel   = (float *)malloc(sizeof(float) * k_size * k_size);

	for (int i = 0; i < src_h; i++) {
		unsigned char *pTmp = pBlur + (padd_size + i)* img_w + padd_size;
		memcpy(pTmp, src + i * src_w, src_w);

		/* col padding */
		for (int j = 0; j < padd_size; j++) {
			if (padd_type == 0) {	// zero padding.
				pTmp[j - padd_size] = 0;
				pTmp[src_w + j] = 0;
			} else {		// reflect padding.
				pTmp[j - padd_size] = pTmp[padd_size - j];
				pTmp[src_w + padd_size - j] = pTmp[src_w - padd_size + j];
			}
		}
	}

	/* rows padding. */
	if (padd_type == 0) {
		memset(pBlur, 0, padd_size * img_w);
		memset(pBlur + img_w * (img_h - 1 - padd_size), 0, padd_size * img_w);
	} else {
		for (int i = 0; i < padd_size; i++) {
			memcpy(pBlur + i * img_w, pBlur + (k_size - 1 - i) * img_w, img_w);
			memcpy(pBlur + (img_h - 1 - i) * img_w, pBlur + (img_h - k_size + i) * img_w, img_w);
		}
	}

	/* Generate filter kernel. */
	/* Eq: 1 / (2 * pi * Sigma^2) * exp( (x - m)^2 + (y - m)^2) */
	kernelGenerate(gaussKernel, k_size, k_sigma);

	/* Gaussian filter. */
	for (int row = padd_size; row < img_h - padd_size; row++) {
		for (int col = padd_size; col < img_w - padd_size; col++) {
			float sum = 0;
			for (int i = -padd_size; i < padd_size; i++) {
				for (int j = -padd_size; j < padd_size; j++) {
					sum += pBlur[(row + i) * img_w + col + j] * gaussKernel[(i + padd_size) * k_size + j + padd_size];
				}
			}
			pDst[(row - padd_size) * src_w + col - padd_size] = sum;
		}
	}

	free(gaussKernel);
	free(pBlur);
}