% Simple classification example for 1D textures
%
% If we fit an autoregressive model to data, then the parameters of the
% model provide information about the data.  Autoregressive models are good
% for describing textures.  This example takes a region from five different
% texture images, and estimates the AR model for each row.  A scatterplot
% of the resulting model coefficients shows that they do tend to cluster
% within a class, so in principle they can be used for classification.
%
% Notes:
% 1. The scatterplot is shown for a second-order AR model, although models
%    with higher order may be more appropriate
% 2. The variance or standard deviation of the driving noise process is
%    actually also a parameter that might describe the texture, but it is 
%    ignored in the plot
% 3. The AR coefficients themselves might actually be bad classification
%    features if a nearest neighbour (Euclidean) measure is used.  AR
%    modeling literature often refers to reflection coefficients, which
%    might be better.  The code below also calculates the estimated
%    spectra, and relevant features could be derived directly from these
% 4. The methods described extend quite easily to two dimensions, in which
%    case one would expect better classification ability since the model
%    also incorporates vertical (and diagonal) structure.

% Load data from five images in texture directory
texdir = 'brodatz';
d = dir(sprintf('%s/*.gif',texdir));
if length(d)<5, error('Not enough images in texture directory'); end
imc = cell(5,1);  sigc = cell(5,1);
for i=1:length(sigc)
  im = imread(sprintf('%s/%s',texdir,d(i).name));
  im = im(1:100,1:100);
  imc{i} = double(im);  % store entire image
  sigc{i} = double(im(1,:));  % store first rows
end

% Fit separate model to each row of each image.  Put results in a cell
% array of structures
p = 2;  % two parameter model to visualise
res = cell(1,length(imc));  % cell array for results (one element per image)
for i=1:length(res)
  im = imc{i};  % current image to analyse
  
  % Storage for each row
  am = zeros(p,size(im,1));  % each column has AR coefficients for one row
  em = zeros(1,size(im,1));  % driving noise variance for each row
  wv = linspace(0,2*pi,200)';
  Hm = zeros(length(wv),size(im,1));  % transfer function values for each row
  
  % Populate storage for current image
  for j=1:size(im,1)
    dat = im(j,:);  dat = dat - mean(dat);
    X0 = flipud(im2col(dat,[1 p+1]));  % data matrix
    X = X0(2:end,:)';  y = X0(1,:)';  % observation vectors and target values
    a = X\y;  % autoregressive coefficients
    w = y-X*a;  % driving noise sequence
    sigsq = var(w);  % driving noise variance
  
    % Store results for row
    am(:,j) = a;
    em(j) = var(w);
    Hm(:,j) = freqz(sqrt(sigsq),[1 -a'],wv);
  end
  
  % Store results for image
  res{i} = struct('am',am,'em',em,'Hm',Hm);
end

% Show scatterplot of AR coefficients for all rows using different colour
% for each image
cspec = {'r','g','b','c','m'}; 
fh = figure;
for i=1:length(res)
  am = res{i}.am;
  sh = scatter(am(1,:),am(2,:),cspec{i});  hold on;
end
hold off;  xlabel('a_1');  ylabel('a_2');
title('Scatterplot of AR coefficients');

