function C = corr(F, epsilon)
% Compute inhibition map among filters
% F: the bank of filters
% epsilon: threshold for correlation
% C: the inhibition matrix

N = size(F, 2);   % N: the number of filters

for (i = 1:N) 
   hi = (size(F{i}, 1)-1)/2;   % half size of filter i
   for (j = 1:N)
      hj = (size(F{j}, 1)-1)/2;    % half size of filter j
      I = zeros(2*(hi+hj)+1)+sqrt(-1)*zeros(2*(hi+hj)+1);   
      I((hj+1):(hj+2*hi+1), (hj+1):(hj+2*hi+1)) = F{i}; % put F{i} at the center
      
      % convolve F{i} by F{j}, 4 real-imaginary combinations
      rr = filter2(real(F{j}), real(I), 'same');
      ri = filter2(real(F{j}), imag(I), 'same');
      ir = filter2(imag(F{j}), real(I), 'same');
      ii = filter2(imag(F{j}), imag(I), 'same');
      C{i, j} = ((rr.*rr + ri.*ri + ir.*ir + ii.*ii)<epsilon)*1.0; 
      % C{i,j} is how filter i at (0,0) inhibites filter j at 
      % locations overlapping with filter i
   end
end

