快速高斯濾波

高斯濾波器是圖像處理中經常用到的濾波器，其濾波核函數為：

為簡單起見，這里省略了歸一化因子。

由的可分離特性：

得：

 

其中為輸入圖像，為輸出圖像，為濾波模板半徑。根據准則，通常使。

由上式可見，我們可以將二維高斯濾波分解為兩次一維高斯濾波。

對於二維高斯濾波，設圖像大小，高斯模板大小，處理每個像素點需要次操作，則算法復雜度。若使用一維高斯核先對圖像逐行濾波，再對中間結果逐列濾波，則處理每個像素需要次操作，算法復雜度，隨着濾波模板尺寸的增大，算法優勢越明顯。

程序：

#include "stdafx.h"
#include<stdlib.h>
#include<math.h>
//邊界處理
int Edge(int i,int x,int Max)
{
    int k=i+x;
    if(k<0)k=-i;
    else if(k>=Max) k=Max-i-1;
    else k=x;
    return k;
}

//二維高斯處理灰度圖像
extern "C" _declspec(dllexport) void GaussFilterGray(unsigned char *A, int nWidth, int nHeight, int Stride,double dSigma)
{
    unsigned char *buffer=(unsigned char*)malloc(Stride*nHeight);
    memcpy(buffer,A,Stride*nHeight);
    int nWindowSize = (int)(1+2*ceil(3*dSigma)); 
    int nCenter = (nWindowSize)/2;
    double* pdKernel = new double[nWindowSize*nWindowSize]; 
    double  dSum = 0.0;
    double scale2X = 0.5/(dSigma*dSigma);
    double dFilter=0.0;
    double ImageData=0.0;
    //生成二維高斯濾波核
    for(int i=0; i<nWindowSize; i++)  
    {  
        for(int j=0; j<nWindowSize; j++)  
        {  
            int nDis_x = i-nCenter;  
            int nDis_y = j-nCenter;  
            pdKernel[j+i*nWindowSize]=exp(-(nDis_x*nDis_x+nDis_y*nDis_y)*scale2X);
            dSum += pdKernel[i*nWindowSize+j];  
        }  
    }  
    //歸一化
    for(int i=0; i<nWindowSize; i++)  
    {  
        for(int j=0; j<nWindowSize; j++)  
        {  
            pdKernel[i*nWindowSize+j] /= dSum;  
        }  
    }  
    //逐像素處理
    for(int inx=0,i=0; i<nHeight; i++)  
    {  
        for(int j=0; j<nWidth; j++,inx++)  
        {  
            dFilter=0;
            //鄰域內加權平均
            for(int n=0,x=-nCenter; x<=nCenter; x++) 
            {  
                int i_x=Edge(i,x,nHeight);
                for(int y=-nCenter; y<=nCenter; y++,n++)
                {  
                    int j_y=Edge(j,y,nWidth);
                    int index=(i+i_x)*Stride+j+j_y;//鄰域內像素在內存中的下標
                    ImageData=buffer[index];
                    dFilter+=ImageData*pdKernel[n];
                }  
            } 
            A[inx]= max(min(255,dFilter),0);
        }  
    }  
    delete[]pdKernel;  
    delete[]buffer;
}

//一維高斯處理灰度圖像
extern "C" _declspec(dllexport) void GaussFilterGray1D(unsigned char *A, int nWidth, int nHeight, int Stride,double dSigma)
{
    unsigned char *buffer=(unsigned char*)malloc(Stride*nHeight);
    memcpy(buffer,A,Stride*nHeight);
    int nWindowSize = (int)(1+2*ceil(3*dSigma)); 
    int nCenter = (nWindowSize)/2;
    double* pdKernel = new double[nWindowSize]; 
    double  dSum = 0.0;
    double scale2X = 0.5/(dSigma*dSigma);
    double dFilter=0.0;
    double ImageData=0.0;    
    //生成一維高斯核
    for(int i=0; i<nWindowSize; i++)  
    {  
        int nDis_x = i-nCenter;  
        pdKernel[i]=exp(-(nDis_x*nDis_x)*scale2X);
        dSum += pdKernel[i];  
    }  
    //歸一化
    for(int i=0; i<nWindowSize; i++)  
    {  
        pdKernel[i] /= dSum;  
    }  
    //橫向濾波
    for(int inx=0,i=0; i<nHeight; i++)  
    {  
        for(int j=0; j<nWidth; j++,inx++)  
        {  
            dFilter=0;
            for(int n=0,x=-nCenter; x<=nCenter; x++,n++) 
            {  
                int j_x=Edge(j,x,nWidth);           
                int index=inx+j_x;
                ImageData=A[index];//從原圖像A中取值
                dFilter+=ImageData*pdKernel[n]; 
            } 
            buffer[inx]= max(min(255,dFilter),0);//中間結果放在buffer中
        }  
    }

    //縱向濾波
    for(int i=0;i<nWidth;i++)
    {
        for(int j=0;j<nHeight;j++)
        {
            dFilter=0;
            for(int n=0,x=-nCenter; x<=nCenter; x++,n++) 
            {  
                int j_x=Edge(j,x,nHeight);           
                int index=(j+j_x)*Stride+i;
                ImageData=buffer[index];//從中間圖像buffer中取值
                dFilter+=ImageData*pdKernel[n]; 
            } 
            A[j*Stride+i]= max(min(255,dFilter),0);
        }
    }
    delete[]pdKernel;  
    delete[]buffer;
}


//二維高斯處理彩色圖像
extern "C" _declspec(dllexport) void GaussFilterColor(unsigned char *A, int nWidth, int nHeight, int Stride,double dSigma)
{
    int Step=3;
    if(Stride==4*nWidth)Step=4;//四通道圖像
    unsigned char *buffer=(unsigned char*)malloc(Stride*nHeight);
    memcpy(buffer,A,Stride*nHeight);
    int nWindowSize = (int)(1+2*ceil(3*dSigma)); 
    int nCenter = (nWindowSize)/2;
    double* pdKernel = new double[nWindowSize*nWindowSize]; 
    double  dSum = 0.0;    
    double scale2X = 0.5/(dSigma*dSigma);
    double dFilterB=0.0;
    double dFilterG=0.0; 
    double dFilterR=0.0;
    int index;
    double ImageData;
    //生成二維高斯濾波核
    for(int i=0; i<nWindowSize; i++)  
    {  
        for(int j=0; j<nWindowSize; j++)  
        {  
            int nDis_x = i-nCenter;  
            int nDis_y = j-nCenter;  
            pdKernel[j+i*nWindowSize]=exp(-(nDis_x*nDis_x+nDis_y*nDis_y)*scale2X);
            dSum += pdKernel[i*nWindowSize+j];  
        }  
    } 
    //歸一化
    for(int i=0; i<nWindowSize; i++)  
    {  
        for(int j=0; j<nWindowSize; j++)  
        {  
            pdKernel[i*nWindowSize+j] /= dSum;  
        }  
    }  
    for(int i=0; i<nHeight; i++)  
    {  
        for(int j=0; j<nWidth;j++)  
        {  
            dFilterB=0;
            dFilterG=0; 
            dFilterR=0; 
            for(int n=0,x=-nCenter; x<=nCenter; x++) 
            {  
                int i_x=Edge(i,x,nHeight);
                for(int y=-nCenter; y<=nCenter; y++,n++)
                {  
                    int j_y=Edge(j,y,nWidth);
                    index=(i+i_x)*Stride+(j+j_y)*Step;
                    //三通道BGR,四通道BGRA
                    ImageData=buffer[index];
                    dFilterB+=ImageData * pdKernel[n];
                    index+=1;
                    ImageData=buffer[index];
                    dFilterG+=ImageData * pdKernel[n];
                    index+=1;
                    ImageData=buffer[index];
                    dFilterR+=ImageData * pdKernel[n];
                }  
            }      
            index=i*Stride+j*Step;
            A[index]=max(min(dFilterB,255),0);
            A[index+1]=max(min(dFilterG,255),0);
            A[index+2]=max(min(dFilterR,255),0);
        }  
    }  
    delete[]pdKernel;  
    delete[]buffer;
}


//一維高斯處理彩色圖像
extern "C" _declspec(dllexport) void GaussFilterColor1D(unsigned char *A, int nWidth, int nHeight, int Stride,double dSigma)
{
    int Step=3;
    if(Stride==4*nWidth)Step=4;//四通道圖像
    unsigned char *buffer=(unsigned char*)malloc(Stride*nHeight);
    memcpy(buffer,A,Stride*nHeight);
    int nWindowSize = (int)(1+2*ceil(3*dSigma)); 
    int nCenter = (nWindowSize)/2;
    double* pdKernel = new double[nWindowSize]; 
    double scale2X = 0.5/(dSigma*dSigma);
    double  dSum = 0.0;
    double dFilterB=0;
    double dFilterG=0; 
    double dFilterR=0; 
    double ImageData;
    int index;
    //一維高斯核
    for(int i=0; i<nWindowSize; i++)  
    {  
        int nDis_x = i-nCenter;  
        pdKernel[i]=exp(-(nDis_x*nDis_x)*scale2X);
        dSum += pdKernel[i];  
    }  
    //歸一化
    for(int i=0; i<nWindowSize; i++)  
    {  
        pdKernel[i] /= dSum;  
    }  
    //橫向濾波
    for(int i=0; i<nHeight; i++)  
    {  
        for(int j=0; j<nWidth;j++)  
        {  
            dFilterB=0;
            dFilterG=0; 
            dFilterR=0; 
            for(int n=0,x=-nCenter; x<=nCenter; x++,n++) 
            {  
                int j_x=Edge(j,x,nWidth);
                index=i*Stride+(j+j_x)*Step;
                ImageData=A[index];//從原圖像A中取值
                dFilterB+=ImageData * pdKernel[n];
                index+=1;
                ImageData=A[index];
                dFilterG+=ImageData * pdKernel[n];
                index+=1;
                ImageData=A[index];
                dFilterR+=ImageData * pdKernel[n];
             } 
            index=i*Stride+j*Step;
            buffer[index]=max(min(dFilterB,255),0);//中間結果放在buffer中
            buffer[index+1]=max(min(dFilterG,255),0);
            buffer[index+2]=max(min(dFilterR,255),0);
        }      
    }  

    //縱向濾波
    for(int i=0;i<nWidth;i++)
    {
        for(int j=0;j<nHeight;j++)
        {
            dFilterB=0;
            dFilterG=0; 
            dFilterR=0; 
            for(int n=0,x=-nCenter; x<=nCenter; x++,n++) 
            {  
                int j_x=Edge(j,x,nHeight);           
                int index=(j+j_x)*Stride+i*Step;
                ImageData=buffer[index];//從中間圖像buffer中取值
                dFilterB+=ImageData * pdKernel[n];
                index+=1;
                ImageData=buffer[index];
                dFilterG+=ImageData * pdKernel[n];
                index+=1;
                ImageData=buffer[index];
                dFilterR+=ImageData * pdKernel[n];
            } 
            index=j*Stride+i*Step;
            A[index]=max(min(dFilterB,255),0);
            A[index+1]=max(min(dFilterG,255),0);
            A[index+2]=max(min(dFilterR,255),0);
        }
    }
    delete[]pdKernel;  
    delete[]buffer;
}

 

演示結果：

 

上面對一幅512*512的彩色圖像，基本的高斯算法耗時1469ms,而快速高斯耗時439ms。選取的，此時的濾波模板大小為。

完整的工程項目下載：快速高斯濾波

參考：

http://blog.csdn.net/markl22222/article/details/10313565

http://m.blog.csdn.net/blog/zxpddfg/45912561