How to create 64 Gabor features at each scale and orientation in the spatial and frequency domain

Question

Normally, a Gabor filter, as its name suggests, is used to filter an image and extract everything that it is oriented in the same direction of the filtering.

In this question, you can see more efficient code than written in this Link

Assume 16 scales of Filters at 4 orientations, so we get 64 gabor filters.

scales=[7:2:37], 7x7 to 37x37 in steps of two pixels, so we have 7x7, 9x9, 11x11, 13x13, 15x15, 17x17, 19x19, 21x21, 23x23, 25x25, 27x27, 29x29, 31x31, 33x33, 35x35 and 37x37.

directions=[0, pi/4, pi/2, 3pi/4].

The equation of Gabor filter in the spatial domain is:

The equation of Gabor filter in the frequency domain is:

Answer 1

In the spatial domain:

function [fSiz,filters,c1OL,numSimpleFilters] = init_gabor(rot, RF_siz)    

    image=imread('xxx.jpg');
    image_gray=rgb2gray(image);
    image_gray=imresize(image_gray, [100 100]);
    image_double=double(image_gray);

    rot = [0 45 90 135]; % we have four orientations
                RF_siz    = [7:2:37]; %we get 16 scales (7x7 to 37x37 in steps of two pixels)
                minFS     = 7; % the minimum receptive field
                maxFS     = 37; % the maximum receptive field
                sigma  = 0.0036*RF_siz.^2 + 0.35*RF_siz + 0.18; %define the equation of effective width
                lambda = sigma/0.8; % it the equation of wavelength (lambda)
                G      = 0.3;   % spatial aspect ratio: 0.23 < gamma < 0.92


                numFilterSizes   = length(RF_siz); % we get 16

                numSimpleFilters = length(rot); % we get 4

                numFilters       = numFilterSizes*numSimpleFilters; % we get 16x4 = 64 filters

                fSiz             = zeros(numFilters,1); % It is a vector of size numFilters where each cell contains the size of the filter (7,7,7,7,9,9,9,9,11,11,11,11,......,37,37,37,37)

                filters          = zeros(max(RF_siz)^2,numFilters); % Matrix of Gabor filters of size %max_fSiz x num_filters, where max_fSiz is the length of the largest filter and num_filters the total number of filters. Column j of filters matrix contains a n_jxn_j filter (reshaped as a column vector and padded with zeros).




            for k = 1:numFilterSizes  
                for r = 1:numSimpleFilters
                    theta     = rot(r)*pi/180; % so we get 0, pi/4, pi/2, 3pi/4
                    filtSize  = RF_siz(k); 
                    center    = ceil(filtSize/2);
                    filtSizeL = center-1;
                    filtSizeR = filtSize-filtSizeL-1;
                    sigmaq    = sigma(k)^2;

                    for i = -filtSizeL:filtSizeR
                        for j = -filtSizeL:filtSizeR

                            if ( sqrt(i^2+j^2)>filtSize/2 )
                                E = 0;
                            else
                                x = i*cos(theta) - j*sin(theta);
                                y = i*sin(theta) + j*cos(theta);
                                E = exp(-(x^2+G^2*y^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
                            end
                            f(j+center,i+center) = E;
                        end
                    end

                    f = f - mean(mean(f));
                    f = f ./ sqrt(sum(sum(f.^2)));
                    p = numSimpleFilters*(k-1) + r;
                    filters(1:filtSize^2,p)=reshape(f,filtSize^2,1);
                    fSiz(p)=filtSize;
                end
            end

            % Rebuild all filters (of all sizes)

            nFilts = length(fSiz);
            for i = 1:nFilts
              sqfilter{i} = reshape(filters(1:(fSiz(i)^2),i),fSiz(i),fSiz(i));

            %if you will use conv2 to convolve an image with this gabor, so you should also add the equation below. But if you will use imfilter instead of conv2, so do not add the equation below.

                    sqfilter{i} = sqfilter{i}(end:-1:1,end:-1:1); %flip in order to use conv2 instead of imfilter (%bug_fix 6/28/2007);

    convv=imfilter(image_double, sqfilter{i}, 'same', 'conv');
    figure;

        imagesc(convv);
        colormap(gray);

                      end 

Answer 2

phi = ij*pi/4; % ij = 0, 1, 2, 3
theta = 3;
sigma = 0.65*theta;
filterSize = 7;   % 7:2:37

G = zeros(filterSize);


for i=(0:filterSize-1)/filterSize
    for j=(0:filterSize-1)/filterSize
        xprime= j*cos(phi);
        yprime= i*sin(phi);
        K = exp(2*pi*theta*sqrt(-1)*(xprime+ yprime));
        G(round((i+1)*filterSize),round((j+1)*filterSize)) =...
           exp(-(i^2+j^2)/(sigma^2))*K;
    end
end

Answer 3

As of R2015b release of the Image Processing Toolbox, you can create a Gabor filter bank using the gabor function in the image processing toolbox, and you can apply it to an image using imgaborfilt.

Answer 4

In the frequency domain:

sigma_u=1/2*pi*sigmaq;
sigma_v=1/2*pi*sigmaq;
u0=2*pi*cos(theta)*lambda(k);
v0=2*pi*sin(theta)*lambda(k);

for u = -filtSizeL:filtSizeR
            for v = -filtSizeL:filtSizeR

                if ( sqrt(u^2+v^2)>filtSize/2 )
                    E = 0;
                else
                    v_theta = u*cos(theta) - v*sin(theta);
                    u_theta = u*sin(theta) + v*cos(theta);

                  E=(1/2*pi*sigma_u*sigma_v)*((exp((-1/2)*(((u_theta-u0)^2/sigma_u^2))+((v_theta-v0)^2/sigma_v^2))) + (exp((-1/2)*(((u_theta+u0)^2/sigma_u^2))+((v_theta+v0)^2/sigma_v^2))));

                end
                f(v+center,u+center) = E;
            end
        end