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


n1=size(mu1,1); % sample size in group 1
n2=size(mu2,1); % sample size in group 2

if size(mu1,3)>1 %3D matrix
    n_vertex=size(mu1,2);  % # of vertices
    d=size(mu1,3);  %dimension
    
    
    for i=1:n_vertex
        
        w1=squeeze(mu1(:,i,:));  %n1 x 3
        w2=squeeze(mu2(:,i,:)); %n2 x 3
        m1=mean(w1,1);   %1x3
        m2=mean(w2,1); %1 x3
        
        z1=w1-kron(ones(n1,1),m1);
        z2=w2-kron(ones(n2,1),m2);
        sigma = z1'*z1 +z2'*z2;
        
        s=size(sigma,1);
        sigmap = sigma + 0.00001*eye(s); %this removes possible singuarity.
        hotel(i)=(m2-m1)*pinv(sigma)*(m2-m1)';
        
    end;
    
else
    n_vertex=1;
    d=size(mu1,2);
    
 
    m1=mean(mu1,1);   %1x2
    m2=mean(mu2,1); %1 x2
    
    z1=mu1-kron(ones(n1,1),m1);
    z2=mu2-kron(ones(n2,1),m2);
    sigma = z1'*z1 +z2'*z2;
    
    s=size(sigma,1);
    sigmap = sigma + 0.00001*eye(s); %this removes possible singuarity.
    hotel=(m2-m1)*inv(sigma)*(m2-m1)';
end;

const=n1*n2*(n1+n2-d-1)/(d*(n1+n2));
h.t=const*hotel;
h.pvalue=1-fcdf(h.t, d, n1+n2-d-1 );

