function  Idef = make_PLFrontEnd(Iparams,fid)

%MAKE_PLFRONTEND -- Set up the front end of the Perceptual Learning Model.
%
%  Idef = make_PLFrontEnd(Iparams)
%  Idef = make_PLFrontEnd(Iparams,fid)
%
%  IPARAMS is a structure generated by PLM_INPUT_PARAMS.
%  FID is a stream for "verbose" messages: 0=none, 1=stdout (default).
%
%  MAKE_PLFRONTEND calculates the front-end filters and everything else
%  necessary for applying the model to real images. It stores all these
%  data structures in the global variable FRONTEND_MEMOIZATION for later
%  access inside PLFRONTEND. The user should not alter this global
%  variable. Instead, MAKE_PLFRONTEND *MUST* be called whenever the
%  IPARAMS have been changed in order to re-create the memoization.
%
%  IDEF is a structure with the same fields as IPARAMS and the following
%  additional fields:
%    * memo_CRC -- for checking the consistency of FRONTEND_MEMOIZATION.
%    * pad_XY -- index [1,1] of the unpadded stimulus maps to PAD_XY+1
%        in the "bakground-padded" image. See also IPARAMS.bkgr_idx.
%
%  See also PLM_FRONTEND, PLM_INPUT_PARAMS.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 2.2 $  $Date: 2003/05/30 10:31 $

global FRONTEND_MEMOIZATION

if (nargin<2)  fid = 1 ; end    % verbose messages to standard output

%-- Overwrite the memoization store
FRONTEND_MEMOIZATION = struct('CRC',now,'Iparams',Iparams) ;
Idef = Iparams ;
Idef.memo_CRC = FRONTEND_MEMOIZATION.CRC ;   % curr.time as unique marker

%-- Call PLM_RFIELDS to calculate the receptive fields
N_orient = Iparams.N_orient ;
N_spfreq = Iparams.N_spfreq ;
if (fid)  fprintf(fid,...
   '\nFront end: %d orientations by %d spatial frequencies.\n',...
   N_orient,N_spfreq) ;  fprintf(fid,...
   'Orientation HWHH: %s deg.\n',mat2str1(Iparams.HW_orient)) ;
   fprintf(fid,'Sp.freq. FWHH: %s octaves.\n',mat2str1(Iparams.octaves));
end

RFields = PLM_RFields(Iparams) ;
if (size(RFields,1)~=N_orient)
   error('IPARAMS.N_orient does not equal length(IPARAMS.D_orient).') ;
end
if (size(RFields,2)~=N_spfreq)
   error('IPARAMS.N_spfreq does not equal length(IPARAMS.D_spfreq).') ;
end


%-- Prepare for zero-padding to avoid wrap-around pollution by FFT
img_size = Iparams.img_size ;
FFT_size = Iparams.FFT_size ;
pad_size = FFT_size - img_size ;
if (any(pad_size<0))
   error('Iparams.FFT_size is less than Iparams.img_size.') ;
end
pad_XY = floor(pad_size./2) ;    % typically [33 33]
Idef.pad_XY = pad_XY ;
padded_RF = zeros(FFT_size) ;    % RField background is always 0

%-- Take the Fourier transform of each RF and store it for filtering
FRONTEND_MEMOIZATION.RFcosFFT = cell(N_orient,N_spfreq) ;
FRONTEND_MEMOIZATION.RFsinFFT = cell(N_orient,N_spfreq) ;
for k = 1:N_orient * N_spfreq
   %- Cosine-phase filter
   G = RFields(k).Wcos ;      % Gabor weights
   padded_RF([1:img_size(1)]+pad_XY(1),[1:img_size(2)]+pad_XY(2)) = G ;
   FRONTEND_MEMOIZATION.RFcosFFT{k} = fft2(fftshift(padded_RF)) ;
   %- Sine-phase filter
   G = RFields(k).Wsin ;      % Gabor weights
   padded_RF([1:img_size(1)]+pad_XY(1),[1:img_size(2)]+pad_XY(2)) = G ;
   FRONTEND_MEMOIZATION.RFsinFFT{k} = fft2(fftshift(padded_RF)) ;
end

%-- Report contrast gain control parameters
if (fid)  fprintf(fid,...
   'Contrast gain control via divisive normalization: %s.\n',...
   Iparams.CGC_mode) ;
end

%-- Calculate weights of the kernels that pool over XY space
sigma_integr = Iparams.FW_integr * Iparams.W2sigma ;
rep_scale = Iparams.rep_scale ;    % scale factor for orientation kernel
centers = Iparams.centers ;        % of integration basins [deg.vis.a]
N_centers = size(centers,1) ;

Xmin = -floor(img_size(2)/2) ; Ymin = -floor(img_size(1)/2) ; %! XY-->JI
X0 = ([1:img_size(2)]+Xmin).* Iparams.deg_pix ;    % pixels-->deg.vis.ang
Y0 = ([1:img_size(1)]+Ymin).* Iparams.deg_pix ;    % Note: Y=first index

for k = 1:N_centers
   Xc = centers(k,1) ; Yc = centers(k,2) ;    % Note: Y=second index
   [X,Y] = meshgrid(X0-Xc,Y0-Yc) ;
   R2 = X.^2 + Y.^2 ;
   % radially Gaussian weights adding up to 1.0*rep_scale
   Gausswgt = exp(R2./(-2*sigma_integr.^2)) ;
   Gausswgt = Gausswgt .* (rep_scale/sum(Gausswgt(:))) ;
   FRONTEND_MEMOIZATION.Gausswgt(k,:) = Gausswgt(:)' ;  % [1x4096]
end

if (fid)  fprintf(fid,...
   'Radial Gaussian integration weights: FW=%.1f deg.vis.ang.\n',...
   Iparams.FW_integr) ; fprintf(fid,'Integration centers:') ;
   fprintf(fid,'  (%.1f,%.1f)',centers) ; fprintf(fid,' deg.vis.ang.\n');
end

%-- Remind the user about the cache
if (fid)
   fprintf(fid,'\nGlobal variable FRONTEND_MEMOIZATION updated.\n\n') ;
end

%--- Return IDEF
%%%%% End of file MAKE_PLFRONTEND.M
