高斯算法实现

高斯算法的原理

首先，高斯滤波算法的一般过程分为两步：

计算掩膜(高斯核)
卷积（即掩膜上每一个位置的值和图像对应位置的像素值的乘积、求和运算）

其次，我们知道高斯分布也叫做正态分布；

在二维空间中，这个公式生成的曲面的等高线是从中心开始呈正态分布的同心圆。分布不为零的像素组成的卷积矩阵与原始图像做变换。每个像素的值都是周围相邻像素值的加权平均。原始像素的值有最大的高斯分布值，所以有最大的权重，相邻像素随着距离原始像素越来越远，其权重也越来越小。

通过高斯分布公式计算出需要使用的高斯核（GaussKern），也就是计算出各个位置的所需的权重大小，使得不同位置的像素值影响度不同，距离中心点的越远的位置，权重越小，反之越近的权重越大。

int* Gauss::buildGaussKern(int winSize, int sigma)
{

int wincenter, x;
float sum = 0.0f;
//计算中心点大小
wincenter = winSize / 2;
//kern用来存储高斯模糊前的数据
//ikern用来存储高斯模糊后与像素值256的乘积值
float *kern = (float*)malloc(winSize*sizeof(float));
int *ikern = (int*)malloc(winSize*sizeof(int));
//计算高斯分布公式的系数
float SQRT_2PI = 2.506628274631f;
float sigmaMul2PI = 1.0f / (sigma * SQRT_2PI);
float divSigmaPow2 = 1.0f / (2.0f * sigma * sigma);

for (x = 0; x < wincenter + 1; x++)
{
    kern[wincenter - x] = kern[wincenter + x] = exp(-(x * x)* divSigmaPow2) * sigmaMul2PI;
    sum += kern[wincenter - x] + ((x != 0) ? kern[wincenter + x] : 0.0);
}
sum = 1.0f / sum;

for (x = 0; x < winSize; x++)
{
    kern[x] *= sum;
    ikern[x] = kern[x] * 256.0f;
}

free(kern);
return ikern;
}

理论上来讲，图像中每点的分布都不为零，这也就是说每个像素的计算都需要包含整幅图像。在实际应用中，在计算高斯函数的离散近似时，在大概3σ距离之外的像素都可以看作不起作用，这些像素的计算也就可以忽略。通常，图像处理程序只需要计算 (6 \sigma + 1) \times (6 \sigma + 1) 的矩阵就可以保证相关像素影响。对于边界上的点，通常采用复制周围的点到另一面再进行加权平均运算。

//忽略3sigma外的像素点
//计算（6*sigma+1）*（6*sigma+1）的矩阵
unsigned int  winsize = (1 + (((int)ceil(3 * sigma)) * 2));

除了圆形对称之外，高斯模糊也可以在二维图像上对两个独立的一维空间分别进行计算，这叫作线性可分。这也就是说，使用二维矩阵变换得到的效果也可以通过在水平方向进行一维高斯矩阵变换加上竖直方向的一维高斯矩阵变换得到。从计算的角度来看，这是一项有用的特性，因为这样只需要O(n * M * N) + O(m * M * N)次计算，而不可分的矩阵则需要O(m * n * M * N)次计算，其中M,N是需要进行滤波的图像的维数，m、n是滤波器的维数。

     for (unsigned int w = 0; w < width; w += channels)
    {
        //存储r g b分量的数据
        unsigned int rowR = 0;
        unsigned int rowG = 0;
        unsigned int rowB = 0;
        unsigned int alpha = 0;
        //将高斯系数赋值给gaussKernPtr
        int * gaussKernPtr = gaussKern;
        int whalfsize = w + width - halfsize;
        //当前位置
        unsigned int  curPos = rowWidth + w;
        for (unsigned int k = 1; k < winsize; k += channels)
        {
            unsigned int  pos = rowWidth + ((k + whalfsize) % width);
            int fkern = *gaussKernPtr++;
            //当前像素值乘以高斯系数，rowR G B这了泛指单通道的当前像素点高斯处理后的数
            rowR += (pixels[pos + MT_RED] * fkern);
            rowG += (pixels[pos + MT_GREEN] * fkern);
            rowB += (pixels[pos + MT_BLUE] * fkern);
            alpha += (pixels[pos + MT_ALPHA] * fkern);
        }

算法实现部分

void Gauss::GaussBlur(unsigned char *pixels, unsigned int width, unsigned int height, unsigned int channels, int sigma)
{
//四个通道。 存储每行数据的宽度
width = 4 * width;

if ((width % 4) != 0)
   width += (4 - (width % 4));
//窗的大小
//忽略3sigma外的像素点
//计算（6*sigma+1）*（6*sigma+1）的矩阵
unsigned int  winsize = (1 + (((int)ceil(3 * sigma)) * 2));
//构建高斯核，计算高斯系数
int *gaussKern = buildGaussKern(winsize, sigma);

winsize *= 3;
//窗的边到中心的距离
unsigned int  halfsize = winsize / 2;
//开辟新的内存存储处理高斯模糊后的数据
unsigned char *tmpBuffer = (unsigned char*)malloc(width * height* sizeof(unsigned char));
//外层循环，图像的高度
for (unsigned int h = 0; h < height; h++)
{
    //当前行的宽度为图像的高度乘以每行图像的数据所占的宽度。
    //因为是按行存储的数组
    unsigned int  rowWidth = h * width;
    
    for (unsigned int w = 0; w < width; w += channels)
    {
        //存储r g b alpha分量的数据
        unsigned int rowR = 0;
        unsigned int rowG = 0;
        unsigned int rowB = 0;
        unsigned int alpha = 0;
        //将高斯系数赋值给gaussKernPtr
        int * gaussKernPtr = gaussKern;
        int whalfsize = w + width - halfsize;
        //当前位置
        unsigned int  curPos = rowWidth + w;
        for (unsigned int k = 1; k < winsize; k += channels)
        {
            unsigned int  pos = rowWidth + ((k + whalfsize) % width);
            int fkern = *gaussKernPtr++;
            //当前像素值乘以高斯系数，rowR G B这了泛指单通道的当前像素点高斯处理后的数
            rowR += (pixels[pos + MT_RED] * fkern);
            rowG += (pixels[pos + MT_GREEN] * fkern);
            rowB += (pixels[pos + MT_BLUE] * fkern);
            alpha += (pixels[pos + MT_ALPHA] * fkern);
        }
        //除以256
        tmpBuffer[curPos + MT_RED] = ((unsigned char)(rowR >> 8));
        tmpBuffer[curPos + MT_GREEN] = ((unsigned char)(rowG >> 8));
        tmpBuffer[curPos + MT_BLUE] = ((unsigned char)(rowB >> 8));
        tmpBuffer[curPos + MT_ALPHA] = ((unsigned char)(alpha >> 8));
    }
}
winsize /= 3;
//窗的边到中心的距离
halfsize = winsize / 2;

for (unsigned int w = 0; w < width; w++)
{
    for (unsigned int h = 0; h < height; h++)
    {
        unsigned int col_all = 0;
        int hhalfsize = h + height - halfsize;
        for (unsigned int k = 0; k < winsize; k++)
        {
            col_all += tmpBuffer[((k + hhalfsize) % height)* width + w] * gaussKern[k];
        }
        pixels[h * width + w] = (unsigned char)(col_all >> 8);
    }
}
free(tmpBuffer);
free(gaussKern);
}