快速高斯模糊

以前写过两篇关于快速高斯滤波的文章，但是代码都没写完整。

算法来源及介绍请参考博主Imageshop的博文http://www.cnblogs.com/Imageshop/archive/2013/01/07/2849782.html。

为了避免需要的朋友在看我的文章时过于浪费时间，所以删除前两篇文章，重写此篇目的为需要的朋友提供有用的信息。

下面给出的是功能完整的 c代码，下面代码无论是格式还是速度都可以进一步优化，请自行处理

/****************************************
* src    : 原始图像数据                 *
* dst    : 模糊后图像数据               *
* width  : 图像宽                       *
* height : 图像高                       *
* stride : 图像每一行字节数                *
* chan   : 图像通道数                   *
* sigma  : 高斯参数                     *
* chan   : 图像通道数                   *
*****************************************/
void IMG_GaussBlur(unsigned char* src, unsigned char* dst, int width, int height, int stride, int chan, float sigma)
{
    int i = 0;
    int h,w;
    int row    = 0;
    int col    = 0;
    int pos    = 0;
    int channel    = 0;
    int n = 0;
    int bufsize = 0;        
    int size = 0;
    int rowstride = 0;
    int itemp0 = 0;
    int itemp1 = 0;    
    float temp = 0;
    float fTab[256] = {0};
    unsigned char* ps;
    float *pIn;

    int blurH = height+3;
    int blurW = width+3;
    for (i=0; i<256; i++)
    {
        fTab[i] = i+1;
    }

    int    channelsize = width*height+(width+height)*6;

    if (width>height)
    {
        bufsize = width+6;
    }
    else
    {
        bufsize = height+6;
    }    

    float* w1    = (float *) malloc (bufsize * sizeof (float));
    float *w2    = (float *) malloc (bufsize * sizeof (float));
    float *in    = (float *) malloc (channelsize * sizeof (float));
    float *out    = (float *) malloc (channelsize * sizeof (float));

//----------------计算高斯核---------------------------------------//
    float  q    = 0; 
    float  q2, q3;    
    double b0;
    double b1;
    double b2;
    double b3;
    double B    = 0;
    int    N    = 3;

    if (sigma >= 2.5)
    {
        q = 0.98711 * sigma - 0.96330;
    }
    else if ((sigma >= 0.5) && (sigma < 2.5))
    {
        q = 3.97156 - 4.14554 * (float) sqrt ((double) 1 - 0.26891 * sigma);
    }
    else
    {
        q = 0.1147705018520355224609375;
    }

    q2 = q * q;
    q3 = q * q2;
    b0 = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
    b1 = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
    b2 = (                   -((1.4281*q2)+(1.26661 *q3)));
    b3 = (                                 (0.422205*q3));
    B = 1.0-((b1+b2+b3)/b0);

    //加速方法 减少循环多次/b0
    b1 /= b0;
    b2 /= b0;
    b3 /= b0;

//----------------计算高斯核结束---------------------------------------//    
    // 处理图像的多个通道
    for (channel = 0; channel < chan; channel++)
    {
        // 获取一个通道的所有像素值
        pIn = in;
        for (h=0; h<height; h++)
        {
            ps = src + h*stride;
            for (w=0;w<width;w++)
            {
                /* 0-255 => 1-256 */
                *pIn++ = fTab[ps[channel]];
                
                if (w==width-1)
                {
                    *pIn++ = fTab[ps[channel]];
                    *pIn++ = fTab[ps[channel]];
                    *pIn++ = fTab[ps[channel]];
                }
                ps+=chan;
            }
        } 
        memcpy(in+(height)*(width+3), in+(height-1)*(width+3), (width+3)*sizeof(float));
        memcpy(in+(height+1)*(width+3), in+(height-1)*(width+3), (width+3)*sizeof(float));
        memcpy(in+(height+2)*(width+3), in+(height-1)*(width+3), (width+3)*sizeof(float));

        //纵向处理
        size        = blurW;    
        bufsize     = size+3;
        size        -= 1;
        for (row=0 ;row < blurH; row++)
        {
            pos      = row * blurW; 
            temp  = (in + pos)[0];  
            w1[0] = temp;
            w1[1] = temp;
            w1[2] = temp;

            for ( i = 0 , n=3; i <= size ; i++, n++)
            {
                w1[n] = (float)(B*(in + pos)[i] +    ((b1*w1[n-1] +     b2*w1[n-2] + b3*w1[n-3] )));
            }

            temp =  w1[size+3];
            w2[size+1]= temp;
            w2[size+2]= temp;
            w2[size+3]= temp;
            for (i = size, n = i; i >= 0; i--, n--)
            {
                 (out + pos)[i] = w2[n] = (float)(B*w1[n] +    ((b1*w2[n+1] +    b2*w2[n+2] + b3*w2[n+3] )));
            }    
        }    
        

        //横向处理
        size        = blurH;
        rowstride   = blurW;
        bufsize     = size+3;
        size        -= 1;
        for (col=0; col < blurW; col++)
        {
            temp  = (out + col)[0];
            w1[0] = temp;
            w1[1] = temp;
            w1[2] = temp;
            for ( i = 0 , n=3; i <= size ; i++, n++)
            {
                w1[n] = (float)(B*(out + col)[i*rowstride] + ((b1*w1[n-1] + b2*w1[n-2] + b3*w1[n-3] )));
            }
            
            temp        = w1[size+3];
            w2[size+1]    = temp;
            w2[size+2]    = temp;
            w2[size+3]    = temp;
            for (i = size, n = i; i >= 0; i--, n--)
            {
                (in + col)[i * rowstride] =w2[n]= (float)(B*w1[n] +  ((b1*w2[n+1] + b2*w2[n+2] + b3*w2[n+3] )));
            }                
        }
    
        //修正偏移的拷贝方法
        for(int y=0; y<height; y++)
        {
            ps = dst+ y*stride;
            itemp1 = (y+3)*blurW;                                // +3
            for (int x=0; x<width; x++)
            {        
                ps[channel] = in[itemp1+x+3]-1;
                ps+=chan;
            }
        }         
    }

    free (w1);
    free (w2);
    free (in);
    free (out);
}

  调用参考

IplImage* src = cvLoadImage("./test.jpg", 1);
IplImage* dst = cvLoadImage("./test.jpg", 1);

IMG_GaussBlur((unsigned char*)src->imageData, 
              (unsigned char*)dst->imageData, 
              src->width, src->height, 
              src->widthStep,
              src->nChannels, 2);