--- a
+++ b/MATLAB/analysis/analyzeBones.m
@@ -0,0 +1,141 @@
+%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';
\ No newline at end of file