function  [E,C,S] = PLM_stim_energy1(stim,Idef)

%PLM_STIM_ENERGY -- Local energy at various orientations and frequencies.
%
%  E = PLM_stim_energy(stim,Idef)
%
%  STIM is a grayscale image (in IJ coordinates when IDEF.flipud=1=deflt).
%  IDEF is a "front-end definition structure" produced by MAKE_PLFRONTEND.
%  E is a 3-D array of size [img_size,N_orient,N_spfreq], eg [64,64,7,5].
%  Each page E(:,:,o,f) is a map of the local stimulus energy at
%  orientation IDEF.D_orient(O) and spatial frequency IDEF.D_spfreq(F).
%  The local energy is calculated in the following way:
%   1. If IDEF.flipud=1 (default), flip the stimulus image upside down.
%      This is typically necessary because MAKE_PLE_STIM uses IJ-style
%      coordinates while all PLModel functions (and PLM_RFields in
%      partucular) use XY-style coordinates. That is, [1,1] denotes the
%      bottom left-hand corner of STIM in PLM_STIM_ENERGY and the top
%      left-hand corner in MAKE_PLE_STIM. See AXIS('XY') vs. AXIS('IJ').
%      Note also that X grows along the 2nd idx of STIM (as in MESHGRID).
%   2. Pad the (flipped) image up to IDEF.FFT_size (see IDEF.bkgr_idx).
%   3. For each orientation 0 and frequency F, convolve the zero-padded
%      image with two Gabor filters at cosine and sine phases.
%   4. Energy_map := cos_phase_map.^2 + sin_phase_map.^2
%   5. Crop ENERGY_MAP back to IDEF.img_size, discarding the padding.
%
%  Reference: Heeger, D.J. (1992). Half-squaring in responses of cat
%    striate cells. Visual Neuroscience, 9, 427-443.
%
%  See also PLFRONTEND, PLM_CONVOLVE, PLM_RFIELDS.

%  Original coding by Alexander Petrov, UC Irvine.
%  $Revision: 2.1 $  $Date: 2002/07/12 20:20 $

global FRONTEND_MEMOIZATION

%-- Make sure that the memoization store is up to date
if (FRONTEND_MEMOIZATION.CRC ~= Idef.memo_CRC)
   error('FRONTEND_MEMOIZATION is obsolete. Call MAKE_PLFRONTEND.') ;
end

%-- Decipher IDEF
N_orient = Idef.N_orient ;
N_spfreq = Idef.N_spfreq ;
img_size = Idef.img_size ;
FFT_size = Idef.FFT_size ;
pad_XY = Idef.pad_XY ;

%-- Normalize the stimulus
s = (stim-Idef.LUT_gray) ./ Idef.LUT_span ;    % now in [-.5 +.5] range
background_color = s(Idef.bkgr_idx) ;
if (Idef.flipud)       % STIM is in IJ- instead of XY-coordinates?
   s = flipud(s) ;     % convert to XY-coors: [1 1]=bottom left
end

%-- Pad with BACKGROUND_COLOR and take the Fourier transform
padded_stim = repmat(background_color,FFT_size) ;
Y_range = [1:img_size(1)]+pad_XY(1) ;      % NB: the first index is Y
X_range = [1:img_size(2)]+pad_XY(2) ;
padded_stim(Y_range,X_range) = s ;
stim_FFT = fft2(padded_stim) ;

%-- Calculate energy using the filters cached in FRONTEND_MEMOIZATION
E = zeros([img_size,N_orient,N_spfreq]) ;   % [[64,64],7,5]
C = zeros(size(E)); S = zeros(size(E)) ;
for k = 1:(N_orient*N_spfreq)     % for each orientation and frequency
   [o,f] = ind2sub([N_orient,N_spfreq],k) ;
   %- convolve with two Gabors: cos-phase and sin-phase
   H = stim_FFT .* FRONTEND_MEMOIZATION.RFcosFFT{o,f} ;
   Ccos = real(ifft2(H)) ;
   H = stim_FFT .* FRONTEND_MEMOIZATION.RFsinFFT{o,f} ;
   Csin = real(ifft2(H)) ;
   %- calculate energy and trim back to IMG_SIZE
   padded_E = Ccos.^2 + Csin.^2 ;
   E(:,:,o,f) = padded_E(Y_range,X_range) ;
   C(:,:,o,f) = Ccos(Y_range,X_range) ;
   S(:,:,o,f) = Csin(Y_range,X_range) ;
end

%--- Return E
%%%%% End of file PLM_STIM_ENERGY.M
