function [h p]=hotelT2(v1,v2)
%
% [h p] =hotelT2(v1,v2)
%
% Computes the Hotelling's T-square statistic for given vector measures v1 and v2
% v1,v2    :  vector measures  
%                For instance if v1 is the size of 2562 x 3 x 47, 
%                v1(:,1,12) is the x-coordinate of the 12-th subject.
%                v1 and v2 do not have to have the same number of subjects.
%
%  h : Hotelling's T2 statistic
%  p : The corresponding p-value
%
% (C) Moo K. Chung 2007-2008
%     Department of Biostatistics and Medical Informatics
%     Waisman Laboratory for Brain Imaging
%     University of Wisconsin-Maison
%  
% email://mkchung@wisc.edu

p=size(v1,2);  % dimension 

n_vertex=size(v1,1);  % # of vertices
n1=size(v1,3); % sample size in group 1
n2=size(v2,3); % sample size in group 2

for i=1:n_vertex
    w1=squeeze(v1(i,:,:));
    w2=squeeze(v2(i,:,:));
    m1=mean(w1,2);
    m2=mean(w2,2);
    
    z1=w1-kron(ones(1,n1),m1);
    z2=w2-kron(ones(1,n2),m2);    
    sigma = z1*z1' +z2*z2';
    
    h(i)=(m2-m1)'*inv(sigma)*(m2-m1);
    
end;

const=n1*n2*(n1+n2-p-1)/(p*(n1+n2));
h=const*h;
p=1-fcdf(h, p, n1+n2-p-1 );

