Study participants

Forty patients with facial asymmetries of mandibular deviations were selected from the Department of Orthodontics, Department of Maxillofacial Surgery, and Department of Prosthetics of the Peking University Stomatology Hospital. The inclusion criteria were as follows: deviation of the chin from the midline of the face perpendicular to the line of the bilateral pupils via the point of the nasal root greater than 3 mm when the patients were in the natural head position [18, 19]. This study was approved by the Ethics Committee of Peking University School and Hospital of Stomatology (PKUSSIRB-202,054,042). The volunteers were fully informed of the content and purpose of the experiment and provided informed consent.

Experimental equipment and software

A 3D facial scanner FaceScan (3D-Shape Corp, Germany) with a scanning speed of 0.2–0.8 s, scanning accuracy of 0.1 mm, range of 270–320° (covering the left ear to right ear), raster scanning principle, and charge-coupled device with 5 million pixels, and approximately 10,000 data points was used for the experiments. The mean distance between facial point clouds was within the range of 0.6–2 mm. The 3D facial data processing software Geomagic Studio 2013 (3D System, Morrisville, NC, USA) was used, the EWPA algorithm program was developed in the Python programming language, and the objective function of the classical PA algorithm was optimized for weighting.\(\text\) is obtained as shown in Eq. (1):

$$} = \mathop \limits_Q \sum\nolimits_^p }_}}}\_}}_}} - } \cdot }\_}}_}}} \right\|}_2}},$$

(1)

where LMK_Org is the original landmark, LMK_Mir is the mirror landmark, LMK_Orgi and LMK_Miri(i = 1,2,…,32) are the corresponding landmarks in the original and mirror point landmarks, respectively, Q is the spatial change matrix, p is the number of landmarks, and Wi(i = 1,2,…,32) is the weight factor of each facial landmark.

Experimental methodsAcquisition and processing of 3D face data

The equipment was calibrated before acquisition to ensure that accurate data were obtained. The patient was seated 135 cm from the facial scanner and guided by the clinician to a natural head position, with the eyes looking forward, keeping the Frankfort horizontal plane (FH plane) parallel to the ground plane, and with a naturally relaxed facial expression. The criteria that could be adopted for 3D facial data were effective in displaying the facial contours, high resolution, and absence of significant movement and occlusion. Using the reverse engineering software Geomagic Studio 2013, the 3D data of the patient’s face were processed as necessary, including removing redundant data for hole repair. The original 3D facial model was adjusted to the natural head position such that the FH plane of the natural head coordinate system coincided with the XZ plane of the global coordinate system, and the sagittal plane coincided with the YZ plane of the global coordinate system. The extraction of important anatomical landmarks of the original facial model (Model_Org) was conducted by a senior expert based on clinical experience. In total, 32 anatomical landmarks (10 in the midline and 11 bilateral) were manually selected by the expert in the full facial region, including the trichion, glabella, superciliary ridge, nasion, pronasale, subnasale, exocanthion, pupil, endocanthion, tragion, zygion, alare, subalare, labiale superius, labiale inferius, sublabiale, pogonion, gnathion, crista philtre, cheilion, and gonion, which were extracted three times at one-week intervals to obtain the mean value of the coordinates as LMK_Org, as shown in Fig. 1; Table 1. The coordinates of the center of gravity of LMK_Org were calculated, and the original model was translated until the center of gravity of LMK_Org coincided with the origin of the global coordinate system and was saved as an .OBJ file.

Fig. 1

The 32 anatomic landmarks (upper facial third: trichion, glabella, superciliary ridge; middle facial third: nasion, pronasale, subnasale, exocanthion, pupil, endocanthion, tragion, zygion, alare, and subalare; lower facial third: labiale superius, labiale inferius, sublabiale, pogonion, gnathion, crista philtre, cheilion, and gonion)

Table 1 Definitions of three-dimensional facial anatomical landmarksQuantitative assessment of the asymmetry of facial landmarks

In this study, the asymmetry of 3D facial anatomical landmarks was quantitatively assessed using the EDMA method, which compares the geometric differences between two individuals based on a matrix of distances between a series of landmarks [20, 21]. For example, if K anatomical landmarks on an organism exist (1, 2, 3, 4 ……K), the distance lines d(1,1), d(1,2), … d(1,k), d(2,1), d(2,2), … d(2,k), … d(k,1), d(k,2), …d(k, k) between any two landmarks can be calculated for a total of k(k-1)/2 line segments. From this a form matrix (FM) of order k × k is formed, which can reflect the common characteristics of a particular data type and the individual differences between the same types. When analyzing the differences between two bodies of the same type, the ratio of the distance of the corresponding landmarks from the line segments in the FM of the two bodies can be formed into a new matrix called the form difference matrix (FDM) of the two bodies.

The specific needs of the quantitative evaluation of anatomical landmark asymmetry of the human face and the set of landmarks on the left and right sides of the face (11 landmarks each) were combined with the set of landmarks on the midline (10 landmarks) to build the “left face morphology” landmark set (L) and the “right face morphology” landmark set (R). “The “left face” and “right face” were treated as mirror images, and their morphological differences were quantified, and a new index based on the EDMA method was established to assess the degree of asymmetry of the landmarks, defined as morphological index E value. Considering landmark set R as an example, the right face morphological matrix was formed by calculating the distance between the midline and right landmark and each right landmark (a total of 10 × 1 + 11 × (11 − 1)/2 = 165 distances), named matrix \(_\). Similarly, we calculated the left face morphology matrix\(_\) (a simplified diagram of the line segments in the matrix is shown in Fig. 2). The FDM was calculated as the ratio of the linear distances between the same name landmarks in the matrix \(_\) with \(_\)(the larger distance is used as the numerator; all ratios are greater than or equal to 1) forms of the left and right faces, with 165 ratio elements, as shown in Eq. (2):

$$_ :=\frac_,_\right)}_,_\right)},$$

(2)

where FDM is the morphological difference matrix of the left and right faces, and MRk and MLk are the kth elements of the matrices MR and ML (k = 1,2,…,165).

Fig. 2

Line segment between landmarks (e.g., exocanthion, zygion, cheilion, pronasale). The yellow line segment is the distance between the midline and bilateral landmarks in the right side of the patient’s face landmark set, the blue line segment is the distance between the midline and bilateral landmarks in the left side of the patient’s face landmark set, and red points inside the blue and yellow dotted lines are the L landmark set and R landmark set respectively)

If the ratio in the matrix is 1, the left and right sides are perfectly symmetrical; if the ratio is greater than 1, a morphological difference (asymmetry) between the two sides is observed. To investigate the symmetry of a single landmark, this study proposes averaging a landmark and all its associated FDM line ratio elements as a morphological indicator of that landmark, defined as the Fi value (i = 1, 2, …,21), where Fi denotes the mean of the line segment ratio elements associated with the ith landmark in the morphological difference matrix. To make the distribution of the definition domain of Fi values as independent variables more consistent with the demand of the value domain of the assignment function, this study established a mapping relationship between the morphological index E values and Fi values, as shown in Eq. (3):

$$_=100\left(_-_\right),$$

(3)

where Ei is the morphological index of the landmark, Fi (i = 1,2,…,21) is the mean of the FDM line segment ratio elements of the landmark, and Fmin is the minimum value of Fi.

The smaller the Ei, the closer it is to 0, and the lower the degree of asymmetry of the landmark, and vice versa, the higher the degree of asymmetry. Two expressions of the assignment function were constructed in this study, namely, the offset power function and the linear function, as shown in Eqs. (4) and (5), respectively:

$$W = - \frac}}}E + 1.$$

(5)

Graphs of the two assignment functions for the newly constructed landmark asymmetry index E are shown in Fig. 3.

In this study, Wi (i = 1,2,…,21) is the weighting factor of each facial landmark, E is the morphological index of the landmark, and Emax is the maximum value of E.

Fig. 3

Function graphs of OP and LP. (a) Offset power function. (b) segmented power function. E is the morphological index of the landmark. Wi is the weight factor for each facial landmark

Facial MSP constructionExperimental group 1: construction of the MSP based on the PA algorithm

For the 40 facial models in this study, the 3D spatial coordinates of the 32 landmarks in LMK_Org were input into the PA algorithm program based on the Python language based on the previous study of our research group [22]. LMK_Mir was obtained by mirroring the original landmark set based on the YZ plane and calculating the overlap effect between the original and mirror landmark set based on the PA algorithm in Python, without weight differences. The transformation matrix of the mirror landmark set was then calculated and loaded onto the LMK_Mir using Geomagic Studio 2013. Finally, the SRP of the facial data for each patient was constructed by taking the union of the original and mirror models (Model Uni_PA) in Geomagic Studio using the function “plane” and “symmetry”, defined as ‘MSP_PA’.

Experimental group 2: construction of the MSP based on the EWPA algorithm

The 3D spatial coordinates of the 32 landmarks in LMK_Org were input into the EWPA algorithm program based on the Python language, and the E value of each landmark was automatically calculated and weighted using the method described in Sect. 2.3.2. Based on the E value offset power function and linear function, Wi had a monotonically decreasing distribution with increasing E value; the highest weight was one, and the lowest was 0. The above algorithms for the weighting and alignment of landmarks based on the offset power function and linear function are defined as EOWPA and ELWPA, respectively.

The original landmark set is based on the YZ plane mirror landmark set, the weighted overlap effect of the paired landmarks between LMK_Org and LMK_Mir is achieved based on the weighted least squares method, and the final weighted optimal overlap position of LMK_Org and LMK_Mir is obtained. The same method in 2.4.1 was used to construct the MSP, defined as the EWPA MSP (MSP_EOWPA, MSP_ELWPA).

Reference group: construction of the reference plane based on the regional ICP algorithm

Based on the above Model_Org and Model_Mir models, a region of good symmetry of the face on the original model was manually selected by a senior expert on the Geomagic Studio software, a mirror model Model_Mir was obtained for Model_Org based on the mirrored YZ plane, and ICP registration was performed based on the region selected by the expert. After overlapping the regions of the original and mirror models, the original-mirror model was obtained, and the MSP was calculated, which was defined as the “reference plane” in this study (MSP_Ref). The final reference plane was constructed for the facial data of the 40 patients included in this study.

The effect of the same 3D face data based on the PA, EOWPA, ELWPA, and regional ICP algorithms on the construction of the MSP is shown in Fig. 4.

Fig. 4

Determining the MSP based on the EOWPA, PA, and regional ICP algorithms in one case. The red plane signifies the MSP of ground truth, the green plane represents the EOWPA algorithm, and the yellow plane represents the PA algorithm

Data analysis

To investigate the intra-observer error of landmark extraction, the senior clinical expert repeated the landmark measurements three times at one-week intervals, and the intra-class coefficient (ICC) was calculated. For the 40 facial models in this study, the angle between the MSP (MSP_PA, MSP_EOWPA, and MSP_ELWPA) and the expert reference plane (MSP_Ref), denoted as Ang_PA, Ang_EOWPA, and Ang_ELWP, respectively, constructed by the PA algorithm and the EWPA algorithm were calculated for each model. The mean and standard deviation of the angle error of each algorithm were calculated.

The Shapiro–Wilk normality test was performed using SPSS software (version 21.0) on the plane angle errors of experimental group 1 (PA algorithm) and experimental group 2 (EWPA algorithm) for the 40 patients. One-way ANOVA with Tukey’s multiple comparison test was used if the samples conformed to a normal distribution and the variances were homogeneous; otherwise, the Kruskal–Wallis H-test was used, and the test levels were all α = 0.05 on both sides as a significant difference.

The reference plane was used to classify the degree of facial asymmetric deformity in the 40 patients, based on the calculation of the distance from the pogonion to the reference plane and to analyze the angle error between the reference planes and MSP_PA, MSP_EOWPA, and MSP_ELWPA for different degrees of deformity in the 3D facial data.