addpath(genpath('gifti-master'));


f = fopen('datasets/testSet.txt');
subs = fscanf(f,'%c');
fclose(f);
subs = strsplit(subs);
subs(end) = [];

shp = [254 126];

DTIsubs = [];
n=0;
for s = 1:length(subs)
    for LR = 'LR'
        n=n+1;
        midthick{n} = ['datasets/HCP1200/results/sub-' subs{s} '/surf_T2w/sub-' subs{s} '_hemi-' LR '_space-T1w_midthickness.surf.gii'];
        
        Vmid{n} = gifti(midthick{n});
        thickness{n} = gifti(strrep(midthick{n},'midthickness.surf','thickness.shape'));
        curvature{n} = gifti(strrep(midthick{n},'midthickness.surf','curvature.shape'));
        gyrification{n} = gifti(strrep(midthick{n},'midthickness.surf','gyrification.shape'));
        
        system(['wb_command -volume-to-surface-mapping ../datasets/HCP1200/results/sub-' subs{s} '/anat/sub-' subs{s} '_desc-preproc_T1w.nii.gz '...
            midthick{n} ' '...
            strrep(midthick{n},'midthickness.surf','T1w.shape') ' -trilinear']);
        T1w{n} = gifti(strrep(midthick{n},'midthickness.surf','T1w.shape'));
        
        system(['wb_command -volume-to-surface-mapping ../datasets/HCP1200/results/sub-' subs{s} '/anat/sub-' subs{s} '_space-T1w_desc-preproc_T2w.nii.gz '...
            midthick{n} ' '...
            strrep(midthick{n},'midthickness.surf','T2w.shape') ' -trilinear']);
        T2w{n} = gifti(strrep(midthick{n},'midthickness.surf','T2w.shape'));
        
        MDfn{s} = ['~/projects/ctb-akhanf/ext-data/hcp1200/deriv/hcp1200_dtifit/results/sub-' subs{s} '/dtifit_MD.nii.gz'];
        if exist(MDfn{s},'file')
            DTIsubsidx = [DTIsubs s];
            
            system(['wb_command -volume-to-surface-mapping ' MDfn{s} ' '...
                midthick{n} ' '...
                strrep(midthick{n},'midthickness.surf','MD.shape') ' -trilinear']);
            MD{n} = gifti(strrep(midthick{n},'midthickness.surf','MD.shape'));
            system(['wb_command -volume-to-surface-mapping ' strrep(MDfn{s},'_MD','_FA') ' '...
                midthick{n} ' '...
                strrep(midthick{n},'midthickness.surf','FA.shape') ' -trilinear']);
            FA{n} = gifti(strrep(midthick{n},'midthickness.surf','FA.shape'));            
        end
    end
end
save('testSet_flatmaps');

%% plot

clear; close all;
load('testSet_flatmaps');
mkdir('flatmaps')

for s=1:length(Vmid)
    V(s,:,:) = Vmid{s}.vertices;
end
V([1:2:size(V,1)],:,1) = -V([1:2:size(V,1)],:,1); %flip left
V = squeeze(mean(V,1));
FV.vertices = V;
FV.faces = Vmid{1}.faces;

load('BigBrain_ManualSubfieldsUnfolded.mat');
subfields_avg = imresize(subfields_avg,shp,'nearest')';
plot_foldunfold(subfields_avg',FV,0,0,'subfields');
custom_colormap;
colormap(Csubfields(1:5,:));
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/subfields.png');

for s=1:length(thickness)
    thick(s,:) = thickness{s}.cdata;
end
thick = reshape(mean(thick,1),shp);
plot_foldunfold(thick,FV,0,1,'mean thickness');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/thickness.png');

for s=1:length(gyrification)
    gyr(s,:) = gyrification{s}.cdata;
end
gyr = reshape(mean(gyr,1),shp);
plot_foldunfold(gyr,FV,0,1,'mean gyrification');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/gyrification.png');

for s=1:length(curvature)
    curv(s,:) = curvature{s}.cdata;
end
curv = reshape(mean(curv,1),shp);
plot_foldunfold(curv,FV,0,1,'mean curvature');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/curvature.png');

for s=1:length(T1w)
    T1(s,:) = T1w{s}.cdata;
end
T1(T1<0) = nan;
for s=1:length(T2w)
    T2(s,:) = T2w{s}.cdata;
end
T2(T2<0) = nan;
T1overT2 = reshape(mean(T1./T2,1),shp);
T1overT2(T1overT2>10) = nan;
plot_foldunfold(T1overT2,FV,0,1,'mean T1/T2w');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/T1T2.png');

n=0;
for s = 1:length(FA)
    if ~isempty(FA{s})
        n=n+1;
        fa(n,:) = FA{s}.cdata;
    end
end
fa = reshape(mean(fa,1),shp);
plot_foldunfold(fa,FV,0,1,'fractional anisotropy');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/FA.png');

n=0;
for s = 1:length(MD)
    if ~isempty(MD{s})
        n=n+1;
        md(n,:) = MD{s}.cdata;
    end
end
md = reshape(mean(md,1),shp);
plot_foldunfold(md,FV,0,1,'mean diffusivity');
hold on; imcontour(subfields_avg,'w');
saveas(gcf,'flatmaps/MD.png');  
        
