corrcoef_directional(t1, t2, doplot)

% [rho,r2,p] = corrcoef_directional(t1, t2 [,doplot])
%  Correlation coefficient for directional data
%
%   t1,t2: The two sets of directional data. t1 and t2 are
%       both Nx2 matrices (the x,y coordinates of N points). 
%   doplot: if 1, then results are plotted (default = 0) on a grid.
%           if grid location is one pair of points: one from t1 (red)
%           and one from t2 (blue).
% 
%   rho is the normalized correlation coefficient (it ranges from 0 to 1, 
%           and has similar properties to Pearson's r).
%   r2 is the unnormalized coefficient (it ranges from 0 to K, 
%           where K = dimensionality of the space)
%   p is the probability that t1 and t2 are drawn from independent 
%           distributions (requires Statistics Toolbox).
%
%  For details, see:
%   Jupp,P.E., & Mardia,K.V. (1980). "A general correlation coefficient for 
%   directional data and related regression problems", Biometrika, 67(1):163-173.
%
%  Example:
%   a = randn(40,2);
%   b = (a(:,1)+i*a(:,2)) * exp(i * 0.3);    % rotate the points
%   b = [real(b) imag(b)] + randn(40,2)*0.25; % and add noise
%   rho = corrcoef_directional(a,b,1)
%
%  Gautam Vallabha, Dec-2002, Florida Atlantic University

function [rho,r2,p] = corrcoef_directional(t1, t2, doplot)

if nargin < 3, doplot = 0; end

% unitize the vector sets. This is equivalent to taking
% angles, and converting to cos,sin
n = size(t1,1);
t1 = t1 ./ repmat(sqrt(sum(t1.^2,2)), 1, 2); 
t2 = t2 ./ repmat(sqrt(sum(t2.^2,2)), 1, 2); 

if doplot,
    figure; hold on; axis equal;
    k = 2*ceil(sqrt(n));
    if doplot==1,
      ang=getangle(t1,t2); 
      [dummy,ii]=sort(abs(ang));
      t11=t1(ii,:); t22=t2(ii,:);
    else
      t11=t1; t22=t2;
    end
    [xx,yy] = meshgrid(1:2:k, 1:2:k);
    pos = [xx(:) yy(:)];
    for i=1:n,
        x = pos(i,1); y = pos(i,2);
        plot([x x+t11(i,1)], [y y+t11(i,2)], 'r');
        plot([x x+t22(i,1)], [y y+t22(i,2)], 'b');
    end
    set(gca,'xtick', 1:2:k, 'ytick', 1:2:k);
    grid on;
end

S = cov([t1 t2]);
s11 = S(1:2,1:2);
s12 = S(1:2,3:4);
s21 = S(3:4,1:2);
s22 = S(3:4,3:4);

r2 = trace(inv(s11) * s12 * inv(s22) * s21);
rs11 = rank(s11); rs22 = rank(s22);

maxr2 = min(rs11,rs22);
rho = sqrt(r2/maxr2); 
% rho is a 0..1 measure that has similar properties to Pearson's r.

if exist('chi2pdf','file')
  df = rs11 * rs22;
  p = 1 - chi2cdf(n*r2,df);
else
  p = nan;
end

%-----------
% gets the angle of vb RELATIVE to va
function ang = getangle(vva,vvb)
% assume vectors are normalized
ang = zeros(size(vva,1));
for i=1:size(vva,1),
  va=vva(i,:); vb=vvb(i,:);
  tc = acos(sum(va .* vb)); % theta, via cosine
  ts = sign(cross([va 0], [vb 0])); % sign of sin(x)
  ang(i) = ts(3) * abs(tc);
end