图像快速高斯滤波实现

1、二维高斯滤波公式：(μ=0,σ）

　　G(x,y) = 1/2πσ²·exp(-(x²+y²)/2σ²)

　　G(x,y)可以分离为：G(x,y) = G(x)·G(y)= 1/2πσ²·exp(-x²/2σ²)·exp(-y²/2σ²)，其中系数1/2πσ²可以不考虑，因为在归一化时又将其消掉了。

　　所以，二维的滤波可以分离为两个一维的滤波。

2、假设灰度图像宽高分别为W、H，二维高斯滤波模板宽高均为m_size，那么二维模板直接卷积计算复杂度为O(W*H*m_size*m_size)；

　　同样采用两次一维（横向、纵向）卷积计算的复杂度为O(W*H*m_szie*2)，显然计算复杂度降低，滤波模板阶数越高计算速度的提升越明显。

　　通常根据3σ原则来计算模板阶数，阶数k = 2*ceil(3σ)+1，因为高斯分布中(μ=0)±3σ几乎涵盖了所有的分布，所以可以根据此规则估算合理的模板阶数（奇数）。

3、代码实现中这里使用了自己创建的二维数组，目的是来比对两种滤波的结果，代码如下：

#include <QCoreApplication>
#include <iostream>
#include <math.h>
using namespace std;

double** CreateGaussMap(int map_size, double sigma);
double* CreateGauss(int map_size, double sigma);
void GaussFilter2(int src[9][9], double dst[9][9], double** map);
void GaussFilter(int src[9][9], double dst[9][9], double* map);

int main(int argc, char *argv[])
{
    QCoreApplication a(argc, argv);

    double** gauss2 = CreateGaussMap(3,0.33);
    for (int i=0;i<3;++i) {
        for (int j=0;j<3;++j) {
            cout<<gauss2[i][j]<< " ";
        }
        cout<<endl;
    }
    cout<<endl;

    double* gauss = CreateGauss(3,0.33);
    for (int i=0;i<3;++i) {
        cout<<gauss[i]<< " ";
    }
    cout<<endl;

    int arr[9][9]={{1,1,1,1,1,1,1,1,1},{2,2,2,2,2,2,2,2,2},
                   {3,3,3,3,3,3,3,3,3},{4,4,4,4,4,4,4,4,4},
                   {5,5,5,5,5,5,5,5,5},{6,6,6,6,6,6,6,6,6},
                   {7,7,7,7,7,7,7,7,7},{8,8,8,8,8,8,8,8,8},
                   {9,9,9,9,9,9,9,9,9}};
    double dst2[9][9];
    double dst[9][9];
    GaussFilter2(arr,dst2,gauss2);

    for (int i=0;i<9;++i) {
        for (int j=0;j<9;++j) {
            cout<<dst2[i][j]<<" ";
        }
        cout<<endl;
    }
    GaussFilter(arr,dst,gauss);
    for (int i=0;i<9;++i) {
        for (int j=0;j<9;++j) {
            cout<<dst[i][j]<<" ";
        }
        cout<<endl;
    }

    return a.exec();
}


double** CreateGaussMap(int map_size, double sigma){
    // g(x,y) = (1/(2*pi*sigma^2))*exp(-(x^2+y^2)/(2*sigma^2))
    int center = (map_size-1)/2;
    //double amplitude = 1/(2*M_PI*sigma*sigma);
    double sigma_sqr = 2*pow(sigma,2);
    double x_sqr = 0, y_sqr = 0;
    double** gauss = new double*[map_size];
    double sum = 0;

//    for (int i=0;i<map_size;++i) {

//    }

    for(int i=0;i<map_size;++i){
        gauss[i] = new double[map_size];
        y_sqr = pow(double(i-center),2);
        for (int j=0;j<map_size;++j) {
            x_sqr = pow(double(j-center),2);
            sum += gauss[i][j] = exp(-(x_sqr + y_sqr)/sigma_sqr);
        }
    }

    if(sum != 0){
        for (int i=0;i<map_size;++i) {
            for (int j=0;j<map_size;++j) {
                gauss[i][j] /= sum;
            }
        }
    }

    return gauss;
}

double* CreateGauss(int map_size, double sigma){
    // g(x,y) = (1/(2*pi*sigma^2))*exp(-(x^2+y^2)/(2*sigma^2))
    int radius = (map_size-1)/2;
    //double amplitude = 1/(2*M_PI*sigma*sigma);
    double sigma_sqr = 2*pow(sigma,2);
    double* gauss = new double [map_size];
    double sum = 0;

    for(int i=0;i<map_size;++i){
        gauss[i] = exp(-(pow(double(i-radius),2))/sigma_sqr);
        sum += gauss[i];
    }

    if(sum != 0){
        for (int i=0;i<map_size;++i) {
            gauss[i] /= sum;
        }
    }

    return gauss;
}

void GaussFilter2(int src[9][9], double dst[9][9], double** map){

    for (int i=0;i<9;++i) {
        for (int j=0;j<9;++j) {
            dst[i][j] = src[i][j];
        }
    }

    int center = 1;
    for (int i=0;i<9-center-1;++i) {
        for (int j=0;j<9-center-1;++j) {
            for (int k=0;k<3;++k) {
                for (int m=0;m<3;++m) {
                    dst[i+center][j+center] += src[i+k][j+m]*map[k][m];
                }
            }
        }
    }
}

void GaussFilter(int src[9][9], double dst[9][9], double* map){
    double tmp[9][9];
    for (int i=0;i<9;++i) {
        for (int j=0;j<9;++j) {
            tmp[i][j] = src[i][j];
        }
    }

    int center = 1;
    int a=0;
    double t=0;
    for (int i=0;i<9-center-1;++i) {
        for (int j=0;j<9-center-1;++j) {
            tmp[i+center][j+center] = 0;
            for (int k=0;k<3;++k) {
                a = src[i+center][j+k];
                t=map[k];
                tmp[i+center][j+center] += (src[i+center][j+k])*map[k];
            }
        }
    }

    for (int i=0;i<9;++i) {
        for (int j=0;j<9;++j) {
            dst[i][j] = tmp[i][j];
        }
    }

    for (int j=0;j<9-center-1;++j) {
        for (int i=0;i<9-center-1;++i) {
            for (int k=0;k<3;++k) {
                dst[i+center][j+center] += tmp[i+k][j+center]*map[k];
            }
        }
    }
}

4、结果中，两中方式滤波结果一致，μ=0，σ=0.33，k = 3。

二维高斯滤波模板：

9.87532e-005  0.00973996  9.87532e-005
0.00973996     0.960645      0.00973996
9.87532e-005  0.00973996  9.87532e-005

对9*9的二维数组滤波结果：
1 1   1   1   1   1   1   1   1
2 4   4   4   4   4   4   4   2
3 6   6   6   6   6   6   6   3
4 8   8   8   8   8   8   8   4
5 10 10 10 10 10 10 10 5
6 12 12 12 12 12 12 12 6
7 14 14 14 14 14 14 14 7
8 16 16 16 16 16 16 16 8
9 9   9   9   9   9   9   9   9

 

一维高斯滤波模板

[0.00993747  0.980125  0.00993747]

两次一维高斯滤波的结果：（注：第一次横向卷积计算后的数据进行第二次纵向卷积计算）
1 1   1   1   1   1   1   1   1
2 4   4   4   4   4   4   4   2
3 6   6   6   6   6   6   6   3
4 8   8   8   8   8   8   8   4
5 10 10 10 10 10 10 10 5
6 12 12 12 12 12 12 12 6
7 14 14 14 14 14 14 14 7
8 16 16 16 16 16 16 16 8
9 9   9   9   9   9   9   9   9