function T = handmadebasis( winsize )
% handmadebasis - create a set of complex-cells
%
% SYNTAX:
% T = handmadebasis( winsize, random, dirs, pos, freqs, length );
%
% PARAMETERS:
% winsize       window size
% random        0: given parameters, n: random parameters (n filters)
% dirs          ordered eg: (0:5)*pi/6, or random eg: [0 pi]
% pos           ordered eg: linspace(-0.4,0.4,8), or random eg: [-0.4 0.4]
% freqs         ordered eg: 10*pi, random eg: [8*pi 12*pi]
% length        RF length/width ratio (one value only please)
%
% OUTPUT:
% T.W           forward transform weights
% T.B           iconic representation of filters
% T.H           complex-cell pooling matrix
% T.Hnz         contains indexes of non-zero elements of T.H
% T.winsize     copy of 'winsize' input parameter
% T.dirs        copy of 'dirs' input parameter
% T.posx        copy of 'pos' input parameter
% T.posy        copy of 'pos' input parameter
% T.freqs       copy of 'freqs' input parameter
% T.length      copy of 'length' input parameter
%

% This is the window size & thickness for the icons (default = 100, 0)
global iconsize iconthick RF

if isempty(iconsize), iconsize = 100; end
if isempty(iconthick), iconthick = 0; end

% Generate the parameter values list
if RF.random<=0,
  
  if RF.random==-1,
    posx = 0;
    posy = RF.pos;
  else
    posx = RF.pos;
    posy = RF.pos;
  end
  
  % Regular lattice
  ind = 0;
  xpos = [];
  ypos = [];
  dirval = [];
  freqval = [];
  
  for x = posx
    for y = posy
      for d = RF.dirs
        for freq = RF.freqs

        ind = ind+1;
        xpos = [xpos x];
        ypos = [ypos y];
        dirval = [dirval RF.dirs];
        freqval = [freqval RF.freqs];
  
        end
      end
    end
  end  
else

  % Random positions
  if size(pos,2)>1,
    xpos = rand(1,RF.random)*(pos(2)-pos(1))+pos(1);
    ypos = rand(1,RF.random)*(pos(2)-pos(1))+pos(1);
  else
    xpos = ones(1,RF.random)*pos;
    ypos = xpos;
  end
    
  % Random dirs
  if size(dirs,2)>1,
    dirval = rand(1,RF.random)*(dirs(2)-dirs(1))+dirs(1);
  else
    dirval = ones(1,RF.random)*dirs;
  end
  
  % Random frequencies
  if size(freqs,2)>1,
    freqval = rand(1,RF.random)*(freqs(2)-freqs(1))+freqs(1);
  else
    freqval = ones(1,RF.random)*freqs;
  end
  
end

totalfilters = 2*size(xpos,2);

% Generate the matrix to hold the filters (and the matrix to hold icons)
A = zeros(winsize^2,totalfilters); 
B = zeros(iconsize^2,totalfilters/2);

% Generate the filters
for ind=1:size(xpos,2)
	
  x = xpos(ind);
  y = ypos(ind);
  d = dirval(ind);
  freq = freqval(ind);
	
  % generate the gabor filters
  win1 = gabor( winsize, x, y, d, freq, 0, 0, RF.length );
  win2 = gabor( winsize, x, y, d, freq, pi/2, 0, RF.length );
    
  % convert into column vectors
  A(:,ind) = reshape( win1, [winsize^2 1] );
  A(:,ind+totalfilters/2) = reshape( win2, [winsize^2 1] );
  
  % draw and save icons
  icon = drawicon( iconsize, x, y, d, freq, RF.length );
  B(:,ind) = reshape( icon, [iconsize^2 1] );
end

% Normalizes to zero mean and unit norm
A = A-ones(size(A,1),1)*mean(A);
A = A./(ones(size(A,1),1)*sqrt(sum(A.^2)));

% Builds the structure to return
T.W = A';
T.B = B;
T.H = [eye(totalfilters/2) eye(totalfilters/2)];

for i=1:(totalfilters/2)
  T.Hnz{i} = find(T.H(i,:));
end

T.winsize = winsize;
T.dirs = RF.dirs;
T.posx = RF.pos;
T.posy = RF.pos;
T.freqs = RF.freqs;
T.length = RF.length;

return

%------------------------------------------------------------------------%


function win = drawicon( winsize, x, y, d, freq, length )
% drawicon - draws the symbol to be used to show activity of complex cell
%

global iconthick

% Initially empty patch
win = zeros(winsize,winsize);

% Dimensions 
longside = 0.32*(length/2)*12*pi/freq;
shortside = 0.14*12*pi/freq;
    
% Draw the oval
if iconthick,
  for coef = linspace(0.95,1.05,10),
    win = drawoval( win, x, y, d, coef*longside, coef*shortside );
  end
else
  win = drawoval( win, x, y, d, longside, shortside );
end

return

%------------------------------------------------------------------------%

function win = drawoval( win, x, y, d, longside, shortside )
% drawoval - draws an ellips
%

% This is the window size
winsize = size(win,2);

% First, calculate bounding oriented rectangle
x1 = x - cos(d)*longside/2 - sin(d)*shortside/2;
y1 = y - sin(d)*longside/2 + cos(d)*shortside/2;
x2 = x + cos(d)*longside/2 - sin(d)*shortside/2;
y2 = y + sin(d)*longside/2 + cos(d)*shortside/2;
x3 = x + cos(d)*longside/2 + sin(d)*shortside/2;
y3 = y + sin(d)*longside/2 - cos(d)*shortside/2;
x4 = x - cos(d)*longside/2 + sin(d)*shortside/2;
y4 = y - sin(d)*longside/2 - cos(d)*shortside/2;

% Convert all coordinates to window coordinates
x1 = x1*winsize + winsize/2;
x2 = x2*winsize + winsize/2;
x3 = x3*winsize + winsize/2;
x4 = x4*winsize + winsize/2;
y1 = y1*winsize + winsize/2;
y2 = y2*winsize + winsize/2;
y3 = y3*winsize + winsize/2;
y4 = y4*winsize + winsize/2;
x = x*winsize + winsize/2;
y = y*winsize + winsize/2;
longside = longside*winsize;
shortside = shortside*winsize;

% Calculate a bounding box (loose!)
left = max(round(min([x1 x2 x3 x4])),1);
right = min(round(max([x1 x2 x3 x4])),winsize);
top = max(round(min([y1 y2 y3 y4])),1);
bottom = min(round(max([y1 y2 y3 y4])),winsize);

% Calculate coefficient matrix
SR = diag([1/(longside/2) 1/(shortside/2)])*[cos(d) sin(d); -sin(d) cos(d)];
A = SR'*SR;

% Step through all rows
for i=top:bottom,

  p = (A(1,2)+A(2,1))/A(1,1)*(i-y);
  q = (A(2,2)*(i-y)^2-1)/A(1,1);
  inroot = p^2/4 - q;
  if inroot>=0,
    p1 = round(-p/2 + sqrt(inroot) + x);
    p2 = round(-p/2 - sqrt(inroot) + x);
    if p1>=1 & p1<=winsize, win(i,p1)=1; end
    if p2>=1 & p2<=winsize, win(i,p2)=1; end
  end
    
end

% Step through all columns
for i=left:right,
  
  p = (A(1,2)+A(2,1))/A(2,2)*(i-x);
  q = (A(1,1)*(i-x)^2-1)/A(2,2);
  inroot = p^2/4 - q;
  if inroot>=0,
    p1 = round(-p/2 + sqrt(inroot) + y);
    p2 = round(-p/2 - sqrt(inroot) + y);
    if p1>=1 & p1<=winsize, win(p1,i)=1; end
    if p2>=1 & p2<=winsize, win(p2,i)=1; end
  end
    
end

return