## 2 Materials and Methods
### Creation of a PCA-Based Statistical Shape Model  
A Statistical Shape Model (SSM) allows the prediction of a general shape as the sum of an average shape and a certain number of \"'transformations\" (the so-called PCs, or modes), tuned by the respective weights:  
\[G_{\text{CNR}}=G_{\text{CNR,avg}}+\sum_{i=1}^{M}W_{i}\cdot PC_{i} \tag{6}\]  
where: \(G_{\text{CNR,avg}}\) is an average geometry, described by a vector of \(N\times 3\) components, where \(N\) is the total number of mesh nodes; \(PC_{i}\) is the \(i\)th principal component (or mode); \(W_{i}\) is a scalar which represents the weight to be assigned to the \(i\)th PC. In this formulation, the average shape provides the mean anatomical information of the considered geometry (bones in the present case); the average shape can be obtained from all geometries belonging to the training set once they have all been scaled to be the same size, centred (i.e., their centres of mass are coincident) and aligned. The PCs represent the main (i.e., the most variable) directions of deviation from the average shape highlighted by the analyzed cohort. The linear combination of these PCs through properly computed weights generates a new shape, belonging to the same \"family\", i.e., to the same statistical shape distribution. The average shape (\(G_{\text{CNR,avg}}\)) can be evaluated through the Generalized Procrustes Analysis (GPA), while the variability model \(W_{i}\)-\(PC_{i}\) is actually defined through PCA. Both procedures require: (1) the definition of landmarks able to describe the 3D shape (correspondence nodes) on all the shapes in the training set; (2) having obtained isotopological meshes through the mesh morphing technique. In the following analyses, the whole node set of each mesh was used as the landmark set. In the following, a 3D geometry \(G_{i}\) will be represented by a \(N\times 3\) matrix, where \(x_{i,j}\), \(y_{i,j}\), and \(z_{i,j}\) are the coordinates of node \(j\) (\(j=1\dots\) N) belonging to shape \(i\) (\(i=1\dots\) M). As a first step, given a set of \(M\) geometries \(G_{i}\) (\(i=1\dots\) M) representing the same shape (either the mandible or the proximal femur in this case), the data were pre-aligned and scaled in order to remove rotation, translation, and to achieve uniform scaling. This was performed using GPA [19], which consisted of finding the transformation \(T_{i}\) so that each shape \(G_{i}\) best aligned with a target mean shape (\(G_{\text{CNR,avg}}\)). The \"best\" alignment was defined in terms of the minimized sum of squared distances over the nodes in the coordinate frame between the aligned shape and the average shape. Data translation could be performed straightforwardly calculating for each \(i\)th shape, the average of \(x\), \(y\), and \(z\) coordinates (that is the centre of gravity of the point cloud) and subtracting its value from the respective node coordinates:  
\[\overline{x_{i}}=\sum_{j=1}^{N}\frac{x_{i,j}}{N};\overline{y_{i}}=\sum_{j=1}^{ N}\frac{y_{i,j}}{N};\overline{z_{i}}=\sum_{j=1}^{N}\frac{z_{i,j}}{N} \tag{7}\]  
\[\left(x_{i,j},y_{i,j},z_{i,j}\right)_{\text{C}}=\left(x_{i,j},y_{i,j},z_{i,j} \right)-\left(\overline{x_{i}},\overline{y_{i}},\overline{z_{i}}\right) \tag{8}\]  
The new node coordinates were listed into a new centered matrix \(G_{i,c}\) whose centroid was located at \((0,0,0)\). As a second step, the matrices were scaled according to a scale factor \(S_{i}\), obtained by dividing each matrix by the respective norm, and multiplying the result by the average norm \(Norm_{\text{avg}}\):  
\[Norm_{\text{avg}}=\sum_{j=1}^{N}\frac{\text{Norm}(G_{i,C})}{N};\qquad\quad G_{ i,CN}=\frac{G_{i,C}}{\text{Norm}(G_{i,C})}Norm_{\text{avg}}=S_{i}\cdot G_{i,C} \tag{9}\]  
Eventually, an iterative procedure was started:  
All matrices were averaged:  
\[G_{\text{CN,avg}}=\sum_{i=1}^{M}\frac{G_{CN}}{M} \tag{10}\]  
in order to provide an initial estimate of the average shape.