Working on Separable Gabor filters in matlab

Question

A filter g is called separable if it can be expressed as the multiplication of two vectors grow and gcol . Employing one dimensional filters decreases the two dimensional filter's computational complexity from O(M^2 N^2) to O(2M N^2) where M and N are the width (and height) of the filter mask and the image respectively.

In this stackoverflow link, I wrote the equation of a Gabor filter in the spatial domain, then I wrote a matlab code which serves to create 64 gabor features.

According to the definition of separable filters, the Gabor filters are parallel to the image axes - theta = k*pi/2 where k=0,1,2,etc.. So if theta=pi/2 ==> the equation in this stackoverflow link can be rewritten as:

The equation above is extracted from this article.

Note: theta can be extented to be equal k*pi/4. By comparing to the equation in this stackoverflow link, we can consider that f= 1 / lambda.

By changing my previous code in this stackoverflow link, I wrote a matlab code to make the Gabor filters separable by using the equation above, but I am sure that my code below is not correct especially when I initialized the gbp and glp equations. That is why I need your help. I will appreciate your help very much.

Let's show now my code:

    function [fSiz,filters1,filters2,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)
filters1          = zeros(max(RF_siz),numFilters);
filters2          = zeros(numFilters,max(RF_siz));

for k = 1:numFilterSizes  
    for r = 1:numSimpleFilters
        theta     = rot(r)*pi/180;
        filtSize  = RF_siz(k);
        center    = ceil(filtSize/2);
        filtSizeL = center-1;
        filtSizeR = filtSize-filtSizeL-1;
        sigmaq    = sigma(k)^2;

        for x = -filtSizeL:filtSizeR

                    fx = exp(-(x^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
                       f1(x+center,1) = fx;
        end
                      for y = -filtSizeL:filtSizeR
                    gy = exp(-(y^2)/(2*sigmaq));
                    f2(1,y+center) = gy;
                      end



        f1 = f1 - mean(mean(f1));
        f1 = f1 ./ sqrt(sum(sum(f1.^2)));
         f2 = f2 - mean(mean(f2));
        f2 = f2 ./ sqrt(sum(sum(f2.^2)));
        p = numSimpleFilters*(k-1) + r;
        filters1(1:filtSize,p)=f1;
        filters2(p,1:filtSize)=f2;

        convv1=imfilter(image_double,  filters1(1:filtSize,p),'conv');
        convv2=imfilter(double(convv1),  filters2(p,1:filtSize),'conv');
        figure
        imagesc(convv2);
        colormap(gray);

    end
end

Answer 1

I think the code is correct provided your previous version of Gabor filter code is correct too. The only thing is that if theta = k * pi/4;, your formula here should be separated to:

fx = exp(-(x^2)/(2*sigmaq))*cos(2*pi*x/lambda(k));
gy = exp(-(G^2 * y^2)/(2*sigmaq));

To be consistent, you may use

f1(1,x+center) = fx;
f2(y+center,1) = gy;

or keep f1 and f2 as it is but transpose your filters1 and filters2 thereafter. Everything else looks good to me.

EDIT

My answer above works for theta = k * pi/4;, with other angles, based on your paper,

x = i*cos(theta) - j*sin(theta);
y = i*sin(theta) + j*cos(theta);
fx = exp(-(x^2)/(2*sigmaq))*exp(sqrt(-1)*x*cos(theta));
gy = exp(-(G^2 * y^2)/(2*sigmaq))*exp(sqrt(-1)*y*sin(theta));

Answer 2

The final code will be:

function [fSiz,filters1,filters2,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];
                            RF_siz    = [7:2:37]; 
                            minFS     = 7; 
                            maxFS     = 37; 
                            sigma  = 0.0036*RF_siz.^2 + 0.35*RF_siz + 0.18; 
                            lambda = sigma/0.8; 
                            G      = 0.3;   


                            numFilterSizes   = length(RF_siz); 

                            numSimpleFilters = length(rot); 

                            numFilters       = numFilterSizes*numSimpleFilters; 

                            fSiz             = zeros(numFilters,1); 
                            filters1          = zeros(max(RF_siz),numFilters);
                            filters2          = zeros(numFilters,max(RF_siz));

        for k = 1:numFilterSizes  
            for r = 1:numSimpleFilters
                theta     = rot(r)*pi/180;
                filtSize  = RF_siz(k);
                center    = ceil(filtSize/2);
                filtSizeL = center-1;
                filtSizeR = filtSize-filtSizeL-1;
                sigmaq    = sigma(k)^2;

                for x = -filtSizeL:filtSizeR

                            fx = exp(-(x^2)/(2*sigmaq))*exp(sqrt(-1)*x*cos(theta));
                               f1(1, x+center) = fx;
                end
                              for y = -filtSizeL:filtSizeR
                            gy=exp(-(y^2)/(2*sigmaq))*exp(sqrt(-1)*y*sin(theta));
                            f2(y+center,1) = gy;
                              end


         f1 = f1 - mean(mean(f1));
            f1 = f1 ./ sqrt(sum(sum(f1.^2)));
             f2 = f2 - mean(mean(f2));
            f2 = f2 ./ sqrt(sum(sum(f2.^2)));
            p = numSimpleFilters*(k-1) + r;

            filters1(1:filtSize,p)=f1;
            filters2(p,1:filtSize)=f2;

            convv1=imfilter(image_double,  filters1(1:filtSize,p),'conv');
            convv2=imfilter(double(convv1),  filters2(p,1:filtSize),'conv');

            figure
            imagesc(imag(convv2));
            colormap(gray);

        end
    end