Statistical shape modelling of the human mandible: 3D shape predictions based on external morphometric features

Abstract

One of the main limitations in subject-centred design approach is represented by getting 3D models of the region of interest. Indeed, 3D reconstruction from imaging data (i.e., computed tomography scans) is expensive and exposes the subject to high radiation doses. Statistical Shape Models (SSMs) are mathematical models able to describe the variability associated to a population and allow predicting new shapes tuning model parameters. These parameters almost never have a physical meaning and so they cannot be directly related to morphometric features. In this study a gender-combined SSM model of the human mandible was setup, using Generalised Procrustes Analysis and Principal Component Analysis on a dataset of fifty mandibles. Twelve morphometric features, able to characterise the mandibular bone and readily collectable during external examinations, were recorded and correlated to SSM parameters by a multiple linear regression approach. Then a cross-validation procedure was performed on a control set to determine the combination of features able to minimise the average deviation between real and predicted shapes. Compactness of the SSM and main modes of deformations have been investigated and results consistent with previous works involving a higher number of shapes were found. A combination of five features was proved to characterise predicted shapes minimising the average error. As completion of the work, a male SSM was developed and performances compared with those of the combined SSM. The features-based model here proposed could represent a useful and easy-to-use tool for the generation of 3D customised models within a virtual interactive design environment.

1 Introduction

Statistical Shape Models (SSMs) represent a valuable mathematical tool to describe and understand the variability that characterises a specific population. SSMs allow identifying the principal modes of deformation, that are the main variations, associated to shapes within a population and, adapting the shape parameters of the model, new realistic shapes can be predicted. This approach is currently having a great diffusion, especially in the medical [1,2,3], ergonomics [4,5,6] and product design [7,8,9] fields, where researches are becoming more and more subject-centred. With reference to the medical field, the classical methodology implemented to customise the design of patient-specific devices or surgical procedures is based on an interactive design approach [10] and involves three main steps [11,12,13]: reverse engineering to obtain a 3D model of the anatomical region of interest and to setup a virtual working environment; 3D modelling within the virtual environment where knowledges belonging to different fields (i.e., medical, engineering, biomechanical) participate in the definition of design parameters; simulation within the interactive environment and, if necessary, prototyping with additive manufacturing techniques. A key point for the implementation of this kind of approach, is the availability of the 3D model replicating the geometry of interest: indeed, precise and detailed information on the investigated anatomical region are mandatory in order to carry out design and simulation within the virtual environment. When the investigated structures are inner body regions (i.e., bones) the classical approach to obtain models is the use of imaging techniques, such as Computed Tomography (CT), Cone Beam CT (CBCT) or Magnetic Resonance Imaging (MRI) [1, 14]. The main shortcoming related to the use of these systems is that they are expensive and/or expose patients to high-radiation doses; for these reasons, generating 3D models from these techniques is not always a feasible option.

Statistical shape models have been used for the mathematical representation of anatomical structures with different applications [15], including image automatic segmentation [16, 17], identification of morphological differences (i.e., discrimination between healthy and pathological structures [18, 19] or sex discrimination [20, 21]), generation of predictive models able to reproduce realistic or real geometries [22,23,24].

In this study an SSM for the prediction of mandibular bone external geometry was setup, basing the mandibular shape prediction on few external morphometric measurements. The approach here described is wanted to overcome difficulties that have to be faced in reverse engineering of this bone (i.e., unavailability of CT scans), providing a numerical tool able to quickly generate an accurate 3D subject-specific model that can be easily integrable within a virtual environment for numerical analyses and customised design purposes. Some previous works have implemented a statistical shape model approach for the analysis of the mandible geometry variability related to gender (male or female) [25] or age (growing mandible changes) [14]. Even when the SSM model was used to predict new mandibular shapes [25,26,27], these models still involve imaging techniques for fitting principal modes weights, using total/partial bone geometries or hard-tissue landmarks to define SSM parameters.

Here a different approach is proposed, having integrated the statistical shape analysis with a morphometric features analysis, following a methodology used for full body statistical shape modelling [6, 9, 28, 29]. This approach is based on identifying a correlation between the deformation modes of the population and selected shape’s features, that are shape’s characteristic with a physical meaning.

A total of fifty mandible’s CT scans were available: forty were used as training set for statistical analyses and ten as control set for model validation and features selection. In order to consider as much mandibles as possible, the training set was a gender-combined dataset, composed by both male and female mandibles. The correspondence problem was here solved using a mesh morphing approach based on Radial Basis Functions (RBFs); this approach allows obtaining iso-topological meshes with the advantage of using all the mesh nodes as correspondence points for statistical analyses. The Generalised Procrustes Analysis (GPA) was applied in order to filtered out translational, rotational and scaling effects from the population and to estimate the average shape geometry of the training set. Main modes of variations of the mandibular shape have been identified through Principal Component Analysis (PCA), extrapolating Principal Components (PCs) and so defining the variability model. Statistical shape model performances have been evaluated using a compactness metric, and main modes of deformation associated to first PCs have been investigated and compared with results obtained in previous studies, where a higher number of mandibles were used.

The SSM model was then modified so that new predicted shapes can be obtained linearly combining PCs on the basis of morphometric features provided as input. In this study, features are represented by mandibular relevant measurements that can be recorded during external examinations with classical instruments (i.e., calliper, goniometer). A detailed literature research was performed in order to identify commonly used facial measurements related to the mandible shape and a set of twelve features was selected. A cross-validation approach was used on the control set to establish the combination of these features that, if provided as input to the SSM, is able to predict the specific mandible geometry of a subject with the minimum deviation from the real one.

The prediction ability of the feature-based SSM model was assessed comparing predicted and real geometries (that are geometries reconstructed from CT scans) in terms of average deviations between shapes. As completion of this study, also a gender-based SSM was created repeating the same statistical analyses and feature selection procedure on a reduced dataset containing only male mandibles, and its performances were compared with those of the combined model. Finally, the combined SSM was also tested for its ability to predict missing or unhealthy bone’s regions in mandibles with a peculiar geometry.

In this work a limited number of mandibles were available and for this reason this research was aimed to provide preliminary results on the statistical shape modelling of the mandibular bone.

2 Materials and methods

The final aim of the approach presented in this study was to generate a ready-to-use virtual tool that requires as input a set of few external morphometric measurements, recorded by the final user on a specific subject, and provides as output the 3D model of the mandibular geometry. The whole process for the statistical shape model generation is summarised in Fig. 1, highlighting main steps described below in this section. Moreover, the integration of this tool within the design process is also shown.

Fig. 1

Statistical shape model generation workflow and its integration within the subject-specific design process

All the analyses presented in the following have been performed in MATLAB (v. 17, The MathWorks, Inc., Natick, MA, USA) through codes written by the author.

2.1 Mandible samples

A total of fifty Computed Tomography (CT) scans of adult skulls have been collected from the hospital ‘Azienda Ospedaliero Universitario Policlinico G. Rodolico-San Marco’ (Catania, Italy); they were obtained during different medical examinations using a Revolution HD scanner (GE Healthcare S.r.l., 0.5 mm slice thickness, 0.43 × 0.43 pixel spacing, 512 × 512 pixels, voltage KVP 120 kV and exposure time 500 ms), and made available for this study after getting approval from patients. All the mandibular bones did not undergo surgical treatments (absence of prosthetic devices or mandibular plates was assessed) and were not affected by pathological features (i.e., bone resorption, partial bone geometries, mandibular asymmetry). Demographics associated to the sample set are reported in Table 1.

Table 1 Demographics of the mandibular dataset

Three-dimensional models of the mandible external geometry have been reconstructed from CT scans following a common procedure for segmentation and post-processing [3, 25, 26, 30,31,32]. 3D model reconstruction was performed with Mimics and 3-Matic softwares (Materialise Inc., Belgium) by an expert operator: this, along with the high resolution of CT scan data, ensured a good level of reconstruction accuracy (sub-voxel mesh accuracy) transforming 2D CT scans into realistic 3D models of exact geometry, as reported in works that have dealt with this topic [33,34,35,36].

Main steps of this procedure are here summarised and corresponding results are shown in Fig. 2:

• a. Automatic segmentation based on thresholding Soft and hard tissues are separated setting the Hounsfield (HU) window for bone between 570:3070 HU [37, 38]; a collection of the pixels of interest (that is a mask) corresponding to hard tissues was identified by the software using a marching square algorithm to threshold and segment the regions of interest [36, 39, 40]. As a result, the whole skull was obtained (Fig. 2a)

• b. Semi-automatic segmentation Bone regions corresponding only to the skull, mandibular bone and teeth can be clearly identified and a ‘mask split’ tool (Mimics) was used to separate them. Portions of regions corresponding to these parts were selected by the operator on CT images and then the software automatically distinguishes between different bodies and assigns each region to different masks. Figure 2b, c show the mandible separated from the skull and the mandible without teeth respectively. Teeth have been removed in order not to consider respective shape’s variability in subsequent statistical analyses

• c. Automatic interpolation The software automatically detects the high-resolution segmentation contours, that are polylines of the mandibular mask, and fills them via an interpolation between CT slices. This procedure was used to remove potential cavities inside the mandibular mask and so obtain from this a manifold geometry (Fig. 2c)

• d. Post-processing operations Consisting in uniform re-mesh and smoothing, in order to improve the external geometry representation and remove potential artifacts deriving from the segmentation process (Fig. 2d). As a result, a mesh with triangular elements uniformly distributed over the geometry surface was generated, in order to avoid regions with denser/sparser nodes distribution [41]; elements’ dimension was established as a trade-off between limiting the average deviation from the original geometry and the number of surface nodes. Average deviations below 0.05 mm have been obtained with 2 mm triangle’s edge length.

• e. Soft tissue segmentation This operation was performed on CT scans in order to obtain the external geometry of the face of each subject (Fig. 2e), required for the identification of external morphometric measurements, as described in Paragraph 2.2; soft tissues have been identified setting the Hounsfield (HU) window to -700:225 HU.

Fig. 2

Segmentation process: a whole skull segmentation; b mandible with teeth geometry; c mandible without teeth and holes filled; d post-processing on mandible geometry; e soft tissue (face) geometry

2.2 Morphometric measurements

A literature research was performed to identify landmarks and morphometric measurements commonly used for mandibular geometry description [25, 42,43,44,45,46,47,48,49,50,51]. Among those available, only measurements that can be taken easily with classical instruments, as calliper and goniometer, during an external examination were considered in this study. As result of this preliminary analysis, ten anatomical landmarks (hard-tissue and soft-tissue landmarks, Table 2) were selected. These landmarks have been identified, on all mandible and face geometries, with 3-Matic software both automatically (‘Extrema Analysis’ tool) and manually based on definitions reported in Table 2 and information provided in literature [50, 51]; from these, twelve morphometric measurements (features) have been computed (Table 3).

Table 2 Hard tissue and soft tissue landmarks description

Table 3 Morphometric measurements description

With reference to the mandibular bone, literature reports many different morphometric measurements defined from hard-tissue landmarks, but far fewer soft-tissue morphometric features can be found. For this reason, some of the selected features reported in Table 3 have been adapted from the respective hard-tissue measurements.

Soft-tissue landmarks have been identified on the face geometry reconstructed from CT scans (Fig. 3), since it was not possible to directly record these measurements on patients. This aspect, along with proved dependency of facial landmarks on BMI (Body Mass Index) [52], might cause an overestimation of the mandible and bicondylar widths as well as of the body length when dealing with over-weight subjects. In these cases, the position of the gonion and lateral condylion soft-tissue landmarks does not correspond to the actual position which, instead, could be established during a direct examination. For these reasons, all those features that depend on these landmarks have been computed referring to the corresponding hard-tissue landmarks, resulting in more accurate feature estimation.

Fig. 3

Hard- and soft-tissue landmarks used for morphometric measurements definition

2.2.1 Morphometric Measurements Analysis

Morphometric measurements reported in Table 3 have been recorded for all fifty subjects in the dataset and mean and standard deviation values were computed for male and female groups. These values have been compared with data reported in literature. Differences in means between male and female mandibles were assessed for each variable using a two-sample t-Test with a significance level evaluated with the Bonferroni correction (p = α/k) as p = 0.05/12 = 0.004. Normality of data distributions was assessed with a Shapiro–Wilk test [53].

Morphometric features have also been analysed in terms of inter-observer reliability: measurements were evaluated on two randomly selected mandibles (one male and one female) by four different operators. All the observers have the same expertise in the use of 3-Matic software and were trained to record morphometric measurements on 3D models. The intraclass correlation coefficient (ICC) was used to assess reliability in measurements taken by different final users [54].

2.3 Mandibular mesh preparation

Before starting the implementation of the statistical shape model, mandibular meshes must be prepared in order to establish correspondence between all the shapes in the dataset; this operation is mandatory to enable the identification of points of correspondence among all the mandibular geometries, required for the subsequent analyses (GPA and PCA). When working with 3D shapes, one solution is the generation of iso-topological surface meshes, which are meshes characterised by the same number of nodes and with nodes representing the same geometrical point [25, 41, 55, 56]. Working with iso-topological meshes facilitate GPA and PCA analyses, reducing the effort of finding corresponding points among shapes, being these all the mesh nodes.

In the present study iso-topological surface meshes have been obtained implementing mesh morphing technique, which consists in non-rigid transformation of geometries; many different approaches exist to perform morphing, and here Radial Basis Function (RBF) method was used. This approach consists in adapting the mesh of a reference mandible (standard mesh), belonging to the dataset, on all the other mandible geometries (target meshes); specific control points are chosen over the standard and target surface meshes, and the whole set of the standard mesh nodes undergo displacements that are obtained by interpolation based on known displacements of these control points. Different type of ‘piecewise smooth’ (i.e., linear, thin plate spline) and ‘infinitely smooth’ (i.e., Gaussian, inverse multiquadric) basis functions exist [57, 58]; among these the RBF able to provide the best morph for the mandibular shape was assessed by some preliminary analyses, testing different formulations (both piecewise and infinitely smooth). The ability of the basis functions to globally morph the standard mesh was assessed qualitatively and the linear RBF was proved to generate the best morph accuracy also in bone regions with the most complex geometry (condylar and coronoid processes): indeed, using different basis functions generated a significant deformation of the morphed mesh and a consequent deviation from the target mesh in these regions, requiring so a higher number of morphing iterations to obtain a good accuracy. A detailed description of this methodology is provided in the previous work by Pascoletti et al. [59]; here the same procedure was followed and main steps are reported:

• a. The reference mandible is arbitrarily selected among all the mandibles in the dataset, as a geometry with no peculiar features (i.e., severe bone resorption or deformities): its mesh becomes the standard mesh

• b. A set of 27 control points have been identified both manually and automatically (3-Matic ‘Extrema Analysis’ tool) over the surface mandible geometry; these points have been selected as relevant anatomical landmarks for the mandible shape. In order to improve the morphing procedure, other 20 ‘auxiliary’ landmarks have been automatically identified in MATLAB using five segmentation planes defined by anatomical landmarks: by the intersection of these planes with the mandible geometry four new landmarks for each plane have been defined. Landmarks were identified for the standard mesh and for all the other target meshes belonging to the dataset, providing up to 47 control points for mesh morphing

• c. A linear RBF was chosen and used to interpolate nodes displacements, deforming the standard mesh based on the shape of the target geometry

• d. After morphing in strict sense, morphed mesh nodes were perpendicularly projected to the closest triangle of the target mesh and then a Laplacian smoothing was applied; these operations allow improving shape reproduction of target geometries

• e. These steps have been repeated iteratively until the average and maximum deviations between the morphed mesh and the original target geometry [60] were below 0.1 mm and 2 mm respectively.

At the end of these process, the dataset was composed by fifty mandibles with 3089 nodes and 6174 triangular elements.

2.4 Statistical shape model

The statistical shape model was generated on a training set composed by 40 mandibles belonging to the dataset (27 males and 13 females), while the remaining 10 mandibles (7 males and 3 females) have been used as control set for features selection and model accuracy evaluation.

A SSM consists in the deformation of an average shape $\overline{\user2{M}}$ through a variability model $\user2{\Phi b}$ to obtain a new predicted shape ${\varvec{M}}_{{\varvec{p}}}$:

$$ {\varvec{M}}_{{\varvec{p}}} = \overline{\user2{M}} + \user2{\Phi b} $$

(1)

The training set is composed by N_t = 40 shapes with N vertices and each shape can be represented by a matrix M $\in {\mathbb{R}}^{N \times 3}$.

First of all, the Generalised Procrustes Analysis was performed, filtering out location, scale and rotation effects through an iterative process: in this way only shape’s information are preserved in order to compute the average shape and perform subsequent statistical variability analysis.

Shapes are reported to a common reference frame by calculating the centre of gravity of each set of nodes and subtracting its coordinates from shapes’ nodes coordinates, resulting in shapes translated and centred in the origin of the new reference frame.

The iterative process starts with a first estimation of the average shape, arbitrarily chosen as the first mandible in the training set. Then mandibles are scaled so that each shape has the same norm as the average shape, and are rotated to be aligned to the average shape through a rotation matrix computed by the singular value decomposition (SVD) method. After these operations a new estimate of the average shape is made and the process starts again computing new scaling factors and rotation matrices. The iterative process ends when the deviation, that is sum of squared distances, between the new estimate of the average shape and the previous one is smaller than 0.001 mm². A more detailed description of all the mathematical steps of generalised Procrustes analysis can be found in the previous work [59]. As results of this analysis the average shape $\overline{\user2{M}} \in {\mathbb{R}}^{N \times 3}$ of the registered training set ${\varvec{M}}_{r} \in {\mathbb{R}}^{N \times 3}$ was obtained.

The variability model was then calculated through principal component analysis applied to the set of registered meshes M_r. Each shape M_r and the average mandible $\overline{\user2{M}}$ were rearranged stacking nodes coordinates of every i-th geometry into column vectors ${\varvec{m}}_{r,i} \in {\mathbb{R}}^{3N \times 1}$ and $\overline{\user2{m}} \in {\mathbb{R}}^{3N \times 1}$:

$$ {\varvec{m}}_{r,i} = \left[ {x_{1} ,y_{1,} z_{1} ,x_{2} ,y_{2,} z_{2} , \ldots ,x_{N} ,y_{N,} z_{N} } \right]^{T} $$

(2)

and new mandible vectors m_r,d were computed as deviation between shape vectors m_r,i and $\overline{\user2{m}}$:

$$ {\varvec{m}}_{rd,i} = {\varvec{m}}_{r,i} - \overline{\user2{m}} $$

(3)

These N_t vectors were then assembled in a matrix ${\varvec{M}}_{rd} \in {\mathbb{R}}^{{3N \times N_{t} }}$:

$$ {\varvec{M}}_{r,d} = \left[ {{\varvec{m}}_{rd,1} ,{\varvec{m}}_{rd,2} , \ldots ,{\varvec{m}}_{{rd,N_{t} }} } \right]^{T} $$

(4)

Principal component analysis was performed computing eigenvectors and eigen values of the covariance matrix of M_r,d. Eigenvectors $\user2{\varphi }_{{\varvec{i}}} \in {\mathbb{R}}^{3N \times 1}$ (i = 1, N_t − 1) are the so-called Principal Components (PCs) and eigenvalues λ_i associated to each of them represent the respective variance. Variances are sorted so that λ₁ > λ₂ > … > λ_Nt − 1 and corresponding PCs are rearranged in a matrix ${\varvec{\varPhi}}\in {\mathbb{R}}^{{3N \times (N_{t} - 1)}}$ accordingly.

These principal components define the main modes of deformation of the mandibular shape and a linear combination of these modes through weight vectors ${\varvec{b}}_{i} \in {\mathbb{R}}^{{(N_{t} - 1) \times 1}}$ allow deforming the average shape and generating new predicted shapes M_p.

2.5 Statistical shape model performance

One of the most used measures for the evaluation of the quality of an SSM is the compactness parameter. Compactness C defines the ability of the model to describe shape instances with as few modes as possible; therefore, a compact model is able to describe a shape from the population with a limited number of PCs. From its definition compactness is a parameter that depends on variances λ_i associated to principal components. Every PCs is able to describe a certain amount of the total variance contained in the dataset, that is the explained variance λ_expl,i:

$$ \lambda_{expl,i} = \frac{{\lambda_{i} }}{{\mathop \sum \nolimits_{i = 1}^{{N_{t} - 1}} \lambda_{i} }} $$

(5)

From this definition the compactness can be computed as:

$$ C\left( n \right) = \mathop \sum \limits_{i = 1}^{n} \lambda_{expl,i} $$

(6)

where n is the number of modes considered and C(n) is the compactness using n modes.

First principal components have been further investigated to identify the main modes of deformation of the mandibular shape; deformation of the average shape associated to these modes have been computed as:

$$ {\varvec{PC}}_{i} = \overline{\user2{M}} \pm 3\sqrt {\lambda_{i} } \cdot \user2{\varphi }_{{\varvec{i}}} $$

(7)

and respective deviations from $\overline{\user2{M}}$ were analysed considering the Euclidean distance between corresponding nodes.

2.6 Feature modification

The statistical shape model defined by Eq. (1) allows generating new shapes by the control of weights values, but it does not provide a direct way to predict shapes based on intuitive morphometric features, like those described in Paragraph 2.2. Indeed, weights almost never have a physical meaning and so SSM parameters have to be modified in order to provide a direct way to explore the range of mandibular geometries. This was here performed relating several variables by learning a linear mapping between selected control features and PCA weights [6, 9, 28, 29].

For each mandible within the training dataset a feature vector ${\varvec{f}}_{{\varvec{i}}} \in {\mathbb{R}}$^(L+1)×1 containing L features (L = 12 in this case) can be created:

$$ {\varvec{f}}_{i} = \left[ {f_{1} ,f_{2} , \ldots ,f_{l} ,1} \right]^{T} $$

(8)

A mapping matrix ${\varvec{K}} \in {\mathbb{R}}$^{(Nt−1)×(L+1)} relating the feature vector of a mandible to its weight vector may be defined as:

$$ {\varvec{Kf}}_{i} = {\varvec{b}}_{i} $$

(9)

Feature vectors can be assembled in a feature matrix ${\varvec{F}} \in {\mathbb{R}}$^(L+1)×Nt and similarly a weight matrix ${\varvec{B}} \in {\mathbb{R}}^{{(N_{t} - 1) \times N_{t} }}$ can be defined; extending Eq. (9) the mapping matrix is obtained solving for K:

$$ {\varvec{KF}} = {\varvec{B}} $$

(10)

$$ {\varvec{K}} = {\varvec{BF}}^{ + } $$

(11)

being F⁺ the pseudoinverse of F. From Eq. (1) weight vectors b_i can be easily computed through a multiple linear regression over the shapes of training dataset and so the mapping matrix can be obtained; using this mapping matrix, new weight vectors b_i,new for the generation of a new shape M_p, characterised by a feature vector f_i,new, can be generated as:

$$ {\varvec{b}}_{i,new} = {\varvec{Kf}}_{i,new} $$

(12)

2.7 Feature selection

Features modification allows defining physical features as direct input parameters of the SSM. A set of twelve measurements, identified as the more significant and easiest to measure external features, has been selected (Paragraph 2.2). Features that are most relevant for mandible shape description were investigated looking for the best combination of these features, were ‘best’ means the combination that minimises errors when new shape are predicted.

A cross-validation analysis was performed on a control set composed by N_c = 10 mandibles [6, 9, 28, 29]: being L = 12 the total number of features, twelve groups l_g have been defined, where every group includes all the possible k combinations of g features. For example, l₂ contains all the possible combinations of two features, for a total of 66 feature vectors; l₆ is composed by 924 feature vectors that are all the combinations of six features among the twelve available, and so on. Therefore, for every mandible belonging to the control set, predicted shapes have been computed considering all these combinations: being g the number of features considered, for each combination of g features belonging to the training set, matrices ${\varvec{F}}_{comb,g} \in {\mathbb{R}}^{{\left( {{\text{g}} + 1} \right) \times N_{t} }}$, were created and Eq. (11) was used to compute in turn mapping matrices ${\varvec{K}}_{comb,g} \in {\mathbb{R}}^{{\left( {N_{t} - 1} \right) \times \left( {{\text{g}} + 1} \right)}}$. Weight vectors can then be evaluated as:

$$ {\varvec{b}}_{comb i,g} = {\varvec{K}}_{comb,g} {\varvec{f}}_{comb i,g} $$

(13)

being f_{comb i,g} the feature vector of the i-th control mandible containing the combination of g features.

Every predicted mandible was compared to the corresponding mandible obtained from the CT scan data; more precisely, predicted shapes have been compared to mandible geometries after mesh preparation, being so the predicted and original meshes iso-topological.

Deviation between shapes was evaluated considering the Euclidean distance between corresponding nodes, being smaller distances associated to smaller deviations. Considering a control mandible i the deviation for the j-th node was:

$$ d_{ji} = \sqrt {\left( {x_{p,ji} - x_{ji} } \right)^{2} + \left( {y_{p,ji} - y_{ji} } \right)^{2} + \left( {z_{p,ji} - z_{ji} } \right)^{2} } $$

(14)

where (x_p,ji, y_p,ji, z_p,ji) and (x_ji, y_ji, z_ji) are respectively node’s coordinates of the predicted and real shapes.

For every feature combination k, an error distribution vector ${\varvec{\varepsilon}}_{k} \in {\mathbb{R}}^{N \times 1}$ has been evaluated over all the mandibles of the control set, being its elements:

$$ {\varvec{\varepsilon}}_{k,j} = \mathop \sum \limits_{i = 1}^{{N_{c} }} d_{ji} $$

(15)

The 90-th percentile value ε_k,90 was retrieved and for every group l_g the combination with the minimum value of this parameter was retained; moreover, the mean and standard deviation of corresponding vectors ε_k over all nodes have been computed:

$$ \mu_{{{\varvec{\varepsilon}}_{k} }} = \frac{1}{N}\mathop \sum \limits_{j = 1}^{N} {\varvec{\varepsilon}}_{k,j} $$

(16)

$$ \sigma_{{{\varvec{\varepsilon}}_{k} }} = \sqrt {\frac{1}{N - 1}\mathop \sum \limits_{j = 1}^{N} \left( {{\varvec{\varepsilon}}_{k,j} - \mu_{{{\varvec{\varepsilon}}_{k} }} } \right)^{2} } $$

(17)

and the minimum of the mean value was considered for identifying the best combination of features k_opt.

As result of the feature selection process, the SSM becomes:

$$ {\varvec{M}}_{{\varvec{p}}} = \overline{\user2{M}} + \user2{\Phi K}_{opt} {\varvec{f}}_{opt} $$

(18)

where f_opt and K_opt are the feature vector and the mapping matrix corresponding to the optimum combination identified.

The variability model in Eq. (18) can also be read as a linear combination of original PCs by weights represented by the element of the feature vector:

$${\varvec{\varPhi}}^{\user2{*}} = \user2{\Phi K}_{opt} $$

(19)

being Φ^*$\in {\mathbb{R}}^{{3N \times k_{opt} }}$^* the new eigenvectors matrix.

3 Results

3.1 Morphometric measurements analysis

Mean and standard deviation values have been computed for all the twelve features for male and female mandibles and results are shown in Fig. 4, where length measurements are in millimetres and angular measurements in degrees. A good agreement was found between these results and data reported in literature [25, 42, 44, 48], confirming also some general trends (i.e., female gonial angles values greater than male). Moreover, from this analysis it can be seen that bilateral measurements are consistent within a group.

Fig. 4

Mean and standard deviation values for male and female features. ** identify features with p < 0.001. Length measurements are in millimetres and angular measurements in degrees

Statistical differences between male and female features have been assessed and significant differences were found for mandible width (mean difference 7.1 mm), right and left ramus height (mean difference 6.3 mm and 6.5 mm respectively), right and left body length (mean difference 7.3 mm and 7.9 mm respectively), with p < 0.001 for all these features. Moreover, for the bicondylar width (mean difference 5 mm) a p = 0.0046 was found, which is very close to the imposed significance level (0.004). These results are in accordance with those reported in [25], where the same feature have been identified as significantly different among gender.

The ICC between different operators has shown a good reliability for the left body width and for the pogonion-infradentale distance (ICC equal to 0.808 and 0.878 respectively), while for all the other measurements an excellent inter-observer reliability was assessed, with ICCs > 0.933.

3.2 Statistical shape model performance

The compactness parameter is shown in Fig. 5 as function of the number of principal components. This analysis pointed out that the first principal component explained the 30% of the total variance and from C(4) on (blue dashed line in Fig. 5), which is equal to 57%, the gradient's changes become less significant. If the 98% of the total variance is wanted to be captured, 29 PCs are required (red dashed line in Fig. 5).

Fig. 5

Compactness curve as function of the number of PCs. Blue and red dashed lines represent C(4) and C(29) respectively (colour figure online)

Shapes representing the first four principal components have been calculated with Eq. (7), analysing deviations from the average shape of ± 3 $\sqrt {\lambda_{i} }$. These shape modes have been superimposed to the average shape geometry and deviations between deformed and average geometries have been computed and reported with a colormap on the average shape (Fig. 6); major deviations are represented by red and yellow areas. This analysis allowed identifying which are the main deformations explained by first four modes. It was pointed out that the first PC is associated to a uniform change of size of all the mandible distributed over the condylar process, the coronoid process and the internal part of the mental protuberance and alveolar process. Second, third and fourth modes generate more localised deformations: PC₂ mainly affects the coronoid process as well as the ramus and gonial regions, resulting in a uniform stretch/shortening in the superior-inferior direction, without changes in the opening angles; PC₃ produces effects on the central part of the mandibular bone (mental protuberance and the alveolar process) and on the gonial region, where deformations in the anterior–posterior and superior-inferior directions can be seen; a more evident widening/narrowing of the mandible shape is associated to PC₄ and deformations in the mylohyoid line and sublingual fossa can be observed.

Fig. 6

First Principal Components shapes for ± 3 $\sqrt {\lambda_{i} }$ and respective deviations from the average shape: a First PC; b Second PC; c Third PC; d Fourth PC. Colormaps are represented on the average shape and red/blue regions represent maximum/minimum deviations respectively (colour figure online)

3.3 Feature Selection

Having considered twelve features, a total of 4095 combinations and the corresponding predicted shapes have been generated for each control mandible following Eq. (18); the eigenvector matrix was limited to the first 29 PCs, in order to reduce the model dimensionality covering the 98% of the variance.

The ε_k,90 parameter has been investigated for all these combinations and Fig. 7 shows the prediction error distribution for the best features combination within a group l_g, that is the one with the minimum 90-th percentile value.

Fig. 7

Average error distribution for the best combination of features for every group l_g

Boxplot upper and lower limits represent the 10-th and 90-th percentile values respectively, while the grey line inside the box is the median; whiskers extrema indicate the minimum and maximum deviation. For every group l_g, the corresponding best combination was retrieved and the corresponding mean and standard deviations have been computed; Table 4 summarises these results showing that the minimum ε_k,90 value is associated to the l₄ combination (3.46 mm). Figure 8 shows colormaps representing the average error distribution ε_k, over the control set, for all the features group in correspondence of the best combinations reported in Table 4.

Table 4 Best features combinations for every group l_g and corresponding error vector parameters (90-th percentile and mean ± standard deviation)

Fig. 8

Error colormaps represented on the average shape surface for best features combinations of every group l_g

Analysing means values μ_εk it can be seen that the minimum is reached for the l₅ combination (2.79 ± 0.57 mm), while from this combination on this parameter increases again, indicating an overfitting trend. This analysis outlined that the most relevant measurements for shape prediction are represented by left body length, mandible and bicondylar widths, right gonial angle and right ramus height; moreover, from l₆ on also the mandibular corpus posterior angle becomes an important feature.

Considering the best features combination (l₅) the average prediction error for every mandible of the control set was retrieved. Results are reported in Fig. 9, where control mandibles M₂, M₃ and M₉ are female. The prediction error appears to be consistent among all the shapes in the control set, with the exception of mandible M₂, which exhibits the highest deviation (5.6 ± 1.97 mm).

Fig. 9

Average errors for all the mandibles of the control set resulting from the deviation between the real geometry and the shape predicted with the l₅

If the geometry of this control mandible is observed in detail, it can be seen that it is characterised by a peculiar morphology of the internal region of the mandibular body at the level of the mylohyoid line, resulting in a small distance between the left and right side. The real and predicted shapes for this mandible are shown in Fig. 10a: as can be seen the predicted shape is not able to reproduce the opening angles of the ramus and coronoid process, resulting in high errors (up to 10 mm), as shown by the deviation colormap in Fig. 10a. On the other side, if the mandible with the minimum deviation (M₃, 1.8 ± 0.81 mm) is considered, the colormap reported in Fig. 10b is obtained.

Fig. 10

Deviation analysis for the control mandibles with the a maximum (M₂) and b minimum (M₃) deviation. Real (grey) and predicted (red) geometries comparison on the left; deviation colormap represented on the real geometry on the right (colour figure online)

Two examples of predicted shapes obtained with the l₅ features combination are reported in Fig. 11a, comparing real (grey mandible) and predicted geometries (red mandible). The prediction of the SSM for one of the mandibles of the control set (M₅) is represented in Fig. 11a, while a mandible with a peculiar morphology (sever bone resorption/asportation on the left side of the mandibular body) was considered in Fig. 11b. This latter shape was out of both training a control set and was here used as benchmark for evaluating model’s prediction performances when only a partial mandibular geometry is available, proving the ability of the model to reproduce also missing or unhealthy bone’s regions.

Fig. 11

Comparison between real and predicted shapes using the l₅ features combination: a mandible belonging to the control dataset (M₅); b mandible with a peculiar morphology

4 Discussion

In this study a statistical shape model was developed, with the main aim of providing a useful and easy-to-use tool to predict mandible geometry based on external measurements readily collectable during classical examinations (i.e., dental visiting procedure).

Twelve external morphometric measurements have been selected to this end; the main aspects that were considered for the selection of these features are the ability to provide information on the mandible morphology, and the fact that they are easy to record by the final user of the tool (i.e., dentist, designer). In order to validate the latter point, the ICC parameter was used to quantify the inter-observer reliability and results have shown that all the selected features exhibit an ICC greater than 0.933, with the exception of two, which shown nevertheless a good reliability.

A mesh morphing technique based on RBF was here implemented to generate iso-topological meshes, guaranteeing so the nodes correspondence among shapes. The reference mesh adapted to all the other target meshes was arbitrarily chosen among those available in the dataset, with the requirement that it should not have any morphological peculiarities. The potential influence of the choice of the standard mesh on the morphing results was evaluated repeating the morphing procedure with a different target mesh: ten randomly selected mandibles were used for this analysis and deviations between the morphed and the target geometries were assessed. This analysis has pointed out that deviations associated to the first and the second standard mesh were comparable, with maximum deviations below 1.5 mm and average deviations not exceeding 0.15 mm, proving that the morphing procedure was not affected by the choice of the reference mandible. Mesh morphing may be influenced by the ability and sensibility of the operator at some extent, but it should be considered that this procedure was performed once to prepare the dataset and so its accuracy does not depend on the final user of the tool: in this study the morphing error obtained for all the mandibles of the dataset is comparable or lower than deviations obtained in other works [55, 60] where the same approach was implemented.

As described in the first part of the work, the 3D mandibular models were all prepared with a uniform mesh. The choice of remeshing the mandibular geometries with a uniformly distributed mesh was twofold: firstly, it allows limiting distortions during the morphing procedure of the template mesh, being itself uniform (Paragraph 2.3); second, no effects on the shape variance captured by a specific principal component (Paragraph 2.4) are introduced, considering that bone regions with a high density of nodes would have a greater weight than regions with a low density in PCs computation. The use of a uniform mesh avoids these issues and guarantees the meaningfulness of statistical analysis results [41].

Another aspect that has been investigated in relation to the mesh morphing, was the selection of anatomical landmarks used as control points. Repeatability of landmarks selection on different mandibles was assessed through an intra-observer reliability analysis: the ICC has been computed for five different mandibles, randomly selected from the dataset, for which the twenty-seven anatomical landmarks have been recorded three times by the same operator at four-day intervals. 3D coordinates of these points have been analysed and for all of them excellent reliability (ICCs > 0.9) was assessed.

The main modes of deformation of the mandibular bone have been assessed as result of the principal component analysis performed over a combined training dataset of forty mandibles, obtained from CT scan 3D reconstruction. This analysis has shown that the first four principal components are able to capture the 57% of the total variance of the dataset, while 29 PCs are required if the 98% of the explained variance is sought. First PCs account for the most of the variation of the mandibular shape, and the main modes of deformation associated to the first four principal components have been investigated varying from − 3 to + 3 standard deviations ($\sqrt {\lambda_{i} }$). A direct comparison between the average shape and the deformed shape based on each principal eigenvectors as well as corresponding colormaps have been generated (Fig. 6): a uniform change of size of all the mandible geometry is associated to the first PC, while more localised variations can be identified for PC₂, PC₃ and PC₄. A physiological asymmetry of the mandibular bone has been proved to characterise also the morphology of healthy bones, as a result, for example, of the masticatory muscles size and action [61,62,63,64], and this aspect seems to be captured by the third and the fourth modes, that exhibited a slight asymmetry in shape’s deformations between left and right side of the mandibular body.

Results obtained from the analysis of principal modes have been compared with results found in similar studies [25, 26], where a higher number of mandibles was available for the statistical analysis (65 and 60 mandibles respectively). Modes of variation identified in the present study, using a combined dataset, are in agreement with those highlighted by the male SSM analysed in [25], where both combined and gender based SSM were generated, revealing that the model here proposed is able to better capture and represent variations related to male mandibles. This result was expected, considering that a predominance of male shapes was used for the statistical analyses (27 male mandibles in the training set). Following these outcomes, as completion of this work a second SSM was built limiting the training and control sets to the male mandibles available in the datasets (27 and 7 geometries respectively). Analysing the first four principal components, the same main modes of deformation as in the combined model were found, with the exception of the fourth PC which appeared to be much more related to stretch and shortening of the coronoid process in the superior-inferior direction. Furthermore, the male SSM model is more compact than the original one: four PCs describe the 61% of the total variance and 20 PCs account for the 98% of the shape variability, in accordance with results obtained by [25].

From the feature selection analysis an optimal combination of five measurements has been identified (left body length, mandible and bicondylar widths, right gonial angle and right ramus height), with an associated average error of 2.79 ± 0.57 mm. A more detailed analysis of the average error distributions (Fig. 8) showed that the coronoid process is the most critical area to reproduce whenever is the number and the combination of features considered; this result stems from the fact that no external features directly related to the coronoid process have been considered, due to the fact that no facial measurements are able to provide useful information on the geometry of this region. If a more detailed prediction of this area is required, the use of X-Ray imaging techniques should be considered in order to obtain hard-tissue landmarks (i.e., the most superior point of the coronoid process, mandibular notch) and measurements able to describe the coronoid process area.

Beyond this, it can be observed that different features combinations generate better predictions in different mandible regions. Groups of features between l₁ and l₄ provide a worst estimation of the condylar and gonial angle regions, while from l₆ on, the coronoid process and the ramus regions are affected by the highest errors. This is the reason why the average errors associated to l₁₂ are comparable to those of the first two combinations: maximum errors on the coronoid process generates higher average errors, even if the deviation between the actual and the predicted shapes in the central area is low and comparable to those provided by other combinations. From this follows that if some specific mandibular regions are of most interest a different feature combination than l₅ might be more appropriate to obtain a more refined prediction. As general rule, the first six combinations provide a better estimation of the condylion and coronoid processes and ramus posterior region, while the last six allow better estimating the gonial angle region and the frontal part of the mandibular body.

If the same analysis is repeated for the male SSM, again the combination l₅ resulted in the minimum average error (2.45 ± 0.57 mm), with an improvement in the shape estimation of about 12%; moreover, this configuration is characterised by a minimum local deviation of 1.04 mm and a maximum local deviation of 4.72 mm, improving the corresponding parameters of the combined SSM of 33.8% and 2.3% respectively. The same general trends in shape prediction as for the combined SSM are confirmed, with the difference that combinations from l₄ to l₈ are able to provide a much better estimation of the most of the mandibular body, ramus regions and condylar processes, while combination with a higher number of features provide worst approximations of the coronoid process, with maximum average errors up to 6 mm (l₁₂). Overall, this allows reducing the estimation error between the original combined and the male SSMs when the same control mandible is considered (male control mandible); this improvement can be appreciated by consulting the figures provided in the Supplementary Material (Online Resource 1), where average error distributions for the different l_g combinations is reported.

A combined training set with both male and female mandibles has been used in this work in order to maximise the number of shapes available for statistical analyses. Despite its combined nature, the features-based SSM model was proved to be able to capture main modes of deformation of the mandibular shape, as emerged from the comparison with previous studies, and to generate new predicted shapes with average errors around 2.8 mm, comparable to those obtained from SSM based on higher number of mandibles (1.5 mm [25]); moreover, the model was able to provide the same level of prediction for both male and female mandibles within the control set (Fig. 9). Results shown in Fig. 9 have also pointed out that one of the mandibles of the control set (female mandible M₂) is characterised by the highest deviation between the real and predicted geometry (average error 5.5 mm) and this was due to the peculiar morphology of this mandible (Fig. 10a). If this mandible is excluded from the control set and the average error distribution analysis is repeated, the same optimal combination of five features is identified but lower average errors are found: with reference to the best combination l₅, the μ_εk parameter becomes 2.48 ± 0.52 mm and it is comparable to the male SSM corresponding value.

Even if the gender-based SSM (male SSM) certainly allows to further reduce the prediction error, the combined SSM here developed could represent a useful tool for mandibular shape prediction, also when a geometry reconstruction of partial bones is required (Fig. 11b).

One of the main critical points when working with statistical shape models is the availability of a high number of shapes; this aspect is much more critical when the objective is to analyse bones geometries, due to difficulties related to retrieving medical data (i.e., CT scans). The greater the number of shapes the more the results obtained from the statistical analysis can be generalised. Despite results here presented may be considered preliminary, due to the limited number of mandible available, this study has allowed setting up a full procedure for the prediction of new bone geometries and relating main shape’s variations to external morphometric features, highlighting potential and criticisms of this geometrical modelling approach.

The assessment of the SSM proposed in this work was based on geometric considerations in terms of accuracy of predicted shapes, that is the deviation between the predicted geometry and the mandible 3D model reconstructed from CT scans. As result of this analysis, it was possible to assess how the use of different combinations of input features affects the shape’s prediction in terms of overall and region-specific geometric accuracy. This validation approach was here considered in view of different possible applications of the 3D models predicted by feature-based SSM. Among these, those of greatest interest are represented by customised design of prostheses or medical devices based on the patient-specific bone geometry; finite elements analyses for evaluating the interaction between medical devices (i.e., mandibular plates, dental implants) and the bone; multibody analyses for the assessment of the mandible functionality (i.e., kinematic analyses for range of motions evaluation). Each of these applications deserve a specific validation analysis, beyond the geometry accuracy assessment; indeed, it is clear that parameters chosen for the model’s assessment depends also on the final purpose of the model itself. In this work the geometrical analysis was chosen in order to provide a general assessment of the proposed methodology. From the results here presented, it is possible to establish the most suitable combination of features to obtain an accurate prediction of the whole mandible or of specific bone regions and these information can be used to guide the creation of the predicted geometry based on its final use; these outcomes will be investigated with different validation analysis for more specific applications in future development of the model.

From the inter-observer reliability analysis performed on the morphometric measurements, it can be stated that the results obtained in this work can be considered as independent of the final user collecting the morphometric features. In order to completely generalise the SSM model for routine clinical use, future activities will involve recording the investigated features directly on real subjects.

5 Conclusions

A comprehensive tool for the prediction of subject-specific 3D models of the human mandible was generated, trying to take one step forward in statistical shape models of the mandibular bone, thanks to the correlation between the statistical variations of the mandible’s shape and morphological external measurements directly used as input to the model. The advantages of this statistical shape approach are twofold: first of all, it allows generating 3D subject-specific geometries without having to resort to expensive and not always available imaging techniques, but only requiring a limited set of morphometric features; in addition to this, the resolution of the mesh correspondence problem with the mesh morphing method makes it possible to generalize patient-specific procedures, thanks to the straightforward identification of the same anatomical points among different geometries.

A gender-combined SSM was developed and results here obtained were promising: predicted mandibles have been found to have low deviations from the original ones, with no relevant differences due to gender. Moreover, this tool could represent a support when pathological or partial mandibular bones are considered, being the model able to predict the whole mandible geometry by providing a statistically reliable estimation of missing parts. This modelling approach can represent an aid for pre-operative planning or for design of prosthetic devices for large bone portions. In this research the attention was focused on the definition of a methodology for the creation of subject-specific models, relating shape variations to morphometric features. In the future activities, this approach will be integrated into a wider interactive design framework, with the aim of replacing 3D model generation from imaging techniques. Moreover, the number of mandibles involved in the statistical analyses will be increased collecting new CT scans for both male and female bones and focusing on the creation of more refined gender-based models.