% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 
%
% This script runs the GLM analysis of human brain fPACT data.
%
% Authors: Shuai Na, Xiaoyun Yuan.
% Contacts: shuai@caltech.edu, yuanxy92@caltech.edu.
%
% Copyright 2019-2020 
% @ Caltech Optical Imaging Lab (COILAB), California Institute of Technology,
% 1200 E California Blvd, Pasadena, CA 91125, USA.
%
% Permission to use, copy, modify, and distribute this software ancdd
% its documentation for any purpose other than its incorporation into
% a commercial product is hereby granted without fee, provided that
% the above copyright notice appear in all copies and that both that       
% copyright notice and this permission notice appear in supporting          
% documentation, and that the name of California Institute of Technology    
% not be used in advertising or publicity pertaining to distribution of the 
% software without specific, written prior permission.
%
% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % 

clear;
close all;
fclose all;

%% filename
taskname = 'FT';      
T1name = 'sub3_T1_anatomy';
dirc = '../fPACT data/';
if strcmp(taskname, 'FT')
    dataname = 'sub3_FT.mat';
    tform_name = 'sub3_tform_FT.mat';
    load([dirc,tform_name]);
elseif strcmp(taskname, 'LP')
    dataname = 'sub3_LP.mat';
    tform_name = 'sub3_tform_LP.mat';
    load([dirc,tform_name]);
end
load([dirc,T1name]);

%% Load
load([dirc,dataname]);

%% GLM
% hrf 
stimulus_hrf = conv(stimulus_sig, hrf);
stimulus_hrf = stimulus_hrf(1:size(stimulus_sig, 1), 1);
stimulus_hrf = stimulus_hrf(1:20:end);
X = zeros(225, 3);
X(:, 1) = stimulus_hrf - mean(stimulus_hrf);
X(:, 2) = t_reg;
X(:, 3) = 1;
% reshape data
vol = single(data);
mean_vol = mean(vol, 4);
clear data;
[h, w, d, t] = size(vol(:, :, :, :));
for i = 1:t
    vol(:, :, :, i) = imgaussfilt3(vol(:, :, :, i), [3.26, 3.26, 3.26]); % spatial smoothing 0.3-mm FWHM
end
Y = double(reshape(vol(:, :, :, :), [h * w * d, t]));
depth = 60; 
vol_mask = zeros(h, w, d);
for r = 1:size(vol_mask, 1)
    for c = 1:size(vol_mask, 2)
        vol_mask(r, c, max(round(surface(r, c) + 0), 1):...
            min(round(surface(r, c) + depth), size(vol_mask, 3))) = 1;
    end
end
valid_inds = find(vol_mask == 1);
Y = Y(valid_inds, :);
mean_vol_vec = mean_vol(valid_inds);
% temporal bandpass filter 
bandwidth = [0.01, 0.2];
fs = 1 / 2;
bandpass_bandwidth = bandwidth / (fs / 2);
[b, a] = butter(4, bandwidth / (fs / 2), 'bandpass');
for ind = 1:size(valid_inds, 1)
    temp = Y(ind, :);
    if sum(isnan(temp)) == 0 && sum(isinf(temp)) == 0
        Y(ind, :) = filtfilt(b, a, double(temp));
    end
end
% solve beta
X_pinv = pinv(X' * X);
beta = X \ Y';
df = size(Y, 2) - rank(X);
% t test
C = [1 0 0]';
error = Y' - X * beta;
s2 = sum(error .^ 2, 1) / df;
t_stat = (C' * beta) ./ sqrt(s2 * (C' * X_pinv * C)); % t score
t_stat(isnan(t_stat)) = 0;
ltp = 1 - tcdf(t_stat, df);                           % p value
% coefficient of determination
ss_tot = Y' - repmat(mean(Y', 1), [t, 1]);
ss_tot = sum(ss_tot.^2, 1);
ss_res = sum(error .^ 2, 1);
rs = 1 - ss_res ./ (ss_tot + 1e-20);
% reshape to 3D volume
t_stat_vol = zeros(h, w, d);
p_val_vol = ones(h, w,  d);
rs_vol = zeros(h, w, d);
t_stat_vol(valid_inds) = t_stat;
p_val_vol(valid_inds) = ltp;
rs_vol(valid_inds) = rs;

%% Results
% position
region = [151, 650, 151, 650];
angle_rot = -90;
% temporal signal
mask = t_stat;
mask(mask <  max(t_stat(:)) * 0.7) = 0;
mask(mask > 0) = 1;
sig = zeros(225, 1);
mean_val = 0;
for ind = 1:size(t_stat, 2)
    if (mask(1, ind) == 1)
        sig = sig + reshape(Y(ind, :), [225, 1]);
        mean_val = mean_val + mean_vol_vec(ind);
    end
end
baseline = mean(sig(1:15));
sig = sig - baseline;
figure('Position', [100 300 1000 200]);hold on;
sticolor = [0.75 0.75 0.75];
for i = 1:7
    rectangle('Position', [round(onsets(i)) * 2, -4 / 100, 30, 12 / 100], ...
        'FaceColor', sticolor, 'EdgeColor', sticolor);
end
xlim([0 450]);
ylim([-4 8] / 100);
plot(0:2:448, sig /  mean_val);
title('temporal signal');

% t map on T1
t_stat_vol2 = imresize3(t_stat_vol, 2, 'nearest');
t_stat_vol3 = single(zeros(800, 800, 200));
t_stat_vol3(region(1):region(2), region(3):region(4), :) = t_stat_vol2;
t_stat_vol3 = permute(t_stat_vol3, [2, 1, 3]);
t_stat_vol3 = imwarp(t_stat_vol3, tform, 'nearest', 'OutputView', imref3d(size(t_stat_vol3)));
th = max(t_stat_vol3(:)) * 0.7;
t_stat_vol3(t_stat_vol3 < th) = th;
t_stat_vol3 = t_stat_vol3 - th;
t_stat_vol3 = t_stat_vol3 ./ max(t_stat_vol3(:));
map_t_stat = max(t_stat_vol3, [], 3);
% color code
fig = figure;
fig.Position = [100 100 800 800];
imshow(map_t_stat, 'border', 'tight');colormap('autumn');
t_stat_rgb = getframe(fig);
t_stat_rgb = t_stat_rgb.cdata;
t_stat_rgb = im2single(t_stat_rgb);
mask_visual = map_t_stat;
mask_visual(mask_visual < 1e-3) = 0;
mask_visual(mask_visual > 1e-3) = 1;
mask_visual = logical(cat(3, mask_visual, mask_visual, mask_visual));
close(fig);
pad_size = 80;
mask_visual2 = padarray(mask_visual, [pad_size, pad_size, 0]);
t_stat_rgb2 = padarray(t_stat_rgb, [pad_size, pad_size, 0]);
anatomy_visual = anatomy_img;
anatomy_visual(mask_visual2) = t_stat_rgb2(mask_visual2);
figure;imshow(imrotate(anatomy_visual, angle_rot)); title('t map, t>0.7t\_max')

% F map on T1
f_stat_val = rs_vol ./ (1 - rs_vol) * 222 / 2;
f_stat_vol2 = imresize3(f_stat_val, 2, 'nearest');
f_stat_vol3 = single(zeros(800, 800, 200));
f_stat_vol3(region(1):region(2), region(3):region(4), :) = f_stat_vol2;
f_stat_vol3 = permute(f_stat_vol3, [2, 1, 3]);
f_stat_vol3 = imwarp(f_stat_vol3, tform, 'nearest', 'OutputView', imref3d(size(f_stat_vol3)));
th_f = max(f_stat_vol3(:)) * 0.5;
f_stat_vol3(f_stat_vol3 < th_f) = th_f;
f_stat_vol3 = f_stat_vol3 - th_f;
f_stat_vol3 = f_stat_vol3 ./ max(f_stat_vol3(:));
map_f_stat = max(f_stat_vol3, [], 3);
% color code
fig = figure;
fig.Position = [100 100 800 800];
imshow(map_f_stat, 'border', 'tight');colormap('cool');
f_stat_rgb = getframe(fig);
f_stat_rgb = f_stat_rgb.cdata;
f_stat_rgb = im2single(f_stat_rgb);
mask_visual_f = map_f_stat;
mask_visual_f(mask_visual_f < 1e-3) = 0;
mask_visual_f(mask_visual_f > 1e-3) = 1;
mask_visual_f = logical(cat(3, mask_visual_f, mask_visual_f, mask_visual_f));
close(fig);
pad_size = 80;
mask_visual_f2 = padarray(mask_visual_f, [pad_size, pad_size, 0]);
f_stat_rgb2 = padarray(f_stat_rgb, [pad_size, pad_size, 0]);
anatomy_visual_f = anatomy_img;
anatomy_visual_f(mask_visual_f2) = f_stat_rgb2(mask_visual_f2);
figure;imshow(imrotate(anatomy_visual_f, angle_rot)); title('F map, F>0.5F\_max')
