how can i understand this code (about Gabor wavelet)

Question

ZhG on 20 Jun 2013

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/79682-how-can-i-understand-this-code-about-gabor-wavelet

Answered: Mohammad Khodajouei on 4 Dec 2014

I download a code at vision lab of UCSB. It is about Gabor wavelet, and there are something that I don't understand in the code. Is there anyone familiar with Gabor?

I listed questions and code below. First, I run the program with the arguments below to generate a Gabor filter. [Gr1, Gi1] = gabor(13, [2, 0], [0.25, 0.5], [2, 45], 0)

    So, questions are:
    First one is about the calculation of removing redundant information.
    In the paper, a = (Uh/Ul)^(1/(s-1)).
    Thus, I don't know why does it need to solve an equation to obtain 'a' in the program.
    Q1. N=13. Is this the size of gabor filter?
    Q2. Later in the codes, there is    side = fix(N/2);
           Is this the Gaussian envelope size, i.e. axis length of Gaussian envelope?
    Then, I confused about the second argument and the fourth one. 
    For the second argument of this function [s, n],
    It still works while I assign 0 to n. I have this question just because of a subscript containing n-1 in the program.
    t1 = cos(pi/orientation*(n-1));
    Does it mean that I can assign several orientations with the fourth argument, such as
    [2, 0:pi/4:pi]
    Q3. What does [s, n] mean at all?
    Q4, And what does [stage, orientation] mean at all?
    Q5, Is there any relationship between these two arguments?
%----------------------------------------------------------------------
% This function generate the spatial domain of the Gabor wavelets
% which are specified by number of scales and orientations and the 
% maximun and minimun center frequency.
% 
%         N : the size of rectangular grid to sample the gabor
%     index : [s,n] specify which gabor filter is selected
%      freq : [Ul,Uh] specify the maximun and minimun center frequency
% partition : [stage,orientation] specify the total number of filters
%      flag : 1 -> remove the dc value of the real part of Gabor
%             0 -> not to remove
%----------------------------------------------------------------------
function [Gr,Gi] = Gabor(N,index,freq,partition,flag)
% get parameters
s = index(1);
n = index(2);
Ul = freq(1);
Uh = freq(2);
stage = partition(1);
orientation = partition(2);
% computer ratio a for generating wavelets
base = Uh/Ul;
C = zeros(1,stage);
C(1) = 1;
C(stage) = -base;
P = abs(roots(C));
a = P(1);
% computer best variance of gaussian envelope
u0 = Uh/(a^(stage-s));
Uvar = ((a-1)*u0)/((a+1)*sqrt(2*log(2)));
z = -2*log(2)*Uvar^2/u0;
Vvar = tan(pi/(2*orientation))*(u0+z)/sqrt(2*log(2)-z*z/(Uvar^2));
% generate the spetial domain of gabor wavelets
j = sqrt(-1);
if (rem(N,2) == 0)
    side = N/2-0.5;
else
    side = fix(N/2);
end;
x = -side:1:side;
l = length(x);
y = x';
X = ones(l,1)*x;
Y = y*ones(1,l);
t1 = cos(pi/orientation*(n-1));
t2 = sin(pi/orientation*(n-1));
XX = X*t1+Y*t2;
YY = -X*t2+Y*t1;
Xvar = 1/(2*pi*Uvar);
Yvar = 1/(2*pi*Vvar);
coef = 1/(2*pi*Xvar*Yvar);
Gr = a^(stage-s)*coef*exp(-0.5*((XX.*XX)./(Xvar^2)+(YY.*YY)./(Yvar^2))).*cos(2*pi*u0*XX);
Gi = a^(stage-s)*coef*exp(-0.5*((XX.*XX)./(Xvar^2)+(YY.*YY)./(Yvar^2))).*sin(2*pi*u0*XX);
% remove the real part mean if flag is 1
if (flag == 1)
   m = sum(sum(Gr))/sum(sum(abs(Gr)));
   Gr = Gr-m*abs(Gr);
end;

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Sign in to answer this question.

Answer 1

Mohammad Khodajouei on 4 Dec 2014

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/79682-how-can-i-understand-this-code-about-gabor-wavelet#answer_161209

the formula you have presented here "a = (Uh/Ul)^(1/(s-1))" would work as well and it generates the same result as this bit of the code:

base = Uh/Ul; C = zeros(1,stage); C(1) = 1; C(stage) = -base; P = abs(roots(C)); a = P(1);

N defines the size of Gabor kernel. So, for example if N = 3 then your kernel size is 3x3.

if you have confused by this bit here you can change this part:

"if (rem(N,2) == 0) side = N/2-0.5; else side = fix(N/2); end; x = -side:1:side; l = length(x); y = x'; X = ones(l,1)*x; Y = y*ones(1,l);"

with this: "[x,y] = meshgrid(fix(-k/2):fix(k/2), fix(-k/2):fix(k/2));" it does the same

The Gabor bank is defined by multiplying the total number of scales and orientations. Here in this code it is specified as "Partition". So if the total number of Stage or Scale is 5 and the total number of orientation is 4 then the Gabor bank is an array {5,4} which contains 20 filters. So here s and n referring index in the Gabor bank.