% [X,sqd,N,S,r,V,H] = imgsqd(imgfile[,S,r])
% Compute feature vector and pairwise square distances between image pixels
%
% "imgfile" contains an image of N = VxH pixels with intensities or RGB
% colours in 0..255. imgsqd computes and optionally saves the following
% information in a file "imgfile".mat:
%
% - X: feature vectors for each pixel. For greyscale images, Nx3 array with
%   rowwise pixel features; each row is (x,y,I) = (positionx in [1,H],
%   positiony in [1,V], intensity in [0,1]). Each feature is further divided
%   by (S1,S2,S3) for relative scaling. For colour images, Nx5 array with
%   columns 3-5 being the CIE L*u*v* values (uniform colour space), divided
%   by (S1,...,S5).
%
% - sqd: sparse symmetric NxN matrix containing the Euclidean square
%   distances between pairs of rows. The distances are not computed (and
%   assumend Inf) for rows whose pixels are farther than r pixels (in the
%   position x,y). sqd is computed only if r > 0.
%
% - S, r, V, H.
%
% The greyscale image can be visualised as:
%   figure(1); imagesc(reshape(X(:,3),V,H)); colormap(gray(256));
%   set(gca,'DataAspectRatio',[1 1 1]); xlabel('H'); ylabel('V');
%   figure(2); plot3(X(:,1),X(:,2),X(:,3),'ko'); xlabel('V'); ylabel('H');
%   set(gca,'DataAspectRatio',[1 1 1]); view(60,30);
% The colour image can be visualised as:
%   figure(1); image(reshape(luv2rgb(X(:,3:5)),V,H,3));
%   set(gca,'DataAspectRatio',[1 1 1]); xlabel('H'); ylabel('V');
%
% In:
%   imgfile: image file name.
%   S: 1xK array, the relative scales of the features. Default: [1 1 0.01]
%      for greyscale images and [1 1 1 1 1] for colour images.
%   r: window size (in pixels). If <= 0, then sqd won't be computed.
%      Default: 5.
% Out:
%   X,sqd,N,S,r,V,H: as described above. If no output arguments are given
%      when calling imgsqd, the outputs are saved in a .mat file.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2013 by Miguel A. Carreira-Perpinan and Max Vladymyrov

function [X,sqd,N,S,r,V,H] = imgsqd(imgfile,S,r)

[F,map] = imread(imgfile);
iscolour = (ndims(F) > 2) | (size(map,2) == 3);

% ---------- Argument defaults ----------
if ~exist('S','var') | isempty(S)
  if iscolour
    S = [1 1 1 1 1];					% Colour image
  else
    S = [1 1 0.01];					% Greyscale image
  end
end;
if ~exist('r','var') | isempty(r) r = 5; end;
% ---------- End of "argument defaults" ----------


% Image features
V = size(F,1); H = size(F,2); [FH,FV] = meshgrid(1:H,1:V);
if iscolour				% Colour image
  if ndims(F) > 2			% RGB
    X = [FV(:) FH(:) rgb2luv(double(reshape(F,H*V,3))/255)]*diag(S.^(-1));
  else					% Indexed
    X = [FV(:) FH(:) rgb2luv(map(1+F(:),:))]*diag(S.^(-1));
  end
else					% Greyscale image
  X = [FV(:) FH(:) double(F(:))/255]*diag(S.^(-1));		% Intensity
% $$$   F = rgb2luv(repmat(double(F(:))/255,1,3));
% $$$   X = [FV(:) FH(:) F(:,1)]*diag(S.^(-1));			% L* value
end
N = size(X,1);
if nargout==1 return; end;


% Matrix of square distances
sqd = sparse(N,N);
if r > 0
  % Build list of pixel indices inside the window
  rr = floor(r);
  [wH,wV] = meshgrid(-rr:0,-rr:rr); wH = wH(1:end-rr-1); wV = wV(1:end-rr-1);
  tmp = find(wH.^2+wV.^2 <= (r^2+100*eps)); w = [wV(tmp)' wH(tmp)'];
  % Pairwise distances
  x2 = sum(X.^2,2);
  for j=1:V
    for i=1:H
      tmp1 = w(:,1)+i; tmp2 = w(:,2)+j;
      tmp = tmp1>0 & tmp1<=V & tmp2>0 & tmp2<=H;
      ii = i+(j-1)*H; jj = tmp1(tmp) + (tmp2(tmp)-1)*H;
      sqd(ii,jj) = x2(ii) + x2(jj) - 2*X(jj,:)*X(ii,:)';
    end
  end
  sqd = sqd + sqd';			% Symmetrise
end


% Save results
if nargout==0
  save([imgfile(1:max(findstr(imgfile,'.'))-1) '.mat'],...
       'X','sqd','N','S','r','V','H');
end

