快速高斯滤波

高斯滤波器是图像处理中经常用到的滤波器，其滤波核函数为：

为简单起见，这里省略了归一化因子。

由的可分离特性：

得：

其中为输入图像，为输出图像，为滤波模板半径。根据准则，通常使。

由上式可见，我们可以将二维高斯滤波分解为两次一维高斯滤波。

对于二维高斯滤波，设图像大小，高斯模板大小，处理每个像素点需要次操作，则算法复杂度。若使用一维高斯核先对图像逐行滤波，再对中间结果逐列滤波，则处理每个像素需要次操作，算法复杂度，随着滤波模板尺寸的增大，算法优势越明显。

程序：

#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