comp_get_next_layer.m

Standard model object recognition matlab code

function y = comp_get_next_layer (x,f,type)

% FUNCTION y = comp_get_next_layer (x,f,type)
%
% This function loops thru different scales, and at each scales,
% calls get_next_layer function to build the next layer, based on
% the operation "type".  This is the main computational function.

y = [];

% read which scales are there in f.
scales = aux_get_filt_info(fieldnames(f));

% Only do scales that are in the filters.
for i = 1:length(scales)


  f1 = getfield(f, ['f1_s' num2str(scales(i))]);
  % If f1 is not defined, use the first one.
  if length(f1) <= 0
    scale_tag = 's1';
  else
    scale_tag = ['s' num2str(scales(i))];
  end

  f1        = getfield(f, ['f1_'    scale_tag]);
  f2        = getfield(f, ['f2_'    scale_tag]);
  i1        = getfield(f, ['i1_'    scale_tag]);
  i2        = getfield(f, ['i2_'    scale_tag]);
  i3        = getfield(f, ['i3_'    scale_tag]);
  siz       = getfield(f, ['size_'  scale_tag]);
  shift     = getfield(f, ['shift_' scale_tag]);

  scale_tag = ['s' num2str(scales(i))];
  if isstruct(x)
    x_tmp = getfield(x, scale_tag);
  else % this is the stimulus (not structure)
    x_tmp = x;
  end

  % x will be either cropped or padded with zeros to ensure proper
  % convolution.  Depending on the operation type, cropping or
  % zeropadding may be more appropriate.
  x_tmp = comp_crop_or_zeropad (x_tmp, siz(1:2), shift, type);

  % Convert subscripts into index, based on the stimulus size.
  idx   = aux_sub2ind(size(x_tmp),i1,i2,i3);
  y_tmp = get_next_layer(x_tmp,f1,f2,siz,idx,shift,type);
  y     = setfield(y,scale_tag,y_tmp);
end

function y = get_next_layer (x,f1,f2,siz,idx,shift,type)

% FUNCTION y = get_next_layer (x,f1,f2,siz,idx,shift,type)
%
% This function computes the response of the next_layer, based on
% the input prev_layer, based on an operation specified by "type",
% which depends on the parameters f1 and f2.  Normally, f1 and f2
% will be the weights in the weighted sum of normalized scalar
% product (f1 is for the numerator, f2 is for the denominator).
% The dimension of each input arguement is:
%   x = x_prev by y_prev by (number of old features),
%   f1 = x_f by y_f by (num of old features) by (num of new features),
%   f2 = x_f by y_f by (num of old features) by (num of new features),
%   y = x_new by y_new by (num of new features)
% This main operations are multiple 3-d convolutions, implemented
% in mex-c files.
%   y = sum(f1.*(x.^p))./(c+(sum(f2.*(x.^q)).^r);
% Depending on the relative magnitude of p, q, and r, we can
% implement max-like invariance operation or gaussian-like
% specificity operation.
% The arguement "shift" determines how many pixels to shift in x-y
% directions during convolution.
% The weights, f1 and f2, are formated in a sparse format as a
% structure, which contains nonzero weights (f1, f2), subscripts
% pointing to those values (i1, i2, i3), and their
% original sizes (siz).  The original filter sizes need to
% be specified, because the size of x and size of filter are not
% necessarily the same.

x = double(x); % if not double, some typecast problem with mex file.

if length(x) <= 0
  % if the size of filter is greater than the input x,
  % for the "cropping" case (ie. nothing to crop), just return zeros.
  y       = zeros(1,1,siz(4));
else
  mn_size = [size(x,1) size(x,2) siz(1) siz(2)];
  x_r     = reshape(x, [prod(size(x)) 1]);

  if idx(end) > 0
    % Avoid infinite loop in the mex file.
    error('Filters do not end in zeros.');
  end

  idx = double(idx-1); % C index starts from 0; Matlab from 1

  switch type
    case ('nullop')
      y = do_sp1(x_r, f1, idx, mn_size, shift, 0);
    case ('gaussf')
      y = do_sp1(x_r, f1, idx, mn_size, shift, 1);
      global gauss_sigma;
      % This makes the Gaussian function to be 0.01 at the origin.
      gauss_sigma = repmat(sum(f1.^2), [size(y,1) 1])+eps;
      gauss_sigma = sqrt(-gauss_sigma/2/log(0.01));
      y = exp(-y/2./gauss_sigma.^2);
    case ('fulmax')
      y = do_sp1(x_r, f1, idx, mn_size, shift, 2);
    case ('dotpro')
      y = do_sp1(x_r, f1, idx, mn_size, shift, 100);
    otherwise
      init_operation_def;
      %tmp1 = do_sp1(x_r.^p, f1, idx, mn_size, shift, 3);
      %tmp2 = do_sp1(x_r.^q, f2, idx, mn_size, shift, 3);
      %y = abs(tmp1)./(repmat(c,[size(tmp1,1) 1])+abs(tmp2).^r);
      [y1, y2] = do_sp2(x_r.^p, x_r.^q, f1, f2, idx, mn_size, shift);

      y = abs(y1)./(repmat(c,[size(y1,1) 1]) + abs(y2.^r));
      y = y./repmat(y_max, [size(y,1) 1]);
      
      % We now rescale (sigmoid) the response in main_model.m.
      %if sigmoid_a > 0
	%y = 1./(1+exp(-sigmoid_a.*(y-sigmoid_b))); 
      %end
  end

  % The input image may not be square.
  size_m = (mn_size(1)-mn_size(3))/shift+1;
  size_n = (mn_size(2)-mn_size(4))/shift+1;
  % Return the output as a 3-D matrix.
  y = reshape(y, [size_m, size_n, siz(4)]);
end