filt_get_s1.m

Standard model object recognition matlab code

function f = filt_get_S1

% FUNCTION f = filt_get_S1
%
% Here we construct orientation-tuned V1 simple cell's receptive
% fields, based on Gabor functions.  We assume that there are
% discrete number of orientations (default = 4) and scales (default
% = 8).  These parameters are more or less fixed and defined
% throughout this function.

f = [];
num_scale = 8; 
for i = 1:num_scale
  scale_tag = ['s' num2str(i)];
  [f1, f2, shift] = get_S1_filters(i);
  f = filt_package (f, f1, f2, shift, scale_tag);
end

function [f1, f2, shift] = get_S1_filters (which_scale)

% FUNCTION  [f1, f2, shift] = get_S1_filters (which_scale)
% 
% This function gets the Gabor filters at different scales, using
% the parameters defined (work of T. Serre) according to the
% physiological data in V1 simple cells.
%
% In the old HMAX (Nature Neurosci 1999) parameters,
%   filter_sizes = [7 9], [11 13 15], [17 19 21], [23 25 27 20];

shift = 1;

div = [4:-.05:3.2];
switch which_scale
  case 1
    filter_sizes = [7 9];
    Div          = div(1:2);
  case 2
    filter_sizes = [11 13];
    Div          = div(3:4);
  case 3
    filter_sizes = [15 17];
    Div          = div(5:6);
  case 4
    filter_sizes = [19 21];
    Div          = div(7:8);
  case 5
    filter_sizes = [23 25];
    Div          = div(9:10);
  case 6
    filter_sizes = [27 29];
    Div          = div(11:12);
  case 7
    filter_sizes = [31 33];
    Div          = div(13:14);
  case 8
    filter_sizes = [35 37 39];
    Div          = div(15:17);
end
num_scale = size(filter_sizes,2);
max_fs = max(filter_sizes);

rot = [90 -45 0 45]; % four default orientations (can be more)
num_orientation = length(rot);

f1 = zeros(max_fs,max_fs,1, num_scale*num_orientation);
f2 = zeros(max_fs,max_fs,1, num_scale*num_orientation);

for i = 1:num_scale
  % embed in zeros.
  half_fs = round((max_fs-filter_sizes(i))/2);
  filt_range = [1:filter_sizes(i)]+half_fs;

  filt1 = zeros(max_fs, max_fs, num_orientation);
  filt1(filt_range,filt_range,:) = get_gabor(rot, filter_sizes(i), Div(i));

  filt2 = zeros(max_fs, max_fs);
  filt2(filt_range,filt_range) = get_circle(filter_sizes(i),1);

  for j = 1:num_orientation
    f1(:,:,1,i+(j-1)*num_scale) = filt1(:,:,j);
    f2(:,:,1,i+(j-1)*num_scale) = filt2;
  end
end

function filters = get_gabor (rot, RF_siz, Div)

% FUNCTION filters = get_gabor (rot, RF_siz, Div)
% 
% Computes and returns Gabor functions.

num_rot = length(rot); % number of rotaion/orientation
lambda  = RF_siz*2./Div;
sigma   = lambda.*0.8;
G       = 0.3;   % spatial aspect ratio: 0.23 < gamma < 0.92

for r = 1:num_rot
  theta     = rot(r)*pi/180;
  filtSize  = RF_siz;
  center    = ceil(filtSize/2);
  filtSizeL = center-1;
  filtSizeR = filtSize-filtSizeL-1;
  sigmaq    = sigma^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);
      end
      f(j+center,i+center) = E;
    end
  end

  % apply circular mask, make it mean-zero, and normalize.
  f      = f.*get_circle(size(f,1),1); 
  f_norm = sum(sum(f))/sum(sum(get_circle(size(f,1),1)));
  f      = f - get_circle(size(f,1),1)*f_norm;
  f      = f/sqrt(sum(sum(f.^2)));
  filters(:,:,r) = f; 
end

function circle_template = get_circle (filter_size, radius)

% FUNCTION  circle_template = get_circle (size, radius)
% This function Returns 2-D circular template 
% of given size and radius.

inc = 2/filter_size;
[x,y] = meshgrid(-1+inc/2:inc:1-inc/2, -1+inc/2:inc:1-inc/2);
circle_template = double(x.^2 + y.^2 <= radius.^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [filt] = gabor (filter_size)

% function [filter] = gabor (filter_size, sig_xy, k, phase)
% The function returns 2-D Gabor functions of 
% the given filter_size at various orientations.
% The 2-D plane is treated as an interval of 2*pi.
% Wave number k determines how many cycles are present.

% Default values
% k = 2.1; sig_xy = [1/3 1/1.8]*2*pi; phi = 0;
%
% THIS IS AN OLDER (M.KOUH's) VERSION OF GABOR FUNCTION.  WE ARE
% USING T.SERRE'S GET_GABOR INSTEAD.

global num_orientation;
global gabor_k;
global gabor_sig_xy;
global gabor_phase;

filt = zeros(filter_size, filter_size, num_orientation);
filt_tmp = zeros(filter_size);

k = gabor_k;
phi = gabor_phase;
Sx = gabor_sig_xy(1)^2;
Sy = gabor_sig_xy(2)^2;

inc = 2*pi/filter_size;
xy_range = [-pi+inc/2:inc:pi-inc/2];
[x,y] = meshgrid(xy_range, xy_range);

% for circular receptive field, use radius = 1.
circ = get_circle(filter_size, 1);
circ_sum = sum(sum(circ));

for i = 1:num_orientation
  % theta = rotation angle.
  theta = pi/num_orientation*(i-1);
  u = x*cos(theta) - y*sin(theta);
  v = x*sin(theta) + y*cos(theta);

  % Gabor function
  filt_tmp = exp(-u.^2/2/Sx-v.^2/2/Sy).*cos(k*u-phi);
  filt_tmp = filt_tmp.*circ;
  filt_tmp_sum = sum(sum(filt_tmp));

  % Normalize the filter.
  filt_tmp = filt_tmp - circ*filt_tmp_sum/circ_sum;
  filt_tmp = filt_tmp / sqrt(sum(sum(filt_tmp.^2)));
  filt(:,:,i) = filt_tmp;
end