function [timeDelay_per_group] = time_delay_groupsPairs(timeDelayMap, computedROIsGroups, sort_idxs)
% Average time delay maps over groups of emitter and receiver ROIs
groupsName = {'Visual','Parietal','Cingulate','Frontal','Temporal','Motor'};
computedROIsGroupsLab = {groupsName{computedROIsGroups}};

allGroups = unique(computedROIsGroupsLab(sort_idxs),'stable'); % groups labels ordered according to hierarchy

tPoints = size(timeDelayMap,3);
dPoints = size(timeDelayMap,4);

% All to all ROIs
groupPairLab = 'all_all';
all_idxs = logical(ones(1,numel(sort_idxs)));
tmpMap = timeDelayMap(all_idxs,all_idxs,:,:);
             
nTmpPairs = (sum(all_idxs)*(sum(all_idxs)-1))/2;
timeDelay_per_group.(groupPairLab).forward = nan(nTmpPairs,tPoints,dPoints);
timeDelay_per_group.(groupPairLab).feedback = nan(nTmpPairs,tPoints,dPoints);

pairCount = 1;
for relatEmitIdx = 1:sum(all_idxs)
    for relatRecIdx = relatEmitIdx+1:sum(all_idxs)
        timeDelay_per_group.(groupPairLab).forward(pairCount,:,:) = tmpMap(relatRecIdx,relatEmitIdx,:,:);
        timeDelay_per_group.(groupPairLab).feedback(pairCount,:,:) = tmpMap(relatEmitIdx,relatRecIdx,:,:);
        pairCount = pairCount+1;
    end
end
timeDelay_per_group.(groupPairLab).forward = squeeze(mean(timeDelay_per_group.(groupPairLab).forward,1));
timeDelay_per_group.(groupPairLab).feedback = squeeze(mean(timeDelay_per_group.(groupPairLab).feedback,1));

% groups pairs
for geIdx = 1:numel(allGroups) % emitter group
    emit_idxs = strcmp(computedROIsGroupsLab(sort_idxs),allGroups{geIdx});
    for grIdx = geIdx:numel(allGroups) % receiver group
        groupPairLab = [allGroups{geIdx},'_',allGroups{grIdx}];
        rec_idxs = strcmp(computedROIsGroupsLab(sort_idxs),allGroups{grIdx});

        all_idxs = emit_idxs | rec_idxs;
        % Squared array with first 2 dimensions equal to communication only
        % within the selected subgroups pair
        tmpMap = timeDelayMap(all_idxs,all_idxs,:,:);
             
        nTmpPairs = (sum(all_idxs)*(sum(all_idxs)-1))/2;
        timeDelay_per_group.(groupPairLab).forward = nan(nTmpPairs,tPoints,dPoints);
        timeDelay_per_group.(groupPairLab).feedback = nan(nTmpPairs,tPoints,dPoints);

        pairCount = 1;
        for relatEmitIdx = 1:sum(all_idxs)
            for relatRecIdx = relatEmitIdx+1:sum(all_idxs)
                timeDelay_per_group.(groupPairLab).forward(pairCount,:,:) = tmpMap(relatRecIdx,relatEmitIdx,:,:);
                timeDelay_per_group.(groupPairLab).feedback(pairCount,:,:) = tmpMap(relatEmitIdx,relatRecIdx,:,:);
                pairCount = pairCount+1;
            end
        end
        timeDelay_per_group.(groupPairLab).forward = squeeze(mean(timeDelay_per_group.(groupPairLab).forward,1));
        timeDelay_per_group.(groupPairLab).feedback = squeeze(mean(timeDelay_per_group.(groupPairLab).feedback,1));
    end
end

end