码迷，mamicode.com
首页 > 其他好文 > 详细

matlab实现gabor滤波器的几种方式

时间：2014-11-05 12:18:08      阅读：283      评论：0      收藏：0      [点我收藏+]

标签：blog   http   io   ar   os   for   sp   strong   div   

转自：http://blog.csdn.net/watkinsong/article/details/7882443

方式一：

 

[csharp] view plaincopy
 
  1. function result = gaborKernel2d( lambda, theta, phi, gamma, bandwidth)  
  2. % GABORKERNEL2D   
  3. % Version: 2012/8/17 by watkins.song  
  4. % Version: 1.0  
  5. %   Fills a (2N+1)*(2N+1) matrix with the values of a 2D Gabor function.   
  6. %   N is computed from SIGMA.  
  7. %  
  8. %   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]  
  9. %   SIGMA - standard deviation of the Gaussian factor [in pixels]  
  10. %   THETA - preferred orientation [in radians]  
  11. %   PHI   - phase offset [in radians] of the cosine factor  
  12. %   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)  
  13. %   BANDWIDTH - spatial frequency bandwidth at half response,  
  14. %       *******************************************************************  
  15. %        
  16. %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH,   
  17. %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.  
  18. %       The actual value of the parameter whose input value is 0 is computed inside the   
  19. %       function from the input vallues of BANDWIDTH and the other parameter.  
  20. %   
  21. %                           pi               -1    x‘^2+gamma^2*y‘^2               
  22. %   G(x,y,theta,f) =  --------------- *exp ([----{-------------------}])*cos(2*pi*f*x‘+phi);  
  23. %                      2*sigma*sigma          2         sigma^2  
  24. %  
  25. %%% x‘ = x*cos(theta)+y*sin(theta);  
  26. %%% y‘ = y*cos(theta)-x*sin(theta);  
  27. %  
  28. % Author: watkins.song  
  29. % Email: watkins.song@gmail.com  
  30.   
  31. % calculation of the ratio sigma/lambda from BANDWIDTH   
  32. % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396  
  33. % note that in Matlab log means ln    
  34. slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );  
  35.   
  36. % calcuate sigma  
  37. sigma = slratio * lambda;  
  38.   
  39. % compute the size of the 2n+1 x 2n+1 matrix to be filled with the values of a Gabor function  
  40. % this size depends on sigma and gamma  
  41. if (gamma <= 1 && gamma > 0)  
  42.     n = ceil(2.5*sigma/gamma);  
  43. else  
  44.     n = ceil(2.5*sigma);  
  45. end  
  46.   
  47. % creation of two (2n+1) x (2n+1) matrices x and y that contain the x- and y-coordinates of  
  48. % a square 2D-mesh; the rows of x and the columns of y are copies of the vector -n:n  
  49. [x,y] = meshgrid(-n:n);  
  50.   
  51. % change direction of y-axis (In Matlab the vertical axis corresponds to the row index  
  52. % of a matrix. If the y-coordinates run from -n to n, the lowest value (-n) comes  
  53. % in the top row of the matrix ycoords and the highest value (n) in the  
  54. % lowest row. This is oposite to the customary rendering of values on the y-axis: lowest value   
  55. % in the bottom, highest on the top. Therefore the y-axis is inverted:  
  56. y = -y;  
  57.   
  58. % rotate x and y  
  59. % xp and yp are the coordinates of a point in a coordinate system rotated by theta.  
  60. % They are the main axes of the elipse of the Gaussian factor of the Gabor function.  
  61. % The wave vector of the Gabor function is along the xp axis.  
  62. xp =  x * cos(theta) + y * sin(theta);  
  63. yp = -x * sin(theta) + y * cos(theta);  
  64.   
  65. % precompute coefficients gamma2=gamma*gamma, b=1/(2*sigma*sigma) and spacial frequency  
  66. % f = 2*pi/lambda to prevent multiple evaluations   
  67. gamma2 = gamma*gamma;  
  68. b = 1 / (2*sigma*sigma);  
  69. a = b / pi;  
  70. f = 2*pi/lambda;  
  71.   
  72. % filling (2n+1) x (2n+1) matrix result with the values of a 2D Gabor function  
  73. result = a*exp(-b*(xp.*xp + gamma2*(yp.*yp))) .* cos(f*xp + phi);  
  74.   
  75. %%%%%%%%  NORMALIZATION  %%%%%%%%%%%%%%%%%%%%  
  76. % NORMALIZATION of positive and negative values to ensure that the integral of the kernel is 0.  
  77. % This is needed when phi is different from pi/2.  
  78. ppos = find(result > 0); %pointer list to indices of elements of result which are positive  
  79. pneg = find(result < 0); %pointer list to indices of elements of result which are negative   
  80.   
  81. pos =     sum(result(ppos));  % sum of the positive elements of result  
  82. neg = abs(sum(result(pneg))); % abs value of sum of the negative elements of result  
  83. meansum = (pos+neg)/2;  
  84. if (meansum > 0)   
  85.     pos = pos / meansum; % normalization coefficient for negative values of result  
  86.     neg = neg / meansum; % normalization coefficient for psoitive values of result  
  87. end  
  88.   
  89. result(pneg) = pos*result(pneg);  
  90. result(ppos) = neg*result(ppos);  
  91.   
  92. end  


方式二：

 

[csharp] view plaincopy
 
  1. function [Efilter, Ofilter, gb] = gaborKernel2d_evenodd( lambda, theta, kx, ky)  
  2. %GABORKERNEL2D_EVENODD Summary of this function goes here  
  3.  % Usage:  
  4.  %  gb =  spatialgabor(im, wavelength, angle, kx, ky, showfilter)  
  5.  % Version: 2012/8/17 by watkins.song  
  6.  % Version: 1.0  
  7.  %  
  8.  % Arguments:  
  9.  %         im         - Image to be processed.  
  10.  %         wavelength - Wavelength in pixels of Gabor filter to construct  
  11.  %         angle      - Angle of filter in degrees.  An angle of 0 gives a  
  12.  %                      filter that responds to vertical features.  
  13.  %         kx, ky     - Scale factors specifying the filter sigma relative  
  14.  %                      to the wavelength of the filter.  This is done so  
  15.  %                      that the shapes of the filters are invariant to the  
  16.  %                      scale.  kx controls the sigma in the x direction  
  17.  %                      which is along the filter, and hence controls the  
  18.  %                      bandwidth of the filter.  ky controls the sigma  
  19.  %                      across the filter and hence controls the  
  20.  %                      orientational selectivity of the filter. A value of  
  21.  %                      0.5 for both kx and ky is a good starting point.  
  22.  % %    lambda = 3;  
  23.     %   theta = 90;  
  24.     %   kx = 0.5;  
  25.     %   ky = 0.5;  
  26.  %   
  27.  %  
  28.  % Author: watkins.song  
  29.  % Email: watkins.song@gmail.com  
  30.   
  31.  % Construct even and odd Gabor filters  
  32. sigmax = lambda*kx;  
  33. sigmay = lambda*ky;  
  34.        
  35. sze = round(3*max(sigmax,sigmay));  
  36. [x,y] = meshgrid(-sze:sze);  
  37.   
  38. evenFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*cos(2*pi*(1/lambda)*x);  
  39.        
  40. % the imaginary part of the gabor filter  
  41. oddFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*sin(2*pi*(1/lambda)*x);      
  42.    
  43. evenFilter = imrotate(evenFilter, theta, ‘bilinear‘,‘crop‘);  
  44. oddFilter = imrotate(oddFilter, theta, ‘bilinear‘,‘crop‘);    
  45.        
  46. gb = evenFilter;  
  47. Efilter = evenFilter;  
  48. Ofilter = oddFilter;  
  49.   
  50. end  


 

方式三：

 

[csharp] view plaincopy
 
  1. function gb = gaborKernel2d_gaborfilter( lambda, theta, phi, gamma, bw)  
  2. %GABORKERNEL2D_GABORFILTER Summary of this function goes here  
  3. % Version: 2012/8/17 by watkins.song  
  4. % Version: 1.0  
  5. %  
  6. %   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]  
  7. %   SIGMA - standard deviation of the Gaussian factor [in pixels]  
  8. %   THETA - preferred orientation [in radians]  
  9. %   PHI   - phase offset [in radians] of the cosine factor  
  10. %   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)  
  11. %   BANDWIDTH - spatial frequency bandwidth at half response,  
  12. %       *******************************************************************  
  13. %        
  14. %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH,   
  15. %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.  
  16. %       The actual value of the parameter whose input value is 0 is computed inside the   
  17. %       function from the input vallues of BANDWIDTH and the other  
  18. %       parameter.  
  19. %                            -1     x‘^2 + y‘^2               
  20. %%% G(x,y,theta,f) =  exp ([----{-----------------})*cos(2*pi*f*x‘+phi);  
  21. %                             2     sigma*sigma  
  22. %%% x‘ = x*cos(theta)+y*sin(theta);  
  23. %%% y‘ = y*cos(theta)-x*sin(theta);  
  24. %  
  25. % Author: watkins.song  
  26. % Email: watkins.song@gmail.com  
  27.   
  28. % bw    = bandwidth, (1)  
  29. % gamma = aspect ratio, (0.5)  
  30. % psi   = phase shift, (0)  
  31. % lambda= wave length, (>=2)  
  32. % theta = angle in rad, [0 pi)  
  33.    
  34. sigma = lambda/pi*sqrt(log(2)/2)*(2^bw+1)/(2^bw-1);  
  35. sigma_x = sigma;  
  36. sigma_y = sigma/gamma;  
  37.   
  38. sz=fix(8*max(sigma_y,sigma_x));  
  39. if mod(sz,2)==0  
  40.     sz=sz+1;  
  41. end  
  42.   
  43. % alternatively, use a fixed size  
  44. % sz = 60;  
  45.    
  46. [x y]=meshgrid(-fix(sz/2):fix(sz/2),fix(sz/2):-1:fix(-sz/2));  
  47. % x (right +)  
  48. % y (up +)  
  49.   
  50. % Rotation   
  51. x_theta = x*cos(theta)+y*sin(theta);  
  52. y_theta = -x*sin(theta)+y*cos(theta);  
  53.    
  54. gb=exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);  
  55.   
  56. end  


 

方式四：

 

[csharp] view plaincopy
 
  1. function gb = gaborKernel2d_wiki( lambda, theta, phi, gamma, bandwidth)  
  2. % GABORKERNEL2D_WIKI 改写的来自wiki的gabor函数  
  3. % Version: 2012/8/17 by watkins.song  
  4. % Version: 1.0  
  5. %  
  6. %   LAMBDA - preferred wavelength (period of the cosine factor) [in pixels]  
  7. %   SIGMA - standard deviation of the Gaussian factor [in pixels]  
  8. %   THETA - preferred orientation [in radians]  
  9. %   PHI   - phase offset [in radians] of the cosine factor  
  10. %   GAMMA - spatial aspect ratio (of the x- and y-axis of the Gaussian elipse)  
  11. %   BANDWIDTH - spatial frequency bandwidth at half response,  
  12. %       *******************************************************************  
  13. %        
  14. %       BANDWIDTH, SIGMA and LAMBDA are interdependent. To use BANDWIDTH,   
  15. %       the input value of one of SIGMA or LAMBDA must be 0. Otherwise BANDWIDTH is ignored.  
  16. %       The actual value of the parameter whose input value is 0 is computed inside the   
  17. %       function from the input vallues of BANDWIDTH and the other  
  18. %       parameter.  
  19. %                            -1     x‘^2 + y‘^2               
  20. %%% G(x,y,theta,f) =  exp ([----{-----------------})*cos(2*pi*f*x‘+phi);  
  21. %                             2     sigma*sigma  
  22. %%% x‘ = x*cos(theta)+y*sin(theta);  
  23. %%% y‘ = y*cos(theta)-x*sin(theta);  
  24. %  
  25. % Author: watkins.song  
  26. % Email: watkins.song@gmail.com  
  27.   
  28.   
  29. % calculation of the ratio sigma/lambda from BANDWIDTH   
  30. % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396  
  31. % note that in Matlab log means ln    
  32. slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );  
  33.   
  34. % calcuate sigma  
  35. sigma = slratio * lambda;  
  36.   
  37. sigma_x = sigma;  
  38. sigma_y = sigma/gamma;  
  39.   
  40. % Bounding box  
  41. nstds = 4;  
  42. xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));  
  43. xmax = ceil(max(1,xmax));  
  44. ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));  
  45. ymax = ceil(max(1,ymax));  
  46. xmin = -xmax; ymin = -ymax;  
  47. [x,y] = meshgrid(xmin:xmax,ymin:ymax);  
  48.   
  49. % Rotation   
  50. x_theta = x*cos(theta) + y*sin(theta);  
  51. y_theta = -x*sin(theta) + y*cos(theta);  
  52.   
  53. % Gabor Function  
  54. gb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);  
  55.   
  56. end  


 

方式五：

 

[csharp] view plaincopy
 
  1. function [GaborReal, GaborImg] = gaborKernel_matlab( GaborH, GaborW, U, V, sigma)  
  2. %GABORKERNEL_MATLAB generate very beautiful gabor filter  
  3. % Version: 2012/8/17 by watkins.song  
  4. % Version: 1.0  
  5. % 用以生成 Gabor   
  6. % GaborReal: 核实部 GaborImg: 虚部  
  7. % GaborH,GaborW: Gabor窗口 高宽.  
  8. % U,V: 方向 大小  
  9. %            ||Ku,v||^2  
  10. % G(Z) = ---------------- exp(-||Ku,v||^2 * Z^2)/(2*sigma*sigma)(exp(i*Ku,v*Z)-exp(-sigma*sigma/2))  
  11. %          sigma*sigma  
  12. %  
  13. % 利用另外一个gabor函数来生成gabor filter, 通过u,v表示方向和尺度.  
  14. % 这里的滤波器模板的大小是不变的，变化的只有滤波器的波长和方向  
  15. % v: 代表波长  
  16. % u: 代表方向  
  17. % 缺省输入前2个参数,后面参数 Kmax=2.5*pi/2, f=sqrt(2), sigma=1.5*pi;  
  18. % GaborH, GaborW, Gabor模板大小  
  19. % U,方向因子{0,1,2,3,4,5,6,7}  
  20. % V,大小因子{0,1,2,3,4}  
  21. % Author: watkins.song  
  22. % Email: watkins.song@gmail.com  
  23.   
  24. HarfH = fix(GaborH/2);  
  25. HarfW = fix(GaborW/2);  
  26.   
  27. Qu = pi*U/8;  
  28. sqsigma = sigma*sigma;  
  29. Kv = 2.5*pi*(2^(-(V+2)/2));  
  30. %Kv = Kmax/(f^V);  
  31. postmean = exp(-sqsigma/2);  
  32.   
  33. for j = -HarfH : HarfH  
  34.     for i =  -HarfW : HarfW  
  35.         
  36.         tmp1 = exp(-(Kv*Kv*(j*j+i*i)/(2*sqsigma)));  
  37.         tmp2 = cos(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - postmean;  
  38.         %tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - exp(-sqsigma/2);  
  39.         tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j);  
  40.          
  41.         GaborReal(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp2/sqsigma;  
  42.         GaborImg(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp3/sqsigma;  
  43.     end  
  44. end  
  45.   
  46. end  


 

最后调用方式都一样：

 

[csharp] view plaincopy
 
    1. % 测试用程序  
    2. theta = [0 pi/8 2*pi/8 3*pi/8 4*pi/8 5*pi/8 6*pi/8 7*pi/8];  
    3. lambda = [4 6 8 10 12];  
    4. phi = 0;  
    5. gamma = 1;  
    6. bw = 0.5;  
    7.   
    8. % 计算每个滤波器  
    9. figure;  
    10. for i = 1:5  
    11.     for j = 1:8  
    12.         gaborFilter=gaborKernel2d(lambda(i), theta(j), phi, gamma, bw);  
    13.         % 查看每一个滤波器  
    14.         %figure;  
    15.         %imshow(real(gaborFilter),[]);  
    16.         % 将所有的滤波器放到一张图像中查看，查看滤波器组  
    17.         subplot(5,8,(i-1)*8+j);  
    18.         imshow(real(gaborFilter),[]);  
    19.     end  
    20. end  

matlab实现gabor滤波器的几种方式

标签：blog   http   io   ar   os   for   sp   strong   div   

原文地址：http://www.cnblogs.com/qxzy/p/4075665.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！