% Calculate the voxel-to-voxel mean iFC (or called Global Correlation, GCOR). We calculated the all voxel-to-voxel mean iFC as follows: 1) normalize the time courses of all the voxels to zero mean and unit variance; 2) calculate the mean signal across the brain (¡°global signal¡±); 3) calculate correlation between this ¡°global signal¡± and all the other voxels (a simple dot product and then divided by n -1); 4) calculate the mean value of the correlation coefficients across brain. This mean correlation coefficient is equivalent to the all voxel-to-voxel mean correlation. This calculation is similar to the L2 norm way recently proposed by Saad et al. (2013), but in a more intuitive form. This script is written for a paper under review.
%__________________________________________________________________________
% Written by YAN Chao-Gan 130515.
% The Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA
% The Phyllis Green and Randolph Cowen Institute for Pediatric Neuroscience, New York University Child Study Center, New York, NY 10016, USA
% ycg.yan@gmail.com
% http://rfmri.org/Scripts_YAN


BrainMaskDir = '/Processing/Masks';

GlobalCorrelation=zeros(AutoDataProcessParameter.SubjectNum,1);

AutoDataProcessParameter.DataProcessDir='/Processing';


AutoDataProcessParameter.StartingDirName = 'FunImgARC';

parfor i=1:AutoDataProcessParameter.SubjectNum
    
    ROIMask_Global={[BrainMaskDir,filesep,AutoDataProcessParameter.SubjectID{i},'_BrainMask_05_91x109x91.nii']};
    
    DataDir= [AutoDataProcessParameter.DataProcessDir,filesep,AutoDataProcessParameter.StartingDirName,filesep,AutoDataProcessParameter.SubjectID{i}];
    
    [AllVolume,VoxelSize,theImgFileList, Header] =y_ReadAll(DataDir);
    
    [nDim1 nDim2 nDim3 nDimTimePoints]=size(AllVolume);
    
    [MaskData,MaskVox,MaskHead]=y_ReadRPI(ROIMask_Global{1}); %Global Mask
    
    % Convert into 2D. NOTE: here the first dimension is voxels,
    % and the second dimension is subjects. This is different from
    % the way used in y_bandpass.
    %AllVolume=reshape(AllVolume,[],nDimTimePoints)';
    AllVolume=reshape(AllVolume,[],nDimTimePoints);
    
    MaskDataOneDim=reshape(MaskData,[],1);
    MaskIndex = find(MaskDataOneDim);
    nVoxels = length(MaskIndex);
    %AllVolume=AllVolume(:,MaskIndex);
    AllVolume=AllVolume(MaskIndex,:);
    
    
    % Zero mean and unit variance for each time series
    AllVolume = (AllVolume - repmat(mean(AllVolume,2),[1,nDimTimePoints])) ./ repmat(std(AllVolume,0,2),[1,nDimTimePoints]);
    
    % Remove the NaN values, those could be caused by dividing by 0 standard deviation.
    AllVolume(find(isnan(AllVolume))) = 0;
    
    % Calculate the global signal
    GlobalMean = mean(AllVolume)';
    
    % Zero mean and unit variance for the global signal
    GlobalMean = (GlobalMean - mean(GlobalMean)) / std(GlobalMean);
    
    %FC Calculation
    FC=AllVolume*GlobalMean/(nDimTimePoints-1);
    
    GlobalCorrelation(i,1) = mean(FC);
    
end



