下面的代码是计算二维高斯核,代码是用纯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;
}