高斯滤波核计算 C版

#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=3;
    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;
}