高斯算法的原理

首先，高斯濾波算法的一般過程分為兩步：

1. 計算掩膜(高斯核)
2. 卷積（即掩膜上每一個位置的值和圖像對應位置的像素值的乘積、求和運算）

其次，我們知道高斯分布也叫做正態分布；

在二維空間中，這個公式生成的曲面的等高線是從中心開始呈正態分布的同心圓。分布不為零的像素組成的卷積矩陣與原始圖像做變換。每個像素的值都是周圍相鄰像素值的加權平均。原始像素的值有最大的高斯分布值，所以有最大的權重，相鄰像素隨著距離原始像素越來越遠，其權重也越來越小。

通過高斯分布公式計算出需要使用的高斯核（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);
}

高斯算法實現