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