% 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 100).
% 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
%
% Be wary of using too many frequency bins. With more bins there will be
% fewer elements per bin, or even none, producing a very noisy result.
% Copyright (c) 2001-2003 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
%
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
%
% The above copyright notice and this permission notice shall be included in
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.
% July 2001 - original version.
% May 2003 - correction of frequency origin for even and odd size images.
% - corrections to allow for non-square images.
% - extra commenting
function [amp, f, slope] = imspect(im, nbins, lowcut)
if nargin < 2
nbins = 100; % Default No of frequency 'bins'
end
if nargin < 3
lowcut = 2; % Ignore lower 2% of curve when finding
end % line of best fit.
amp = zeros(1, nbins); % preallocate some memory
fcount = ones(1, nbins);
mag = fftshift(abs(fft2(double(im)))); % Amplitude spectrum
% Generate a matrix 'radius' every element of which has a value
% given by its distance from the centre. This is used to index
% the frequency values in the spectrum.
[rows, cols] = size(im);
[x,y] = meshgrid([1:cols],[1:rows]);
% The following fiddles the origin to the correct position
% depending on whether we have and even or odd size.
% In addition the values of x and y are normalised to +- 0.5
if mod(cols,2) == 0
x = (x-cols/2-1)/cols;
else
x = (x-(cols+1)/2)/(cols-1);
end
if mod(rows,2) == 0
y = (y-rows/2-1)/rows;
else
y = (y-(rows+1)/2)/(rows-1);
end
radius = sqrt(x.^2 + y.^2);
% Quantise radius to the desired number of frequency bins
radius = round(radius/max(max(radius))*(nbins-1)) + 1;
% Now go through the spectrum and build the histogram of amplitudes
% vs frequences.
for r = 1:rows
for c = 1:cols
ind = radius(r,c);
amp(ind) = amp(ind)+mag(r,c);
fcount(ind) = fcount(ind)+1;
end
end
% Average the amplitude at each frequency bin. We also add 'eps'
% to avoid potential problems later on in taking logs.
amp = amp./fcount + eps;
% Generate corrected frequency scale for each quantised frequency bin.
% Note that the maximum frequency is sqrt(0.5)
f = [1:nbins]/nbins*sqrt(.5);
% Plots
figure(1), clf
plot(f,amp), xlabel('frequency'), ylabel('amplitude');
title('Histogram of amplitude vs frequency');
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('Histogram of log amplitude vs log frequency');
slope = p(1);