%ANALYZEBONES checks & analyzes the surface models of the VSD bones.
%
% AUTHOR: Maximilian C. M. Fischer
% COPYRIGHT (C) 2023 Maximilian C. M. Fischer
% LICENSE: EUPL v1.2
%

clearvars; close all

addpath(genpath('..\src'))
VSD_addPathes('..\..\..\..\..\')

% Load subjects & meta data
subjectXLSX = '..\res\VSD_Subjects.xlsx';
Subjects = readtable(subjectXLSX);
NoS = size(Subjects, 1);

load(['..\..\Bones\' Subjects.ID{1} '.mat'], 'B')
NoB = length(B);
boneNames = {B.name};
clear B

%% Mesh sanity checks
NoCC = nan(NoS, NoB);
for s=1:NoS
    % Import the bones
    load(['..\..\Bones\' Subjects.ID{s} '.mat'], 'B')
    for b=1:length(B)
        if ~isempty(B(b).mesh)
            stats = statistics(B(b).mesh.vertices, B(b).mesh.faces, 'MinArea', 1e-8);
            sanityStats = rmfield(stats, {'num_faces', 'num_vertices', 'num_edges', ...
                'num_connected_components', 'num_handles', 'euler_characteristic'});
            if ~all(full(cell2mat(struct2cell(sanityStats)))==0)
                warning(['Failed mesh check for subject ' ...
                    Subjects.ID{s} ' (' num2str(s) '): ' B(b).name ' (' num2str(b) ')'])
            end
            NoCC(s,b) = stats.num_connected_components;
        end
    end
end

save('NumberOfConnComp.mat','NoCC')

% Check for intersections of adjacent bones
for s=1:NoS
    % Import the bones
    load(['..\..\Bones\' Subjects.ID{s} '.mat'], 'B')
    mesh = concatenateMeshes([B(1:end).mesh]);
    stats = statistics(mesh.vertices, mesh.faces);
    if stats.num_selfintersecting_pairs ~= 0
        error([Subjects.ID{s} ' has intersections between adjacent bones!'])
    end
end

%% Mesh volume
% Remove incomplete or inconsistent subjects
Subjects = Subjects(cellfun(@(x) isempty(strfind(lower(x),'cut off')), Subjects.Comment),:); %#ok<STREMP>
Subjects = Subjects(cellfun(@(x) isempty(strfind(lower(x),'total knee')), Subjects.Comment),:); %#ok<STREMP>
Subjects = Subjects(cellfun(@(x) isempty(strfind(lower(x),'gender')), Subjects.Comment),:); %#ok<STREMP>
Subjects = Subjects(cellfun(@(x) isempty(strfind(lower(x),'conflicting')), Subjects.Comment),:); %#ok<STREMP>
Subjects = Subjects(~isnan(Subjects.Weight),:);
% Calculate volume of the bone models
NoS = size(Subjects, 1);
volume = nan(NoS, NoB);
for s=1:NoS
    % Import the bones
    load(['..\..\Bones\' Subjects.ID{s} '.mat'], 'B')
    for b=1:length(B)
        volume(s,b) = VolumeIntegrate(B(b).mesh.vertices,B(b).mesh.faces);
        % volume(s,b) = sum(arrayfun(@(x) meshVolume(x), splitMesh(B(b).mesh)));
    end
end

cmm3toccm3 = 0.001;
volumeTable = cell2table(cell(NoB,13), 'VariableNames', {...
    'Bone name', 'Min.', 'Mean', 'Median', 'Max.',...
    'M_Min.', 'M_Mean', 'M_Median', 'M_Max.',...
    'F_Min.', 'F_Mean', 'F_Median', 'F_Max.'});
[minIdx, maxIdx] = deal(nan(NoS, 1));
mIdx = strcmp(Subjects.Sex, 'M');
assert(sum(mIdx) == 10)
fIdx = strcmp(Subjects.Sex, 'F');
assert(sum(fIdx) == 10)
assert(sum(fIdx & mIdx) == 0)
for b=1:NoB
    % All
    volumeTable(b,1) = boneNames(b);
    [minValue, minIdx(b)] = min(volume(:,b)*cmm3toccm3);
    volumeTable{b,2} = {minValue};
    volumeTable(b,3) = meanStats(volume(:,b)*cmm3toccm3,'% 1.0f','format','short');
    volumeTable(b,4) = medianStats(volume(:,b)*cmm3toccm3,'% 1.0f','format','short');
    [maxValue, maxIdx(b)] = max(volume(:,b)*cmm3toccm3);
    volumeTable{b,5} = {maxValue};
    % Male
    [minValue, minIdx(b)] = min(volume(mIdx,b)*cmm3toccm3);
    volumeTable{b,6} = {minValue};
    volumeTable(b,7) = meanStats(volume(mIdx,b)*cmm3toccm3,'% 1.0f','format','short');
    volumeTable(b,8) = medianStats(volume(mIdx,b)*cmm3toccm3,'% 1.0f','format','short');
    [maxValue, maxIdx(b)] = max(volume(mIdx,b)*cmm3toccm3);
    volumeTable{b,9} = {maxValue};
    % Female
    [minValue, minIdx(b)] = min(volume(fIdx,b)*cmm3toccm3);
    volumeTable{b,10} = {minValue};
    volumeTable(b,11) = meanStats(volume(fIdx,b)*cmm3toccm3,'% 1.0f','format','short');
    volumeTable(b,12) = medianStats(volume(fIdx,b)*cmm3toccm3,'% 1.0f','format','short');
    [maxValue, maxIdx(b)] = max(volume(fIdx,b)*cmm3toccm3);
    volumeTable{b,13} = {maxValue};
end

writetable(volumeTable, 'volumeResults.xlsx', 'Sheet','Volume', 'Range','B5',...
    'WriteVariableNames',0, 'WriteRowNames',0)

figure('color','w','numbertitle','off','name','Bone Volume')
bph = boxplot(volume*cmm3toccm3,boneNames,'LabelOrientation','inline');
set(findobj(get(bph(1), 'parent'), 'type', 'text'), 'interpreter','tex');
ylabel('[cm³]')

%% Age, weight, height
SubjectStats(1,1) = meanStats(Subjects.Age,'% 1.0f','format','short');
SubjectStats(1,2) = meanStats(Subjects.Weight,'% 1.0f','format','short');
SubjectStats(1,3) = meanStats(Subjects.Height,'% 1.1f','format','short');
SubjectStats(2,1) = medianStats(Subjects.Age,'% 1.0f','format','short');
SubjectStats(2,2) = medianStats(Subjects.Weight,'% 1.0f','format','short');
SubjectStats(2,3) = medianStats(Subjects.Height,'% 1.1f','format','short');

%% Number of vertices
NoV = nan(NoS, NoB);
for s=1:NoS
    % Import the bones
    load(['..\..\Bones\' Subjects.ID{s} '.mat'], 'B')
    for b=1:length(B)
        NoV(s,b) = size(B(b).mesh.vertices,1);
    end
end

figure('color','w','numbertitle','off','name','Number of mesh vertices')
bph = boxplot(NoV,boneNames,'LabelOrientation','inline');
set(findobj(get(bph(1), 'parent'), 'type', 'text'), 'interpreter','tex');

% [List.f, List.p] = matlab.codetools.requiredFilesAndProducts([mfilename '.m']);
% List.f = List.f'; List.p = List.p';