高斯模糊算法

项目上对图像处理需要用到点高斯算法，网上找到一篇对原理及部分问题分析讲解的还不错的文章，分享一下，后付自己的一段代码

 

- 理论 -

高斯分布函数可表示为一个一维的函数G(x)

 或者一个二维的函数G(x,y)

在这些函数中， X和Y代表了相对于原始中心点（center tap）像素的偏移（pixel offsets）值。也就是说，他们距离中心多少像素。这里的center tap，通常翻译为“中心抽头”，它在电学中的概念是：在整个次级线圈的中心拉出的一段导线上，它相对于另外两边的抽头电压居中，而为0，两边的电压就是一正一负。在这里，我们也可以做相似的理解。即，它表示，以某个像素为中心进行取样，假设它的坐标为（x,y），那么周围四个点的坐标就是(x-1,y)(x+1,y)(x,y-1)(x,y+1)。对于整个分布来说，我们也可以用一个平均参数来实施位移，从而取代center tap的作用。但这样的模糊不符合我们的要求，我们希望我们的分布是基于中心的。

函数中的σ，即“西格玛”。在统计学中，它常用来表示“标准差”。高斯模糊的标准差，表示模糊的延伸距离。它的缺省值一般设为1，你可以提高它，来获得更强的效果。

函数中的e是欧拉数（Euler's number），它的值一般取2.7。你可以到http://blog.sina.com.cn/s/blog_54ab9583010007u5.html看看精确到小数点后10000位的e。

这个函数的结果，就是以x（在一个方向上），或者(x,y)（在两个方向上），为中心的加权（weight）值，或者说该点在多大程度上影响模糊后的像素。

那么，它能获得什么样的结果呢？当σ 值取1时，我们可以得到：

G(0) = 0.3989422804

G(1) = G(-1) = 0.2419707245

G(2) = G(-2) = 0.0539909665

G(3) = G(-3) = 0.0044318484

......

这些结果有什么意义呢？

当x = 0时，表示这一点就是原始中心点（center tap）, 这个点需要保留0.39（39%）的颜色

当x = 1以及x = -1时, 即有一个像素的位移时, 这一点需要保留0.24的颜色值

当x = 2以及x = -2时, 即有两个像素的位移时, 这一点需要保留0.05的颜色值

......

我们可以进行任意次这样的计算。 采样越高，精度越高。

相信现在你已经大致了解高斯模糊的原理了。

 

- 实际应用 -

当然，理论只是基于理想环境下的计算。在一次Bloom shader的测试中，我发现了两个问题：

1、模糊后的画面，其发光度和亮度大都比原始图像要低，画面显得很暗
2、上述的采样点（tap）的标准差是随意的吗？或者说，它有什么要求（或限制）？

以下各段，我将以一维高斯分布为例
首先解决明度问题。假设你有一个全白的单通道采样贴图，那么它所有的像素值都为1。经过第一次（水平方向）的高斯模糊运算之后，每个像素的值的总和，是所有高斯采样的加权总合。对于一个5 × 5的分布，标准差为1时，有：

ΣHBlur = G(0) + G(-1) + G(1) + G(-2) + G(2)

即

ΣHBlur = 0.9908656625

看起来似乎是正确的，其实不然。为了使明度不变，这个结果必须绝对是1。

基于此结果，在第二次运算后，它的值进一步下降为：ΣHVBlur = 0.9818147611

这就是说，我们必须增强加权值，通过确定的一个商，以取得正确的结果。答案非常简单，缺少了这个商，我们当然不能得到1：ΣHBlur^-1.

AugFactor = 1 / ΣHBlur

AugFactor = 1.009218543

如果我们再用这个增量，乘以该点的加权以及两个方向（水平、垂直）的运算，ΣHVBlur的值就是1了！

其次是标准差的问题。我使用的解决方案并不太“官方”，但是它确实有效。我们假设一个5 × 5的模糊必须有一个其值为1的标准差。现在我们知道，在一个5 × 5的分布中，距离中心最远的采样点，其加权值最小。其它尺寸的采样分布（比如7× 7），计算出各点的加权值也有这一规律（函数已经有了，你可以简单证明一下）。

要做的是，只调整标准差，以获取最后的采样点加权值，该值越接近这一结果越好。

σ = 1 -> G(2)aug = 0.054488685

σ = 1.745 -> G(3)aug = 0.054441945

σ = 2.7 -> G(4)aug = 0.054409483

这些都是我反复调试出的最接近的结果。所有这些权重增加，他们的和总是1。因此，寻找最佳标准差，最主要的是：尽可能多提取出模糊的采样数，而不产生偏差。

- 一维复合运算 | 二维单一运算 -

为什么会有一维和二维这两个函数呢？因为我们有两种选择：

1、使用一维函数的两次;一次横向，然后（以横向模糊的结果）再做一次垂直运算 ；

2、用二维函数，只需要一步就能完成。

你也许会问，既然有了二维函数，用它不就行了，何必要用一维函数进行两次计算，那多麻烦呀！我的第一反应是，这也许是因为早期Shader模型施加的一些限制。在做了更多的研究之后，我发现有一个更好的解释，来说明为什么要分开运算。

假设n是模糊的线性尺寸（或者说在水平/垂直方向的采样点数量），而p是整个画面的像素总数，那么：

在一维函数中，你要做2np次纹理采样

在一维函数中，你要做n²p次纹理采样

这意味着，如果你对一个512×512尺寸的图像做7×7采样的模糊，一位维函数需要做367万次搜索，而如果用二维函数，那么要完成这幅图像需要1284万次,是前者的3.5倍！

把一个二维坐标的点分解为两个一维的线性矢量是完全可能的，因为高斯模糊是一个可分割的卷积运算。具体的内容你可以到这里看看：http://nuttybar.drama.uga.edu/pipermail/dirgames-l/2002-December/020874.html

所以，为了追求更快的速度、更高的效率，我宁可放弃二维函数。

以上内容引自『听阿泽哥哥讲道理 http://blog.sina.com.cn/s/blog_4b97ab670100aa3a.html 』  

以下是自己C语言实现的一个gaussian函数

001

/**

002

* 高斯模糊算法

003

*

004	* @Author : Lei Guoting

005	* @Date   : 2012-07-04

006

*

007	* @Param inPixels 原始图片Pixels

008	* @Param outPixels 处理后输出图片Pixels

009	* @Circle 自定义结构体，定义一个需要被高斯模糊处理的园的位置大小

010	* @AndroidBitmapInfo 原始图片Bitmap信息

011

*/

012	void gaussian(void** inPixels,void** outPixels,Circle tGassCir,AndroidBitmapInfo bitmapInfo)

013

{

014	uint32_t x,y;

015	uint32_t tx,ty;

016	uint8_t *tInPixel, *tOutPixel;

017	double tR = 0.0,tG = 0.0,tB = 0.0,tCount = 0.0;

018	double ar = 8,dif,gauss;

019	int range = ar * 3;

020	int startx,endx,starty,endy;

021

int m;

022	int sx,sy;

023	uint32_t top,left,right,bottom;

024	LOGD("filter @gaussian center.x=%d,center.y=%d",tGassCir.center.x,tGassCir.center.y);

025	top = tGassCir.center.y - tGassCir.r;

026	bottom = tGassCir.center.y + tGassCir.r;

027	left = tGassCir.center.x - tGassCir.r;

028	right = tGassCir.center.x + tGassCir.r;

029	if(top < 0){

030

top = 0;

031

}

032	if(bottom > bitmapInfo.height - 1){

033	bottom = bitmapInfo.height - 1;

034

}

035	if(left < 0){

036

left = 0;

037

}

038	if(right > bitmapInfo.width - 1){

039	right = bitmapInfo.width - 1;

040

}

041	//The dimension of X

042	LOGD("filter @gaussian begin..");

043	LOGD("top=%d,bottom=%d,left=%d,right=%d",top,bottom,left,right);

044	Point tDo[(bottom - top +1) * (right - left +1)];

045	uint32_t tDoLen = 0;

046	for(y = top; y <= bottom; y ++){

047	for(x = left; x <= right; x ++){

048

/*

049	if(sqrt(pow((x - tGassCir.center.x),2) + pow((y - tGassCir.center.y),2)) > tGassCir.r){

050

continue;

051

}

052

*/

053

054	if(sqrt((x - tGassCir.center.x) * (x - tGassCir.center.x) + (y - tGassCir.center.y) * (y - tGassCir.center.y)) > tGassCir.r){

055

continue;

056

}

057	Point p = {x,y};

058	tDo[tDoLen] = p;

059	tDoLen ++;

060

061	tOutPixel = (uint8_t* )*outPixels + (y * bitmapInfo.width + x) * 4;

062	startx = (x < range ? -x : -range);

063	endx = ((bitmapInfo.width <= x + range) ? (bitmapInfo.width - x - 1) : range);

064

065	//x dimension

066	tR = tG = tB = tCount = 0.0;

067	for(m = startx; m <= endx; m ++){

068	tInPixel = (uint8_t* )*inPixels + (y * bitmapInfo.width + (x + m)) * 4;

069	dif = m*m;

070	gauss = exp(-dif/(2*ar*ar));

071	tR += gauss * *(tInPixel + 0);

072	tG += gauss * *(tInPixel + 1);

073	tB += gauss * *(tInPixel + 2);

074	tCount += gauss;

075

}

076	*(tOutPixel + 0) = tR / tCount;

077	*(tOutPixel + 1) = tG / tCount;

078	*(tOutPixel + 2) = tB / tCount;

079

}

080

}

081

082	//y dimension

083	uint32_t ind = 0;

084

Point p;

085	for(ind = 0; ind < tDoLen; ind ++){

086	p = tDo[ind];

087	tOutPixel = (uint8_t* )*outPixels + (p.y * bitmapInfo.width + p.x) * 4;

088	starty = (p.y < range ? -p.y : -range);

089	endy = ((bitmapInfo.height <= p.y + range) ? (bitmapInfo.height - p.y - 1) : range);

090	tR = tG = tB = tCount = 0.0;

091	for(m = starty; m <= endy; m ++){

092	tInPixel = (uint8_t* )*outPixels + ((p.y + m) * bitmapInfo.width + p.x) * 4;

093	dif = m*m;

094	gauss = exp(-dif/(2*ar*ar));

095	tR += gauss * *(tInPixel + 0);

096	tG += gauss * *(tInPixel + 1);

097	tB += gauss * *(tInPixel + 2);

098	tCount += gauss;

099

}

100	*(tOutPixel + 0) = tR / tCount;

101	*(tOutPixel + 1) = tG / tCount;

102	*(tOutPixel + 2) = tB / tCount;

103

}

104	LOGD("filter @gaussian end!");

105

}

1	typedef struct{

2	uint32_t x;

3	uint32_t y;

4

}Point;

5

6	typedef struct{

7	Point center;

8

int r;

9

}Circle;

gaussian加权的时候为了尽量快一点时间，稍微省略了点。

效果图: