% VMHC.m

% Note: requires Tools for NIfTI and ANALYZE image
% (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)
% Requires DPABI toolbox (http://rfmri.org/dpabi) for cluster inference

main_path = '';

useGlobalSignal = false;

%% CALCULATE HOMOTOPIC VOXEL-WISE CORRELATIONS:

% Load Mask and Get Voxel-wies Coordinates:
% Middle 2 voxels: 138-139
% Struct: .0703 x .0703 mm
% Funct: .1406 x .5 mm (2 x 7.12 struct voxels equivalent)
mask = load_nii(fullfile(main_path, 'mirrored_groupwise_template.nii'));
ind = find(mask.img(:)>0); [x,y,z] = ind2sub(size(mask.img),ind);
ind = x<139; numvox = sum(ind);
x = x(ind); y = y(ind); z = z(ind);
r_x = zeros(size(x)); for i = 1:length(x); r_x(i) = 139+(138-x(i)); end
nscan = 19;

% Perform denoising:
Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2;
intercept = ones(175,1); homotopy_img = zeros(276,256,59,nscan); % voxel-wise homotopy data
h_wait = waitbar(0,sprintf('Loading functional series %02g...',i));
signals = cell(1,nscan); globalSignal = zeros(175,19);
funct_path_outer = ''; 
for i = 1:nscan
    % Check which physio data is usable:
    if all(physio_params{i}(:,2)==mode(physio_params{i}(:,2))) % iso is a constant, therefore unnecessary
        physio_use = physio_params{i}(:,1); 
    else
        physio_use = physio_params{i}; 
    end
    
    % Load Functional:
    waitbar(0, h_wait,sprintf('Loading functional series %02g...',i));
    funct_path = fullfile(funct_path_outer,sprintf('funct_%02g',i));
    funct4D = zeros(276,256,59,175);
    for j = 6:180
        waitbar(j/180,h_wait);
        funct = load_nii(fullfile(funct_path,sprintf('s_R_mirror_funct3D_bird_%02g_time_%03g.nii',[i,j])));
        globalSignal(j-5, i) = nanmean(funct.img(mask.img(:)>0));
        funct4D(:,:,:,j-5) = funct.img;
    end
    
    % Calculate correlations:
    waitbar(0,h_wait,sprintf('Calculating homotopic correlations for series %02g...',i));
    signals{i}{1} = zeros(175,numvox); signals{i}{2} = zeros(175,numvox); 
    for j = 1:numvox
        if rem(j,10000)==0; waitbar(j/numvox,h_wait); end
        signals{i}{1}(:,j) = squeeze(funct4D(x(j),y(j),z(j),:));
        signals{i}{2}(:,j) = squeeze(funct4D(r_x(j),y(j),z(j),:));
        if useGlobalSignal
            [~,~,signals{i}{1}(:,j)] = regress(signals{i}{1}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use,globalSignal(:,i)]);
            [~,~,signals{i}{2}(:,j)] = regress(signals{i}{2}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use,globalSignal(:,i)]);
        else
            [~,~,signals{i}{1}(:,j)] = regress(signals{i}{1}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use]);
            [~,~,signals{i}{2}(:,j)] = regress(signals{i}{2}(:,j),[intercept,time,time_sq,motion_params{i},CSF{i},physio_use]);
        end
    end
    [signals{i}{1},~,~] = bst_bandpass_filtfilt(signals{i}{1}',Fs,.008,.1,0,'iir');
    [signals{i}{2},~,~] = bst_bandpass_filtfilt(signals{i}{2}',Fs,.008,.1,0,'iir');
    signals{i}{1} = signals{i}{1}'; signals{i}{2} = signals{i}{2}';
    r = fast_corr(signals{i}{1},signals{i}{2});
    for j = 1:numvox
        homotopy_img(x(j),y(j),z(j),i) = r(j); 
        homotopy_img(r_x(j),y(j),z(j),i) = r(j);
    end
end; close(h_wait)

%% t-statistics across subjects

controlCovariates = true;
covariates = [zscore(age), zscore(temps), zscore(iso)];

% Calculate t-stats summarizing voxel-wise homotopy across birds/scans:
homotopy_t = zeros(size(mask.img));
age_t = homotopy_t;
h_wait = waitbar(0,'Calculating t-statistics...');
for i = 1:numvox
    if rem(i,1000)==0; waitbar(i/numvox, h_wait); end
    zvec = fisherz(squeeze(homotopy_img(x(i),y(i),z(i),:)));
    if controlCovariates
        lm = fitlm(covariates, zvec);
        homotopy_t(x(i),y(i),z(i)) = lm.Coefficients.tStat(1);
        homotopy_t(r_x(i),y(i),z(i)) = lm.Coefficients.tStat(1);
        age_t(x(i),y(i),z(i)) = lm.Coefficients.tStat(2);
        age_t(r_x(i),y(i),z(i)) = lm.Coefficients.tStat(2);
    else
        [~,~,~,stats] = ttest(zvec);
        homotopy_t(x(i),y(i),z(i)) = stats.tstat;
        homotopy_t(r_x(i),y(i),z(i)) = stats.tstat;
    end
end; close(h_wait)

% Save t-Stats Images:
% Middle 2 voxels: 138-139
img = mask; img.img = homotopy_t; img.img(137:140,:,:) = 0; % Discard 2 voxels on each side of midline
if controlCovariates
    save_nii(img,fullfile(main_path,'t_stat_homotopy_controlled.nii'))
    img2 = mask; img2.img = age_t;
    save_nii(img2, fullfile(main_path,'t_stat_age_controlled.nii'))
else
    save_nii(img, fullfile(main_path,'t_stat_homotopy.nii'))
end

img.img(135:142,:,:) = 0; % Discard 4 voxels on each side of midline
if controlCovariates
    save_nii(img,fullfile(main_path,'t_stat_homotopy_8voxgap_controlled.nii'))  
    img2.img(136:141,:,:) = 0; 
    save_nii(img2, fullfile(main_path,'t_stat_age_8voxgap_controlled.nii'))
else
    save_nii(img,fullfile(main_path,'t_stat_homotopy_8voxgap.nii'))
end

%% Run cluster inference: Homotopy
% Extract the significant clusters of homotopy using cluster-extent method
% on only 1 hemisphere (hemispheres are mirror images)

if controlCovariates
    StatsImgFile = fullfile(main_path,'t_stat_homotopy_8voxgap_controlled.nii'); % more conservative midline deletion
    OutputName = fullfile(main_path,'T_homotopy_clusters_controlled.nii');
else
    StatsImgFile = fullfile(main_path,'t_stat_homotopy_8voxgap.nii'); % more conservative midline deletion
    OutputName = fullfile(main_path,'T_homotopy_clusters.nii');
end

MaskFile = fullfile(main_path, 'mirrored_groupwise_template.nii');
IsTwoTailed = false; VoxelPThreshold = .001; ClusterPThreshold = .05;
Flag = 'T'; Df1 = 18; Df2 = []; VoxelSize = mask.hdr.dime.pixdim(2:4);
Header = []; dLh = []; 
[Data_Corrected, ClusterSize, Header] = y_GRF_Threshold(StatsImgFile, ...
    VoxelPThreshold,IsTwoTailed,ClusterPThreshold,OutputName,MaskFile,...
    Flag,Df1,Df2,VoxelSize,Header, dLh);
y_ClusterReport(Data_Corrected, Header, 18) % 18 is the cluster (voxel neighbor) connectivity criterion

% Save bilateral version of cluster-thresholded t-stats:
if controlCovariates
    clusterT = load_nii(fullfile(main_path,'ClusterThresholded_T_homotopy_clusters_controlled.nii'));
    clusterT.img(143:276,:,:) = clusterT.img(134:-1:1,:,:); % 135:142 voxels were discarded, so start from 143:(134+142)
    save_nii(clusterT,fullfile(main_path,'T_homotopy_clusters_bilat_controlled.nii'));
else
    clusterT = load_nii(fullfile(main_path,'ClusterThresholded_T_homotopy_clusters.nii'));
    clusterT.img(143:276,:,:) = clusterT.img(134:-1:1,:,:); % 135:142 voxels were discarded, so start from 143:(134+142)
    save_nii(clusterT,fullfile(main_path,'T_homotopy_clusters_bilat.nii'));
end