高斯模糊實現小結

http://blog.csdn.net/zddblog/article/details/7450033

注：博客中圖表的大小難以調整，導致閱讀不便，這裡有其pdf版本：高斯模糊實現小結.pdf

 

高斯模糊是一種圖像濾波器，它使用正態分布(高斯函數)計算模糊模板，並使用該模板與原圖像做卷積運算，達到模糊圖像的目的。

N維空間正態分布方程為：

其中，σ是正態分布的標准差，σ值越大，圖像越模糊(平滑)。r為模糊半徑，模糊半徑是指模板元素到模板中心的距離。如二維模板大小為m*n，則模板上的元素(x,y)對應的高斯計算公式為：

在二維空間中，這個公式生成的曲面的等高線是從中心開始呈正態分布的同心圓。分布不為零的像素組成的卷積矩陣與原始圖像做變換。每個像素的值都是周圍相鄰像素值的加權平均。原始像素的值有最大的高斯分布值，所以有最大的權重，相鄰像素隨著距離原始像素越來越遠，其權重也越來越小。這樣進行模糊處理比其它的均衡模糊濾波器更高地保留了邊緣效果。

理論上來講，圖像中每點的分布都不為零，這也就是說每個像素的計算都需要包含整幅圖像。在實際應用中，在計算高斯函數的離散近似時，在大概3σ距離之外的像素都可以看作不起作用，這些像素的計算也就可以忽略。通常，圖像處理程序只需要計算(6σ+1)*(6σ+1)的矩陣就可以保證相關像素影響。

1、使用給定高斯模板平滑圖像函數

σ=0.84089642的7行7列高斯模糊矩陣為：

現使用該模板對源圖像做模糊處理，其函數如下：

/高斯平滑  
//未使用sigma，邊緣無處理  
void GaussianTemplateSmooth(const Mat &src, Mat &dst, double sigma)  
{  
    //高斯模板(7*7)，sigma = 0.84089642，歸一化後得到  
    static const double gaussianTemplate[7][7] =   
    {  
        {0.00000067, 0.00002292, 0.00019117, 0.00038771, 0.00019117, 0.00002292, 0.00000067},  
        {0.00002292, 0.00078633, 0.00655965, 0.01330373, 0.00655965, 0.00078633, 0.00002292},  
        {0.00019117, 0.00655965, 0.05472157, 0.11098164, 0.05472157, 0.00655965, 0.00019117},  
        {0.00038771, 0.01330373, 0.11098164, 0.22508352, 0.11098164, 0.01330373, 0.00038771},  
        {0.00019117, 0.00655965, 0.05472157, 0.11098164, 0.05472157, 0.00655965, 0.00019117},  
        {0.00002292, 0.00078633, 0.00655965, 0.01330373, 0.00655965, 0.00078633, 0.00002292},  
        {0.00000067, 0.00002292, 0.00019117, 0.00038771, 0.00019117, 0.00002292, 0.00000067}  
    };  
  
    dst.create(src.size(), src.type());  
    uchar* srcData = src.data;  
    uchar* dstData = dst.data;  
  
    for(int j = 0; j < src.cols-7; j++)  
    {  
        for(int i = 0; i < src.rows-7; i++)  
        {  
            double acc = 0;  
            double accb = 0, accg = 0, accr = 0;   
            for(int m = 0; m < 7; m++)  
            {  
                for(int n = 0; n < 7; n++)  
                {  
                    if(src.channels() == 1)  
                        acc += *(srcData + src.step * (i+n) + src.channels() * (j+m)) * gaussianTemplate[m][n];  
                    else  
                    {  
                        accb += *(srcData + src.step * (i+n) + src.channels() * (j+m) + 0) * gaussianTemplate[m][n];  
                        accg += *(srcData + src.step * (i+n) + src.channels() * (j+m) + 1) * gaussianTemplate[m][n];  
                        accr += *(srcData + src.step * (i+n) + src.channels() * (j+m) + 2) * gaussianTemplate[m][n];  
                    }  
                }  
            }  
            if(src.channels() == 1)  
                *(dstData + dst.step * (i+3) + dst.channels() * (j+3))=(int)acc;  
            else  
            {  
                *(dstData + dst.step * (i+3) + dst.channels() * (j+3) + 0)=(int)accb;  
                *(dstData + dst.step * (i+3) + dst.channels() * (j+3) + 1)=(int)accg;  
                *(dstData + dst.step * (i+3) + dst.channels() * (j+3) + 2)=(int)accr;  
            }  
        }  
    }  
      
}  

其效果如圖1所示，7*7的高斯模板與源圖像做卷積運算時，會產生半徑為3的邊緣，在不精確的圖像處理中，可用源圖像像素填充，或者去掉邊緣。

2、二維高斯模糊函數

上述的例子中，如何求得高斯模糊矩陣是高斯模糊的關鍵。根據高斯函數的性質，圖像處理程序只需要計算(6σ+1)*(6σ+1)的矩陣就可以保證相關像素影響。因此，可根據σ的值確定高斯模糊矩陣的大小。高斯矩陣可利用公式(1-2)計算，並歸一化得到。歸一化是保證高斯矩陣的值在[0,1]之間，

其處理函數如下：

void GaussianSmooth2D(const Mat &src, Mat &dst, double sigma)  
{  
    if(src.channels() != 1)  
        return;  
  
    //確保sigma為正數   
    sigma = sigma > 0 ? sigma : 0;  
    //高斯核矩陣的大小為(6*sigma+1)*(6*sigma+1)  
    //ksize為奇數  
    int ksize = cvRound(sigma * 3) * 2 + 1;  
      
//  dst.create(src.size(), src.type());  
    if(ksize == 1)  
    {  
        src.copyTo(dst);      
        return;  
    }  
  
    dst.create(src.size(), src.type());  
  
    //計算高斯核矩陣  
    double *kernel = new double[ksize*ksize];  
  
    double scale = -0.5/(sigma*sigma);  
    const double PI = 3.141592653;  
    double cons = -scale/PI;  
  
    double sum = 0;  
  
    for(int i = 0; i < ksize; i++)  
    {  
        for(int j = 0; j < ksize; j++)  
        {  
            int x = i-(ksize-1)/2;  
            int y = j-(ksize-1)/2;  
            kernel[i*ksize + j] = cons * exp(scale * (x*x + y*y));  
  
            sum += kernel[i*ksize+j];  
//          cout << " " << kernel[i*ksize + j];  
        }  
//      cout <<endl;  
    }  
    //歸一化  
    for(int i = ksize*ksize-1; i >=0; i--)  
    {  
        *(kernel+i) /= sum;  
    }  
    //圖像卷積運算，無邊緣處理  
    for(int j = 0; j < src.cols-ksize; j++)  
    {  
        for(int i = 0; i < src.rows-ksize; i++)  
        {  
            double acc = 0;  
  
            for(int m = 0; m < ksize; m++)  
            {  
                for(int n = 0; n < ksize; n++)  
                {  
                    acc += *(srcData + src.step * (i+n) + src.channels() * (j+m)) * kernel[m*ksize+n];   
                }  
            }  
  
          
            *(dstData + dst.step * (i + (ksize - 1)/2) + (j + (ksize -1)/2)) = (int)acc;  
        }  
    }  
    delete []kernel;  
}  

利用該函數，取σ=0.84089642即可得到上例中7*7的模板，該模板數據存在kernel中。利用該函數計算的歸一化後的3*3，5*5階高斯模板如表2,3所示：

由上表可以看出，高斯模板是中心對稱的。

模糊效果如圖2所示。

 

對圖2中邊緣的處理：

int center = (ksize-1) /2;  
//圖像卷積運算,處理邊緣  
for(int j = 0; j < src.cols; j++)  
{  
    for(int i = 0; i < src.rows; i++)  
    {  
        double acc = 0;  
  
        for(int m = -center, c = 0; m <= center; m++, c++)  
        {  
            for(int n = -center, r = 0; n <= center; n++, r++)  
            {  
                if((i+n) >=0 && (i+n) < src.rows && (j+m) >=0 && (j+m) < src.cols)  
                {  
                      
                    acc += *(srcData + src.step * (i+n) + src.channels() * (j+m)) * kernel[r*ksize+c];   
              
                }  
            }  
        }  
  
  
        *(dstData + dst.step * (i) + (j)) = (int)acc;  
    }  
}  

結果圖為：

如上圖所示，邊緣明顯變窄，但是存在黑邊。

3、改進的高斯模糊函數

上述的二維高斯模糊函數GaussianSmooth2D達到了高斯模糊圖像的目的，但是如圖2所示，會因模板的關系而造成邊緣圖像缺失，σ越大，缺失像素越多，額外的邊緣處理會增加計算量。並且當σ變大時，高斯模板(高斯核)和卷積運算量將大幅度提高。根據高斯函數的可分離性，可對二維高斯模糊函數進行改進。

高斯函數的可分離性是指使用二維矩陣變換得到的效果也可以通過在水平方向進行一維高斯矩陣變換加上豎直方向的一維高斯矩陣變換得到。從計算的角度來看，這是一項有用的特性，因為這樣只需要O(n*M*n)+O(m*M*N)次計算，而二維不可分的矩陣則需要O(m*n*M*n)次計算，其中，m,n為高斯矩陣的維數，M,N為二維圖像的維數。

另外，兩次一維的高斯卷積將消除二維高斯矩陣所產生的邊緣。

(關於消除邊緣的論述如下圖2.4所示， 對用模板矩陣超出邊界的部分——虛線框，將不做卷積計算。如圖2.4中x方向的第一個模板1*5，將退化成1*3的模板，只在圖像之內的部分做卷積。)

改進的高斯模糊函數如下：

void GaussianSmooth(const Mat &src, Mat &dst, double sigma)  
{  
    if(src.channels() != 1 && src.channels() != 3)  
        return;  
  
    //  
    sigma = sigma > 0 ? sigma : -sigma;  
    //高斯核矩陣的大小為(6*sigma+1)*(6*sigma+1)  
    //ksize為奇數  
    int ksize = ceil(sigma * 3) * 2 + 1;  
  
    //cout << "ksize=" <<ksize<<endl;  
    //  dst.create(src.size(), src.type());  
    if(ksize == 1)  
    {  
        src.copyTo(dst);      
        return;  
    }  
  
    //計算一維高斯核  
    double *kernel = new double[ksize];  
  
    double scale = -0.5/(sigma*sigma);  
    const double PI = 3.141592653;  
    double cons = 1/sqrt(-scale / PI);  
  
    double sum = 0;  
    int kcenter = ksize/2;  
    int i = 0, j = 0;  
    for(i = 0; i < ksize; i++)  
    {  
        int x = i - kcenter;  
        *(kernel+i) = cons * exp(x * x * scale);//一維高斯函數  
        sum += *(kernel+i);  
  
//      cout << " " << *(kernel+i);  
    }  
//  cout << endl;  
    //歸一化,確保高斯權值在[0,1]之間  
    for(i = 0; i < ksize; i++)  
    {  
        *(kernel+i) /= sum;  
//      cout << " " << *(kernel+i);  
    }  
//  cout << endl;  
  
    dst.create(src.size(), src.type());  
    Mat temp;  
    temp.create(src.size(), src.type());  
  
    uchar* srcData = src.data;  
    uchar* dstData = dst.data;  
    uchar* tempData = temp.data;  
  
    //x方向一維高斯模糊  
    for(int y = 0; y < src.rows; y++)  
    {  
        for(int x = 0; x < src.cols; x++)  
        {  
            double mul = 0;  
            sum = 0;  
            double bmul = 0, gmul = 0, rmul = 0;  
            for(i = -kcenter; i <= kcenter; i++)  
            {  
                if((x+i) >= 0 && (x+i) < src.cols)  
                {  
                    if(src.channels() == 1)  
                    {  
                        mul += *(srcData+y*src.step+(x+i))*(*(kernel+kcenter+i));  
                    }  
                    else   
                    {  
                        bmul += *(srcData+y*src.step+(x+i)*src.channels() + 0)*(*(kernel+kcenter+i));  
                        gmul += *(srcData+y*src.step+(x+i)*src.channels() + 1)*(*(kernel+kcenter+i));  
                        rmul += *(srcData+y*src.step+(x+i)*src.channels() + 2)*(*(kernel+kcenter+i));  
                    }  
                    sum += (*(kernel+kcenter+i));  
                }  
            }  
            if(src.channels() == 1)  
            {  
                *(tempData+y*temp.step+x) = mul/sum;  
            }  
            else  
            {  
                *(tempData+y*temp.step+x*temp.channels()+0) = bmul/sum;  
                *(tempData+y*temp.step+x*temp.channels()+1) = gmul/sum;  
                *(tempData+y*temp.step+x*temp.channels()+2) = rmul/sum;  
            }  
        }  
    }  
  
      
    //y方向一維高斯模糊  
    for(int x = 0; x < temp.cols; x++)  
    {  
        for(int y = 0; y < temp.rows; y++)  
        {  
            double mul = 0;  
            sum = 0;  
            double bmul = 0, gmul = 0, rmul = 0;  
            for(i = -kcenter; i <= kcenter; i++)  
            {  
                if((y+i) >= 0 && (y+i) < temp.rows)  
                {  
                    if(temp.channels() == 1)  
                    {  
                        mul += *(tempData+(y+i)*temp.step+x)*(*(kernel+kcenter+i));  
                    }  
                    else  
                    {  
                        bmul += *(tempData+(y+i)*temp.step+x*temp.channels() + 0)*(*(kernel+kcenter+i));  
                        gmul += *(tempData+(y+i)*temp.step+x*temp.channels() + 1)*(*(kernel+kcenter+i));  
                        rmul += *(tempData+(y+i)*temp.step+x*temp.channels() + 2)*(*(kernel+kcenter+i));  
                    }  
                    sum += (*(kernel+kcenter+i));  
                }  
            }  
            if(temp.channels() == 1)  
            {  
                *(dstData+y*dst.step+x) = mul/sum;  
            }  
            else  
            {  
                *(dstData+y*dst.step+x*dst.channels()+0) = bmul/sum;  
                *(dstData+y*dst.step+x*dst.channels()+1) = gmul/sum;  
                *(dstData+y*dst.step+x*dst.channels()+2) = rmul/sum;  
            }  
          
        }  
    }  
      
    delete[] kernel;  
}  

該函數中使用的水平方向和垂直方向的高斯矩陣為同一矩陣，實際計算時可根據需要取不同。

模糊效果如圖3所示：

比較

使用GetTickCount()進行比較，GetTickCount()函數的精度為1ms。

以下表格中的數據均為作者機器上的某兩次運行結果取均值，編程環境為vs2010+opencv2.2

上表中Debug版本的GaussianTemplateSmooth竟然比GaussianSmooth2D運行時間長，難道是二維數組比不上一維指針，或者是Debug版本的問題？實驗結果確實如上。

將本文所寫函數與opencv2.2提供的高斯模糊函數GaussianBlur一起進行比較。

結論

如上表4,5所示，對GaussianSmooth2D的改進函數GaussianSmooth,當越大時，提速效果越明顯，這種速度的改進在Debug模式下尤為明顯。無論是在Debug,還是在Release模式下，Opencv2.2提供的GaussianBlur完勝本文所用的函數。建議在學習算法時參考本文內容，實際項目中使用GaussianBlur。

 

本例代碼鍵連接：http://download.csdn.net/detail/zddmail/4217704