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

转自: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. % calculation of the ratio sigma/lambda from BANDWIDTH
  31. % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396
  32. % note that in Matlab log means ln
  33. slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );
  34. % calcuate sigma
  35. sigma = slratio * lambda;
  36. % compute the size of the 2n+1 x 2n+1 matrix to be filled with the values of a Gabor function
  37. % this size depends on sigma and gamma
  38. if (gamma <= 1 && gamma > 0)
  39. n = ceil(2.5*sigma/gamma);
  40. else
  41. n = ceil(2.5*sigma);
  42. end
  43. % creation of two (2n+1) x (2n+1) matrices x and y that contain the x- and y-coordinates of
  44. % a square 2D-mesh; the rows of x and the columns of y are copies of the vector -n:n
  45. [x,y] = meshgrid(-n:n);
  46. % change direction of y-axis (In Matlab the vertical axis corresponds to the row index
  47. % of a matrix. If the y-coordinates run from -n to n, the lowest value (-n) comes
  48. % in the top row of the matrix ycoords and the highest value (n) in the
  49. % lowest row. This is oposite to the customary rendering of values on the y-axis: lowest value
  50. % in the bottom, highest on the top. Therefore the y-axis is inverted:
  51. y = -y;
  52. % rotate x and y
  53. % xp and yp are the coordinates of a point in a coordinate system rotated by theta.
  54. % They are the main axes of the elipse of the Gaussian factor of the Gabor function.
  55. % The wave vector of the Gabor function is along the xp axis.
  56. xp = x * cos(theta) + y * sin(theta);
  57. yp = -x * sin(theta) + y * cos(theta);
  58. % precompute coefficients gamma2=gamma*gamma, b=1/(2*sigma*sigma) and spacial frequency
  59. % f = 2*pi/lambda to prevent multiple evaluations
  60. gamma2 = gamma*gamma;
  61. b = 1 / (2*sigma*sigma);
  62. a = b / pi;
  63. f = 2*pi/lambda;
  64. % filling (2n+1) x (2n+1) matrix result with the values of a 2D Gabor function
  65. result = a*exp(-b*(xp.*xp + gamma2*(yp.*yp))) .* cos(f*xp + phi);
  66. %%%%%%%% NORMALIZATION %%%%%%%%%%%%%%%%%%%%
  67. % NORMALIZATION of positive and negative values to ensure that the integral of the kernel is 0.
  68. % This is needed when phi is different from pi/2.
  69. ppos = find(result > 0); %pointer list to indices of elements of result which are positive
  70. pneg = find(result < 0); %pointer list to indices of elements of result which are negative
  71. pos = sum(result(ppos)); % sum of the positive elements of result
  72. neg = abs(sum(result(pneg))); % abs value of sum of the negative elements of result
  73. meansum = (pos+neg)/2;
  74. if (meansum > 0)
  75. pos = pos / meansum; % normalization coefficient for negative values of result
  76. neg = neg / meansum; % normalization coefficient for psoitive values of result
  77. end
  78. result(pneg) = pos*result(pneg);
  79. result(ppos) = neg*result(ppos);
  80. 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. % Construct even and odd Gabor filters
  31. sigmax = lambda*kx;
  32. sigmay = lambda*ky;
  33. sze = round(3*max(sigmax,sigmay));
  34. [x,y] = meshgrid(-sze:sze);
  35. evenFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*cos(2*pi*(1/lambda)*x);
  36. % the imaginary part of the gabor filter
  37. oddFilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2).*sin(2*pi*(1/lambda)*x);
  38. evenFilter = imrotate(evenFilter, theta, \'bilinear\',\'crop\');
  39. oddFilter = imrotate(oddFilter, theta, \'bilinear\',\'crop\');
  40. gb = evenFilter;
  41. Efilter = evenFilter;
  42. Ofilter = oddFilter;
  43. 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. % bw = bandwidth, (1)
  28. % gamma = aspect ratio, (0.5)
  29. % psi = phase shift, (0)
  30. % lambda= wave length, (>=2)
  31. % theta = angle in rad, [0 pi)
  32. sigma = lambda/pi*sqrt(log(2)/2)*(2^bw+1)/(2^bw-1);
  33. sigma_x = sigma;
  34. sigma_y = sigma/gamma;
  35. sz=fix(8*max(sigma_y,sigma_x));
  36. if mod(sz,2)==0
  37. sz=sz+1;
  38. end
  39. % alternatively, use a fixed size
  40. % sz = 60;
  41. [x y]=meshgrid(-fix(sz/2):fix(sz/2),fix(sz/2):-1:fix(-sz/2));
  42. % x (right +)
  43. % y (up +)
  44. % Rotation
  45. x_theta = x*cos(theta)+y*sin(theta);
  46. y_theta = -x*sin(theta)+y*cos(theta);
  47. gb=exp(-0.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);
  48. 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. % calculation of the ratio sigma/lambda from BANDWIDTH
  28. % according to Kruizinga and Petkov, 1999 IEEE Trans on Image Processing 8 (10) p.1396
  29. % note that in Matlab log means ln
  30. slratio = (1/pi) * sqrt( (log(2)/2) ) * ( (2^bandwidth + 1) / (2^bandwidth - 1) );
  31. % calcuate sigma
  32. sigma = slratio * lambda;
  33. sigma_x = sigma;
  34. sigma_y = sigma/gamma;
  35. % Bounding box
  36. nstds = 4;
  37. xmax = max(abs(nstds*sigma_x*cos(theta)),abs(nstds*sigma_y*sin(theta)));
  38. xmax = ceil(max(1,xmax));
  39. ymax = max(abs(nstds*sigma_x*sin(theta)),abs(nstds*sigma_y*cos(theta)));
  40. ymax = ceil(max(1,ymax));
  41. xmin = -xmax; ymin = -ymax;
  42. [x,y] = meshgrid(xmin:xmax,ymin:ymax);
  43. % Rotation
  44. x_theta = x*cos(theta) + y*sin(theta);
  45. y_theta = -x*sin(theta) + y*cos(theta);
  46. % Gabor Function
  47. gb= exp(-.5*(x_theta.^2/sigma_x^2+y_theta.^2/sigma_y^2)).*cos(2*pi/lambda*x_theta+phi);
  48. 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. HarfH = fix(GaborH/2);
  24. HarfW = fix(GaborW/2);
  25. Qu = pi*U/8;
  26. sqsigma = sigma*sigma;
  27. Kv = 2.5*pi*(2^(-(V+2)/2));
  28. %Kv = Kmax/(f^V);
  29. postmean = exp(-sqsigma/2);
  30. for j = -HarfH : HarfH
  31. for i = -HarfW : HarfW
  32. tmp1 = exp(-(Kv*Kv*(j*j+i*i)/(2*sqsigma)));
  33. tmp2 = cos(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - postmean;
  34. %tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j) - exp(-sqsigma/2);
  35. tmp3 = sin(Kv*cos(Qu)*i+Kv*sin(Qu)*j);
  36. GaborReal(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp2/sqsigma;
  37. GaborImg(j+HarfH+1, i+HarfW+1) = Kv*Kv*tmp1*tmp3/sqsigma;
  38. end
  39. end
  40. 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. figure;
    9. for i = 1:5
    10. for j = 1:8
    11. gaborFilter=gaborKernel2d(lambda(i), theta(j), phi, gamma, bw);
    12. % 查看每一个滤波器
    13. %figure;
    14. %imshow(real(gaborFilter),[]);
    15. % 将所有的滤波器放到一张图像中查看,查看滤波器组
    16. subplot(5,8,(i-1)*8+j);
    17. imshow(real(gaborFilter),[]);
    18. end
    19. end