Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Ann Biomed Eng. Author manuscript; available in PMC 2017 May 1.
Published in final edited form as:
PMCID: PMC4833683

An eFace-Template Method for Efficiently Generating Patient-Specific Anatomically-Detailed Facial Soft Tissue FE Models for Craniomaxillofacial Surgery Simulation


Accurate surgical planning and prediction of craniomaxillofacial surgery outcome requires simulation of soft-tissue changes following osteotomy. This can only be accomplished on an anatomically-detailed facial soft tissue model. However, current anatomically-detailed facial soft tissue model generation is not appropriate for clinical applications due to the time intensive nature of manual segmentation and volumetric mesh generation. This paper presents a novel semi-automatic approach, named eFace-template method, for efficiently and accurately generating a patient-specific facial soft tissue model. Our novel approach is based on the volumetric deformation of an anatomically-detailed template to be fitted to the shape of each individual patient. The adaptation of the template is achieved by using a hybrid landmark-based morphing and dense surface fitting approach followed by a Thin-Plate Spline interpolation. This methodology was validated using 4 visible human datasets (regarded as gold standards) and 30 patient models. The results indicated that our approach can accurately preserve the internal anatomical correspondence (i.e., muscles) for finite element modeling (FEM). Additionally, our hybrid approach was able to achieve an optimal balance among the patient shape fitting accuracy, anatomical correspondence and mesh quality. Furthermore, the statistical analysis showed that our hybrid approach was superior to two previously published methods: mesh-matching and landmark-based transformation. Ultimately, our eFace-template method can be directly and effectively used clinically to simulate the facial soft tissue changes in the clinical application.

Keywords: template deformation, soft-tissue-change simulation, finite element modeling, surgical planning, CMF surgery, visible human, surface matching

1. Introduction

Craniomaxillofacial (CMF) deformities affect a human’s head and facial appearance. CMF surgery is a reconstructive procedure to correct such head deformities. The success of CMF surgery requires extensive pre-surgical planning. A good surgical plan involves accuate simulation of not only the bony segment movements (called osteotomy), but also facial soft tissue changes, following the osteotomy. However, only the osteotomy can currently be accurately simulated33. Accurately predicting facial soft tissue changes, following the osteotomy, still poses a major challenge to clinicians, due to the complex nature of the facial soft tissue anatomy.

Finite Element Modelling (FEM) is a widely accepted technique for virtually simulating facial soft tissue changes following virtual osteotomies under complex loading conditions. FEM allows incorporating highly complex and non-linear biomechanical tissues behavior9,11,22,31, and it is based on a volumetric discretization of the tissue structure through the definition of 3D meshes. The success of FEM depends on the geometrical fidelity of the FE model that represents the anatomy of interest and the accuracy of constitutive material properties used to characterize the tissue mechanical behavior20,32. There are two main challenges in the use of FEM. The first challenge is that the segmentation of facial soft tissue anatomy structure can be very labor-intense and time-consuming. The second challenge is the process of constructing a valid 3D mesh, with either tetrahedral or hexahedral finite elements. Compared to tetrahedrons, hexahedrons have advantages in convergence, error estimation and computation time3,38. However, it commonly takes hours, even days, to generate a good quality 3D hexahedral mesh for an object with complex geometry like a human face17. This time constraint significantly impacts the usability of such technology in a clinical setting. Besides FEM, mass spring models (MSM)27 and mass tensor models (MTM)7 are two alternatives used to simulate facial soft tissue changes following osteotomy22. MTM also requires anatomically-detailed structure information similar to FEM. Although MSM does not require detailed anatomic structures, its accuracy is limited22. Additionally, statistical and stochastic methods have also been used to predict soft tissue change16,25. The statistical model is trained by examining the pre- and post-operative data of a selective set of real patients. The stochastic approach uses probability of features to model the facial soft tissue changes based on input face features. The examination of facial soft tissue changes can then be achieved using purely imaging techniques1. However, such shape based statistical techniques are limited in their ability to model surgical procedures and changes to the underlying tissue. The accuracy of the prediction depends on the number of training patient samples. Most importantly, both, the statistical and stochastic methods, are inherently limited in extrapolating to patient cases beyond the datasets that are used to generate them.

The before mentioned limitations in generating patient-specific FE models may be overcome by building an anatomically-detailed template FE model and morphing it to an individual patient, thus enabling the CMF surgical simulation. Although this approach is time efficient, it also creates new challenges during the adaptation (morphing) of the template model. The first new challenge is shape accuracy. The morphing algorithm must accurately preserve all the anatomical correspondence, while fitting the shape of the patient. Another new challenge is the element distortion caused by superimposing the template shape to an individual patient with its own specifics. One solution is to deform a FEM template model based on skin and bone surface registration. Previously published methods include automatic elastic registration5 and landmark-based surface registration12,19. The elastic registration (mesh-matching algorithm)5 registers skin and skull surfaces automatically, depending on only the local topology of the mesh without feature or landmark constraints. The deformation of the generic model is then conducted by using the parameters of the surface registration. However, the correspondence of anatomical structure in the deformation has limited accuracy. Although this approach may be applicable for facial animations, it is not for clinical purpose of CMF surgical planning10. The landmark-based transformation method19 uses about 120 landmarks on the face skin surface (including anatomical landmarks, mathematical landmarks and pseudo landmarks) and 16 landmarks on the skull surface. It is difficult to digitize the mathematical landmarks and pseudo landmarks manually due to the lack of clear landmark definition. Based on the manually digitized landmarks, more sparse landmarks are generated to guide the surface registration. The registration accuracy depends largely on the number of manually digitized landmarks. The work by Huang et al.12 is incapable of preserving anatomical correspondence for large deformation using landmark-only “non-linear” affine transformation, thus this method can be efficient when the model deformation is limited. However, it cannot be used for our application where the deformation of the facial soft tissue model is quite large. Instead of directly deforming a template model, others reported on using image registration in retrieving the mapping field between the FEM template and the patient13,30. However, accurate image registration still faces challenges, e.g., feature extraction, similarity measures, computational cost, and interference of artifacts, which can be significant in head computed tomography (CT) images.

The objective of our study was to develop a novel strategy to significantly improve patient-specific anatomically-detailed facial soft tissue FE modeling. Furthermore, or goal was to validate our modeling approach on a fixed set of human subject data before starting a clinical trial. The clinical contribution of this project is that the resulting patient-specific soft tissue FE model can be directly and effectively used to simulate soft tissue changes following the osteotomy. Our key technical contribution is that our hybrid template morphing approach can achieve an optimal balance among the patient shape fitting accuracy, anatomical correspondence and mesh quality. The second contribution is that our approach can accurately preserve the internal anatomical correspondence (i.e., muscles) for future FE modeling of surgical procedures.

2. Materials and Methods

2.1 The template model

CT data and color gross anatomical cross-sectional images of a Chinese Visible Female36 were used to generate the template model. The CT data contained axial scans of 512×512 of scanning matrix, 300 mm field of view and 1.0 mm of thickness (voxel size: 0.586×0.586×1.0 mm3) slice. The scan parameters were set at 140 kV and 171 mAs. The resolution of the cross-sectional images for gross anatomy were 3072×2048, 0.17 mm of pixel size, and 0.25 mm of image intervals. Eleven muscles (buccinator, depressor anguli oris, depressor labii, levator anguli oris, levator labii, levator labii alaeque nasi, mentalis, orbicularis oris, zygomaticus major, zygomaticus minor and masseter)24 were segmented from the color images by 2 CMF surgeons (Z.T. and J.J.X.) together. Skin, mucosa and skull were segmented from CT data. All the segmented anatomical structures from CT and color cross-sectional images were registered, resulting in an anatomically-detailed soft tissue template model. This model was composed of following tissues: skin, mucosa (in close contact with skull), muscles and filler (as fat). The thickness of the skin and mucosa was fixed at 2 mm2. The template was transformed into a hexahedral mesh using Mimics 17 (Materialise NV, Leuven, Belgium). It contained 613,458 elements and 678,014 vertices. Within the hexahedral mesh of the template, the upper and lower lips were modeled separately. This was necessary to improve the simulation accuracy in the lip area, and account for open/closed lip variation in the image data sets. Finally, the triangulated skin surface (7586 vertices) and skull surface (15,444 vertices) were also constructed for surface registration between the template and patient dataset.

2.2 Novel eFace-template method

The framework of our eFace-template method is shown in Figure 1. The template deformation was divided into skin and skull (mucosa) surface registration and volume interpolation. For preserving the anatomical correspondence, a group of easily identifiable anatomical landmarks (45 on the skin, 26 on the maxilla and 22 on the mandible) were manually digitized to guide the surface registration. The landmarks on skin, maxilla, and mandible are shown in Figure 2. All but 4 of these landmarks (landmarks marked in red in Figure 2) were anatomical landmarks that were commonly used by clinicians8,23,26, and were associated with high repeatability and reproducibility in being identified18,28,29,35. The 4 landmarks around the lip area (marked in green in Figure 2) were mathematically derived landmarks - evenly dividing the right and left lips into 3 equal sections. The rationale for adding these 4 landmarks was to increase the accuracy of registration around the lip area. The 4 green landmarks, along with red landmark in the middle of the lip (called stomion), were used to separate the upper and lower lips. In case the patients were scanned in an open mouth position, the 5 landmarks were correspondingly digitized twice on both upper and lower lips, in pairs.

Figure 1
Framework of the proposed method. For each new patient, the skin and skull surfaces and “built-in” landmarks of the template are deformed to fit the shape of the patient through surface registration. Furthermore, the volumetric mesh of ...
Figure 2
Landmarks on skin (Top), maxilla (Middle) and mandible (Bottom). Landmarks in red are the anatomical landmarks, and landmarks in green on the lip are derived landmarks between the stomion (red landmark located in the middle of the lip) and cheilion (two ...

During the surface registration, landmark-based Thin-Plate Splines (TPS) were used, aligning the template skin and bone surfaces to the patient facial surfaces. TPS was chosen, because it produced smooth non-rigid deformation simultaneously based on multiple landmarks. However, the deformation of a vertex away from the landmarks must be interpolated based on its distance to the landmarks, which might increase distance errors. Therefore, we developed a dense surface fitting registration method to further refine the initial result from TPS. Based on the deformations of skin and bone surfaces, the volume deformation of the template was also interpolated by using TPS. This strategy works well even for large deformations in morphing from the template to the patient surfaces.

2.2.1 Landmark-based TPS registration

Given two sets of corresponding points X = {xi} and Y ={yi}, i = 1, 2, … , N , the TPS registration is to fit a mapping function by minimizing an energy function6

(Eq. 1)

For a 3D point (D = 3), xi = (xix, xiy, xiz, 1). By minimizing the first term in (Eq. 1), the point set X is mapped to Y. The second term is a smoothness constraint controlled by γ. The setting of γ directly affects the accuracy of mapping and smoothness. The function f minimizing the energy function (Eq. 1) has the form

(Eq. 2)

where d is an affine transformation matrix in dimension (D + 1) × (D + 1) and wi is a coefficient for the linear combination of a series of thin-plate radial basis functions [var phi], [var phi](a) = ||a||2 log ||a||. By substituting f in (Eq. 2) to the energy function (Eq. 1), we get

(Eq. 3)

where W is the coefficient matrix in dimension N × (D+1) and Ø is a (N × N) matrix formed from [var phi] (xixj). After determining d and W by minimizing (Eq. 3) based on the landmarks, the deformation for other vertices is interpolated using (Eq. 2).

2.2.2 Smoothness-constrained dense surface registration

The landmark-based TPS is very limited for generating a dense mapped mesh. In order to obtain an accurate mapping, the deformed template surface mesh is further projected to the target patient surface mesh. The projection results in a good geometrical match. However, at the same time, it might also yield a poor quality mesh because of intersecting triangles and distorted shapes. To overcome this problem, a smoothing term is used to relocate each vertex of the mesh to the centroid of its neighboring vertices. The smoothing definition is as follows:

(Eq. 4)

where iC is the index set of neighborhood vertices to the vertex vi in the vertex set V = {vi, i = 1, … , M}, and ni is the number of neighborhood vertices to the vertex vi. In the smoothing and projection, the correspondence of landmarks must to be preserved. Therefore, an optimization function is used to balance the smoothing, surface projection and landmark preservation. It is formulized as follows:

(Eq. 5)

The first term represents the dense surface projection. B is the set of indices of the rest of vertices in V except landmark vertices, B = {p1, …, pNB}. vp is the projected vertex of vp. The projected vertex on the target patient surface is found by calculating the closest point that has the minimum distance between a vertex and the target surface. The second term accounts for the landmark correspondence preservation. C is the set of indices of landmark vertices, C = {k1, … kMNB}. The third term is the smoothing constraint. α is a scale constant to control the degree of smoothness. The minima of (Eq. 5) can be determined by taking its partial derivatives, which results in a linear system

(Eq. 6)

where A=(FBFCαL), and b=(bBbC0). L is a M × M matrix, Lij={1i=j-1/nijiC0otherwise.

FB is a NB × M matrix. Each row of it contains only one non-zero element to constrain the position of the projected vertices. Elements of F are:

Eq. 7)

FC is a (MNB) × M matrix to constrain the position of landmarks, with elements as:

(Eq. 8)

V′ is a M × 1 vector corresponding to one component of coordinates of vertices. bB is a NB × 1 vector, and it represents the corresponding coordinates of the projected target vertices for the vertices in B. bc is a (MNB) × 1 vector, and it represents the corresponding coordinates of target landmark vertices. The linear system in (Eq. 6) can be solved in a sense of least square as V′ (ATA)−1ATb.

2.2.3 Volumetric mesh interpolation

Based on the deformation of skin and skull surface vertices, the deformation of hexahedral element vertices can be interpolated according to their distance to surface vertices. In our approach, the interpolation was achieved by using TPS. According to the correspondences between the original template skin and skull surface vertices (Hts and Htb), and the deformed skin and skull surface vertices ( Hts and Htb), the TPS deformation parameters d and W were calculated by minimizing the energy function (Eq. 3). The deformation of the volumetric mesh was then interpolated using (Eq. 2).

The summary of our template deformation algorithm is:

Input: Template hexahedral volume mesh Mt, template skin and skull surface meshes Hts, Htb, and patient skin and skull surface meshes Hps, Hpb
Step 1. Digitize corresponding landmarks manually on Hts and Hps, and Htb and Hpb
Step 2. Deform template surface meshes Hts and Htb based on landmarks using TPS, results in [H with macron]ts and [H with macron]tb
Step 3. Project the vertices on [H with macron]ts and [H with macron]tb to the patient skin and skull surfaces, respectively
Step 4. Optimize smoothing and projection using (Eq. 5), results in H′ts and H′tb
Step 5. Repeat Step 3 to Step 4 until the update of H′ts and H′tb is within a prior tolerance (10−3 in average).
Step 6. Interpolate the deformation of hexahedral volume mesh by using TPS, results in M′t.
Output: Deformed hexahedral volume mesh M′t

In the volumetric mesh interpolation, bad surface triangle vertices were eliminated to avoid creating invalid elements. The bad surface triangles were defined as these with a triangle angle smaller than 10° or intersected with other triangles. They were mainly distributed on the boundary of the surface, the connection areas of upper and lower lips, the eyes and the connection areas between teeth. We verified that eliminating these vertices did not affect the accuracy of the surface fitting.

In normal subjects and some patients (depending on the types of the deformities), the reposed upper and lower lips are naturally closed and the anterior teeth are covered by the lips. However, in several other patients, the anterior teeth are exposed when the lips are at naturally reposed. During the volumetric interpolation using the skin and skull surface deformations for a patient with exposed anterior teeth, the inner layer of the upper lip is deformed and stretched to be visible. This will cause interpolation errors in the lip area. To avoid this problem, the anterior teeth exposed outside of the lips are marked and excluded during the volumetric interpolation.

2.3 Evaluate the effectiveness of our eFace-template method

In this experiment, data of 30 patients with Class I (4 patients), Class II (5 patients), and Class III (21 patients) dentofacial deformity were used. The protocol was approved by our Institutional Review Board prior to the study (IRB0413-0045). For each patient, an axial CT scan of head (512×512 of scanning matrix, 1.25 mm of thickness slice, and 250 mm filed of view, captured in 120 kV and 250 mAs) and a 3D facial soft tissue photo (3dMD, Atlanta, GA) were used. The 3dMD camera was operated by a clinician to ensure each patient’s facial expression was neutral and lips were reposed. 3D CT soft tissue was replaced by the 3dMD soft tissue by registering the 3dMD photo to the CT model. This was done to prevent any possible facial soft tissue strain during the CT scanning. The criterion of the registration of 3dMD facial surface to the CT model was to minimize the error between them. This was done by using surface based registration in Mimics software (Materialise NV, Leuven, Belgium)21,34. The effectiveness of our eFace-template method was evaluated in three aspects: surface registration, volumetric interpolation and the comparison with 2 previously published methods, mesh-matching (M-M) method5 and landmark-based transformation (LbT) method19. The segmentation of the patient’s CT data, 3dMD image registration and landmark digitization were completed by 2 CMF surgeons (Z.T. and C.-M. C.) together.

2.3.1 Surface registration evaluation

The accuracy of surface registration (registration error) was measured using the deviation map between the deformed template surface and the original patient surface. The registration error was measured as the shortest distance between each given vertex on the deformed template surface to the original patient surface. We then calculated the root mean square (RMS) value of all distances and used the resulting value for evaluation.

2.3.2 Volumetric interpolation

In volumetric interpolation, if a higher density of surface mesh is used to interpolate the volumetric mesh, the surface of the template can be more accurately aligned to that of the patient’s model. Ho ever, the quality of hexahedrons may be reduced as there are more stringent constraints on the surface points, and the computational time will be increased accordingly. In order to determine an optimal density of surface mesh for the volumetric interpolation, the changes of mesh quality, registration error of interpolated surface to target surface, and computational time based on the mesh densities were first analyzed. The mesh densities were controlled by selecting vertex samples in different space resolutions.

The hexahedron mesh quality was measured by the hexahedron shape skew metric15 and Scaled Jacobian37. The shape skew metric was in range [0, 1], in which 1 was for a rectangular brick and 0 for a degenerate element. Scaled Jacobian was evaluated for each node cell of a hexahedron, and it was in range [−1, 1]. The higher the metric the least distorted the hexahedron. Negative values of the Scaled Jacobian signified that the elements were invalid (invalid elements). The average of the 8 Scaled Jacobian values of a hexahedral element was used to describe the quality of this element.

2.3.3 Comparison ours with M-M and LbT methods

To evaluate the advantage of our eFace-template method, the efficacy of our method was compared to M-M and LbT methods. Our manually digitized landmarks were used as the input landmarks required for the LbT method. The same γ was used to control the smoothness in LbT method. Three transformations were conducted using LbT. The first transformation was based on the manually digitized landmarks, the second transformation was using automatically chosen landmarks by dividing surface into cells 10 mm × 10 mm × 10 mm, and the landmarks for the third transformation were chosen by dividing surface into cells 3 mm × 3 mm × 3 mm. In our volumetric interpolation, we also used 3 mm resolution surface vertex samples to interpolate the deformation of hexahedral mesh. The M-M method used the rigid, affine, and local spline registration parameters of skin surface to deform the volume mesh. Only nodes, those were physically in contact with the skull, were deformed based on the registration parameters of the skull surface.

The efficacy of each method was evaluated using the following paramters: absolute surface registration error (RE, in mm), landmark correspondence errors (CE, in mm), mesh quality, and mean computational time. The mesh quality was measured by shape skew (SS)15, scaled Jacobian (SJ)37, and number of invalid elements (InE). The mean values of SS and SJ and InE from all the testing patients were used for final comparison.

Since our patient samples included three different types of dentofacial deformities: Class I, II, and III, we determined through statistical analysis if there was a statistically significant difference in the efficacy of our method among the Class I, II and III patients. We used F-tests to determine if the variance of the dependent outcome measures (RE, CE, SS, SJ and InE) were comparable between the patient groups. Furthermore, we used T-test to determine if there were any statistical differences in means. The calculation was only done for our eFace-template method.

Based on the result of above statistical analysis, a two-factor repeated measures analysis of variance (ANOVA) and post-hoc (Tukey’s) tests were used to determine whether there was a statistically significant difference among the different modeling methods. Between-subject factor was the modeling methods: M-M, LbT and our eFace-template methods, and within-subject factor was the outcome measure, which included all 5 parameters (RE, CE, SS, SJ and InE).

2.4 Methods for validation of internal corresponding structures

Validation of our eFace-template method was carried out through analyzing the accuracy of internal corresponding structures (i.e. muscles). Four anatomically-detailed soft tissue models, including skin, mucosa and individual muscles, were manually segmented and generated using the following datasets: NLM Visible Male (NLM-M), NLM Visible Female (NLM-F), Chinese Visible Human 3 (CVH3) and Korea Visible Male (KVM). The segmentation was completed by two CMF surgeons (Z.T. and J.J.X.) together. The segmented muscles served as ground truth during the validation. The template model was then deformed and adapted to the facial soft tissue model of each visible human model using our hybrid method. Two overlap ratios (measured by recall and precision) were calculated to determine the overlap between the mapped template muscles and the manually segmented ground truth muscles. For the volume of mapped muscles A and the volume of ground truth muscles B, recall=ABB, and precision=ABA. Higher recall and precision values imply a better overlap between the resulted and ground truth muscles.

Our eFace-template method was compared to M-M and LbT methods. Finally, a two-factor repeated measures ANOVA was used to detect whether there was a statistically significant difference among the three modeling methods. Within-subject factor was the 2 measures: recall and precision. Between-subject factor was the modeling methods: M-M, LbT and our eFace-template methods. If there was a statistically significant difference, post-hoc (Tukey’s) tests were conducted to determine if any method was statistically better. The correlation of overlap ratios from the three methods also was studied.

3 Results

Based on experiments, smoothness constraint γ = 0.1 worked best for our application (achieving mapping error within 0.02 mm in (Eq. 1)), and the constant scale α = 1 produced a good balance between mesh quality and mapping accuracy. The algorithms were implemented in MATLAB using a PC with a regular office personal computer (3.4 GHz i7 central processing unit (CPU) and 16 GB random-access memory (RAM)). In the dense surface registration, it usually took 10–15 iterations for the mapping procedure to converge.

3.1 Effectiveness of our eFace-template method

All patient specific models were successfully generated using our eFace-template method despite the types of dentofacial deformities. Figure 3 shows four representative examples. The effectiveness of our method was further quantitatively assessed in the following three aspects: surface registration, volumetric interpolation, and the comparison with the M-M and LbT methods.

Figure 3
Generated soft tissue models for 4 representative deformities ((a) Class III; (b) Class III with anterior open bite; (c) Class II; and (d) Class II with strained lip and anterior open-bite) using our method. Each patient shows an original 3dMD skin surface, ...

3.1.1 Surface registration evaluation

The surface registration results of a patient are shown in Figure 4 (e) (f). The deformed skin and skull surfaces were visually identical to the original surfaces of the patient in Figure 4 (c) (d). The color-coded deviation map and RMS of the registration errors for the skin and skull surfaces are shown in Figure 4 (g) (h). They indicate that the registration errors in the majority of the regions were smaller than 0.15 mm. Only small regions on the skin surface around the lips, and on the skull surface around teeth, temporal and zygomatic process had a registration error of 0.3 mm or larger. One reason of a slightly larger registration error on these regions might be that the original patient mesh density was higher than the template mesh (Figure 5 (a–c)). Another reason might be that the proposed surface projection method in the surface registration might reduce the mesh density during the projection from a smoothing surface to a more curved surface (Figure 5 (d)), resulting in smoothing out some of the details on the template model. Finally, registration errors larger than 1 mm were only found in small regions around the eyes on the skin surface and the boundary of the skull surface (Figure 4 (g) (h)). This was due to the different eye status (closed versus open) between the template and the patient. Nonetheless, it was clinically insignificant as the eyes were not the region of interest for facial soft-tissue-change simulation following the osteotomies.

Figure 4
Surface registration result after the template is deformed to a given patient. The registration error with value larger than 1 mm is marked in red. (a) Skin surface of the template; (b) Skull surface of the template; (c) Original patient’s skin ...
Figure 5
The skull surface in (a) is from the patient in Figure 4. (b) and (c) are close-up views of the corresponding region in the red box in (a). (b) is original mesh of the patient and (c) is mesh resulted from the template deformation. (d) is the surface ...

3.1.2 Volumetric interpolation

Figure 6 shows the changes of the mesh quality, registration error of interpolated surface to target surface, and computational time based on mesh densities (vertex samples in resolutions from 1 mm to 5 mm). The plotted curves indicated that the registration error of interpolated surface would be increased while the number of surface points used for interpolation was reduced. However, the mesh quality of hexahedrons was increased (shape skew metric and Scaled Jacobian were increased slightly and invalid elements were reduced distinctively) and the computation time was reduced largerly while reducing the number of vertex samples from 1 mm to 5 mm resolutions. To balance among surface error, mesh quality and computation time, we concluded that 3 mm might be the optimal resolution of the mesh density. Therefore, the 3 mm resolution was also used in our experiments.

Figure 6
Influence of interpolation vertex samples. 1 mm resolution has 12390 vertices, 2 mm has 9869 vertices, 3 mm has 7269 vertices, 4 mm has 5127 vertices, 5 mm has 3691 vertices. The calculated Shape Skew, Scaled Jacobian, registration error and calculation ...

The invalid elements were mainly distributed around the eyes, the jaw, and sharp corners close to the skull surface (Figure 7 (a)). This was because that the eyes were closed in the template, while the eyes of the patient 3dMD photographs were mostly open. Therefore the elements around the eyes on the template were forced to deform largely. This might cause large distortion of elements, which in turn might render them as invalid. The density of the jaw meshes was often reduced after surface registration due to the sharp curve of the jaw (see Figure 5 (d)). Therefore, elements around this area were also easily distorted. Similarly, sharp corners close to the skull surfaces were a source of distorted and invalid elements. Among elements in the original template corresponding to these invalid elements, about 96% of them had a minimum Scaled Jacobian less than 0.4 (Figure 7 (b)). The distortion of most of these elements could be avoided in a future template with conformal elements. In addition, the invalid elements could be repaired by a post processing step using the method by Chabanas et al.5.

Figure 7
(a) Invalid elements (in green) of a generated patient FE model. (b) The distribution of the minimum Scaled Jacobian of the elements in the original template corresponding to these invalid elements.

3.1.3 Comparison ours with M-M and LbT methods

The performance (RE, CE, SS, SJ, and InE) of our eFace-template method was calculated as mean ± standard deviation and tabulated (Table 1) based on the different types of the dentofacial deformities. The result of the F-test for variance showed that there was no statistically significant difference among the Class I and II patients in all 5 measures (P > 0.05). The subsequent T-test for means also showed no statistical difference among Class I and II patients in all 5 measures (P > 0.05). Therefore, the patients with Class I and II types of the deformities were pooled and analyzed as a single group. The statistical comparison of the patients with Class III type of deformities with the remaining patients only showed significantly difference in means of RE and CE. Although there was a statistical significance between Class III and remaining patients in means of RE and CE, the absolute difference is only around 0.02mm and 0.03mm, which was not clinical significance. Therefore, in the following comparison to the M-M and LbT methods, the three groups of patients were pooled together.

Table 1
Performance of our method in three patient groups (mean ± standard deviation).

The results of the repeated measures ANOVA showed there was a statistically significant difference among the 3 methods [F (2, 210) = 681.44, P [double less-than sign] 0.01]. The differences calculated in a post-hoc Tukey’s test showed that our eFace-template method had a significantly smaller RE and significantly larger SS and SJ than LbT method. In addition, our method had significantly smaller CE and InE than M-M method. This statistical analysis showed that our method was able to combine the benefits of the M-M and the LbT methods.

The results of comparing our eFace-template method with the M-M and LbT methods showed that using our method, both RE (0.2 mm) and CE (0.3 mm) were small and clinically acceptable (within 1 mm is the clinical consensuses) (Table 2). It indicated that our method could accurately morph the template model to individual patient’s shape while preserving the anatomical correspondence. M-M method was able to fit the shape of patient with a smaller surface RE (0.1 mm). However, the anatomical correspondence was clinically unacceptable (CE was 12 mm). It used a rigid registration plus an affine transformation to initially align the template to the patient. For large deformation, the alignment was not accurate, and it may only match part of the patient shape (Figure 8 (b)). LbT method had a relatively large RE (1.6 mm) while the CE was the smallest (< 0.1 mm). Unlike the smooth skin generated using our method (Figure 8 (d)), the skin surface generated by LbT was extremely coarse and also clinically unacceptable.

Figure 8
Surface registration results. (a) Original, (b) M-M, (c) LbT, (d) Ours. Problematic regions are marked in red circles. In the result by M-M method, the chin is missed. This is because that the initial alignment by a rigid registration plus an affine transformation ...
Table 2
Comparison with two related methods (30 subjects, *mean ± standard deviation).

The result of the mesh quality evaluation showed that using our eFace-template method, InE was the smallest while SS and SJ were comparable to the other two methods (see Table 2 for details). As shown in Figure 7, most of these InE could be avoided by using a template with conformal elements. Finally, total CPU time of our method (for both registration and interpolation) was similar with the other two methods. The main cost of our CPU time was the volumetric mesh interpolation. The volumetric mesh interpolation cost about 3.5 minutes while the surface registration cost 0.5 minutes. The high CPU cost by volumetric mesh interpolation was due to the computation of thin-plate radial basis functions [var phi] in (Eq. 2) for large number of hexahedron elements (613,458 elements). For further time improvement, we may need to find an optimized distance calculation method for calculating [var phi]. We also believe that by converting our MATLAB code to executable Microsoft Foundation Class C/C++ code, there will be a significant speed improvement - expect to be within a minute based on our programming experience.

3.2 Validation of internal corresponding structures

Our anatomically detailed template was successfully registered to the gold standard models of the four visible human subjects using all three methods (M-M, LbT and our eFace-template methods). Recall and precision values for the four visible human subjects were calculated to determine the overlap between the mapped template muscles and the manually segmented ground truth muscles. The overlap values of muscles by the three methods are shown in Figure 9. Our eFace-template method achieved a recall value of 0.71 in average with standard deviation 0.06, and a precision value of 0.77 in average with standard deviation 0.08. This demonstrated that the generated patient-specific muscles were more than 70% correct. The comparison among the 3 methods showed that the muscles mapped using our eFace-template method had the largest overlaps (Figure 9).

Figure 9
Bar plot of recall and precision values for comparing accuracy of muscle matching by the three methods: M-M, LbT, and our eFace-template methods. Our eFace-template method had the largest recall and precision values on all the four visible human subjects ...

The result of the repeated measures ANOVA showed that our eFace-template method was statistically significant different (F (2, 18) = 13.97, P = 2×10−4) from the other two methods in volumetric overlap. The results of the post-hoc (Tukey’s) tests showed that the overlap of our method was significantly better than M-M method (P < 0.01 for both recall and precision). In the comparison with LbT method, our method had a significantly larger precision value than LbT (P = 0.01), but the difference of recall was not that distinct (P = 0.08). The statistical test generally concluded that our method better preserved the correspondence of anatomical structure than the other two methods for the facial soft tissue model generation. The correlation analysis showed that there was no correlation (P > 0.05 of the correlation measure) between the three methods in the muscle overlap measure.

4. Discussion

This paper presents an efficient hybrid approach, the eFace-template, to generate accurate patient-specific anatomically-detailed facial soft tissue FE models. It is based on the volumetric deformation of a template model to the individual shape of each patient’s facial soft tissue. The deformation of the template is achieved by using a hybrid landmark-based morphing and dense surface fitting approach, followed by a TPS interpolation.

Our eFace-template method has the combined benefits from M-M and LbT methods, the two most common approaches for template model mapping. The surface registration of M-M method only depends on the local topology of the mesh without feature- or landmark-based constraints. This results in a less accurate anatomical correspondence. LbT method achieves surface registration based on first manually digitized and later automatically generated landmarks. The registration accuracy depends on the number of manually digitized landmarks. Therefore, the resulting model after the registration using only 93 landmarks is coarse with larger surface deviation errors. In contrast, our method combines landmark-based surface morphing with dense surface projection to accurately fit the template model to the shape of patients while preserving anatomical correspondence. In addition, a smoothing term is used to avoid the reduction of the mesh quality in the dense surface fitting, and this helps reduce the possibility to generate invalid elements after interpolating the volumetric deformation based on the surface deformations. Our method could potentially generate the patient-specific facial soft tissue model with accurate internal corresponding structures (overlap of muscles larger than 70%), accurate surface matching (surface fitting error of 0.2 mm, corresponding error of important anatomic points of 0.3 mm), and acceptable hexahedral mesh quality (only less than 0.3% invalid elements). It is clinically significant as the resulted patient-specific soft tissue FE model can be directly and effectively used to simulate soft tissue changes following the osteotomy.

We also observe some areas that could be improved in the future. They include template model, landmark digitization and errors in lip area. The template model contains a number of distorted hexahedral elements, which contribute as the main source to the InE in our experiment. We are currently testing an approach of using voxel elements for the template. With this approach, the number of distorted hexahedral elements tends to be reduced to less than 0.1%.

Landmark constraint is necessary and effective in the preservation of anatomical correspondence. However, we recognized the landmark digitization error. Some landmarks, e.g. the landmarks on the cheek, may have large variations due to the smooth nature of the cheek – the landmarks lack of clear anatomical definition. These landmarks are important to define the facial shape. In addition, although we select landmarks which are easy to be identified anatomically, it still takes around 20 minutes to complete the landmark digitization. An automatic process for digitizing landmarks may significantly improve the accuracy and efficiency. We currently are working on a separate project to achieve automatic landmark digitization.

Finally, in a real clinical case, some areas of the facial soft tissues, e.g., lips and cheeks, are not attached to bone or teeth. Currently there is no accurate model to describe such sliding movement of the facial soft tissue. In this paper, the soft tissue with lip was assumed to be attached on the bone or teeth. For most of the region, this assumption does not create a larger error as the relative positions of tissue and bone is similar between patients and template. But the difference of relative positions of lip and teeth between patients is much more distinctive and complex than other mouth regions. This may cause problems in the lip area while interpolating the deformation of the soft tissue based on the deformation of skin and bone surfaces. We are currently improving the template model by including the lip sliding deformation.

In the future, we will continue to improve the mesh quality during the template deformation. One potential framework is to first generate the patient-specific surface soft tissue model with skin, skull and muscles by deforming a template model, and then generate high quality hexahedrons mesh models by developing an effective hexahedral mesh generation algorithm. Another future work is to use our generated patient-specific FE model to conduct parameter optimization to characterize the tissue mechanical behavior, and then achieve tissue deformation simulation and compare with clinical outcome.

5. Conclusions

This paper presents a novel patient-specific facial soft tissue mesh generation algorithm, the eFace-template method. It is based on the volumetric deformation of a template model with predefined anatomic structure to fit the shape of the patient. The deformation of the template is achieved by using landmark-constrained surface deformation followed by a TPS interpolation. Experimental results show the effectiveness in the generation of facial soft tissue model while preserving the anatomic correspondence. The presented method has significant clinical meaningfulness as it tremendously saves the time in the generation of subject-specific facial soft tissue mesh for facial soft-tissue-change simulation following the virtual osteotomies.


This work is funded by NIH/NIDCR grants R01DE022676 and R01DE021863. Dr. Tang was sponsored by China Scholarship Council and Dr. Chang was sponsored by Taichung Veterans General Hospital while they were working at the Surgical Planning Laboratory, Department of Oral and Maxillofacial Surgery, Houston Methodist Research Institute, Houston, TX, USA.

The authors also acknowledge that Chinese Visual Human data were provided by Institute of Basic Medical Science, Third Military Medical University, Chongqing, China.


Analysis of Variance
Correspondence error
Central processing unit
Computed tomography
Chinese Visible Human 3
Finite Element
Finite Element Modeling
Invalid elements
Korea Visible Male
Landmark-based transformation
Mass spring model
Mass tensor model
National Library of Medicine - Visible Female
National Library of Medicine - Visible Male
Registration error
Scaled Jacobian
Root Mean Square
Shape skew
Thin-Plate Splines


1. Baik HS, Kim SY. Facial soft-tissue changes in skeletal Class III orthognathic surgery patients analyzed with 3-dimensional laser scanning. Am J Orthod Dentofacial Orthop. 2010;138(2):167–178. [PubMed]
2. Barbarino GG, Jabareen M, Trzewik J, Nkengne A, Stamatas G, Mazza E. Development and validation of a three-dimensional finite element model of the face. J Biomec Eng. 2009;131(4):041006. [PubMed]
3. Benzley SE, Perry E, Merkley K, Clark B. A comparison of all hexagonal and all tetrahedral finite element meshes for elastic and elasto-plastic analysis. 4th International Meshing Roundtable; Sandia National Laboratories; 1995. pp. 179–191.
4. Bucki M, Lobos C, Payan Y, Hitschfeld N. Jacobian-based repair method for finite element meshes after registration. Engineering with Computers. 2011;27(3):285–297.
5. Chabanas M, Luboz V, Payan Y. Patient specific finite element model of the face soft tissues for computer-assisted maxillofacial surgery. Med Image Anal. 2003;7:131–151. [PubMed]
6. Chui H, Rangarajan A. CVPR. Vol. 2. Hilton Head Island; South Carolina: 2000. A new algorithm for non-rigid point matching; pp. 44–51.
7. Cotin S, Delingette H, Ayache N. A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation. Visual Comput. 2000;16(8):437–452.
8. Farkas LG. Anthropometry of the head and face in clinical practice. 2. New York: Raven Press; 1994.
9. Freutel M, Schmidt H, Dürselen L, Ignatius A, Galbusera F. Finite element modeling of soft tissues: material models, tissue interaction and challenges. Clinical Biomechanics. 2014;29:363–372. [PubMed]
10. Gladilin E, Zachow S, Deuflhard P, Hege HC. Anatomy- and physics-based facial animation for craniofacial surgery simulations. Med Biol Eng Comput. 2004;42(2):167–170. [PubMed]
11. Hambli R, Allaoui S. A robust 3d finite element simulation of human proximal femur progressive fracture under stance load with experimental validation. Annals of Biomedical Engineering. 2013;41(12):2515–2527. [PubMed]
12. Hung APL, Wu T, Hunter P, Mithraratne K. A framework for generating anatomically detailed subject-specific human facial models for biomechanical simulations. Visual Comput. 2015;31(5):527–539.
13. Ji S, Ford JC, Greenwald RM, Beckwith JG, Paulsen KD, Flashman LA, McAllister TW. Automated subject-specific, hexahedral mesh generation via image registration. Finite Elem Anal Des. 2011;47(10):1178–1185. [PMC free article] [PubMed]
14. Joe B. Technical report. Zhou Computing Services Inc; 2008. Shape measures for quadrilaterals, pyramids, wedges, and hexahedra.
15. Knupp PM. Algebraic mesh quality metrics for unstructured initial meshes. Finite Elem Anal Des. 2003;39:217–241.
16. Lee HJ, Suha HY, Leea YS, Leeb SJ, Donatellic RE, Dolced C, Wheelere TT. A better statistical method of predicting postsurgery soft tissue response in Class II patients. Angle Orthodontist. 2014;84(2):322–328. [PubMed]
17. Li Y, Liu Y, Xu W, Wang W, Guo B. All-hex meshing using singularity-restricted field. ACM Trans on Gr (SIGGRAPH) 2012;31(6):177.
18. Lin H, Zhu P, Lin Y, Zheng Y, Xu Y. Reliability and reproducibility of landmarks on three-dimensional soft-tissue cephalometrics using different placement methods. Plast Reconstr Surg. 2014;134(1):120e–110e. [PubMed]
19. Lou HD, Chen S, Chen G, Xu TM, Rong QG. Patient-specific modeling of facial soft tissue based on radial basis functions transformations of a standard three-dimensional finite element model. Chin Med J. 2012;125(22):4066–4071. [PubMed]
20. Luboz V, Promayon E, Payan Y. Linear elastic properties of the facial soft tissues using an aspiration device: towards patient specific characterization. Annals of Biomedical Engineering. 2014;42(11):2369–2378. [PubMed]
21. Maal TJ, Plooij JM, Rangel FA, Mollemans W, Schutyser FAC, Bergé SJ. The accuracy of matching three-dimensional photographs with skin surfaces derived from cone-beam computed tomography. Int J Oral Maxillofac Surg. 2008;37:641–646. [PubMed]
22. Mollemans W, Schutyser F, Nadjmi N, Maes F, Suetens P. Predicting soft tissue deformations for a maxillofacial surgery planning system: from computational strategies to a complete clinical validation. Med Image Anal. 2007;11:282–301. [PubMed]
23. Netter FH. Anatomy of the head and neck. Summit, New Jersey: Ciba-Geigy Corporation; 1993.
24. Pan B, Xia JJ, Yuan P, Gateno J, Lp HH, Lee PK, Chow B, Zhou X. MICCAI. Vol. 7510. France: 2012. Incremental kernel ridge regression for the prediction of soft tissue deformations; pp. 99–106. [PMC free article] [PubMed]
25. Sub HY, Lee SJ, Lee YS, Donatelli RE, Wheeler TT, Kim SH, Eo SH, Seo BM. A more accurate method of predicting soft tissue changes after mandibular setback surgery. J Oral Maxillofac Surg. 2012;70(10):553–562. [PubMed]
26. Swennen GRJ, Schutyser F, Hausamen JE. Three-dimensional cephalometry. Heidelberg: Springer Berlin Heidelberg; 2006.
27. Teschner M. PhD Thesis. Shaker Verlag; Aachen, Germany: 2001. Direct computation of soft-tissue deformation in craniofacial surgery simulation.
28. Titiz I, Laubinger M, Keller T, Hertrich K, Hirschfelder U. Repeatability and reproducibility of landmarks-a three-dimensional computed tomography study. European Journal of Orthodontics. 2012;34:276–286. [PubMed]
29. Toma AM, Zhurov A, Playle R, Ong E, Richmond S. Reproducibility of facial soft tissue landmarks on 3D laser-scanned facial images. Orthod Craniofac Res. 2009;12:33–42. [PubMed]
30. Tong C. PhD Thesis. National University of Singapore; Singapore: 2013. Generation of patient-specific finite-element mesh from 3d medical images.
31. Wu TS-H. PhD thesis. The University of Auckland; 2014. A computational framework for modeling the biomechanics of human facial expressions.
32. Wang Q, Sun W. Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans. Annals of Biomedical Engineering. 2012;41(1):142–153. [PubMed]
33. Xia JJ, Gateno J, Teichgraeber JF. A new paradigm for complex midface reconstruction - a reversed approach. J Oral Maxillofac Surg. 2009;67(3):693–703. [PMC free article] [PubMed]
34. Xin P, Yu H, Cheng H, Shen S, Shen SGF. Image fusion in craniofacial virtual reality modeling based on CT and 3dMD photogrammetry. The Journal of Craniofacial Surgery. 2013;24(5):1573–1576. [PubMed]
35. Zamora N, Llamas JM, Cibrián R, Gandia JL, Paredes V. A study on the reproducibility of cephalometric landmarks when undertaking a three-dimensional (3D) cephalometric analysis. Med Oral Patol Oral Cir Bucal. 2012;17(4):e678–88. [PMC free article] [PubMed]
36. Zhang SX, Heng PA, Liu ZJ, et al. The Chinese visible human (CVH) datasets incorporate technical and imaging advances on earlier digital humans. J Anat. 2004;204(3):165–173. [PubMed]
37. Zhang Y, Bajaj C. Adaptive and quality quadrilateral/hexahedral meshing from volumetric data. Comput Methods Appl Mech Eng. 2006;195(9):942–960. [PMC free article] [PubMed]
38. Zienkiewicz OC, Taylor RL. The finite element method: basic formulation and linear problems. Maidenhead: MacGraw-Hill College; 1987.