首页 > 代码库 > matlab实现gabor滤波器的几种方式

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

2024-07-31 07:31:45   229人阅读

转自：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滤波器的几种方式


联系
我们