imspect.m

MATLAB Functions for Multiple View Geometry

% 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);