% IMSPECT - Plots image amplitude spectrum averaged over all orientations.
%
% Usage:  [amp, f, slope] = imspect(im, nbins, lowcut)
%                                          \     /
%                                          optional
% Arguments:
%               im     - Image to be analysed.
%               nbins  - No of frequency bins to use (defaults to 400).
%               lowcut - Percentage of lower frequencies to ignore when
%                        finding line of best fit. This avoids problems
%                        with amplitude spikes at low frequencies. 
%                        (defaults to 2%) .
%
% Returns:
%               amp    - 1D array of amplitudes.
%               f      - Corresponding array of frequencies.
%               slope  - Slope of line of best fit to the log-log data

% Peter Kovesi  
% Department of Computer Science & Software Engineering
% The University of Western Australia
% pk@cs.uwa.edu.au    www.cs.uwa.edu.au/~pk
%
% July 2001

function [amp, f, slope] = imspect(im, nbins, lowcut)
    
    if nargin < 2
	nbins = 400;            % Default No of frequency 'bins'
    end
    if nargin < 3
	lowcut = 2;             % Ignore lower 2% of curve when finding
                                % line of best fit.
    end

    amp = zeros(1, nbins);
    fcount = ones(1, nbins);
    
    mag = fftshift(abs(fft2(im)));

    [rows, cols] = size(im);    
    [x,y] = meshgrid([1:cols],[1:rows]);
    x = x-cols/2;
    y = y-rows/2;
    freq = sqrt(x.^2 + y.^2);
    freq = round(freq/max(max(freq))*(nbins-1)) + 1;
    
    for r = 1:rows
	for c = 1:cols
	    ind = freq(r,c); 
	    amp(ind) = amp(ind)+mag(r,c);
	    fcount(ind) = fcount(ind)+1;
	end
    end
    
    amp = amp./fcount + eps;             % Add eps to avoid potential
    f = [1:nbins]/nbins*sqrt(.5) + eps;  % problems in taking logs.

    figure(2), clf
    loglog(f,amp), xlabel('log frequency'), ylabel('log amplitude');
    
    % Find line of best fit (ignoring specified fraction of low frequency values) 
    fst = round(nbins*lowcut/100);
    p = polyfit(log(f(fst:end)), log(amp(fst:end)),1);
    y = exp(p(1)*log(f) + p(2));        % log(y) = p(1)*log(f) + p(2)
    hold on, loglog(f, y,'Color',[1 0 0])

    n = round(nbins/10);
    text(f(n), amp(n), sprintf('Slope = %f.2',p(1)));
    title('Log amplitude vs log frequency');

    slope = p(1);
