% Example of 2D parametric spectral analysis using autoregressive process.
%
% Uses a random image from a specified directory (with little error
% checking) to demonstrate autoregressive spectral estimation in 2D.  A 2D
% AR model of the specified order is fitted to the data, and it is shown
% how the model can be used for synthesis and spectrum estimation.
%
% Notes:
% 1. Matlab ordering of elements is used throughout - it is admittedly
%    complicated to keep corresponding elements in the right places
% 2. The AR model used here is not separable, so it should be able to
%    capture diagonal structure

% Load random image from texture directory
texdir = 'brodatz';
d = dir(sprintf('%s/*.gif',texdir));
ii = ceil(length(d)*rand);  if ii==0, ii=1; end
im = imread(sprintf('%s/%s',texdir,d(ii).name));
im = double(im(1:100,1:100));  im = im - mean(im(:));

% Get samples for linear prediction.  To understand the function im2dat
% used below, look carefully at the inputs and outputs of the commands
%   im = reshape(1:25,5,5);  im2col(im,[3 4]);
% Getting samples into the form of a data matrix makes it easy to do linear
% algebra on the problem (which matlab does well)
bs = [3 4];  % model order (vertical and horizontal) <<<=== big effect?
X0 = im2col(im,bs);  % data matrix
X = X0(1:end-1,:)';  y = X0(e nd,:)';  % observation vectors and target values

% Solve for model
a = X\y;  % AR parameter estimates
w = y-X*a;  % innovation sequence
sigsq = var(w);  % sample variance of (white?) innovation sequence
wi = reshape(w,size(im)-bs+1);  % driving image
fh = figure;  imagesc(wi);  colormap(gray);  
title('Estimated white noise image from training data (should be white)');

% Create synthetic image, assuming zero valued samples for initialisation.  
% Matlab's filter2 function only does FIR, and this process is hard to
% vectorise, so the implementation is slow (try to never iterate over
% pixels in matlab)
ak = reshape([a; 0],bs);  % reorder AR parameters to filter kernel
wim = sqrt(sigsq)*randn(100,100);  % white noise image for synthesis
sim = zeros(100,100);  % initially blank synthetic image
for r=bs(1):size(sim,1)
  for c=bs(2):size(sim,2)
    % Generate output pixel (r,c) from model and innovation
    sim(r,c) = sum(sum(ak.*sim(r-bs(1)+1:r,c-bs(2)+1:c))) + wim(r,c);
  end
end
figure;  subplot(1,2,1);  imagesc(im);  axis equal;  axis tight;  colormap(gray);  title('Input data');
subplot(1,2,2);  imagesc(sim);  axis equal;  axis tight;  colormap(gray);
title(sprintf('Synthetic data: bs=[%i %i]\n',bs(1),bs(2)));

% Spectrum estimate (is the DC gain correct?)
h = -ak;  h(end,end) = 1;
[Hi,Fx,Fy] = freqz2(h,50,50);  % frequency response of inverse filter
H = sqrt(sigsq)./Hi;  % actual frequency response
fh = figure;  imagesc(Fx,Fy,abs(H));  colormap(gray);
xlabel('\omega_1');  ylabel('\omega_2');
title('Magnitude response of estimated AR filter');
