IIR递归高斯滤波总结

本文主要是介绍IIR递归高斯滤波总结，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

本文主要总结递归高斯算法的基本原理，分析，实现，指令优化实现等

原始高斯滤波算法计算量很大，网上资料很多，不赘述了，自行search

二维高斯滤波在图像处理中非常非常重要，最常见的是高斯模糊，做模糊滤镜，其本质是空域的低通滤波器，在本人另外的两篇博客同态滤波和MSRCR中有提及。对于同态滤波，需要对对数变换后的结果进行高通滤波，常规是变换到频域进行处理，也可以不进行频域变换，直接进行空域高斯滤波实现低通滤波，然后再用减法实现高通滤波。对于MSRCR算法，每个颜色通道需要进行三个强度的高斯滤波，复杂度可想而知。而且常规的高斯滤波需要考虑滤波窗大小，可以自己参考opencv里面的高斯滤波，可以自己搜一下复杂度分析，对于视频等实时处理，原始算法显然没什么实用性。

后面有人提出先进行水平一维滤波，然后再进行竖直方向一维滤波，对于一个点的滤波，计算数据量由NxN减少为2*N，但还是和窗口大小N有关

有看到一种递归减半滤波平面大小的递归实现，原理不太清楚，可参考如下资料，文中fastFilter函数的定义：

http://blog.csdn.net/smallstones/article/details/41787079

然后是本文主要介绍的IIR递归高斯滤波，这个滤波器能近似接近高斯滤波器，也是分成两次一维滤波，关键是不需要设定滤波窗口大小，复杂度跟窗口大小无关，对于一维滤波先进行一次正向滤波，然后进行一次逆向滤波，所以对于大小为NXN大小的二维矩阵，计算复杂度为O(NXN)

计算公式参考：

http://www.cnblogs.com/ImageVision/archive/2012/06/11/2545555.html

实现代码：

GIMP中有代码：http://gimp.sourcearchive.com/documentation/2.6.1/contrast-retinex_8c-source.html

[html]  view plain copy  
 typedef struct  
 {  
   gint     scale;  
   gint     nscales;  
   gint     scales_mode;  
   gfloat   cvar;  
 } RetinexParams;  
   
   
 /*  
  * 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;  
   
 /*  
   g_printerr ("q %f\n", q);  
   g_printerr ("q2 %f\n", q2);  
   g_printerr ("q3 %f\n", q3);  
   g_printerr ("c->b[0] %f\n", c->b[0]);  
   g_printerr ("c->b[1] %f\n", c->b[1]);  
   g_printerr ("c->b[2] %f\n", c->b[2]);  
   g_printerr ("c->b[3] %f\n", c->b[3]);  
   g_printerr ("c->B %f\n", c->B);  
   g_printerr ("c->sigma %f\n", c->sigma);  
   g_printerr ("c->N %d\n", c->N);  
 */  
 }  
   
 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);  
 }  

[html]  view plain copy  
   

[cpp]  view plain copy  
 /* 
  @function: 1channel 2D gauss filter 
  */  
 void gausssmooth1ch(float *src, float *dst, int width, int height, gauss3_coefs *c)  
 {  
     for (int y=0; y<height; y++) {  
         gausssmooth(src+y*width, dst+y*width, width, 1, c);  
     }  
       
     for (int x=0; x<width; x++) {  
         gausssmooth(dst+x, dst+x, height, width, c);  
     }  
       
 }  