IIR型高斯滤波的原理及实现

二、实现

GIMP中有IIR型高斯滤波的实现，代码位于contrast-retinex.c中，读者可自行查看。下面给出本人实现的核心代码：

typedef struct
{
	float B;
	float b[4];
} gauss_coefs;

//参数计算
void compute_coefs3(gauss_coefs *c,float sigma)
{
	float q, q2, q3;
	if (sigma >= 2.5)
	{
		q = 0.98711 * sigma - 0.96330;
	}
	else if ((sigma >= 0.5) && (sigma < 2.5))
	{
		q = 3.97156 - 4.14554 * (float) sqrt ((float) 1 - 0.26891 * sigma);
	}
	else
	{
		q = 0.1147705018520355224609375;
	}
	q2 = q * q;
	q3 = q * q2;
	c->b[0] = 1/(1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
	c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3)) *c->b[0];
	c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3)))*c->b[0];
	c->b[3] = (                                 (0.422205*q3)) *c->b[0];
	c->B = 1.0-(c->b[1]+c->b[2]+c->b[3]);

}

//IIR型高斯滤波
//**************************************************************//
//参考文献：Recursive implementation of the Gaussian filter
//Src:输入图像
//Dest:输出图像
//sigma:高斯标准差
//**************************************************************//
IS_RET  IIRGaussianFilter(TMatrix *Src,TMatrix *Dest,float sigma)
{
	int X,Y,Width=Src->Width,Height=Src->Height,Stride=Src->WidthStep,Channel=Src->Channel;
	gauss_coefs  c;
	compute_coefs3(&c,sigma);
	unsigned char *LinePS,*LinePD;
	float *Buff,*BPointer,*w1,*w2;	
	Buff=(float*)malloc(sizeof(float)*Height*Width);//缓存
	if(Buff==NULL){return IS_RET_ERR_OUTOFMEMORY;}
	for(int i=0;i<Channel;i++)
	{
		LinePS=Src->Data+i;
		LinePD=Dest->Data+i;
		//拷贝原图至缓存Buff
		BPointer=Buff;
		for (Y=0;Y<Height;Y++)
		{
			for (X=0;X<Width;X++)
			{
				BPointer[0]=float(LinePS[0]);
				BPointer++;
				LinePS+=Channel;
			}
			LinePS+=Stride-Width*Channel;
		}
		//横向滤波
		BPointer=Buff;
		w1=(float*)malloc(sizeof(float)*(Width+3));
		if(w1==NULL){return IS_RET_ERR_OUTOFMEMORY;}
		w2=(float*)malloc(sizeof(float)*(Width+3));
		if(w2==NULL){return IS_RET_ERR_OUTOFMEMORY;}
		for(Y=0;Y<Height;Y++)
		{
			//前向滤波
			w1[0]=w1[1]=w1[2]=BPointer[0];
			for(int n=3,i=0;i<Width;n++,i++)
			{
				w1[n]=c.B*BPointer[i]+(c.b[1]*w1[n-1]+c.b[2]*w1[n-2]+c.b[3]*w1[n-3]);
			}
			//后向滤波
			w2[Width]=w2[Width+1]=w2[Width+2]=w1[Width+2];
			for(int n=Width-1;n>=0;n--)
			{
				BPointer[n]=w2[n]=c.B*w1[n+3]+(c.b[1]*w2[n+1]+c.b[2]*w2[n+2]+c.b[3]*w2[n+3]);
			}
			BPointer+=Width;
		}

		//纵向滤波
		BPointer=Buff;
		w1=(float*)realloc(w1,sizeof(float)*(Height+3));
		if(w1==NULL){return IS_RET_ERR_OUTOFMEMORY;}
		w2=(float*)realloc(w2,sizeof(float)*(Height+3));
		if(w2==NULL){return IS_RET_ERR_OUTOFMEMORY;}
		for (X=0;X<Width;X++)
		{
			//前向滤波
			w1[0]=w1[1]=w1[2]=BPointer[0];
			for(int n=3,i=0;i<Height;n++,i++)
			{
				w1[n]=c.B*BPointer[i*Width]+(c.b[1]*w1[n-1]+c.b[2]*w1[n-2]+c.b[3]*w1[n-3]);
			}
			//后向滤波
			w2[Height]=w2[Height+1]=w2[Height+2]=w1[Height+2];
			for(int n=Height-1;n>=0;n--)
			{
				BPointer[n*Width]=w2[n]=c.B*w1[n+3]+(c.b[1]*w2[n+1]+c.b[2]*w2[n+2]+c.b[3]*w2[n+3]);
			}
			BPointer++;
		}
		//拷贝缓存至结果
		BPointer=Buff;
		for(Y=0;Y<Height;Y++)
		{
			for(X=0;X<Width;X++)
			{
				LinePD[0]=BPointer[0];
				if(BPointer[0]>255)LinePD[0]=255;
				if(BPointer[0]<0)LinePD[0]=0;	
				BPointer++;
				LinePD+=Channel;
			}
			LinePD+=Stride-Width*Channel;
		}
		free(w1);
		free(w2);
	}
	return IS_RET_OK;
}

实验结果：对一幅1024*1024的彩色图像，算法耗时175ms。

参考文献：

Young I T, Van Vliet L J. Recursive implementation of the Gaussian filter[J]. Signal processing, 1995, 44(2): 139-151.