高斯滤波

高斯滤波器 (Gaussian Filter)

• 最有用的滤波器 (尽管不是最快的)。 高斯滤波是将输入数组的每一个像素点与 高斯内核 卷积将卷积和当作输出像素值。

• 还记得1维高斯函数的样子吗?

  

  假设图像是1维的,那么观察上图，不难发现中间像素的加权系数是最大的， 周边像素的加权系数随着它们远离中间像素的距离增大而逐渐减小。

Note

2维高斯函数可以表达为 :

其中 为均值 (峰值对应位置)， 代表标准差 (变量 和 变量 各有一个均值，也各有一个标准差)。

 

高斯滤波的实现比较简单，opencv和matlab均有函数实现，参看图像处理之滤波器

这里有一个简单的例子对比：

#include <opencv2/opencv.hpp>
#include <vector>
using namespace std;
using namespace cv;
void OpenCv_Gauss_Test()
{
    Mat src = imread("lena.jpg", 0);

    int scale = 1;
    resize(src, src, Size(src.cols * scale, src.rows*scale));
    Mat dst = src.clone();
    cout << src.cols << "*" << src.rows << " " << src.rows * src.cols / 10000 << "M" << endl;
    int ksize = 11;
    GaussianBlur(src, dst, Size(ksize, ksize), 0, 0);
    namedWindow("dstopencv");
    imshow("dstopencv",dst);
}

void mygausstest()
{
    Mat src = imread("lena.jpg", 0);

    int scale = 1;
    resize(src, src, Size(src.cols * scale, src.rows*scale));
    Mat dst = src.clone();
    cout << src.cols << "*" << src.rows << " " << src.rows * src.cols / 10000 << "M" << endl;
    int ksize = 11;
    int step = ksize / 2;
    int m = src.rows - step * 2 - 1;
    int n = src.cols - step * 2 - 1;
    float sigma = 0.3*((ksize - 1)*0.5 - 1) + 0.8;
    Mat  gausskernel = getGaussianKernel(ksize, sigma, CV_64F);
    Mat showkernel = gausskernel*gausskernel.t();
    for (int i = step + 1; i < m; ++i)
    {
        for (int j = step + 1; j < n; ++j)
        {
            Mat sub =src(Range(i - step, i + step + 1), Range(j - step, j + step + 1));
            sub.convertTo(sub, CV_64F);
            dst.at<uchar>(i, j) = cv::sum(sub.mul(showkernel))[0];
        }
    }
    namedWindow("src");
    namedWindow("dst");
    imshow("src", src);
    imshow("dst", dst);
    waitKey(0);
}

void main()
{
    OpenCv_Gauss_Test();
    mygausstest();
}

这里主要讨论高斯滤波的优化加速问题。

 

参考：

http://blog.csdn.net/app_12062011/article/details/27082049

http://www.opencv.org.cn/opencvdoc/2.3.2/html/doc/tutorials/imgproc/gausian_median_blur_bilateral_filter/gausian_median_blur_bilateral_filter.html

http://homepages.inf.ed.ac.uk/rbf/HIPR2/gsmooth.htm

http://blog.csdn.net/app_12062011/article/details/27082049

IIR高斯递归滤波器

http://blog.csdn.net/bluecol/article/details/45584065 (IIR高斯滤波器)

http://q.cnblogs.com/q/54532/

http://www.cnblogs.com/Imageshop/archive/2013/01/07/2849782.html

gimp源码

/*
 * Calculate the coefficients for the recursive filter algorithm
 * Fast Computation of gaussian blurring.
 */
static void
compute_coefs3 (gauss3_coefs *c, gfloat sigma)
{
  /*
   * Papers:  "Recursive Implementation of the gaussian filter.",
   *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
   * formula: 11b       computation of q
   *          8c        computation of b0..b1
   *          10        alpha is normalization constant B
   */
  gfloat q, q2, q3;

  q = 0;

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

  q2 = q * q;
  q3 = q * q2;
  c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));
  c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3));
  c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3)));
  c->b[3] = (                                 (0.422205*q3));
  c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);
  c->sigma = sigma;
  c->N = 3;
}


static void
gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)
{
  /*
   * Papers:  "Recursive Implementation of the gaussian filter.",
   *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.
   * formula: 9a        forward filter
   *          9b        backward filter
   *          fig7      algorithm
   */
  gint i,n, bufsize;
  gfloat *w1,*w2;

  /* forward pass */
  bufsize = size+3;
  size -= 1;
  w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
  w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));
  w1[0] = in[0];
  w1[1] = in[0];
  w1[2] = in[0];
  for ( i = 0 , n=3; i <= size ; i++, n++)
    {
      w1[n] = (gfloat)(c->B*in[i*rowstride] +
                       ((c->b[1]*w1[n-1] +
                         c->b[2]*w1[n-2] +
                         c->b[3]*w1[n-3] ) / c->b[0]));
    }

  /* backward pass */
  w2[size+1]= w1[size+3];
  w2[size+2]= w1[size+3];
  w2[size+3]= w1[size+3];
  for (i = size, n = i; i >= 0; i--, n--)
    {
      w2[n]= out[i * rowstride] = (gfloat)(c->B*w1[n] +
                                           ((c->b[1]*w2[n+1] +
                                             c->b[2]*w2[n+2] +
                                             c->b[3]*w2[n+3] ) / c->b[0]));
    }

  g_free (w1);
  g_free (w2);
}

opencv下面实现IIR高斯滤波器

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

#include "Timer.h"
Timer timer;        //    timer
double elapsedTime;    //    time in millisecond
double  avgtime = 0;
double testnum = 0;
void timebegin()
{
    timer.start();
}

void timeend(string str)
{
    elapsedTime = timer.getElapsedTimeInMilliSec();
    avgtime += elapsedTime;
    str = str + "time is: ";// %lf ms
";
    cout << str << elapsedTime << "ms" << endl;
}

/****************************************
*    src        : 原始图像数据                *
*    dst        : 模糊后图像数据            *
*    width    : 图像宽                    *
*    height    : 图像高                    *
*    sigma    : 高斯参数                    *
*    chan    : 图像通道数                *
*****************************************/
void IMG_GaussBlur(unsigned char* src, unsigned char*& dst, int width, int height, float sigma, int chan)
{
    int    i = 0;
    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    channelsize = width*height;

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

    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));

    //----------------计算高斯核---------------------------------------//
    //compute_coefs3(coef, sigma);    
    float  q = 0;
    float  q2, q3;
    double b[4] = { 0 };
    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;
    b[0] = (1.57825 + (2.44413*q) + (1.4281 *q2) + (0.422205*q3));
    b[1] = ((2.44413*q) + (2.85619*q2) + (1.26661 *q3));
    b[2] = (-((1.4281*q2) + (1.26661 *q3)));
    b[3] = ((0.422205*q3));
    B = 1.0 - ((b[1] + b[2] + b[3]) / b[0]);
    //----------------计算高斯核结束---------------------------------------//

    // 处理图像的多个通道
    for (channel = 0; channel < chan; channel++)
    {
        // 获取一个通道的所有像素值
        for (i = 0, pos = channel; i < channelsize; i++, pos += chan)
        {
            /* 0-255 => 1-256 */
            in[i] = (float)(src[pos] + 1.0);
        }

        //纵向处理
        for (row = 0; row < height; row++)
        {
            pos = row * width;
            //gausssmooth(in + pos, out + pos, width, 1, &coef);
            size = width;
            rowstride = 1;
            bufsize = size + 3;
            size -= 1;

            w1[0] = (in + pos)[0];
            w1[1] = (in + pos)[0];
            w1[2] = (in + pos)[0];

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

            /** backward pass */
            w2[size + 1] = w1[size + 3];
            w2[size + 2] = w1[size + 3];
            w2[size + 3] = w1[size + 3];
            for (i = size, n = i; i >= 0; i--, n--)
            {
                w2[n] = (out + pos)[i * rowstride] = (float)(B*w1[n] +
                    ((b[1] * w2[n + 1] +
                    b[2] * w2[n + 2] +
                    b[3] * w2[n + 3]) / b[0]));
            }
        }

        //下面的横向处理直接将数据 out 与 in 对调 省去memcpy
        //memcpy(in,  out, channelsize * sizeof(float));
        //memset(out, 0  , channelsize * sizeof(float));

        //横向处理
        for (col = 0; col < width; col++)
        {
            //gausssmooth(out + col, in + col, height, width, &coef);            
            size = height;
            rowstride = width;
            bufsize = size + 3;
            size -= 1;

            w1[0] = (out + col)[0];
            w1[1] = (out + col)[0];
            w1[2] = (out + col)[0];
            for (i = 0, n = 3; i <= size; i++, n++)
            {
                w1[n] = (float)(B*(out + col)[i*rowstride] +
                    ((b[1] * w1[n - 1] +
                    b[2] * w1[n - 2] +
                    b[3] * w1[n - 3]) / b[0]));
            }

            w2[size + 1] = w1[size + 3];
            w2[size + 2] = w1[size + 3];
            w2[size + 3] = w1[size + 3];
            for (i = size, n = i; i >= 0; i--, n--)
            {
                w2[n] = (in + col)[i * rowstride] = (float)(B*w1[n] +
                    ((b[1] * w2[n + 1] +
                    b[2] * w2[n + 2] +
                    b[3] * w2[n + 3]) / b[0]));
            }
        }

        //拷贝结果到函数输出指针
        for (i = 0, pos = channel; i < channelsize; i++, pos += chan)
        {
            dst[pos] = in[i] - 1;
        }
    }

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

void main()
{
    Mat src = imread("src.jpg");
    Mat dst1 = src.clone();
    Mat dst2 = src.clone();
    float sigma = 15.0;
    timebegin();
    IMG_GaussBlur(src.data, dst1.data, src.cols, src.rows, sigma, 3);
    timeend("IIR_GAUSS");

    int ksize = 71;
    timebegin();
    GaussianBlur(src, dst2, Size(ksize,ksize),0, 0);
    timeend("OPNECV_GAUSS");

    //namedWindow("src",0);
    namedWindow("dst1",0);
    namedWindow("dst2",0);
    //imshow("src", src);
    imshow("dst1", dst1);
    imshow("dst2", dst2);
    waitKey(0);
}

附加 测试时间的函数代码，这个代码可以在linux以及windows下面通用，也可以参看博文：
Time.cpp

//////////////////////////////////////////////////////////////////////////////
// Timer.cpp
// =========
// High Resolution Timer.
// This timer is able to measure the elapsed time with 1 micro-second accuracy
// in both Windows, Linux and Unix system 
//
//  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
// CREATED: 2003-01-13
// UPDATED: 2006-01-13
//
// Copyright (c) 2003 Song Ho Ahn
//////////////////////////////////////////////////////////////////////////////

#include "Timer.h"
#include <stdlib.h>

///////////////////////////////////////////////////////////////////////////////
// constructor
///////////////////////////////////////////////////////////////////////////////
Timer::Timer()
{
#ifdef WIN32
    QueryPerformanceFrequency(&frequency);
    startCount.QuadPart = 0;
    endCount.QuadPart = 0;
#else
    startCount.tv_sec = startCount.tv_usec = 0;
    endCount.tv_sec = endCount.tv_usec = 0;
#endif

    stopped = 0;
    startTimeInMicroSec = 0;
    endTimeInMicroSec = 0;
}



///////////////////////////////////////////////////////////////////////////////
// distructor
///////////////////////////////////////////////////////////////////////////////
Timer::~Timer()
{
}



///////////////////////////////////////////////////////////////////////////////
// start timer.
// startCount will be set at this point.
///////////////////////////////////////////////////////////////////////////////
void Timer::start()
{
    stopped = 0; // reset stop flag
#ifdef WIN32
    QueryPerformanceCounter(&startCount);
#else
    gettimeofday(&startCount, NULL);
#endif
}



///////////////////////////////////////////////////////////////////////////////
// stop the timer.
// endCount will be set at this point.
///////////////////////////////////////////////////////////////////////////////
void Timer::stop()
{
    stopped = 1; // set timer stopped flag

#ifdef WIN32
    QueryPerformanceCounter(&endCount);
#else
    gettimeofday(&endCount, NULL);
#endif
}



///////////////////////////////////////////////////////////////////////////////
// compute elapsed time in micro-second resolution.
// other getElapsedTime will call this first, then convert to correspond resolution.
///////////////////////////////////////////////////////////////////////////////
double Timer::getElapsedTimeInMicroSec()
{
#ifdef WIN32
    if(!stopped)
        QueryPerformanceCounter(&endCount);

    startTimeInMicroSec = startCount.QuadPart * (1000000.0 / frequency.QuadPart);
    endTimeInMicroSec = endCount.QuadPart * (1000000.0 / frequency.QuadPart);
#else
    if(!stopped)
        gettimeofday(&endCount, NULL);

    startTimeInMicroSec = (startCount.tv_sec * 1000000.0) + startCount.tv_usec;
    endTimeInMicroSec = (endCount.tv_sec * 1000000.0) + endCount.tv_usec;
#endif

    return endTimeInMicroSec - startTimeInMicroSec;
}



///////////////////////////////////////////////////////////////////////////////
// divide elapsedTimeInMicroSec by 1000
///////////////////////////////////////////////////////////////////////////////
double Timer::getElapsedTimeInMilliSec()
{
    return this->getElapsedTimeInMicroSec() * 0.001;
}



///////////////////////////////////////////////////////////////////////////////
// divide elapsedTimeInMicroSec by 1000000
///////////////////////////////////////////////////////////////////////////////
double Timer::getElapsedTimeInSec()
{
    return this->getElapsedTimeInMicroSec() * 0.000001;
}



///////////////////////////////////////////////////////////////////////////////
// same as getElapsedTimeInSec()
///////////////////////////////////////////////////////////////////////////////
double Timer::getElapsedTime()
{
    return this->getElapsedTimeInSec();
}

//////////////////////////////////////////////////////////////////////////////
// Timer.h
// =======
// High Resolution Timer.
// This timer is able to measure the elapsed time with 1 micro-second accuracy
// in both Windows, Linux and Unix system 
//
//  AUTHOR: Song Ho Ahn (song.ahn@gmail.com)
// CREATED: 2003-01-13
// UPDATED: 2006-01-13
//
// Copyright (c) 2003 Song Ho Ahn
//////////////////////////////////////////////////////////////////////////////

#ifndef TIMER_H_DEF
#define TIMER_H_DEF

#ifdef WIN32   // Windows system specific
#include <windows.h>
#else          // Unix based system specific
#include <sys/time.h>
#endif


class Timer
{
public:
    Timer();                                    // default constructor
    ~Timer();                                   // default destructor

    void   start();                             // start timer
    void   stop();                              // stop the timer
    double getElapsedTime();                    // get elapsed time in second
    double getElapsedTimeInSec();               // get elapsed time in second (same as getElapsedTime)
    double getElapsedTimeInMilliSec();          // get elapsed time in milli-second
    double getElapsedTimeInMicroSec();          // get elapsed time in micro-second


protected:


private:
    double startTimeInMicroSec;                 // starting time in micro-second
    double endTimeInMicroSec;                   // ending time in micro-second
    int    stopped;                             // stop flag 
#ifdef WIN32
    LARGE_INTEGER frequency;                    // ticks per second
    LARGE_INTEGER startCount;                   //
    LARGE_INTEGER endCount;                     //
#else
    timeval startCount;                         //
    timeval endCount;                           //
#endif
};

#endif // TIMER_H_DEF