从GIMP的Retinex算法里发现了一种高斯模糊的快速实现方法【开发记录】。

[laviewpbt]

这段时间在研究Retinex 技术（Retinex技术的难点其实还是个高斯模糊） ，看例程代码时翻到了GIMP的源代码，结果却找到了一种简单而又快速的高斯模糊的实现方式。

这种高斯模糊的实现同GIMP内嵌的高斯模糊算法也有所不同，并且速度上还有一定的优势，具体的代码可以参考GIMP下的contrast-retinex.c里面的代码。

GIMP自带的高斯模糊的代码在 blur-gauss里。

不过使用Retinex里的高斯模糊的代码会有一个小问题，就是多次模糊会发现图像像素整体向右下角或某个方向偏移，这个问题的解决很简单，有朋友遇到的时候在来问问，这里先卖个关子。

具体的算法论文可以再百度搜索 ： Recursive Implementation of the gaussian filter.

贴一段核心代码：

/*
 * 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;

}

View Code 
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);
}

有兴趣的朋友可以去看这两个源代码，C的算法代码转换为C#基本上很容易。g_try_malloc 之类的函数可以用Marshal.AllocHGlobal实现。

这段代码的优化还有很大的空间， 并且高斯模糊算法很容易并行化。

我用单线程处理 在i3的CPU上处理一副4000*3000的数码彩色照片的耗时基本在725ms（任意模糊半径理论耗时都是一样的，耗时仅为算法本身的，不包括图像数据读取和显示），说明这个算法很实用。

[rayxue001]

我知道matlab里面有自带的，但是我现在做的是硬件实现，高斯函数里面带指数肯定不行的，所以就找到了这种逼近的方法，想在MATLAB里面仿真一下。

[laviewpbt]

哦，有点出入是和谁比啊？

[isKanghua]

。。。。。那种像素偏移的问题我遇到了 楼主怎么解决的啊

[laviewpbt]

好像是在起点的位置处补上三个边缘像素，然后再处理的。

[rayxue001]

以下是我的matlab代码，其中delta即C代码中的sigma，N的大小为C代码中的width或height，这里我模拟了一个冲激函数in作为输入。这样理想的结果应该得到一个中心为N/2的高斯函数，但是，绘出来的图形有偏差，楼主能否帮我看一下问题出在哪呢？

delta=54; N=256; in=zeros(1,N); in(N/2)=1; if delta>=2.5 q=0.98711*delta-0.96330; else if (delta>=0.5)&&(delta<=2.5) q=3.97156-4.14554*sqrt(1-0.26891*delta); else q = 0.1147705018520355224609375; end end b0=1.57825+2.44413*q+1.4281*q*q+0.422205*q*q*q; b1=2.44413*q+2.85619*q*q+1.26661*q*q*q; b2=-1.4281*q*q-1.26661*q*q*q; b3=0.422205*q*q*q; B=1-(b1+b2+b3)/b0; w1=zeros(1,N+3); w1(1)=in(1); w1(2)=in(1); w1(3)=in(1); for i=4:N+3 w1(i)=B*in(i-3)+(b1*w1(i-1)+b2*w1(i-2)+b3*w1(i-3))/b0; end figure,plot(w1),title('w1'); w2=zeros(1,N+3); w2(N+1)=w1(N+3); w2(N+2)=w1(N+3); w2(N+3)=w1(N+3); for i=N:-1:1 w2(i)=B*w1(i)+(b1*w2(i+1)+b2*w2(i+2)+b3*w2(i+3))/b0; end out=w2(1:N); figure,plot(out),title('out');

[isKanghua]

为什么是三个？对任意设置的sigma的值都适用的吗？

[向立天]

过来支持一下 朋友这些年来一直钻研这个领域啊 精神可嘉

[isKanghua]

楼主能具体说下修改哪个地方吗？小弟初接触opencv，还望楼主指点下

[卧_槽]

楼主，能改成C#不？C++改起来很累啊！！

[laviewpbt]

一些中间变量展开，每行和每列都可以并行处理等等，都是一些语法上结合各自使用的语言来优化了。 都不是算法层面的处理的。

[isKanghua]

还是我，问下对于彩图是RGB通道三个通道分别处理，我试过转到YUV只处理Y这么一个通道但效果不行，楼主有什么办法对于彩图只处理一个通道达到同样效果吗？

[laviewpbt]

NO WAY.

[isKanghua]

还是谢谢楼主能这么快的给我回复，楼主好人一生平安

[isKanghua]

楼主你好，看你博文提到了说这段代码有很大的优化空间，不知道还可不可以在速度上作一些优化，希望楼主能给个方向

[wxfred]

这个算法能够进行大半径的滤波吗？

[laviewpbt]

引用 61 楼 wxfred 的回复:

这个算法能够进行大半径的滤波吗？

任意打半径都可以。

[laviewpbt]

// PS半径为r的高斯模糊效果约等于其3次半径同为r的方框模糊 引起这个的原因是因为高斯模糊是的半径实际上是高斯函数的方差值δ, 不是真正物理意义的半径。 而高斯函数的值大约在[-3δ,3δ]区间之外衰减为0。 这个回溯算法本身是O(1)算法，就是算法本身没啥好优化的，所能做的优化主要从编码层次上了，这个层次熟手和菜鸟也能达到1个数量级的区别。

[wxfred]

感谢楼主分享的代码。 不过我这台E6700的机子处理40MP的彩色图像（不包括图像的读取和存储），居然要耗18s！简直达不到要求啊。 PS的高斯滤波才2-3s，方框滤波0.6s左右。有人说PS的高斯是3次方框的结果，根据耗时来说还真有可能。不过方框模糊的半径只能调整数，但PS的高斯却可以用小数（比如半径可以用0.1），这太奇怪了。 楼主有木有研究过用3次方框模糊（box blur）来模拟高斯模糊，根据我目前的研究， 当半径r较大时，PS半径为r的高斯模糊效果约等于其3次半径同为r的方框模糊，但还是有细微的差异。当r较小时，这种差异就很大。

[wxfred]

嗯 我把compute_coefs3里w1,w2内存重复申请释放放到主函数里了，一来只申请 w1 = new float[w+3]; w2 = new float[w+3]; h1 = new float[h+3]; h2 = new float[h+3]; 一次。 40MP三通道图耗时减到16s了，还是太慢啊。 你的i3是3代以上？

[ahpong]

如何解决模糊后图像右下方偏移？