function  RF = PLM_RFields(Iparams)

%PLM_RFIELDS -- Receptive fields of "V1 neurons" in the PercLearn Model
%
%  RF = PLM_RFields(Iparams)
%
%  This function makes a set of receptive fields that define the first
%  linear-filtering stage of the front-end of the PercLearning model.
%  All RFs are Gabors and their support is the same: IPARAMS.img_size.
%  IPARAMS is an "input parameter structure" produced by PLM_INPUT_PARAMS.
%  The following fields are relevant: 'img_size', 'D_orient', 'D_spfreq',
%  'HW_orient', 'octaves', 'deg_pix' (see PLM_INPUT_PARAMS for details).
%
%  The fields 'D_orient' and 'D_spfreq' define the orientations [in deg
%  from the vertical] and spatial frequencies [in cycles/deg.vis.ang].
%  The result is a structure array RF with N_orient rows and N_spfreq
%  columns. Each structure has the following fields:
%   1. orient -- orientation of the sinusoidal grating [deg]
%   2. spfreq -- spatial frequency of the grating [cycles/deg.v.a.]
%   3. S_perp -- std.deviation along the perpendicular direction
%   4. S_prll -- std.dev. along the parallel dir. (see RFIELD_SIZE)
%   5. Wcos -- weights of the kernel with phase 0, [64x64]=img_size
%   6. Wsin -- weights of the kernel with phase pi/2, [64x64]
%
%  Example:
%    RF=PLM_RFields(PLM_input_params); for o=1:7 for k=1:5 ...
%    subplot(7,5,5*o+k-5);imagesc(RF(o,k).Wsin);axis square;end;end
%
%  See also MAKE_PLFRONTEND, PLM_INPUT_PARAMS, RFIELD_SIZE, GABOR.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 2.1 $  $Date: 2005/03/11 12:00 $

%-- Decipher Iparams
orient = Iparams.D_orient.*(pi/180) ;    % degrees -->radians
N_orient = length(orient) ;
spfreq = Iparams.D_spfreq ;
N_spfreq = length(spfreq) ;

%-- Prepare the mesh
img_sz = Iparams.img_size ;
deg_pix = Iparams.deg_pix ;
Xmin = -floor(img_sz(2)/2) ; Ymin = -floor(img_sz(1)/2) ;
x = ([0:img_sz(2)-1]+Xmin) .* deg_pix ;
y = ([0:img_sz(1)-1]'+Ymin) .* deg_pix ;
[X,Y] = meshgrid(x,y) ;

%-- Standard deviations of the Gabor envelopes
[all_S_perp,all_S_prll] = RField_size(Iparams) ;  % [N_orient,N_spfreq]

%-- Main loop
for r = 1:N_spfreq
   omega = 2*pi*spfreq(r) ;
   for c = 1:N_orient
      theta0 = orient(c) ;
      S_perp = all_S_perp(c,r) ; % envelope stdev perpendicular to grtg
      S_prll = all_S_prll(c,r) ; % envelope std.dev parallel to grating
      %- Rotate so that the bars are parallel to the Y axis
      sin_theta0 = sin(theta0) ;
      cos_theta0 = cos(theta0) ;
      rotX = X.*cos_theta0 - Y.*sin_theta0 ;
      rotY = Y.*cos_theta0 + X.*sin_theta0 ;
      %- Sinusoidal gratings (with ORIENT=0 in rotated coordinates)
      Gcos = cos(rotX.*omega) ;
      Gsin = sin(rotX.*omega) ;
      %- Gaussian envelope (rotX and rotY are independent)
      gaussX = exp(-0.5 .* (rotX./S_perp).^2) ;
      gaussY = exp(-0.5 .* (rotY./S_prll).^2) ;
      envelope = gaussX .* gaussY ;
      envelope = envelope ./ (sum(envelope(:))) ;   % normalize
      %- Put it all together to produce weights
      Wcos = Gcos .* envelope ;
      Wsin = Gsin .* envelope ;
      %- Pack and return
      s.orient = theta0 ;
      s.spfreq = spfreq(r) ;
      s.S_perp = S_perp ;
      s.S_prll = S_prll ;
      s.Wcos = Wcos ;
      s.Wsin = Wsin ;
      RF(c,r) = s ;
   end  % for c = 1:N_orient
end  % for r = 1:N_spfreq

%--- Return RF
%%%%% End of file PLM_RFIELDS.M
