function  [N,nE] = PLM_gain_control(E,Idef)

%PLM_GAIN_CONTROL -- Contrast gain control via divisive normalization.
%
%  N = PLM_gain_control(E,Idef)
%  [N,nE] = PLM_gain_control(E,Idef)
%
%  E is a 4-D array of local stimulus energies made by PLM_STIM_ENERGY.
%  IDEF is a "front-end definition structure" made by MAKE_PLFRONTEND.
%  N is a 3-D normalization matrix of size [img_size,N_spfreq]. The
%  orientation dimension (the third dimension) is omitted because
%  the normalization pools across all orientations. The normalization
%  is divisive and local in space. It is done independently for each
%  location [x,y] according to the following algorithm:
%   1. Within a given frequency band F (and for each [x,y] location),
%      sum the energies across all orientations.
%   2. For each frequency F, pool across nearby frequency bands as
%      prescribed by the weight matrix IDEF.CGC_Fpool. The result forms
%      the "normalization signal" N(:,:,f) for that frequency band.
%   3. Add the constant IDEF.CGC_const to all normalization signals.
%      The resulting 3-D matrix N is the main output of the function.
%   4. Optionally, PLM_GAIN_CONTROL can also return the normalized-
%      energy matrix NE. It is the same size as the raw energy E:
%      nE = E ./ repmat(N,<each_orientation>)
%
%  Reference: Heeger, D.J. (1992). Normalization of cell responses in
%   cat striate cortex, Visual Neuroscience, 9, 181-197. [Equation 4]
%
%  See also PLM_STIM_ENERGY, PLFRONTEND.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 2.0 $  $Date: 2002/06/20 16:00 $

%-- Decipher IDEF
sigmasq = Idef.CGC_const ;
Fpool = Idef.CGC_Fpool ;

sz = size(E) ;     % typically [64,64,37,4]
img_size = sz([1,2]) ;
N_orient = sz(3) ;
N_spfreq = sz(4) ;

%-- Pool across all orientations
E_sumOrient = reshape(sum(E,3),[img_size,N_spfreq]) ;

%-- Pool across neighboring frequency bands
E_sumFreq = zeros([img_size,N_spfreq]) ;
total_sum = zeros(1,N_spfreq) ;
for k = 1:N_spfreq
   W = repmat(reshape(Fpool(k,:),[1,1,N_spfreq]),[img_size,1]) ;
   %E_sumFreq(:,:,k) = sum(E_sumOrient.*W,3) ;  % W1*f1+W2*F2+...+Wn*Fn
   foo = sum(E_sumOrient.*W,3) ;
   total_sum(k) = mean(foo(:)) ;
end

%-- Add the semisaturation constant and return N
%N = E_sumFreq + sigmasq ;
N = reshape(repmat(total_sum,64*64,1),[img_size,N_spfreq]) ;

%-- If normalized energy is also requested, replicate N and divide
if (nargout==2)
   nE = reshape(N,[img_size,1,N_spfreq]) ;      % [64x64x 1x4]
   nE = E ./ repmat(nE,[1,1,N_orient,1]) ;      % [64x64x37x4]
end

%--- Return N and possibly NE
%%%%% End of file PLM_GAIN_CONTROL.M
