传统高斯模糊与优化算法(附完整C++代码)

高斯模糊（英语：Gaussian Blur），也叫高斯平滑，是在Adobe Photoshop、GIMP以及Paint.NET等图像处理软件中广泛使用的处理效果，通常用它来减少图像噪声以及降低细节层次。这种模糊技术生成的图像，其视觉效果就像是经过一个半透明屏幕在观察图像，这与镜头焦外成像效果散景以及普通照明阴影中的效果都明显不同。高斯平滑也用于计算机视觉算法中的预先处理阶段，以增强图像在不同比例大小下的图像效果（参见尺度空间表示以及尺度空间实现）。 从数学的角度来看，图像的高斯模糊过程就是图像与正态分布做卷积。由于正态分布又叫作高斯分布，所以这项技术就叫作高斯模糊。图像与圆形方框模糊做卷积将会生成更加精确的焦外成像效果。由于高斯函数的傅立叶变换是另外一个高斯函数，所以高斯模糊对于图像来说就是一个低通滤波器。

高斯模糊是一种图像模糊滤波器，它用正态分布计算图像中每个像素的变换。N维空间正态分布方程为

在二维空间定义为

其中r是模糊半径 （），σ是正态分布的标准偏差。在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。这样进行模糊处理比其它的均衡模糊滤波器更高地保留了边缘效果，参见尺度空间实现。

理论上来讲，图像中每点的分布都不为零，这也就是说每个像素的计算都需要包含整幅图像。在实际应用中，在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。通常，图像处理程序只需要计算的矩阵就可以保证相关像素影响。对于边界上的点，通常采用复制周围的点到另一面再进行加权平均运算。

除了圆形对称之外，高斯模糊也可以在二维图像上对两个独立的一维空间分别进行计算，这叫作线性可分。这也就是说，使用二维矩阵变换得到的效果也可以通过在水平方向进行一维高斯矩阵变换加上竖直方向的一维高斯矩阵变换得到。从计算的角度来看，这是一项有用的特性，因为这样只需要次计算，而不可分的矩阵则需要次计算，其中,是需要进行滤波的图像的维数，、是滤波器的维数。

对一幅图像进行多次连续高斯模糊的效果与一次更大的高斯模糊可以产生同样的效果，大的高斯模糊的半径是所用多个高斯模糊半径平方和的平方根。例如，使用半径分别为6和8的两次高斯模糊变换得到的效果等同于一次半径为10的高斯模糊效果，。根据这个关系，使用多个连续较小的高斯模糊处理不会比单个高斯较大处理时间要少。

在减小图像尺寸的场合经常使用高斯模糊。在进行欠采样的时候，通常在采样之前对图像进行低通滤波处理。这样就可以保证在采样图像中不会出现虚假的高频信息。高斯模糊有很好的特性，如没有明显的边界，这样就不会在滤波图像中形成震荡。

以上资料摘自维基百科(高斯模糊词条)：

https://zh.wikipedia.org/wiki/%E9%AB%98%E6%96%AF%E6%A8%A1%E7%B3%8A

那么具体如何实现呢？

代码献上:

inline int* buildGaussKern(int winSize, int sigma)
{
	int wincenter, x;
	float   sum = 0.0f;
	wincenter = winSize / 2;
	float *kern = (float*)malloc(winSize*sizeof(float));
	int *ikern = (int*)malloc(winSize*sizeof(int));
	float SQRT_2PI = 2.506628274631f;
	float sigmaMul2PI = 1.0f / (sigma * SQRT_2PI);
	float divSigmaPow2 = 1.0f / (2.0f * sigma * sigma);
	for (x = 0; x < wincenter + 1; x++)
	{
		kern[wincenter - x] = kern[wincenter + x] = exp(-(x * x)* divSigmaPow2) * sigmaMul2PI;
		sum += kern[wincenter - x] + ((x != 0) ? kern[wincenter + x] : 0.0);
	}
	sum = 1.0f / sum;
	for (x = 0; x < winSize; x++)
	{
		kern[x] *= sum;
		ikern[x] = kern[x] * 256.0f;
	}
	free(kern);
	return ikern;
}

void GaussBlur(unsigned char*  pixels, unsigned int    width, unsigned int  height, unsigned  int channels, int sigma)
{
	width = 3 * width;
	if ((width % 4) != 0) width += (4 - (width % 4));

	unsigned int  winsize = (1 + (((int)ceil(3 * sigma)) * 2));
	int *gaussKern = buildGaussKern(winsize, sigma);
	winsize *= 3;
	unsigned int  halfsize = winsize / 2;

	unsigned char *tmpBuffer = (unsigned char*)malloc(width * height* sizeof(unsigned char));

	for (unsigned int h = 0; h < height; h++)
	{
		unsigned int  rowWidth = h * width;

		for (unsigned int w = 0; w < width; w += channels)
		{
			unsigned int rowR = 0;
			unsigned int rowG = 0;
			unsigned int rowB = 0;
			int * gaussKernPtr = gaussKern;
			int whalfsize = w + width - halfsize;
			unsigned int  curPos = rowWidth + w;
			for (unsigned int k = 1; k < winsize; k += channels)
			{
				unsigned int  pos = rowWidth + ((k + whalfsize) % width);
				int fkern = *gaussKernPtr++;
				rowR += (pixels[pos] * fkern);
				rowG += (pixels[pos + 1] * fkern);
				rowB += (pixels[pos + 2] * fkern);
			}

			tmpBuffer[curPos] = ((unsigned char)(rowR >> 8));
			tmpBuffer[curPos + 1] = ((unsigned char)(rowG >> 8));
			tmpBuffer[curPos + 2] = ((unsigned char)(rowB >> 8));

		}
	}
	winsize /= 3;
	halfsize = winsize / 2;
	for (unsigned int w = 0; w < width; w++)
	{
		for (unsigned int h = 0; h < height; h++)
		{
			unsigned    int col_all = 0;
			int hhalfsize = h + height - halfsize;
			for (unsigned int k = 0; k < winsize; k++)
			{
				col_all += tmpBuffer[((k + hhalfsize) % height)* width + w] * gaussKern[k];
			}
			pixels[h * width + w] = (unsigned char)(col_all >> 8);
		}
	}
	free(tmpBuffer);
	free(gaussKern);
}

备注:

之于原始算法，我做了一些小改动，主要是为了考虑一点点性能上的问题。

有时会写太多注释反而显得啰嗦，所以将就着看哈。

这份代码，实测速度非常糟糕，处理一张5000x3000在半径大小5左右都要耗时十来秒至几十秒不等，实在难以接受。

由于速度的问题，网上就有不少优化算法的实现。

之前我也发过一篇《快速高斯模糊算法》，在同等条件下，这个算法已经比如上算法快上十几倍。

由于这份代码实在难以阅读学习，所以，我对其进行了进一步的调整和优化。

void GaussianBlur(unsigned char* img,  unsigned int width, unsigned int height, unsigned int channels, unsigned int radius)
{
	radius = min(max(1, radius), 248);
	unsigned int kernelSize = 1 + radius * 2;
	unsigned int* kernel = (unsigned int*)malloc(kernelSize* sizeof(unsigned int));
	memset(kernel, 0, kernelSize* sizeof(unsigned int));
	int(*mult)[256] = (int(*)[256])malloc(kernelSize * 256 * sizeof(int));
	memset(mult, 0, kernelSize * 256 * sizeof(int));

	int	xStart = 0;
	int	yStart = 0;
	width = xStart + width - max(0, (xStart + width) - width);
	height = yStart + height - max(0, (yStart + height) - height);
	int imageSize = width*height;
	int widthstep = width*channels;
	if (channels == 3 || channels == 4)
	{
		unsigned char *    CacheImg = nullptr;
		CacheImg = (unsigned char *)malloc(sizeof(unsigned char) * imageSize * 6);
		if (CacheImg == nullptr) return;
		unsigned char *    rCache = CacheImg;
		unsigned char *    gCache = CacheImg + imageSize;
		unsigned char *    bCache = CacheImg + imageSize * 2;
		unsigned char *    r2Cache = CacheImg + imageSize * 3;
		unsigned char *    g2Cache = CacheImg + imageSize * 4;
		unsigned char *    b2Cache = CacheImg + imageSize * 5;
		int sum = 0;
		for (int K = 1; K < radius; K++){
			unsigned int szi = radius - K;
			kernel[radius + K] = kernel[szi] = szi*szi;
			sum += kernel[szi] + kernel[szi];
			for (int j = 0; j < 256; j++){
				mult[radius + K][j] = mult[szi][j] = kernel[szi] * j;
			}
		}
		kernel[radius] = radius*radius;
		sum += kernel[radius];
		for (int j = 0; j < 256; j++){
			mult[radius][j] = kernel[radius] * j;
		}
		for (int Y = 0; Y < height; ++Y) {
			unsigned char*     LinePS = img + Y*widthstep;
			unsigned char*     LinePR = rCache + Y*width;
			unsigned char*     LinePG = gCache + Y*width;
			unsigned char*     LinePB = bCache + Y*width;
			for (int X = 0; X < width; ++X) {
				int     p2 = X*channels;
				LinePR[X] = LinePS[p2];
				LinePG[X] = LinePS[p2 + 1];
				LinePB[X] = LinePS[p2 + 2];
			}
		}
		int kernelsum = 0;
		for (int K = 0; K < kernelSize; K++){
			kernelsum += kernel[K];
		}
		float fkernelsum = 1.0f / kernelsum;
		for (int Y = yStart; Y < height; Y++){
			int heightStep = Y * width;
			unsigned char*     LinePR = rCache + heightStep;
			unsigned char*     LinePG = gCache + heightStep;
			unsigned char*     LinePB = bCache + heightStep;
			for (int X = xStart; X < width; X++){
				int cb = 0;
				int cg = 0;
				int cr = 0;
				for (int K = 0; K < kernelSize; K++){
					unsigned    int     readPos = ((X - radius + K + width) % width);
					int * pmult = mult[K];
					cr += pmult[LinePR[readPos]];
					cg += pmult[LinePG[readPos]];
					cb += pmult[LinePB[readPos]];
				}
				unsigned int p = heightStep + X;
				r2Cache[p] = cr* fkernelsum;
				g2Cache[p] = cg* fkernelsum;
				b2Cache[p] = cb* fkernelsum;
			}
		}
		for (int X = xStart; X < width; X++){
			int WidthComp = X*channels;
			int WidthStep = width*channels;
			unsigned char*     LinePS = img + X*channels;
			unsigned char*     LinePR = r2Cache + X;
			unsigned char*     LinePG = g2Cache + X;
			unsigned char*     LinePB = b2Cache + X;
			for (int Y = yStart; Y < height; Y++){
				int cb = 0;
				int cg = 0;
				int cr = 0;
				for (int K = 0; K < kernelSize; K++){
					unsigned int   readPos = ((Y - radius + K + height) % height) * width;
					int * pmult = mult[K];
					cr += pmult[LinePR[readPos]];
					cg += pmult[LinePG[readPos]];
					cb += pmult[LinePB[readPos]];
				}
				int    p = Y*WidthStep;
				LinePS[p] = (unsigned char)(cr * fkernelsum);
				LinePS[p + 1] = (unsigned char)(cg * fkernelsum);
				LinePS[p + 2] = (unsigned char)(cb* fkernelsum);

			}
		}
		free(CacheImg);
	}
	else if (channels == 1)
	{
		unsigned char *    CacheImg = nullptr;
		CacheImg = (unsigned char *)malloc(sizeof(unsigned char) * imageSize * 2);
		if (CacheImg == nullptr) return;
		unsigned char *    rCache = CacheImg;
		unsigned char *    r2Cache = CacheImg + imageSize;

		int sum = 0;
		for (int K = 1; K < radius; K++){
			unsigned int szi = radius - K;
			kernel[radius + K] = kernel[szi] = szi*szi;
			sum += kernel[szi] + kernel[szi];
			for (int j = 0; j < 256; j++){
				mult[radius + K][j] = mult[szi][j] = kernel[szi] * j;
			}
		}
		kernel[radius] = radius*radius;
		sum += kernel[radius];
		for (int j = 0; j < 256; j++){
			mult[radius][j] = kernel[radius] * j;
		}
		for (int Y = 0; Y < height; ++Y) {
			unsigned char*     LinePS = img + Y*widthstep;
			unsigned char*     LinePR = rCache + Y*width;
			for (int X = 0; X < width; ++X) {
				LinePR[X] = LinePS[X];
			}
		}
		int kernelsum = 0;
		for (int K = 0; K < kernelSize; K++){
			kernelsum += kernel[K];
		}
		float fkernelsum = 1.0f / kernelsum;
		for (int Y = yStart; Y < height; Y++){
			int heightStep = Y * width;
			unsigned char*     LinePR = rCache + heightStep;
			for (int X = xStart; X < width; X++){
				int cb = 0;
				int cg = 0;
				int cr = 0;
				for (int K = 0; K < kernelSize; K++){
					unsigned    int     readPos = ( (X - radius + K+width)%width);
					int * pmult = mult[K];
					cr += pmult[LinePR[readPos]];
				}
				unsigned int p = heightStep + X;
				r2Cache[p] = cr * fkernelsum;
			}
		}
		for (int X = xStart; X < width; X++){
			int WidthComp = X*channels;
			int WidthStep = width*channels;
			unsigned char*     LinePS = img + X*channels;
			unsigned char*     LinePR = r2Cache + X;
			for (int Y = yStart; Y < height; Y++){
				int cb = 0;
				int cg = 0;
				int cr = 0;
				for (int K = 0; K < kernelSize; K++){
					unsigned int   readPos = ((Y - radius + K+height)%height) * width;
					int * pmult = mult[K];
					cr += pmult[LinePR[readPos]];
				}
				int    p = Y*WidthStep;
				LinePS[p] = (unsigned char)(cr* fkernelsum);
			}
		}
		free(CacheImg);
	}
	free(kernel);
	free(mult);
}

　　其中有部分算法优化技巧，想来也能起到一点抛砖引玉的作用。

贴个效果图:

本文只是抛砖引玉一下，若有其他相关问题或者需求也可以邮件联系我探讨。

邮箱地址是:

[email protected]

时间： 2024-08-02 18:22:38

传统高斯模糊与优化算法(附完整C++代码)的相关文章

mser 最大稳定极值区域(文字区域定位)算法 附完整C代码

mser 的全称:Maximally Stable Extremal Regions 第一次听说这个算法时,是来自当时部门的一个同事, 提及到他的项目用它来做文字区域的定位,对这个算法做了一些优化. 也就是中文车牌识别开源项目EasyPR的作者liuruoze,刘兄. 自那时起就有一块石头没放下,想要找个时间好好理理这个算法. 学习一些它的一些思路. 因为一般我学习算法的思路:3个做法, 第一步,编写demo示例. 第二步,进行算法移植或效果改进. 第三步,进行算法性能优化. 然后在这三个过程中

音频降噪算法 附完整C代码

降噪是音频图像算法中的必不可少的. 目的肯定是让图片或语音 更加自然平滑,简而言之,美化. 图像算法和音频算法 都有其共通点. 图像是偏向 空间 处理,例如图片中的某个区域. 图像很多时候是以二维数据为主,矩形数据分布. 音频更偏向 时间 处理,例如语音中的某短时长. 音频一般是一维数据为主,单声道波长. 处理方式也是差不多,要不单通道处理,然后合并,或者直接多通道处理. 只是处理时候数据参考系维度不一而已. 一般而言, 图像偏向于多通道处理,音频偏向于单通道处理. 而从数字信号的角度来看,也可

自动曝光修复算法 附完整C代码

众所周知, 图像方面的3A算法有: AF自动对焦(Automatic Focus)自动对焦即调节摄像头焦距自动得到清晰的图像的过程 AE自动曝光(Automatic Exposure)自动曝光的是为了使感光器件获得合适的曝光量 AW自动白平衡(Automatic White Balance)白平衡的本质是使白色物体在任何光源下都显示白色 前面的文章也有提及过,在刚开始做图像算法的时候,我是先攻克的自动白平衡算法. 后来攻克自动曝光的时候,傻啦吧唧的,踩了不少坑. 我相信一定不止我一个,一开始的时

音频自动增益 与 静音检测 算法 附完整C代码

前面分享过一个算法<音频增益响度分析 ReplayGain 附完整C代码示例> 主要用于评估一定长度音频的音量强度, 而分析之后,很多类似的需求,肯定是做音频增益,提高音量诸如此类做法. 不过在项目实测的时候,其实真的很难定标准, 到底在什么样的环境下,要增大音量,还是降低. 在通讯行业一般的做法就是采用静音检测, 一旦检测为静音或者噪音,则不做处理,反之通过一定的策略进行处理. 这里就涉及到两个算法,一个是静音检测,一个是音频增益. 增益其实没什么好说的,类似于数据归一化拉伸的做法. 静音检

音频自动增益 与 静音检测 算法 附完整C代码【转】

转自:https://www.cnblogs.com/cpuimage/p/8908551.html 前面分享过一个算法<音频增益响度分析 ReplayGain 附完整C代码示例> 主要用于评估一定长度音频的音量强度, 而分析之后,很多类似的需求,肯定是做音频增益,提高音量诸如此类做法. 不过在项目实测的时候,其实真的很难定标准, 到底在什么样的环境下,要增大音量,还是降低. 在通讯行业一般的做法就是采用静音检测, 一旦检测为静音或者噪音,则不做处理,反之通过一定的策略进行处理. 这里就涉及到

基于傅里叶变换的音频重采样算法 (附完整c代码)

前面有提到音频采样算法: WebRTC 音频采样算法 附完整C++示例代码 简洁明了的插值音频重采样算法例子 (附完整C代码) 近段时间有不少朋友给我写过邮件,说了一些他们使用的情况和问题. 坦白讲,我精力有限,但一般都会抽空回复一下. 大多数情况,阅读一下代码就能解决的问题, 也是要尝试一下的. 没准,你就解决了呢? WebRtc的采样算法本身就考虑到它的自身应用场景, 所以它会有一些局限性,例如不支持任意采样率等等. 而简洁插值的这个算法, 我个人也一直在使用,因为简洁明了,简单粗暴. 我自

图片文档倾斜矫正算法 附完整c代码

2年前在学习图像算法的时候看到一个文档倾斜矫正的算法. 也就是说能将一些文档图像进行旋转矫正, 当然这个算法一般用于一些文档扫描软件做后处理 或者用于ocr 文字识别做前处理. 相关的关键词: 抗倾斜 反倾斜  Deskew 等等. 最简单算法实现思路,采用 霍夫变换(Hough Transform)进行直线检测, 当然也可以用霍夫变换检测圆. 在倾斜矫正算法中,自然就是检测直线. 通过对检测出来的直线进行角度判断, 一般取 认可度最高的几条直线进行计算, 最后求取均衡后的角度值. 进行图像角度

【数据结构与算法】内部排序之一：插入排序和希尔排序的N中实现（不断优化，附完整源码）

转载请注明出处:http://blog.csdn.net/ns_code/article/details/20043459   前言 本来想将所有的内部排序总结为一篇博文,但是随着研究的深入,还是放弃了这个念头,斟前酌后,还是觉得分开来写比较好,具体原因,看完本篇博文也就自然明了了. 本篇文章主要探讨插入排序和希尔排序,之所将二者放在一起,很明显,是因为希尔排序是建立在插入排序的基础之上的.     注:以下各排序算法的N种实现方法大部分都是我根据算法思想,自己写出来的,或者是参考其本身的经典实

音频算法之小黄人变声 附完整C代码

前面提及到<大话音频变声原理 附简单示例代码>与<声音变调算法PitchShift(模拟汤姆猫) 附完整C++算法实现代码> 都稍微讲过变声的原理和具体实现. 大家都知道,算法从实现到最后工程应用,中间的环节和问题特别多. 尤其是编码的架构设计,好的数据结构和代码逻辑封装肯定是可复用,组件化的. 前几天写完<音频识别算法思考与阶段性小结>的时候, 我也提及到了. 会做一些算法编码优化相关的分享. 而有时候我总觉得文字表达很苍白, 所以我尽可能地把代码写得简洁易懂, 一方