function [X Qa Qb L] =crossmodularity(A,B,alpha,beta,T)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Given two input (symmetrical) matrices A and B, this function
    % calculates the crossmodularity index X
    % Input's description:
    % A and B are squared matrices of N*N (typically connectivity matrices), being N the number of ROIs
    % alpha and beta are two artibitrary thersholds to binarize the two  matrices (necessary for the similarity calculation)
    % T is the label vector: each element vector is defined as an integer corresponding to the module that ROI belongs to
    % Output's description:
    % X crossmodularity
    % Qa and Qb are modularities of inA and inB associatted to partition T
    % L is the similarity between A and B
    % Ibai Diez. Biocruces Health Research Institute, Spain
    % Please, cite this reference if you use it
    % Diez, I. et al. A novel brain partition highlights the modular skeleton shared by structure and function.
    % Sci. Rep. 5, 10532; doi: 10.1038/srep10532 (2015).
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Get the different labels of the modules
    labels=unique(T);
    
    %For each module compute sorensen index
    sorensen=zeros(1,length(labels));
    indx_m=[];
    for m=1:length(labels)
     %Select the rois of each module and binarizes the resulting matrices using alpha and betha
     indx=find(T==labels(m));
     bin_A=A(indx,indx)>alpha;
     bin_B=B(indx,indx)>beta;
     if length(indx)==1
       continue;
     end
     sorensen(m)=sum(2*(bin_A(:).*bin_B(:)))/(sum(bin_A(:))+sum(bin_B(:)));
     indx_m=[indx_m m];
    end
    %The total similarity is the mean similarity of all the modules
    L=mean(sorensen(indx_m));
    
    
    %Compute the modularity index
    Qa=modularity_index(abs(A),T);
    Qb=modularity_index(abs(B),T);
    
    %Compute the cross modularity
    X=(Qa.*Qb.*L).^(1/3);
end