matlab中fspecial中生成高斯模板的解释以及c语言实现

下面这个函数是模拟fspecial中的生成高斯模板的函数。

1. function h=makefilter(size,sigma)  
2. %size为模板大小  
3. %sigma为标准差  
4.   
5.   
6.   
7. %下面的代码其实是从fspecial中摘录出来的，我做了一些更改放到自己写的函数里面便于解释  
8. %计算高斯模板的中心位置  
9. siz    = ([size size]-1)/2;  
10. sig    = sigma;  
11.   
12.   
13. %用meshgrid是为了加速，不用for循环  
14. [x y]  = meshgrid(-siz(2):siz(2),-siz(1):siz(1));  
15.   
16.   
17. %计算exp(-(x^2+y^2)/(2*sig^2))  
18. %我想你肯定有一个疑问，那就是为什么不除以2*pi*sig^2  
19. %因为不除也没有关系，因为最后归一化之后会约掉  
20. arg    = -(x.*x+y.*y)./(2*sig*sig);  
21. h      = exp(arg);  
22.   
23.   
24.   
25.   
26. %不知道它为什么要这样，忘懂得人解释一下  
27. h(h<eps*max(h(:))) = 0;  
28.   
29.   
30. %求和，用来归一化  
31. sumh   = sum(h(:));  
32.   
33.   
34. %防止求和之后出现为0的情况，然后再归一化一下使高斯，模板为小数  
35. if sumh ~=0  
36.     h=h/sumh;  
37. end  
38.   
39.   
40.   
41.   
42.   
43.   
44. end  

运行结果：

1. >> makefilter(3,0.5)  
2.   
3. ans =  
4.   
5.     0.0113    0.0838    0.0113  
6.     0.0838    0.6193    0.0838  
7.     0.0113    0.0838    0.0113  

下面给出C的代码

1. #include <conio.h>  
2. #include <stdio.h>  
3. #include <math.h>  
4.    
5.   double** makeGaussianFilter(int iSize, double sigma);  
6.   double** malloc2DArray(int iRow, int iCol);  
7.   void free2DArray(double** p,int iRow);  
8.   int _tmain(int argc, _TCHAR* argv[])  
9.   {  
10.       double** filter = makeGaussianFilter(3,0.5);  
11.       for (int i=0; i<3; i++)  
12.       {  
13.           for (int j=0; j<3; j++)  
14.           {  
15.               printf("%lf\t",filter[i][j]);  
16.           }  
17.           printf("\n");  
18.       }  
19.       _getch();  
20.       free2DArray(filter,3);  
21.       return 0;  
22.   }  
23.   
24.   double** makeGaussianFilter(int iSize, double sigma)  
25.   {  
26.       double** filter = malloc2DArray(iSize,iSize);  
27.       int center = (iSize-1)/2;  
28.       double sum = 0;  
29.       double x2 = 0;  
30.       double y2 = 0;  
31.       for (int i=0; i<iSize; i++)  
32.       {  
33.           x2 = pow(double(i-center),2);  
34.           for (int j=0; j<iSize; j++)  
35.           {  
36.               y2 = pow(double(j-center),2);  
37.               sum += filter[i][j] = exp(-(x2+y2)/(2*sigma*sigma));  
38.           }  
39.       }  
40.   
41.   
42.       if (sum!=0)  
43.       {  
44.           //归一化  
45.           for (int i=0; i<iSize; i++)  
46.           {  
47.               for (int j=0; j<iSize; j++)  
48.               {  
49.                   filter[i][j] /= sum;  
50.               }  
51.           }  
52.       }  
53.   
54.       return filter;  
55.   }  
56.   
57.   
58.     
59.   double** malloc2DArray(int iRow, int iCol)  
60.   {  
61.       double **filter = new double*[iRow];  
62.       for (int i=0; i<iRow; i++)  
63.       {  
64.           filter[i] = new double[iCol];  
65.       }  
66.       return filter;  
67.   }  
68.   void free2DArray(double** p,int iRow)  
69.   {  
70.       for (int i=0; i<iRow; i++)  
71.       {  
72.           delete []p[i];  
73.       }  
74.       delete []p;  
75.   }  

结果：

总结：有时候不一定按照公式的算法去做，可以简化计算。

http://blog.csdn.net/hdanbang/article/details/45335253