学步园

下面的代码是计算二维高斯核，代码是用纯C写的，可移植性应该不错。计算得到的二维高斯核可以用来对图像进行高斯滤波。由于高斯滤波器是可分的，所以可以用一维高斯核得到二维高斯核（其实中间涉及到一些证明，这里直接借用结论），如代码所示。

//====================================================================      
// 作者   : quarryman      
// 邮箱   : quarrying{at}qq.com      
// 主页   : http://blog.csdn.net/quarryman      
// 日期   : 2013年12月22日      
// 描述   : 计算二维高斯核    
//====================================================================   
#include <stdio.h>
#include <memory.h>
#include <stdlib.h>
#include <math.h>
const double eps=1e-6;;

// 高斯函数
double guassian(double x,double sigma)
{
	return exp(-(x*x)/(2*sigma*sigma+eps));
}

// 计算高斯核
void guassiankernel(double* kernel,int size, double sigma)
{
	double sum=0;
	double* data=kernel;
	int i;
	for(i=0;i<size;++i)
	{
		double index=(size>>1)-i;
		if(size&1) // size为奇数
		{
			*(data+i)=guassian(index,sigma);
		}
		else // size为偶数
		{
			index-=0.5;
			*(data+i)=guassian(index,sigma);
		}
		sum+=*(data+i);
	}
	// 归一化
	for(i=0;i<size;++i)
	{
			*(data+i)/=sum;
	}
}

void matproduct(double a[],double b[],double c[],int m,int n,int p)
{
	int i,j,k;
	for(i=0;i<m;++i)
	{
		for(j=0;j<p;++j)
		{
			double sum=0;
			for(k=0;k<n;++k)
			{
				sum+=a[i*n+k]*b[k*p+j];
			}
			c[i*p+j]=sum;
		}
	}
}

void matshow(double a[],int m,int n)
{
	int i,j;
	for(i=0;i<m;++i)
	{
		for(j=0;j<n;++j)
		{
			printf("%lf\t",a[i*n+j]);
		}
		printf("\n");;
	}
}

// 计算二维高斯核
void guassiankernel2d(double* kernel,int sizeX,int sizeY,double sigmaX,double sigmaY)
{
	double *matX=(double *)malloc(sizeX*sizeof(double));
	double *matY=(double *)malloc(sizeY*sizeof(double));
	guassiankernel(matX,sizeX,sigmaX);
	guassiankernel(matY,sizeY,sigmaY);
	matproduct(matX,matY,kernel,sizeX,1,sizeY);
	free(matX);
	free(matY);
}

int main()
{
	int sizeX=4;
	int sizeY=3;
	double *kernel=(double *)malloc(sizeX*sizeY*sizeof(double));
	guassiankernel2d(kernel,sizeX,sizeY,1,1);
	matshow(kernel,sizeX,sizeY);
	free(kernel);
	return 0;
}