% MakeGaborFilters.m
% Make the Gabor filters for VisNet
% Copyright Edmund Rolls 29 June 2021 https://www.oxcns.org Oxford Centre for Computational Neuroscience
clear; close all hidden; format compact
% List of frequencies and angles to iterate over
freqList = [1 2 4 8]; % from 28 Jan 2013
% freqList = [1 2 4 8 16 32 64 128]; % This makes large filters with low spatial frequencies
asymAngList = [ 0 45 90 135 180 225 270 315]; % the asymmetrical filters
% Note that for the symmetric filters, only angles of 0 45 90 and 135 are needed for VisNet 
N_ANGLES = 8;
TH0 = pi / N_ANGLES * 2.0      % the angle between the filters which range from 0 to 360 deg 
dd = 1. / sqrt(2. * pi);
Display = 1; % 0 for no display. 1 for display
fig = 1;

% Make the gabor filters in Matlab
for freq = 1 : length(freqList) % the default is 4 for VisNet 
     rr(freq) = (2 * power(2, freq)); % the filter size in pixels
%      auxgaba = power(2, -(freq-1)); % this can cause problems if freq==1
%      auxgabaa = power(2, -(freq-1));
     auxgaba = power(2, -(freq-0.99));
     auxgabaa = power(2, -(freq-0.99));
     gsize = 2 * rr(freq); % the size of the filter
     % for each of the orientations at each frequency
     for (ang = 1 : N_ANGLES)
         auxgabc = cos((ang-1) * TH0);
         auxgabs = sin((ang-1) * TH0);

         %build the filter
         for i = 1 : gsize
             for j = 1 : gsize
                 temp_x = auxgaba * i - auxgabaa * rr(freq);
                 temp_y = auxgaba * j - auxgabaa * rr(freq);
                 x = temp_x * auxgabc + temp_y * auxgabs;
                 y = temp_y * auxgabc - temp_x * auxgabs;
                 temp = dd * auxgaba * exp(-(4 * x * x + y * y) / 8); % the 4 influences the width of the filter, and the 4 the length
                 gabor_r(i, j) = temp * (cos(pi * x) - exp(-pi * pi / 2));
                 gabor_i(i, j) = temp * sin(pi * x);
             end
         end
         filename = strcat('gabor', '.', int2str(freqList(freq)), '.', int2str(asymAngList(ang)),'.dat')
         fid = fopen(filename,'w');
         fwrite(fid, gabor_r, 'single'); 
         fclose(fid);
         filename = strcat('gabor', 'a.', int2str(freqList(freq)), '.', int2str(asymAngList(ang)),'.dat')
         fid = fopen(filename,'w');
         fwrite(fid, gabor_i, 'single'); 
         fclose(fid);

         if Display >= 0
             figure(fig)
             imagesc(gabor_r);
             fig = fig + 1;
             figure(fig)
             imagesc(gabor_i);
             fig = fig + 1;
         end                 
     end
end
   
% The following is just a check
for i = asymAngList
    figure(100+i/45+1);
    freq=1;
    ang=i;
    strcat('gabor', 'a.', int2str(freq), '.', int2str(ang),'.dat')
    fid = fopen( strcat('gabor', 'a.', int2str(freq), '.', int2str(ang),'.dat') );
    image= fread(fid, [8*freq 8*freq], 'single'); 
    fclose(fid);
    imagesc(image);
end

%FilterConv4('ramp.dat','',1) % A line like this can be used to test convolultion with these filters