%% Berron et al manual subfields in unfolded space
addpath(genpath([getenv('AUTOTOP_DIR') '/tools']));

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

FSlbls = [203:210 226];
FSlbls2 = [203:210 214 226]; % includes SRLM

mri_convert = 'singularity exec dependencies/pwighton_freesurfer_7.1.0.sif mri_convert';
APres = 256; PDres = 128;
unfoldedSubfields = nan([APres PDres,length(subs)]);
n=0;
warning off;
for s = 1:length(subs)
    for LR = 'LR'
        try
        if LR=='L'; lr = 'l'; else lr = 'r'; end
        sub = subs{s};
        AP = load_untouch_nii(['datasets/HCP1200/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-AP.nii.gz']);
        PD = load_untouch_nii(['datasets/HCP1200/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-PD.nii.gz']);
                
        manlbl = ['datasets/FS7/HCPtest/sub-' sub '/mri/' lr 'h.hippoAmygLabels-T1-*.v21.FS60.FSvoxelSpace.mgz'];
        manlbl_resamp = ['datasets/HCP1200/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-FSsubfields_dseg.nii.gz'];
        system([mri_convert ' ' manlbl ' ' manlbl_resamp]);
        aff = ['datasets/HCP1200/work/sub-' sub '/anat/sub-' sub '_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt'];
        ref = ['datasets/HCP1200/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-subfields_dseg.nii.gz'];
        system(['antsApplyTransforms -d 3 -i ' manlbl_resamp ' -t ' aff ' -o ' manlbl_resamp ' -r ' ref ' -n MultiLabel']);
        FSlbl = load_untouch_nii(manlbl_resamp);
        AP = round(AP.img * APres);
        PD = round(PD.img * PDres);
        n = n+1;
        for ap = 1:APres
            for pd = 1:PDres
                lbls = FSlbl.img(AP==ap & PD==pd);
                lbls = lbls(ismember(lbls,FSlbls));
                unfoldedSubfields(ap,pd,n) = mode(lbls);
            end
        end
        
        %% also get volumes (in voxels)
        for l = 1:length(FSlbls2)
            volumeFS(l,n) = length(find(FSlbl.img==FSlbls2(l)));
        end
        hippUnfold = load_untouch_nii(['datasets/HCP1200/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-subfields_dseg.nii.gz']);
        for l = 1:8
            volumeHippUnfold(l,n) = length(find(hippUnfold.img==l));
        end
        catch
            disp(['failed: ' sub '_hemi-' LR]);
        end
    end
end
save('FS_unfoldedSubfields');


%% plot
clear; close all;
load('FS_unfoldedSubfields.mat');

tmp = unfoldedSubfields;
tmp(unfoldedSubfields==203) = 1; % parasub
tmp(unfoldedSubfields==204) = 2; % presub
tmp(unfoldedSubfields==205) = 3; % sub
tmp(unfoldedSubfields==206) = 4; % ca1
tmp(unfoldedSubfields==208) = 5; % ca2&3
tmp(unfoldedSubfields==209) = 6; % Cca4
tmp(unfoldedSubfields==210) = 7; % dg
tmp(unfoldedSubfields==226) = 8; % tail
unfoldedSubfields=tmp;

custom_colormap;
Csubfields(4,:) = []; %no CA3
Csubfields = [Csubfields(1,:)+20/255; Csubfields]; % prosub
Csubfields = [Csubfields(1,:)+20/255; Csubfields]; % parasub
Csubfields(end-1:end,:) = []; %no cyst or SRLM
Csubfields = [Csubfields; [.5 .5 .5]]; %tail

% figure('units','normalized','outerposition',[0 0 1 1]);
for s = 1:size(unfoldedSubfields,3)
    unfoldedSubfields(:,:,s) = fillmissing(unfoldedSubfields(:,:,s),'nearest');
%     subplot(10,15,s); imagesc(squeeze(unfoldedSubfields(:,:,s))); axis equal tight;
%     colormap(Csubfields);
end

% plot all boundaries
figure('units','normalized','outerposition',[0 0 1 1]);
imagesc(mode(unfoldedSubfields,3)); axis equal tight;
colormap(Csubfields);
saveas(gcf,'FS_unfoldedMean.png');
% set(gca,'YDir','normal');
hold on;
alpha = 0.05;
msz = 4;
se = zeros([3 3]);
se(2,:) = 1; se(:,2) = 1;
for s = 1:size(unfoldedSubfields,3)
    img = squeeze(unfoldedSubfields(:,:,s));
    img = imresize(img,2,'nearest');
    lbls = unique(img(:));
    for i = 1:length(lbls)
        [x,y] = ind2sub(size(img),find(imdilate(img==lbls(i),se) & img~=lbls(i)));
        scatter(y/2,x/2,msz,'MarkerFaceColor','k','MarkerEdgeColor','none', 'MarkerFaceAlpha',alpha,'MarkerEdgeAlpha',alpha);
    end
end
saveas(gcf,'FS_unfoldedManualSubfieldsAlltogether.png');

%% plant-altman plots
addpath('BlandAltman');

totvolumeFS = sum(volumeFS,1) * 0.3^3; % this should include SRLM and hippocampal fissure and hippocampal tail
% totvolumeFS(totvolumeFS<1000) = nan;
totvolumeHippunfold = sum(volumeHippUnfold,1) * 0.3^3;

totvolumeFSLR = [totvolumeFS(1:2:end); totvolumeFS(2:2:end)];
totvolumeHippunfoldLR = [totvolumeHippunfold(1:2:end); totvolumeHippunfold(2:2:end)];

[c,h] = BlandAltman(totvolumeFSLR',totvolumeHippunfoldLR',{'FS','HippUnfold'});
saveas(gcf,'FS_wholehippBlantAltman.png');

%% reconstruct borders

FS_unfoldedavg = mode(unfoldedSubfields,3);

for s = 1:length(subs)
    for LR = 'LR'
        if LR=='L'; lr = 'left'; else lr = 'right'; end
        sub = subs{s};
        AP = load_untouch_nii(['datasets/HCP1200/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-AP.nii.gz']);
        PD = load_untouch_nii(['datasets/HCP1200/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-PD.nii.gz']);
        origheader = AP;
        mask = find(AP.img>0 & PD.img>0);
        qpts = [AP.img(mask) PD.img(mask)];
        qpts(:,1) = ceil(qpts(:,1)*256);
        qpts(:,2) = ceil(qpts(:,2)*128);
        labelled_pts = [];
        for n = 1:length
            
            
        origheader.img(mask) = labelled_pts;
        save_untouch_nii(origheader,['datasets/HCP1200/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-FSUnfoldedSubfields_dseg.nii.gz']);
        disp([sub '_hemi-' LR ' done']);
    end
end

