%% energyhist.m
% 2D energy histograms
% compute and display distribution of energy by radial frequency

%  S. Allen Broughton - 1 Oct 10

%% get and display image and shifted DFT
 load mandrill
% load earth
% load earth; X = X-mean(X(:)); % eliminate the DC coefficient 
%  X = zeros(128,128); X(61:80,41:80)=1; % rectangle
fX = fft2(X); 
sfX = fftshift(fX);

figure(1)
subplot(1,2,1)
imagesc(X), colormap(gray), axis equal, axis tight 
title('original image');
subplot(1,2,2)
imagesc(log(1+abs(sfX))), colormap(gray), axis equal, axis tight 
title('shifted log(1+|dft|)');

%% get image data and set up shifted frequency matrices
[m,n] = size(X);
N = m*n;
hm = round(m/2);
hn = round(n/2);
maxfreq = sqrt(hm^2+hn^2);

V = ((0:(m-1))-hm)'*ones(1,n);   % vertical frequencies
H = ones(m,1)*((0:(n-1))-hn);    % horizontal frequencies
R = sqrt(V.^2+H.^2);             % radial frequencies

%% display frequency maps
figure(2)
subplot(1,3,1)
imagesc([-hm,hm],[-hm,hm],V), colormap(gray), axis equal, axis tight 
title('vertical frequencies');
subplot(1,3,2)
imagesc([-hm,hm],[-hm,hm],H), colormap(gray), axis equal, axis tight 
title('horizontal frequencies')
subplot(1,3,3)
imagesc([-hm,hm],[-hm,hm],R), colormap(gray), axis equal, axis tight 
title('radial frequencies')


 %% set up freqency data and calculate energy histogram
nbins = 50;
cEhist = [ ]; 
df = maxfreq/nbins;
E = energy(X) % display this to indicate calulation is proceeding

for j = 1: nbins,  
    J = find (R <= j*df);
    Y = zeros(size(sfX)); 
    Y(J) = sfX(J); 
    cEhist = [cEhist, energy(Y)/N];    
end
Ehist = cEhist -[0,cEhist(1:nbins-1)]; % show energies


%% display results 
figure(3)
freqs = (1: nbins)*df;
subplot(2,1,1)
stem(freqs,Ehist),
xlabel('frequency')
ylabel('energy')
title('frequency-energy distribution')
subplot(2,1,2)
plot(freqs,cEhist,'o',freqs,E*ones(size(freqs)),'k-',freqs,0,'k-') % put in zero to include axis
xlabel('frequency')
ylabel('energy')
title('frequency-energy cumulative distribution')