% 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 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.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2005 by Miguel A. Carreira-Perpinan

function imgsqd(imgfile,S,r)

[F,map] = imread(imgfile);
iscolour = (ndims(F) > 2) | (size(map,2) == 3);
isimg = 1;				% Image flag

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


% 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
  for i=2:N
    % Select window indices inside the image
    [tmp1 tmp2] = ind2sub([V H],i); tmp1 = w(:,1)+tmp1; tmp2 = w(:,2)+tmp2;
    tmp = tmp1>0 & tmp1<=V & tmp2>0 & tmp2<=H;
    for j=sub2ind([V H],tmp1(tmp),tmp2(tmp))'
      sqd(i,j) = sum((X(i,:)-X(j,:)).^2);
    end
  end
  sqd = sqd + sqd';			% Symmetrise
end


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