%%%% get AHD and DICE metrics for segmentations in relation to the 7TAMI
%%%% ground-truth.
%%%%% this scripts assumes you already have all native-space segmentations from
%%%% the three methods ready as .nii files.

%%%%% Add SPM-----------------------------------
spm_path = 'C:\Program Files\spm12';addpath(spm_path)
spm('Defaults', 'fMRI');
spm_jobman('initcfg');
%%%%----------------------------------------
%%%% DEFINE THE FOLLOWING
%%%%---------------------------------------
threshold = 0.25;
STAMI_directory = 'C:\Users\gargy\OneDrive\Desktop\7TAMI\mni_masks';
analysis_output_directory = 'C:\Users\gargy\OneDrive - University of Stirling\HIPS_THOMAS_Manuscript';
%%%%-----------------------------------------------------
% %%%% for HCP35
directory = 'D:\HCP35\unzipped'; % for HCP35
subfolder = 'anat\T1'; % for HCP35
probabilistic_maps_dir = 'D:\HCP35_probabilistic_maps'; % for HCP35
thresh_prob_maps_dir = 'D:\HCP35_thresh_prob_maps'; % for HCP35
MNI_vs_STAMI_analysis_name = 'HCP35_MNI_vs_STAMI.xlsx'; %for HCP35
Dice_AHD_Native_FileName='DICE_AHD_STAMI_native_HCP35.csv';
% 
% % % % HCP35 folder list
folder_list = {
'mgh_1001',...
'mgh_1002',...
'mgh_1003',...
'mgh_1004',...
'mgh_1005',...
'mgh_1006',...
'mgh_1007',...
'mgh_1008',...
'mgh_1009',...
'mgh_1010',...
'mgh_1011',...
'mgh_1012',...
'mgh_1013',...
'mgh_1014',...
'mgh_1015',...
'mgh_1016',...
'mgh_1017',...
'mgh_1018',...
'mgh_1019',...
'mgh_1021',...
'mgh_1022',...
'mgh_1023',...
'mgh_1024',...
'mgh_1025',...
'mgh_1026',...
'mgh_1027',...
'mgh_1028',...
'mgh_1029',...
'mgh_1030',...
'mgh_1031',...
'mgh_1032',...
'mgh_1033',...
'mgh_1034'};

% % %%%%%%---------------------------------------
[~, num] = size(folder_list);
cd(directory);
% % % %-----------------------------------------------------

% % %%%% 6. MAKE PROBABILISTIC MAPS for 7T-AMI's lateral group
% % %%%%--------------------------------------------
%%%%% List of ROI names (same names across subjects)
roi_names = {
   'resliced_native_AV_L_FS_DTI_MNI',...
   'resliced_native_AV_L_FS_T1_MNI',...
   'resliced_native_AV_L_HT_New_MNI',...
   'resliced_native_AV_R_FS_DTI_MNI',...
   'resliced_native_AV_R_FS_T1_MNI',...
   'resliced_native_AV_R_HT_New_MNI',...
   'resliced_native_CL_L_FS_DTI_MNI',...
   'resliced_native_CL_L_FS_T1_MNI',...
   'resliced_native_CL_L_HT_New_MNI',...
   'resliced_native_CL_R_FS_DTI_MNI',...
   'resliced_native_CL_R_FS_T1_MNI',...
   'resliced_native_CL_R_HT_New_MNI',...
   'resliced_native_CM_L_FS_DTI_MNI',...
   'resliced_native_CM_L_FS_T1_MNI',...
   'resliced_native_CM_L_HT_New_MNI',...
   'resliced_native_CM_R_FS_DTI_MNI',...
   'resliced_native_CM_R_FS_T1_MNI',...
   'resliced_native_CM_R_HT_New_MNI',...
   'resliced_native_CeM_L_FS_DTI_MNI',...
   'resliced_native_CeM_L_FS_T1_MNI',...
   'resliced_native_CeM_R_FS_DTI_MNI',...
   'resliced_native_CeM_R_FS_T1_MNI',...
   'resliced_native_LD_L_FS_DTI_MNI',...
   'resliced_native_LD_L_FS_T1_MNI',...
   'resliced_native_LD_R_FS_DTI_MNI',...
   'resliced_native_LD_R_FS_T1_MNI',...
   'resliced_native_LGN_L_FS_DTI_MNI',...
   'resliced_native_LGN_L_FS_T1_MNI',...
   'resliced_native_LGN_L_HT_New_MNI',...
   'resliced_native_LGN_R_FS_DTI_MNI',...
   'resliced_native_LGN_R_FS_T1_MNI',...
   'resliced_native_LGN_R_HT_New_MNI',...
   'resliced_native_LP_L_FS_DTI_MNI',...
   'resliced_native_LP_L_FS_T1_MNI',...
   'resliced_native_LP_R_FS_DTI_MNI',...
   'resliced_native_LP_R_FS_T1_MNI',...
   'resliced_native_LSg_L_FS_DTI_MNI',...
   'resliced_native_LSg_L_FS_T1_MNI',...
   'resliced_native_LSg_R_FS_DTI_MNI',...
   'resliced_native_LSg_R_FS_T1_MNI',...
   'resliced_native_MD_Pf_L_FS_DTI_comb_MNI',...
   'resliced_native_MD_Pf_L_FS_T1_comb_MNI',...
   'resliced_native_MD_Pf_L_HT_New_MNI',...
   'resliced_native_MD_Pf_R_FS_DTI_comb_MNI',...
   'resliced_native_MD_Pf_R_FS_T1_comb_MNI',...
   'resliced_native_MD_Pf_R_HT_New_MNI',...
   'resliced_native_MDl_L_FS_DTI_MNI',...
   'resliced_native_MDl_L_FS_T1_MNI',...
   'resliced_native_MDl_R_FS_DTI_MNI',...
   'resliced_native_MDl_R_FS_T1_MNI',...
   'resliced_native_MDm_L_FS_DTI_MNI',...
   'resliced_native_MDm_L_FS_T1_MNI',...
   'resliced_native_MDm_R_FS_DTI_MNI',...
   'resliced_native_MDm_R_FS_T1_MNI',...
   'resliced_native_MGN_L_FS_DTI_MNI',...
   'resliced_native_MGN_L_FS_T1_MNI',...
   'resliced_native_MGN_L_HT_New_MNI',...
   'resliced_native_MGN_R_FS_DTI_MNI',...
   'resliced_native_MGN_R_FS_T1_MNI',...
   'resliced_native_MGN_R_HT_New_MNI',...
   'resliced_native_MVRe_L_FS_DTI_MNI',...
   'resliced_native_MVRe_L_FS_T1_MNI',...
   'resliced_native_MVRe_R_FS_DTI_MNI',...
   'resliced_native_MVRe_R_FS_T1_MNI',...
   'resliced_native_Pf_L_FS_DTI_MNI',...
   'resliced_native_Pf_L_FS_T1_MNI',...
   'resliced_native_Pf_R_FS_DTI_MNI',...
   'resliced_native_Pf_R_FS_T1_MNI',...
   'resliced_native_PuA_L_FS_DTI_MNI',...
   'resliced_native_PuA_L_FS_T1_MNI',...
   'resliced_native_PuA_R_FS_DTI_MNI',...
   'resliced_native_PuA_R_FS_T1_MNI',...
   'resliced_native_PuI_L_FS_DTI_MNI',...
   'resliced_native_PuI_L_FS_T1_MNI',...
   'resliced_native_PuI_R_FS_DTI_MNI',...
   'resliced_native_PuI_R_FS_T1_MNI',...
   'resliced_native_PuL_L_FS_DTI_MNI',...
   'resliced_native_PuL_L_FS_T1_MNI',...
   'resliced_native_PuL_R_FS_DTI_MNI',...
   'resliced_native_PuL_R_FS_T1_MNI',...
   'resliced_native_PuM_L_FS_DTI_comb_MNI',...
   'resliced_native_PuM_L_FS_T1_MNI',...
   'resliced_native_PuM_R_FS_DTI_comb_MNI',...
   'resliced_native_PuM_R_FS_T1_MNI',...
   'resliced_native_PuMl_L_FS_DTI_MNI',...
   'resliced_native_PuMl_R_FS_DTI_MNI',...
   'resliced_native_PuMm_L_FS_DTI_MNI',...
   'resliced_native_PuMm_R_FS_DTI_MNI',...
   'resliced_native_Pu_L_FS_DTI_comb_MNI',...
   'resliced_native_Pu_L_FS_T1_comb_MNI',...
   'resliced_native_Pu_R_FS_DTI_comb_MNI',...
   'resliced_native_Pu_R_FS_T1_comb_MNI',...
   'resliced_native_Pul_L_HT_New_MNI',...
   'resliced_native_Pul_R_HT_New_MNI',...
   'resliced_native_VA_L_FS_DTI_MNI',...
   'resliced_native_VA_L_FS_DTI_comb_MNI',...
   'resliced_native_VA_L_FS_T1_MNI',...
   'resliced_native_VA_L_FS_T1_comb_MNI',...
   'resliced_native_VA_L_HT_New_MNI',...
   'resliced_native_VA_R_FS_DTI_MNI',...
   'resliced_native_VA_R_FS_DTI_comb_MNI',...
   'resliced_native_VA_R_FS_T1_MNI',...
   'resliced_native_VA_R_FS_T1_comb_MNI',...
   'resliced_native_VA_R_HT_New_MNI',...
   'resliced_native_VAmc_L_FS_DTI_MNI',...
   'resliced_native_VAmc_L_FS_T1_MNI',...
   'resliced_native_VAmc_R_FS_DTI_MNI',...
   'resliced_native_VAmc_R_FS_T1_MNI',...
   'resliced_native_VLP_L_HT_New_MNI',...
   'resliced_native_VLP_R_HT_New_MNI',...
   'resliced_native_VL_L_HT_New_MNI',...
   'resliced_native_VL_R_HT_New_MNI',...
   'resliced_native_VLa_L_FS_DTI_MNI',...
   'resliced_native_VLa_L_FS_T1_MNI',...
   'resliced_native_VLa_L_HT_New_MNI',...
   'resliced_native_VLa_R_FS_DTI_MNI',...
   'resliced_native_VLa_R_FS_T1_MNI',...
   'resliced_native_VLa_R_HT_New_MNI',...
   'resliced_native_VLp_L_FS_DTI_MNI',...
   'resliced_native_VLp_L_FS_T1_MNI',...
   'resliced_native_VLp_R_FS_DTI_MNI',...
   'resliced_native_VLp_R_FS_T1_MNI',...
   'resliced_native_VPL_L_FS_DTI_MNI',...
   'resliced_native_VPL_L_FS_T1_MNI',...
   'resliced_native_VPL_L_HT_New_MNI',...
   'resliced_native_VPL_R_FS_DTI_MNI',...
   'resliced_native_VPL_R_FS_T1_MNI',...
   'resliced_native_VPL_R_HT_New_MNI',...
   'resliced_native_VLPd_L_HT_New_MNI',...
   'resliced_native_VLPd_R_HT_New_MNI',...
   'resliced_native_VLPv_L_HT_New_MNI',...
   'resliced_native_VLPv_R_HT_New_MNI',...
   'resliced_native_V_comb_R_FS_DTI_MNI',...
   'resliced_native_V_comb_R_FS_T1_MNI',...
   'resliced_native_V_comb_R_HT_New_MNI',...
   'resliced_native_V_comb_L_FS_DTI_MNI',...
   'resliced_native_V_comb_L_FS_T1_MNI',...
   'resliced_native_V_comb_L_HT_New_MNI'
   };
% 
n_roi = numel(roi_names);
% List of subject folders
n_subj = numel(folder_list);

%%%% % Load one example file to get dimensions and header info
example_file = [directory '\' folder_list{1} '\' subfolder '\' roi_names{1} '.nii'];
V_ref = spm_vol(example_file);
dims = V_ref.dim;

%%% % Loop over each ROI
for r = 1:n_roi
    roi = [roi_names{r} '.nii'];
    roi_count = zeros(dims);  % Accumulate voxel-wise counts

    % Loop over each subject
    for s = 1:n_subj
        subj = folder_list{s};
        roi_file = [directory '\' subj '\' subfolder '\' roi];

        if ~isfile(roi_file)
            warning('Missing: %s', roi_file);
            continue;
        end

        V = spm_vol(roi_file);
        img = spm_read_vols(V);
        img = img > 0.5;  % Ensure binary

        roi_count = roi_count + img;
    end

    % Compute proportion map
    roi_prop = roi_count / n_subj;

    % Write result
    V_out = V_ref;
    V_out.fname = fullfile(probabilistic_maps_dir, sprintf('proportion_%s', roi));    
    spm_write_vol(V_out, roi_prop);

    fprintf('Saved: %s\n', V_out.fname);
end

fprintf('Done: all ROI proportion maps saved.\n');

for r = 1:numel(roi_names)
    roi = roi_names{r};
    prop_file = fullfile(probabilistic_maps_dir, ['proportion_', roi, '.nii']);

    if ~isfile(prop_file)
        warning('File not found: %s', prop_file);
        continue;
    end

    % Load proportion image
    V = spm_vol(prop_file);
    img = spm_read_vols(V);

    % Threshold and binarize
    binary_mask = img >= threshold;

    % Save binary ROI
    V_out = V;
    V_out.fname = fullfile(thresh_prob_maps_dir, ['thresh_', roi, '.nii']);
    spm_write_vol(V_out, binary_mask);

    fprintf('Saved thresholded probabilistic ROI: %s\n', V_out.fname);
end

fprintf('All thresholded binary ROIs saved.\n');

% %-----------------------------------------------
% % % % 7. Get Dice and AHD metrics for STAMI masks vs segmentations in MNI space
% % % %-------------------------------------------
% % % Get file lists
% 
STAMI_files_left = {
'r7TAMI_L_AV_MNI.nii',...
'r7TAMI_L_CL_MNI.nii',...
'r7TAMI_L_CM_MNI.nii',...
'r7TAMI_L_LGN_MNI.nii',...
'r7TAMI_L_MD_MNI.nii',...
'r7TAMI_L_MGN_MNI.nii',...
'r7TAMI_L_Pu_MNI.nii',...
'r7TAMI_L_V_MNI.nii'
};

STAMI_files_right = {
'r7TAMI_R_AV_MNI.nii',...
'r7TAMI_R_CL_MNI.nii',...
'r7TAMI_R_CM_MNI.nii',...
'r7TAMI_R_LGN_MNI.nii',...
'r7TAMI_R_MD_MNI.nii',...
'r7TAMI_R_MGN_MNI.nii',...
'r7TAMI_R_Pu_MNI.nii',...
'r7TAMI_R_V_MNI.nii'
};

fsdtifiles_left = {
'thresh_resliced_native_AV_L_FS_DTI_MNI.nii',...
'thresh_resliced_native_CL_L_FS_DTI_MNI.nii',...
'thresh_resliced_native_CM_L_FS_DTI_MNI.nii',...
'thresh_resliced_native_LGN_L_FS_DTI_MNI.nii',...
'thresh_resliced_native_MD_Pf_L_FS_DTI_comb_MNI.nii',...
'thresh_resliced_native_MGN_L_FS_DTI_MNI.nii',...
'thresh_resliced_native_Pu_L_FS_DTI_comb_MNI.nii',...
'thresh_resliced_native_V_comb_L_FS_DTI_MNI.nii'
};

fsdtifiles_right = {
    'thresh_resliced_native_AV_R_FS_DTI_MNI.nii',...
'thresh_resliced_native_CL_R_FS_DTI_MNI.nii',...
'thresh_resliced_native_CM_R_FS_DTI_MNI.nii',...
'thresh_resliced_native_LGN_R_FS_DTI_MNI.nii',...
'thresh_resliced_native_MD_Pf_R_FS_DTI_comb_MNI.nii',...
'thresh_resliced_native_MGN_R_FS_DTI_MNI.nii',...
'thresh_resliced_native_Pu_R_FS_DTI_comb_MNI.nii',...
'thresh_resliced_native_V_comb_R_FS_DTI_MNI.nii'
};


fst1files_left = {
'thresh_resliced_native_AV_L_FS_T1_MNI.nii',...
'thresh_resliced_native_CL_L_FS_T1_MNI.nii',...
'thresh_resliced_native_CM_L_FS_T1_MNI.nii',...
'thresh_resliced_native_LGN_L_FS_T1_MNI.nii',...
'thresh_resliced_native_MD_Pf_L_FS_T1_comb_MNI.nii',...
'thresh_resliced_native_MGN_L_FS_T1_MNI.nii',...
'thresh_resliced_native_Pu_L_FS_T1_comb_MNI.nii',...
'thresh_resliced_native_V_comb_L_FS_T1_MNI.nii'
};

fst1files_right = {
'thresh_resliced_native_AV_R_FS_T1_MNI.nii',...
'thresh_resliced_native_CL_R_FS_T1_MNI.nii',...
'thresh_resliced_native_CM_R_FS_T1_MNI.nii',...
'thresh_resliced_native_LGN_R_FS_T1_MNI.nii',...
'thresh_resliced_native_MD_Pf_R_FS_T1_comb_MNI.nii',...
'thresh_resliced_native_MGN_R_FS_T1_MNI.nii',...
'thresh_resliced_native_Pu_R_FS_T1_comb_MNI.nii',...
'thresh_resliced_native_V_comb_R_FS_T1_MNI.nii'
};


htfiles_left = {
'thresh_resliced_native_AV_L_HT_New_MNI.nii',...
'thresh_resliced_native_CL_L_HT_New_MNI.nii',...
'thresh_resliced_native_CM_L_HT_New_MNI.nii',...
'thresh_resliced_native_LGN_L_HT_New_MNI.nii',...
'thresh_resliced_native_MD_Pf_L_HT_New_MNI.nii',...
'thresh_resliced_native_MGN_L_HT_New_MNI.nii',...
'thresh_resliced_native_Pul_L_HT_New_MNI.nii',...
'thresh_resliced_native_V_comb_L_HT_New_MNI.nii'
    };

htfiles_right = {
'thresh_resliced_native_AV_R_HT_New_MNI.nii',...
'thresh_resliced_native_CL_R_HT_New_MNI.nii',...
'thresh_resliced_native_CM_R_HT_New_MNI.nii',...
'thresh_resliced_native_LGN_R_HT_New_MNI.nii',...
'thresh_resliced_native_MD_Pf_R_HT_New_MNI.nii',...
'thresh_resliced_native_MGN_R_HT_New_MNI.nii',...
'thresh_resliced_native_Pul_R_HT_New_MNI.nii',...
'thresh_resliced_native_V_comb_R_HT_New_MNI.nii'
    };

roiFiles1set = {STAMI_files_left, STAMI_files_right};
roiFiles2set = {fsdtifiles_left, fsdtifiles_right, fst1files_left, fst1files_right, htfiles_left, htfiles_right};

roiFiles1suffix = {'STAMI_L', 'STAMI_R'};
roiFiles2suffix = {'fsdti_L', 'fsdti_R', 'fst1_L', 'fst1_R', 'ht_L', 'ht_R'};

for m = 1:length(roiFiles1set)
    for l = 1:length(roiFiles2set)
        roiFiles1 = cellfun(@(f) fullfile(STAMI_directory, f), roiFiles1set{m}, 'UniformOutput', false);
        roiFiles2 = cellfun(@(f) fullfile(thresh_prob_maps_dir, f), roiFiles2set{l}, 'UniformOutput', false);

        n1 = numel(roiFiles1);
        n2 = numel(roiFiles2);

        diceMatrix = zeros(n1, n2);
        hausdorffMatrix = zeros(n1, n2);

[~, labels1, ext1] = cellfun(@fileparts, roiFiles1, 'UniformOutput', false);
labels1 = strcat(labels1, ext1);  % get full file names like 'L_AV.nii'


[~, labels2, ext2] = cellfun(@fileparts, roiFiles2, 'UniformOutput', false);
labels2 = strcat(labels2, ext2);

        for i = 1:n1
            V1 = spm_vol(roiFiles1{i});
            img1 = spm_read_vols(V1) > 0;

            for j = 1:n2
                fprintf('Processing: %s <-> %s\n', roiFiles1{i}, roiFiles2{j});
                V2 = spm_vol(roiFiles2{j});
                img2 = spm_read_vols(V2) > 0;

                % %%%Reslice img2 to match img1 if needed
                if any(abs(V1.mat(:) - V2.mat(:)) > 1e-4) || any(size(img1) ~= size(img2))
                    img2 = reslice_and_binarize(V1, V2);
                end

                % %%%Dice coefficient
                intersection = sum(img1(:) & img2(:));
                diceMatrix(i,j) = 2 * intersection / (sum(img1(:)) + sum(img2(:)));

                %%%% Hausdorff (average)

%%%Step 1: Get voxel coordinates of non-zero (mask) voxels
                idx1 = find(img1); %%%%% linear indices of non-zero voxels in img1
                idx2 = find(img2); %%%%%% linear indices of non-zero voxels in img2
                if isempty(idx1) || isempty(idx2)
                    hausdorffMatrix(i,j) = NaN;
                    continue;
                end

%%%%%%%%%%Step 2: Convert to subscripts (x,y,z coordinates)
                [x1,y1,z1] = ind2sub(size(img1), idx1);
                [x2,y2,z2] = ind2sub(size(img2), idx2);
                coords1 = [x1,y1,z1];
                coords2 = [x2,y2,z2];

              %%%  % Compute distances
                d12 = zeros(size(coords1,1), 1);
%%%%%Step 3: Compute average directed Hausdorff distance from 1 → 2
                for k = 1:size(coords1,1)
                    diffs = coords2 - coords1(k,:);
                    dists = sqrt(sum(diffs.^2, 2));
                    d12(k) = min(dists); %%%For each voxel in ROI1, find the closest voxel in ROI2 and store the distance.
                end
%%%%%%%%%Step 4: Same for 2 → 1
                d21 = zeros(size(coords2,1), 1);
                for k = 1:size(coords2,1)
                    diffs = coords1 - coords2(k,:);
                    dists = sqrt(sum(diffs.^2, 2));
                    d21(k) = min(dists);
                end
%%%%%%%%%%%Step 5: Hausdorff distance
%%%%%%%%hausdorffMatrix(i,j) stores the AHD between roiFiles1{i} and roiFiles2{j}.
                hausdorffMatrix(i,j) = max(mean(d12),mean(d21));

            end
        end

        % Define unique Excel file names for each pair
        xlsxname_dice = [analysis_output_directory '\' roiFiles1suffix{m} '_' roiFiles2suffix{l} '_Dice_' MNI_vs_STAMI_analysis_name];
        xlsxname_hausdorff = [analysis_output_directory '\' roiFiles1suffix{m} '_' roiFiles2suffix{l} '_AHD_' MNI_vs_STAMI_analysis_name];


        % Write Dice results to dice_results.xlsx
        xlswrite(xlsxname_dice, [{'ROI1/ROI2'}, labels2], 'Dice', 'A1');
        xlswrite(xlsxname_dice, [labels1', num2cell(diceMatrix)], 'Dice', 'A2');

        % Write Hausdorff results to hausdorff_results.xlsx
        xlswrite(xlsxname_hausdorff, [{'ROI1/ROI2'}, labels2], 'Hausdorff', 'A1');
        xlswrite(xlsxname_hausdorff, [labels1', num2cell(hausdorffMatrix)], 'Hausdorff', 'A2');

        disp(['Saved to ' xlsxname_dice 'and' xlsxname_hausdorff]);

    end
end
% % % % %%% ---------------------------------------------------------

%%%----------------------------------
% 10. Calculate and save Dice and AHD in native space
%%---------------------------------------------
%%-%%-%%-%%-%%-Define a list of file pairs (each pair is a 1x2 cell array)
filePairs = {
  {'r7TAMI_L_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_L_FS_DTI.nii'};
  {'r7TAMI_L_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_L_FS_T1.nii'};
  {'r7TAMI_L_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_L_HT_New.nii'};
  {'r7TAMI_L_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_L_FS_DTI.nii'};
  {'r7TAMI_L_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_L_FS_T1.nii'};
  {'r7TAMI_L_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_L_HT_New.nii'};
  {'r7TAMI_L_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_L_FS_DTI.nii'};
  {'r7TAMI_L_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_L_FS_T1.nii'};
  {'r7TAMI_L_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_L_HT_New.nii'};
  {'r7TAMI_L_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_L_FS_DTI.nii'};
  {'r7TAMI_L_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_L_FS_T1.nii'};
  {'r7TAMI_L_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_L_HT_New.nii'};
  {'r7TAMI_L_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_L_FS_DTI.nii'};
  {'r7TAMI_L_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_L_FS_T1.nii'};
  {'r7TAMI_L_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_L_HT_New.nii'};
  {'r7TAMI_L_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_L_FS_DTI_comb.nii'};
  {'r7TAMI_L_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_L_FS_T1_comb.nii'};
  {'r7TAMI_L_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_L_HT_New.nii'};
  {'r7TAMI_L_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pu_L_FS_DTI_comb.nii'};
  {'r7TAMI_L_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pu_L_FS_T1_comb.nii'};
  {'r7TAMI_L_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pul_L_HT_New.nii'};
  {'r7TAMI_L_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_L_HT_New.nii'};
  {'r7TAMI_L_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_L_FS_T1.nii'};
  {'r7TAMI_L_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_L_FS_DTI.nii'};
  {'r7TAMI_R_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_R_HT_New.nii'};
  {'r7TAMI_R_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_R_FS_T1.nii'}; 
  {'r7TAMI_R_V_MNI_7TAMI_native.nii', 'resliced_native_V_comb_R_FS_DTI.nii'};
  {'r7TAMI_R_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_R_FS_DTI.nii'};
  {'r7TAMI_R_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_R_FS_T1.nii'};
  {'r7TAMI_R_AV_MNI_7TAMI_native.nii', 'resliced_native_AV_R_HT_New.nii'};
  {'r7TAMI_R_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_R_FS_DTI.nii'};
  {'r7TAMI_R_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_R_FS_T1.nii'};
  {'r7TAMI_R_CL_MNI_7TAMI_native.nii', 'resliced_native_CL_R_HT_New.nii'};
  {'r7TAMI_R_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_R_FS_DTI.nii'};
  {'r7TAMI_R_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_R_FS_T1.nii'};
  {'r7TAMI_R_CM_MNI_7TAMI_native.nii', 'resliced_native_CM_R_HT_New.nii'};
  {'r7TAMI_R_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_R_FS_DTI.nii'};
  {'r7TAMI_R_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_R_FS_T1.nii'};
  {'r7TAMI_R_LGN_MNI_7TAMI_native.nii', 'resliced_native_LGN_R_HT_New.nii'};
  {'r7TAMI_R_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_R_FS_DTI.nii'};
  {'r7TAMI_R_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_R_FS_T1.nii'};
  {'r7TAMI_R_MGN_MNI_7TAMI_native.nii', 'resliced_native_MGN_R_HT_New.nii'};
  {'r7TAMI_R_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_R_FS_DTI_comb.nii'};
  {'r7TAMI_R_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_R_FS_T1_comb.nii'};
  {'r7TAMI_R_MD_MNI_7TAMI_native.nii', 'resliced_native_MD_Pf_R_HT_New.nii'};
  {'r7TAMI_R_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pu_R_FS_DTI_comb.nii'};
  {'r7TAMI_R_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pu_R_FS_T1_comb.nii'};
  {'r7TAMI_R_Pu_MNI_7TAMI_native.nii', 'resliced_native_Pul_R_HT_New.nii'};
 };

% Loop through each pair
for i = 1:length(filePairs)
    % Define the names of the two ROI files
    roiName1 = filePairs{i}{1};
    roiName2 = filePairs{i}{2};
    % Example processing
    fprintf('Processing pair %d:\n', i);
    fprintf('  File 1: %s\n', roiName1);
    fprintf('  File 2: %s\n', roiName2);
    % Initialize a cell array to store the results
    results = {};

    for j = 1:length(folder_list)
        folderName = folder_list{j};
        folderPath1 = [directory '\' folder_list{j} '\' subfolder]; % the subject's normalised rois
        folderPath2 = [directory '\' folder_list{j} '\' subfolder]; % the location of the STAMI Atlas ROIs (resliced in native space)
        cd(folderPath1); % go there, in case it otherwise writes stuff elsewhere.
        %---------------------------------------------------------------------------------------------
        % get forward and backward deformation fields
        fprintf('processing %s\n', folder_list{j});

        % Define the full paths to the two ROI files
        roiFilePath1 = fullfile(folderPath1, roiName1);
        roiFilePath2 = fullfile(folderPath2, roiName2);

        % Check if both ROI files exist
        if exist(roiFilePath1, 'file') && exist(roiFilePath2, 'file')
            % Read the .nii files using SPM
            V1 = spm_vol(roiFilePath1);
            Y1 = spm_read_vols(V1);

            V2 = spm_vol(roiFilePath2);
            Y2 = spm_read_vols(V2);

            % Change data type (example: float32)
            V1.dt = [16 0];  % 16 = float32, 0 = no scaling
            V2.dt = [16 0];  % 16 = float32, 0 = no scaling
            % Change scaling to identity, so that '''Intensity: Y = 3.05185e-05 X' changes to 'Intensity: Y = 1 X'
            V1.pinfo(1:2,:) = [1; 0];  % scl_slope = 1, scl_inter = 0 
            V2.pinfo(1:2,:) = [1; 0];  % scl_slope = 1, scl_inter = 0 
            % Binarize in case of partial volume values with anything below
            % 0.5
            Y1 = Y1 > 0.5;
            Y2 = Y2 > 0.5;
            % Save the new binary mask
            V1.fname = roiName1; % <-- rename with same mask name.
            V2.fname = roiName2; % <-- rename with same mask name.
            spm_write_vol(V1, Y1);
            spm_write_vol(V2, Y2);
            % Calculate the number of voxels in each ROI
            numVoxels1 = sum(Y1(:));
            numVoxels2 = sum(Y2(:));

            % Calculate the Dice coefficient
            intersection = sum((Y1(:) > 0) & (Y2(:) > 0));
            diceCoefficient = (2 * intersection) / (numVoxels1 + numVoxels2);

            % Get voxel coordinates for each ROI
            [x1, y1, z1] = ind2sub(size(Y1), find(Y1));
            [x2, y2, z2] = ind2sub(size(Y2), find(Y2));
            % Convert to real-world (mm) coordinates
            coords1 = V1.mat * [x1'; y1'; z1'; ones(1, numel(x1))];
            coords2 = V2.mat * [x2'; y2'; z2'; ones(1, numel(x2))];

            coords1 = coords1(1:3, :)';  % Nx3
            coords2 = coords2(1:3, :)';  % Mx3

            % Manual pairwise Euclidean distances between coords1 (N×3) and coords2 (M×3)
            N = size(coords1, 1);
            M = size(coords2, 1);

            D = zeros(N, M);
            for j = 1:N
                diffs = coords2 - coords1(j, :);         % M×3
                D(j, :) = sqrt(sum(diffs.^2, 2))';       % Row vector: distances from coords1(i,:) to all coords2
            end

            %%%%% Average of minimal distances (both directions)
            d1 = mean(min(D, [], 2));  % From ROI1 to ROI2
            d2 = mean(min(D, [], 1));  % From ROI2 to ROI1

            %%%% AHD
            hausdorff = max(d1,d2);

            % clean up names
            roiName1_clean = erase(roiName1, {'resliced_native_w', '.nii'});
            roiName2_clean = erase(roiName2, {'resliced_native_', '.nii'});

            % Store the results in the results array
            results{end+1, 1} = folderName;
            results{end, 2} = roiName1_clean;
            results{end, 3} = roiName2_clean;
            results{end, 4} = diceCoefficient;
            results{end, 5} = hausdorff;
            results{end, 6} = numVoxels1;
            results{end, 7} = numVoxels2;
        else
            disp(['ROI files not found for subject: ', folderName]);
        end
    end

    % Convert the results cell array to a table
    resultsTable = cell2table(results, 'VariableNames', {'Subject', 'STAMI_Atlas_ROI', 'Segmentation', 'DICE', 'AHD', 'NumVoxelsROI1', 'NumVoxelsROI2'});

    % Define the output CSV file path
    outputCsvFile = fullfile(analysis_output_directory, Dice_AHD_Native_FileName);

    % Check if file exists to determine whether to write header
    writeHeader = ~isfile(outputCsvFile);

    % Write or append table
    writetable(resultsTable, outputCsvFile, 'WriteMode', 'append', 'WriteVariableNames', writeHeader);

disp(['Results saved to ', outputCsvFile]);
end

% -------- Helper function --------
function resliced = reslice_and_binarize(Vref, Vmov)
    flags = struct('interp', 0, 'wrap', [0 0 0], 'mask', 0, 'mean', 0, 'which', 1, 'prefix', 'r');
    spm_reslice_quiet({Vref.fname, Vmov.fname}, flags);
    [~, name, ext] = fileparts(Vmov.fname);
    reslicedVol = spm_vol(['r' name ext]);
    resliced = spm_read_vols(reslicedVol) > 0;
    delete(['r' name ext]);  % clean up
end