Apply Gabor filter to input image using Matlab

[Liszt2015]

The ventral stream is described by the following route:

v1--->v2--->v4--->IT

Units of V1 (in S1 layer of) implement Gabor filters. In this layer, an input grayscale image (120x120 or 160x160) is densely filtered by a battery of gabor filters at each scale and orientation. Therefore, at each pixel of the input image, filters of each size and orientation are centered.

Note: The filters come in 4 orientations and 16 scales (so 16x4 = 64 maps) that are arranged in 8 bands as you show in the image below:


After applying the equation below to the input image, the result of this layer will be a set of images filtered with the different sizes and orientations.





For Now, i am very interested to implement this layer by using matlab code. In my first step, i will use only a gabor filter of size 7x7(i.e lambda=3.5 and bandwidth=2.8) and orientation = 90 degree.

Let us consider the input image "Best_Friends.jpg".



By the way, i wrote the matlab code for this layer as shown below, but i am not sure that it's true, especially when i replace the variables x and y by the pixels of image. That's why i need your help in order to optimize this code and to correct it properly. In my future work, i will apply the same code (after being corrected) to all sizes and orientations. PLEASE it's so urgent, i need your appreciated help.

%Read the original RGB input image
        image=imread('Best_Friends.jpg');
    %convert it to gray scale
        image_gray=rgb2gray(image);
    %resize the image to 160x160 pixels
        image_resize=imresize(image_gray, [160 160]);
    %apply im2double
        image_resize=im2double(image_resize);
    %show the image
        figure(1);
        imshow(image_resize);
        title('Input Image');
        
        %Gabor filter size 7x7 and orientation 90 degree
    %declare the variables
        gamma=0.3; %aspect ratio
        psi=0; %phase
        theta=90; %orientation
        bw=2.8; %bandwidth or effective width
        lambda=3.5; % wavelength
        pi=180;
        
        for x=1:160
            for y=1:160
                
        x_theta=image_resize(x,y)*cos(theta)+image_resize(x,y)*sin(theta);
        y_theta=-image_resize(x,y)*sin(theta)+image_resize(x,y)*cos(theta);
        
        gb(x,y)= exp(-(x_theta.^2/2*bw^2+ gamma^2*y_theta.^2/2*bw^2))*cos(2*pi/lambda*x_theta+psi);
        
            end
        end
        
        figure(2);
        imshow(gb);
        title('filtered image');