function  canvas = spectral_diag(img,Fparams)

%SPECTRAL_DIAG -- Diagram of a Fourier power spectrum.
%
%  spectral_diag(img)
%  canvas = spectral_diag(img,Fparams)
%
%  IMG is the image whose spectrum is being diagrammed.
%  FPARAMS is a parameter structure produced by SPECTRAL_DIAG_PARAMS.
%  If not provided, FPARAMS = spectral_diag_params(size(img)).
%  Note that all spectral frequencies are plotted in log2 scale.
%
% Example:
%  P=PLExp2_params; F=spectral_diag_params(P.size);
%  spectral_diag(P.Gabor(:,:,1,1),F);   % colormap(flipud(gray))
%
%  See also SPECTRAL_DIAG_PARAMS, FFT2, HINTON_DIAG.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2004/12/02 14:55 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

if (nargin==1)  Fparams = spectral_diag_params(sizeof(img)) ; end

sz = Fparams.size ;   % [64 64] or [256 256]
if (any(size(img)~=sz))  error('Incompatible image size.') ; end

%-- Main job: take the Fourier power spectrum
sp = abs(fftshift(fft2(img))) ;    % [64 64] or [256 256]
max_power = Fparams.max_power ;

%- The spectrum is mirror symmetric. Discard the half with Fx<0:
sp(:,1:sz(1)/2) = [] ;    % [64 32] or [256 128]

%- SP is rotated 90* relative to IMG. Rotate it back:
sp = flipud(sp') ;        % [32 64] or [128 256]

%-- Take square root if specified
if (Fparams.sqrtp)
   sp = sqrt(sp) ;
end

%-- Reconstruct the FFT mesh
f_horiz = Fparams.f_horiz ;           % [ 1 64] or [1 256]
f_vert = Fparams.f_vert ;             % [32  1] or [128 1]
[Fy,Fx] = meshgrid(f_horiz,f_vert) ;  % [32 64] or [128 256]

%-- Transformed from Cartesian to polar coordinates:
THi = Fparams.THi ;  % [1 121]
Ri = Fparams.Ri ;    % [50 1]
psp = cart2pol2(Fy,Fx,sp,THi,Ri,'*cubic') ;
%imagesc(90-THi.*(180/pi),Ri,psp,[0 100]) ; axis xy ;

%-- Log2-transform the R axis. Arbitrarily map R=1cyc/deg to LOG_ORIGIN.
lRi = log2(Ri) + Fparams.log_origin ;

%-- (Part of) the log2-transformed spectrum is then projected
%   back onto the canvas, in Cartesian coordinates.
u = Fparams.u ;          % [1 337]
v = Fparams.v ;          % [169 1]
[U,V] = meshgrid(u,v) ;  % [169 337]
canvas = pol2cart2(THi,lRi,psp,u,v,'cubic') ;   % [169 337]
canvas = max(0,min(max_power,canvas)) ;
%imagesc(u,v,canvas,[0 max_power]) ; axis image xy ;

%-- Eliminate the very low frequencies (they are spurious) and
% the very high frequencies, restricting the canvas to an annulus:
R2 = U.^2 + V.^2 ;
canvas(find(R2<Fparams.ann_rng(1)^2)) = -1 ;   % < R_min
canvas(find(R2>Fparams.ann_rng(2)^2)) = -1 ;   % > R_max
%imagesc(u,v,canvas,[-1 max_power]) ; axis image xy ;

%-- Generate a cobweb
spokes = Fparams.spokes ;
rings = Fparams.rings ;
cob = cobweb(U,V,spokes,rings,Fparams.linewidth,Fparams.cob_rng) ;

%-- Superimpose the cobweb on the spectrum
cob_color = max_power + 1 ;   % 51 under current defaults
canvas(find(cob==0)) = cob_color ;
%imagesc(u,v,canvas,[0 cob_color]) ; axis image xy ;

%-- Finally, produce the plot and return
imagesc(u,v,canvas) ; axis image xy ;
set(gca,'ytick',[],'box','off') ;   % meaningless in polar coordinates
set(gca,'xtick',Fparams.xtick,'xticklabel',Fparams.xticklabel) ;
xlabel(Fparams.xlabel) ;
colormap(Fparams.colormap) ;

%--- Return CANVAS
%%%%% End of file HINTON_DIAG.M
