% homotopy_denoising.m

% Within-subjects:
% 1. Linear trends
% 2. Quadratic trends
% 3. 6 rigid-body motion parameters
% 4. 5 principal components of CSF signals
% 5. Body temperature
% 6. Isoflurane dose
% 7. Global signal regression (optional)
% 8. Bandpass filtering

% Between-subjects:
% 1. Age
% 2. Mean Body Temperature
% 3. Mean Isoflurane Dose
% 4. ROI pair Euclidean distance
% 5. ROI volume

% Note: requires Tools for NIfTI and ANALYZE image
% (https://www.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image)

% Use smoothed or unsmoothed data for ROI extraction?
useSmoothedData = false;

% Use homotopic-SC as well as homotopic-nSC?
useSC = true;

% Regress out global signal?
useGlobalSignal = false;

%% LOAD DATA

% Collect file paths to normalized functional images:
functPath = '';
nscan = 19;
timeSeries = cell(1,nscan); % 5 normals, 4 scans each (except normal #5, 3 scans)
for scan = 1:nscan
    timeSeries{scan} = cell(1,175);
    for timepoint = 6:180 % skip first 5 volumes (T1 equilibration)
        if useSmoothedData
            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
                sprintf('s_w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
        else
            timeSeries{scan}{timepoint-5} = fullfile(functPath,sprintf('funct_%02g',scan),...
                sprintf('w_funct3D_bird_%02g_time_%03g.nii',[scan, timepoint]));
        end
    end
end

% Load ROIs:
ROI_path = '';
SongSystem = fullfile(ROI_path,'SongSystemROIs.nii'); % n = 20
mask = load_nii(fullfile(ROI_path,'new_temp_brain_mask.nii'));
nROI = 20;

% Get ROI Subscripts/Coordinates:
rois = load_nii(SongSystem);
roi_subscripts = cell(1,nROI);
for i = 1:nROI
    [x, y, z] = ind2sub(size(rois.img),find(rois.img(:)==i));
    roi_subscripts{i} = [x,y,z];
end

% Get CSF Subscripts/Coordinates:
CSF_mask = load_nii(fullfile(ROI_path, 'CSF_mask.nii'));
[CSF_x,CSF_y,CSF_z] = ind2sub(size(CSF_mask.img), find(CSF_mask.img(:)>0));
n_CSF = length(CSF_x); 

% Extract ROI Data & CSF Signals:
ROI_dat = cell(1,nscan); CSF = cell(1,nscan); globalSignal = zeros(175,nscan); % Initialize
hWait = waitbar(0,'Extracting functional data...');
for scan = 1:nscan
    waitbar(scan/nscan, hWait);
    
    ROI_dat{scan} = zeros(175, nROI); CSF_dat = zeros(175, n_CSF); % Initialize
    for j = 1:175 % iterate through time (j)
        
        img = load_nii(timeSeries{scan}{j});
        globalSignal(j, scan) = mean(img.img(mask.img(:)>0)); % Average whole brain signal for time point j
        
        for k = 1:nROI % iterate through ROIs (k)
            nvox = size(roi_subscripts{k}, 1); % # voxels in ROI
            roi_holder = zeros(1, nvox);
            for l = 1:nvox % iterate through voxels of ROI
                roi_holder(l) = img.img(roi_subscripts{k}(l,1), roi_subscripts{k}(l,2), roi_subscripts{k}(l,3));
            end
            ROI_dat{scan}(j,k) = mean(roi_holder); % average voxels of ROI
        end    
        
        for k = 1:n_CSF % iterate through CSF voxels (k)
            CSF_dat(j,k) = img.img(CSF_x(k), CSF_y(k), CSF_z(k));
        end  
    end
    
    [~, PCs] = pca(CSF_dat); % pca across voxels of CSF mask
    CSF{scan} = PCs(:,1:5); % save PC's 1-5
    
end; delete(hWait)

% Load Motion Data:
motionPath = '';
motion_params = cell(1,nscan);
for i = 1:nscan
    load(fullfile(motionPath, sprintf('motion_scan_%02g.mat',i)))
    motion_params{i} = motion;
end

% Load Physiology Data (Temperature (degrees Celsius), Isoflurane (%)):
physio_path = '';
physio_params = cell(1,nscan);
physio_times = [60:60:540,576]; % 60 seconds to 576 seconds
TR = 1.0671*3; scan_times = TR:TR:TR*175;
for i = 1:nscan
    tbl = readtable(fullfile(physio_path,'All_Normal_Condition_Physiology_Records.xlsx'),'Sheet',i);
    pp_temp = spline(physio_times,tbl.Temperature); % piecewise polynomial of cubic spline
    pred_temp = ppval(pp_temp,scan_times); % evaluate PP at x
    pp_iso = spline(physio_times,tbl.Isoflurane); % piecewise polynomial of cubic spline
    pred_iso = ppval(pp_iso,scan_times); % evaluate PP at x
    physio_params{i} = [pred_temp',pred_iso'];
end

% Labels:
SongSystem_labels = {'L Area X','R Area X','L LMAN','R LMAN','L RA','R RA',...
    'L HVC','R HVC','L Aud. Forebrain','R Aud. Forebrain','L Field L',...
    'R Field L','L Diencephalon','R Diencephalon','L Medulla','R Medulla',...
    'L Uva','R Uva','L DM','R DM'};

ages = [25, 45, 65, 89, 24, 45, 65, 90, 25, 46, 67, 89, 24, 46, 66, 90, 24, 45, 67];

% Get ROI-level tSNR Data:
tSNR_mat = zeros(nscan, nROI); 
for i = 1:nscan
    for j = 1:nROI
        mean_sig = mean(ROI_dat{i}(:,j));
        noise_SD = std(ROI_dat{i}(:,j));
        tSNR_mat(i,j) = mean_sig ./ noise_SD;
    end
end

% Get Mean Physiological Parameters:
temps = zeros(19,1); iso = zeros(19,1);
for i = 1:19
    temps(i) = mean(physio_params{i}(:,1));
    iso(i) = mean(physio_params{i}(:,2));
end

% Get Mean Volumes and Euclidean Distances:
vox_dim = rois.hdr.dime.pixdim(2:4); 
vox_volume = prod(vox_dim);
m_volumes = zeros(nROI); euclideans = zeros(nROI);

% Get ROI Centroids:
centroids = zeros(nROI,3);
for i = 1:nROI
    [x,y,z] = ind2sub(size(rois.img),find(rois.img==i));
    centroids(i,:) = [mean(x),mean(y),mean(z)].*vox_dim;
end
for i = 1:nROI
    for j = 1:nROI
        m_volumes(i,j) = mean([sum(rois.img(:)==i)*vox_volume,sum(rois.img(:)==j)*vox_volume]);
        euclideans(i,j) = sqrt(sum((centroids(i,:)-centroids(j,:)).^2));
    end
end

%% RUN DENOISING

% Labels for FC:
labels = {'L_Area_X','R_Area_X','L_LMAN','R_LMAN','L_RA','R_RA',...
    'L_HVC','R_HVC','L_Aud_Forebrain','R_Aud_Forebrain','L_Field_L',...
    'R_Field_L','L_Diencephalon','R_Diencephalon','L_Medulla','R_Medulla',...
    'L_Uva','R_Uva','L_DM','R_DM'};

useScans = 1:19;

if useSC
    nROI = 20;
else
    nROI = 12;
end

% Preprocess ROI Data / Within-Subjects Nuisance Regression:
Fs = .3124; TR = 1.0671*3; time = (TR:TR:TR*175)'; time_sq = time.^2; intercept = ones(175,1);
hWait = waitbar(0,'Denoising...');
for scan = useScans
    if ishandle(hWait); waitbar(scan/nscan, hWait); end
    
    % Get nuisance regression covariates:
    if all(physio_params{scan}(:,2)==mode(physio_params{scan}(:,2))) % iso is a constant, so don't need
        if useGlobalSignal
            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1), globalSignal(:,scan)];
        else
            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}(:,1)];
        end
    else
        if useGlobalSignal
            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}, globalSignal(:,scan)];
        else
            nuisance = [intercept, time, time_sq, motion_params{scan}, CSF{scan}, physio_params{scan}];
        end
    end
    
    % Perform nuisance regression, retain residuals:
    for j = 1:nROI
        [~, ~, ROI_dat{scan}(:,j)] = regress(ROI_dat{scan}(:,j), nuisance);
    end
    
    % Perform bandpass filtering: 
    [ROI_dat{scan},~,~] = bst_bandpass_filtfilt(ROI_dat{scan}', Fs, .008, .1, 0, 'iir');
    ROI_dat{scan} = ROI_dat{scan}';

end
if ishandle(hWait); delete(hWait); end 

% Calculate Correlation Matrix:
zmats = zeros(nROI, nROI, nscan);
for scan = useScans
    zmats(:,:,scan) = fisherz(corr(ROI_dat{scan}(:,1:nROI)));
end

% Between-Subjects Nuisance Regression:
nscan = length(useScans);
numconn = nROI*(nROI-1)/2*nscan;
birds = [1,1,1,1,2,2,2,2,3,3,3,3,4,4,4,4,5,5,5];
scan = zeros(numconn,1); bird = zeros(numconn,1);  
FC = zeros(numconn,1); age = zeros(numconn,1); 
connection_type = zeros(numconn,1); euclidean = zeros(numconn,1); 
m_volume = zeros(numconn,1); mean_temp = zeros(numconn,1); 
mean_iso = zeros(numconn,1); conn_label = cell(numconn,1); count = 0; 
for i = 1:nROI-1
    for j = i+1:nROI
        FC(count+1:count+nscan) = squeeze(zmats(i,j,:));
        conn_label(count+1:count+nscan) = {[labels{i},'_',labels{j}]};
        age(count+1:count+nscan) = ages(useScans);
        bird(count+1:count+nscan) = birds(useScans);
        scan(count+1:count+nscan) = useScans;
        euclidean(count+1:count+nscan) = repmat(euclideans(i,j),nscan,1);
        m_volume(count+1:count+nscan) = repmat(m_volumes(i,j),nscan,1);
        if ((mod(i,2) && ~mod(j,2)) || (~mod(i,2) && mod(j,2))) && ~(mod(i,2) && j==(i+1)) % 4: heterotopic (contralateral)
            connection_type(count+1:count+nscan) = 4*ones(nscan,1);
        elseif mod(i,2) && j==(i+1)
            if i>12 || j>12
                connection_type(count+1:count+nscan) = 3*ones(nscan,1); % 3: homotopic-SC
            else
                connection_type(count+1:count+nscan) = 2*ones(nscan,1); % 2: homotopic-nSC
            end
        elseif (mod(i,2) && mod(j,2)) || (~mod(i,2) && ~mod(j,2)) % 1: ipsilateral
            connection_type(count+1:count+nscan) = ones(nscan,1);
        end
        mean_temp(count+1:count+nscan) = temps(useScans);
        mean_iso(count+1:count+nscan) = iso(useScans);
        count = count+nscan;
    end
end
bird = categorical(bird); scan = categorical(scan);
if useSC
    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic-nSC','Homotopic-SC','Heterotopic'});
else
    connection_type = categorical(connection_type,1:4,{'Ipsilateral','Homotopic','Homotopic','Heterotopic'});
end

% Mixed Model:
FC_data = table(FC, bird, scan, conn_label, connection_type, age, mean_temp, ...
    mean_iso, euclidean, m_volume, 'VariableNames', {'FC', 'Bird', 'Scan', ...
    'Connection', 'Type', 'Age', 'Temp', 'Iso', 'Euclidean', 'Volume'});

% Select Between-subjects nuisance parameters for formula:
FC_formula = 'FC ~ Type*Age + Type*Iso + Temp + Euclidean + Volume + (1|Bird)';

% Run LME random intercept model:
lme = fitlme(FC_data, FC_formula,  'StartMethod','random',  'verbose',true,  'CheckHessian',true); 
lme.display

% Get random effect coefficients:
[B, Bnames, stats] = randomEffects(lme) %#ok

% Check residuals:
figure; lme.plotResiduals

% Get marginal (fixed-effects only) vs. conditional (w/ random) R2
lme.Rsquared

resid_marginal = residuals(lme, 'Conditional',false); % fixed-effects only
resid_conditional = residuals(lme, 'Conditional',true); % fixed-effects only
R2_marginal = 1 - sum(resid_marginal.^2) / lme.SST
R2_conditional = 1 - sum(resid_conditional.^2) / lme.SST