IIR递归高斯滤波

高斯滤波在项目里很常用，尤其是SIFT特征点提取的时候，PCA也要用。

但是原始的高斯滤波是一个二维的卷积，速度很慢。即使采用优化后的分离高斯滤波（先在x方向滤波，然后在y方向滤波），依然不快。

查阅了很多国内外的文献，在项目中实现了递归高斯滤波， 已经量产运行，效果很不错，运行时间是分离高斯滤波的三分之一到四分之一， 也可以运行到MCU里了！！

递归滤波器能近似模拟高斯滤波器，也是分成两次一维滤波，但是不需要设定滤波窗口大小，复杂度跟窗口大小无关，对于一维滤波先进行一次正向滤波，然后进行一次逆向滤波，每次滤波的结果迭代更新。

公式如下：

Forward：

　　Backward:

代码实现：

1. typedef struct  
2. {  
3.   gint     scale;  
4.   gint     nscales;  
5.   gint     scales_mode;  
6.   gfloat   cvar;  
7. } RetinexParams;  
8.   
9.   
10. /*  
11.  * Calculate the coefficients for the recursive filter algorithm  
12.  * Fast Computation of gaussian blurring.  
13.  */  
14. static void  
15. compute_coefs3 (gauss3_coefs *c, gfloat sigma)  
16. {  
17.   /*  
18.    * Papers:  "Recursive Implementation of the gaussian filter.",  
19.    *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.  
20.    * formula: 11b       computation of q  
21.    *          8c        computation of b0..b1  
22.    *          10        alpha is normalization constant B  
23.    */  
24.   gfloat q, q2, q3;  
25.   
26.   q = 0;  
27.   
28.   if (sigma >= 2.5)  
29.     {  
30.       q = 0.98711 * sigma - 0.96330;  
31.     }  
32.   else if ((sigma >= 0.5) && (sigma < 2.5))  
33.     {  
34.       q = 3.97156 - 4.14554 * (gfloat) sqrt ((double) 1 - 0.26891 * sigma);  
35.     }  
36.   else  
37.     {  
38.       q = 0.1147705018520355224609375;  
39.     }  
40.   
41.   q2 = q * q;  
42.   q3 = q * q2;  
43.   c->b[0] = (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));  
44.   c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3));  
45.   c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3)));  
46.   c->b[3] = (                                 (0.422205*q3));  
47.   c->B = 1.0-((c->b[1]+c->b[2]+c->b[3])/c->b[0]);  
48.   c->sigma = sigma;  
49.   c->N = 3;  
50.   
51. /*  
52.   g_printerr ("q %f\n", q);  
53.   g_printerr ("q2 %f\n", q2);  
54.   g_printerr ("q3 %f\n", q3);  
55.   g_printerr ("c->b[0] %f\n", c->b[0]);  
56.   g_printerr ("c->b[1] %f\n", c->b[1]);  
57.   g_printerr ("c->b[2] %f\n", c->b[2]);  
58.   g_printerr ("c->b[3] %f\n", c->b[3]);  
59.   g_printerr ("c->B %f\n", c->B);  
60.   g_printerr ("c->sigma %f\n", c->sigma);  
61.   g_printerr ("c->N %d\n", c->N);  
62. */  
63. }  
64.   
65. static void  
66. gausssmooth (gfloat *in, gfloat *out, gint size, gint rowstride, gauss3_coefs *c)  
67. {  
68.   /*  
69.    * Papers:  "Recursive Implementation of the gaussian filter.",  
70.    *          Ian T. Young , Lucas J. Van Vliet, Signal Processing 44, Elsevier 1995.  
71.    * formula: 9a        forward filter  
72.    *          9b        backward filter  
73.    *          fig7      algorithm  
74.    */  
75.   gint i,n, bufsize;  
76.   gfloat *w1,*w2;  
77.   
78.   /* forward pass */  
79.   bufsize = size+3;  
80.   size -= 1;  
81.   w1 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));  
82.   w2 = (gfloat *) g_try_malloc (bufsize * sizeof (gfloat));  
83.   w1[0] = in[0];  
84.   w1[1] = in[0];  
85.   w1[2] = in[0];  
86.   for ( i = 0 , n=3; i <= size ; i++, n++)  
87.     {  
88.       w1[n] = (gfloat)(c->B*in[i*rowstride] +  
89.                        ((c->b[1]*w1[n-1] +  
90.                          c->b[2]*w1[n-2] +  
91.                          c->b[3]*w1[n-3] ) / c->b[0]));  
92.     }  
93.   
94.   /* backward pass */  
95.   w2[size+1]= w1[size+3];  
96.   w2[size+2]= w1[size+3];  
97.   w2[size+3]= w1[size+3];  
98.   for (i = size, n = i; i >= 0; i--, n--)  
99.     {  
100.       w2[n]= out[i * rowstride] = (gfloat)(c->B*w1[n] +  
101.                                            ((c->b[1]*w2[n+1] +  
102.                                              c->b[2]*w2[n+2] +  
103.                                              c->b[3]*w2[n+3] ) / c->b[0]));  
104.     }  
105.   
106.   g_free (w1);  
107.   g_free (w2);  
108. }  

 

注意边界的处理，这对结果很重要。