%% Berron et al manual subfields in unfolded space
addpath(genpath([getenv('AUTOTOP_DIR') '/tools']));
APres = 256; PDres = 128;
unfoldedSubfields = nan([APres PDres,35]);
n=0;
for s = 0:34
    sub = sprintf('%03d',s);
    for LR = 'LR'
        if LR=='L'; lr = 'left'; else lr = 'right'; end
        AP = load_untouch_nii(['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-AP.nii.gz']);
        PD = load_untouch_nii(['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-PD.nii.gz']);
        
        manlbl = ['datasets/ASHSatlas_7T-TSE/train/train' sub '/seg_' lr '.nii.gz'];
        manlbl_resamp = ['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-BerronManualSubfield_dseg.nii.gz'];
        aff = ['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/anat/sub-' sub '_desc-affine_from-T1w_to-corobl_type-itk_xfm.txt'];
        ref = ['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-subfields_dseg.nii.gz'];
        system(['antsApplyTransforms -d 3 -i ' manlbl ' -t ' aff ' -o ' manlbl_resamp ' -r ' ref ' -n MultiLabel']);
        Manuallbl = 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 = Manuallbl.img(AP==ap & PD==pd);
                lbls(lbls==0 | lbls >8) = []; % background label or surrounding labels
                unfoldedSubfields(ap,pd,n) = mode(lbls);
            end
        end
        
        %% also get volumes (in voxels)
        lbls = unique(Manuallbl.img(:));
        for l = 1:length(lbls)
            volumeBerron(l,n) = length(find(Manuallbl.img==lbls(l)));
        end
        hippUnfold = load_untouch_nii(['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-subfields_dseg.nii.gz']);
        lbls = unique(hippUnfold.img(:));
        for l = 1:length(lbls)
            volumeHippUnfold(l,n) = length(find(hippUnfold.img==lbls(l)));
        end
    end
end
save('Berron_unfoldedManualSubfields');


%% plot
load('Berron_unfoldedManualSubfields');
custom_colormap;
close all;

unfoldedSubfields(unfoldedSubfields==0) = nan;

% remap
u = unfoldedSubfields;
u(unfoldedSubfields==8) = 1;
u(unfoldedSubfields==1) = 2;
u(unfoldedSubfields==2) = 3;
u(unfoldedSubfields==4) = 4;
u(unfoldedSubfields==3) = 5;
u(unfoldedSubfields==5) = 6;
unfoldedSubfields = u;
Csubfields(6,:) = [.5 .5 .5];

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

% plot all boundaries
figure('units','normalized','outerposition',[0 0 1 1]);
imagesc(mode(unfoldedSubfields,3)); axis equal tight;
colormap(Csubfields(1:6,:));
saveas(gcf,'Berron_unfoldedMean.png');
% set(gca,'YDir','normal');
hold on;
alpha = 0.1;
msz = 4;
se = zeros([3 3]);
se(2,:) = 1; se(:,2) = 1;
for s = 1:70
    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,'Berron_unfoldedManualSubfieldsAlltogether.png');

%% plant-altman plots

addpath('BlandAltman');
load('Berron_unfoldedManualSubfields');

totvolumeBerron = sum(volumeBerron(2:7,:),1) * 0.3^3;
totvolumeHippunfold = sum(volumeHippUnfold(2:end,:),1) * 0.3^3;

totvolumeBerronLR = [totvolumeBerron(1:2:end); totvolumeBerron(2:2:end)];
totvolumeHippunfoldLR = [totvolumeHippunfold(1:2:end); totvolumeHippunfold(2:2:end)];

[c,h] = BlandAltman(totvolumeBerronLR',totvolumeHippunfoldLR',{'Berron','HippUnfold'});
saveas(gcf,'Berron_wholehippBlantAltman.png');

%% reconstruct borders

Berron_unfoldedavg = mode(unfoldedSubfields,3);

for s = 0:34
    sub = sprintf('%03d',s);
    for LR = 'LR'
        if LR=='L'; lr = 'left'; else lr = 'right'; end
        AP = load_untouch_nii(['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/sub-' sub '_hemi-' LR '_space-corobl_desc-cropped_modality-T2w_autotop/coords-AP.nii.gz']);
        PD = load_untouch_nii(['datasets/ASHSatlas_7T-TSE/HippUnfold/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(qpts)
            labelled_pts(n) = Berron_unfoldedavg(qpts(n,1),qpts(n,2));
        end
        origheader.img(mask) = labelled_pts;
        save_untouch_nii(origheader,['datasets/ASHSatlas_7T-TSE/HippUnfold/work/sub-' sub '/seg_T2w/sub-' sub '_hemi-' LR '_space-corobl_desc-BerronUnfoldedSubfields_dseg.nii.gz']);
    end
    disp(s);
end

