快速高斯滤波函数[修正完善版]

原文地址：http://blog.csdn.net/markl22222/article/details/10313565

进行了修正和变量优化。原来作者的函数只支持2次方图片，这次修正了（windows的bitmap行宽是4字节对齐的）。

基本完善了，但是在某些条件下，Y方向的底边还是会出现偏差，一时找不到原因，暂且发表，希望有人能提醒一下。

函数结构我规整了一下，很清晰，很好阅读。

int gauss_blur(
	byte_t* image,	//位图数据
	int linebytes,	//位图行字节数，BMP数据在windows中是4字节对齐的。否则在处理非二次幂的图像时会有偏差
	int width,		//位图宽度
	int height,		//位图高度
	int cbyte,		//颜色通道数量
	float sigma		//高斯系数
	)
{
	int x = 0, y = 0, n = 0;
	int channel = 0;
	int srcline = 0, dstline = 0;
	int channelsize = width*height;
	int bufsize = width > height ? width + 4 : height + 4;
	float *w1 = NULL, *w2 = NULL, *imgbuf = NULL;
	int time = 0;

	#if defined(_INC_WINDOWS)
		time = GetTickCount();
	#elif defined(_CLOCK_T)
		time  = clock();
	#endif

	w1 = (float*)malloc(bufsize * sizeof(float));
	if(!w1)
	{
		return -1;
	}
	w2 = (float*)malloc(bufsize * sizeof(float));
	if(!w2)
	{
		free(w1);
		return -1;
	}
	imgbuf = (float*)malloc(channelsize * sizeof(float));
	if(!imgbuf)
	{
		free(w1);
		free(w2);
		return -1;
	}

//----------------计算高斯核---------------------------------------//
	float q  = 0;
	float q2 = 0, q3 = 0;
	float b0 = 0, b1 = 0, b2 = 0, b3 = 0;
	float B  = 0;

	if (sigma >= 2.5f)
	{
		q = 0.98711f * sigma - 0.96330f;
	}
	else if ((sigma >= 0.5f) && (sigma < 2.5f))
	{
		q = 3.97156f - 4.14554f * (float) sqrt (1.0f - 0.26891f * sigma);
	}
	else
	{
		q = 0.1147705018520355224609375f;
	}

	q2 = q * q;
	q3 = q * q2;
	b0 = (1.57825+ (2.44413f*q)+(1.4281f *q2)+(0.422205f*q3));
	b1 = (         (2.44413f*q)+(2.85619f*q2)+(1.26661f* q3));
	b2 = (                     -((1.4281f*q2)+(1.26661f* q3)));
	b3 = (                                    (0.422205f*q3));
	B = 1.0-((b1+b2+b3)/b0);

	b1 /= b0;
	b2 /= b0;
	b3 /= b0;
//----------------计算高斯核结束---------------------------------------//

<span style="white-space:pre">	</span>// 处理图像的多个通道
	for (channel = 0; channel < cbyte; ++channel)
	{
		// 获取一个通道的所有像素值，并预处理
		for(y=0; y<height; ++y)
		{
			srcline = y*linebytes;
			dstline = y*width;
			for(x=0, n=channel; x<width; ++x, n+=cbyte)
			{
				(imgbuf+dstline)[x] = float((image+srcline)[n]) + 1.0f;
			}
		}

		for (int x=0; x<width; ++x)
		{//横向处理

			w1[0] = (imgbuf + x)[0];
			w1[1] = (imgbuf + x)[0];
			w1[2] = (imgbuf + x)[0];

			for (y=0; y<height; ++y)
			{
				w1[y+3] = B*(imgbuf + x)[y*width] + (b1*w1[y+2] + b2*w1[y+1] + b3*w1[y+0]);
			}

			w2[width+0]= w1[width+2];
			w2[width+1]= w1[width+1];
			w2[width+2]= w1[width+0];

			for (int y=height-1; y>=0; --y)
			{
				(imgbuf + x)[y*width] = w2[y] = B*w1[y+3] + (b1*w2[y+1] + b2*w2[y+2] + b3*w2[y+3]);
			}
		}//横向处理

		for (y=0 ; y<height; ++y)
		{//纵向处理
			srcline = y * width;
			dstline = y * linebytes;

			//取当前行数据
			w1[0] = (imgbuf + srcline)[0];
			w1[1] = (imgbuf + srcline)[0];
			w1[2] = (imgbuf + srcline)[0];

			//正方向横向处理3个点的数据
			for (x=0; x<width ; ++x)
			{
				w1[x+3] = B*(imgbuf + srcline)[x] + (b1*w1[x+2] + b2*w1[x+1] + b3*w1[x+0]);
			}

			w2[width+0]= w1[width+2];
			w2[width+1]= w1[width+1];
			w2[width+2]= w1[width+0];

			//反方向处理
			for (x=width-1; x>=0; --x)
			{
				//(imgbuf + dstline)[x] = w2[x] = B*w1[x+3] + (b1*w2[x+1] + b2*w2[x+2] + b3*w2[x+3]);
				w2[x] = B*w1[x+3] + (b1*w2[x+1] + b2*w2[x+2] + b3*w2[x+3]);

				//存储返回数据
				(image + dstline)[x * cbyte + channel] = w2[x]-1;
			}

		}//纵向处理

		/*
		//存储处理完毕的通道
		for(int y=0; y<height; y++)
		{
			int dstline = y*linebytes;
			int srcline = y*width;
			for (int x=0; x<width; x++)
			{
				//(image + dstline)[x * cbyte + channel] = (imgbuf + srcline)[x]-1;
					//byte_comp((imgbuf + srcline)[x]-1);

			}
		}//存储循环
		//*/
	}//通道循环

	free (w1);
	free (w2);
	free(imgbuf);

	#if defined(_INC_WINDOWS)
		return GetTickCount() - time;
	#elif defined(_CLOCK_T)
		return clock() - time;
	#else
		return 0;
	#endif
}

快速高斯滤波函数[修正完善版]