A surgical simulation system for predicting facial soft tissue deformation

In the field of cranio-maxillofacial (CMF) surgery, surgical simulation is becoming a very powerful tool to plan surgery and simulate surgical results before actually performing a CMF surgical procedure. Reliable prediction of facial soft tissue changes is in particular essential for better preparation and to shorten the time taken for the operation. This paper presents a surgical simulation system to predict facial soft tissue changes caused by the movement of bone segments during CMF surgery. Two experiments were designed to test the feasibility of this simulation system. The test results demonstrate the feasibility of fast and good prediction of post-operative facial appearance, with texture. Our surgical simulation system is applicable to computer-assisted CMF surgery.


Introduction
With the development of techniques in computerassisted orthognathic surgery, research into craniofacial morphology and models of bone and soft tissue have played an important role in maxillo-craniofacial therapy [1]. 3D reconstruction of bone and soft tissues, quantitative measurements, and virtual surgical simulation results provide significant underpinnings to optimize surgical programs. Prediction of new facial appearance is particularly important. Without any risk to patients, surgeons can operate on a patient virtually until an aesthetically satisfactory result is obtained. Therefore, building a precise prediction system for facial soft tissue deformation is very important and useful.
In this paper, we describe a simulation system for facial soft tissue prediction which is based on the linear finite element method. The simulation system can be used by surgeons to obtain an estimate of facial appearance. Our work utilizes individual patient data for both skull and skin surfaces to build the deformation model. The data needed can be acquired from 3D photography and 3D computed tomography (CT). To represent facial soft tissue deformation given the movement of bone structures, we build a correspondence between a face mesh and a bone mesh. Our facial models are composed of tetrahedral elements to permit finite element analysis. Facial tissue changes caused by bone movement can be well simulated.
The remainder of this paper is organized as follows. In Section 2, we present related research and underlying technologies related to our proposed method. We briefly review the principle of the linear elasticity finite element model (FEM) algorithm in Section 3. We then describe in detail in Section 4 our method for building the deformation model and the functions of each module in the deformation prediction pipeline. In Section 5, we report test results obtained with our system, showing that our approach is feasible. Finally, in Section 6, we conclude our work and discuss possible future improvements.

Related work
In the early 1980s, the prediction of post-operative facial appearance was initiated. Some research [2] presented the idea of predicting soft tissue changes due to bone movement by combining bone structure models together with soft tissue models. Mass spring models (MSM) were originally used to predict the soft tissue behavior. They originated in computer animation and could achieve real-time effects, but are considered to be rather inaccurate when precise modeling of physiological behavior is required. Finite element models (FEM) provide another popular method to simulate facial soft tissue deformation. This is an accurate approach to simulating the mechanical and physiological characteristics of soft tissue, but it is more computationally expensive [3]. Keeve et al. [4] compared these two most popular soft tissue models for maxillofacial surgery. They concluded that FEM could be superior to MSM for biomechanically accuracy, but inferior in terms of speed and sample architecture.
Koch et al. [5] proposed an approach to simulate soft tissue changes arising because of realignment of underlying bone structure. However, their methods were not integrated into a surgical simulation system. Maciel et al. [6] designed a system for surgical simulation. However, their attempts focused on the interactions between the user and the computer. Due to the high efficiency requirements, most commercial software products for predicting facial soft tissue deformations for maxillofacial surgery planning are based on statistical strategies. However, the prediction accuracy depends on the number of patient samples involved. Among these strategies, physically based methods can simulate soft tissue deformation more accurately and realistically.

Basic principles of linear FEM
In elasticity theory, the finite element method (FEM) is the most commonly used method to describe the continuous mechanical behavior of soft tissues [7]. This method is based on a volumetric discretization of the structure, using a 3D mesh. For 3D modeling, the most commonly used elements are hexahedra and tetrahedra. In our paper, we use tetrahedra as the mesh elements.
After the FEM discretization, the motion of a deformable tissue can be described by the Euler-Lagrange equation, which is a second order system of ordinary differential equations: where M ∈ R 3n×3n denotes the mass matrix, U ∈ R 3n denotes the displacement of the mesh points, U denotes the first derivative of the displacement vector, andÜ denotes its second derivativer. f ext denotes the external force (e.g., a user imposed force, or the force produced by the movement of bone). f in denotes the internal force or volume force. The mass matrix M can be pre-computed from the object's mesh density distribution in the rest state. The damping forces D(U ,U ) are independent of time, and can be represented as where α and β are Rayleigh damping constants, which are the positive real-valued parameters, which are variable in this system. We use K to denote the tangent stiffness matrix, which can be represented by the Jacobian matrix of the internal forces f in . For each mesh element, f denotes the internal force acting on it. This satisfies the mechanical equilibrium equations: where σ σ σ is the stress tensor. We use u to describe the displacement of the element. From u we can calculate the strain tensor ε ε ε, which describes the deformation per unit area. There are two popular choices in the field of computer graphics, the Green-Lagrange nonlinear strain tensor ε ε ε G : ε ε ε G = 1 2 (∇u + ∇u T + ∇u × ∇u T ), and the Cauchy linear strain tensor ε ε ε C : ε ε ε C = 1 2 (∇u + ∇u T ), where ∇u denotes the gradient of the displacement vector u. In a large deformation setup, the Green-Lagrange strain tensor ε ε ε G can provide a better result than the Cauchy strain tensor ε ε ε C , while when the deformation displacements are relatively small, the Cauchy strain tensor ε ε ε C suffices to calculate the result accurately. Using the Cauchy linear strain tensor can reduce the computation required. In this paper, we use the Cauchy strain tensor ε ε ε C . Hooke's Law states that the stress-strain relationship is linear, σ σ σ = Dε ε ε, where: Here λ and µ are Lamé constants directly proportional to the well-known material properties of Young's modulus E and Poisson's ratio v: Figure 1 shows our deformation prediction pipeline. It includes the following four parts.
• Data preparation. The pre-operative CT scan taken to construct the model of bone structure and the pre-operative 3D photography taken to obtain a textured mesh model of facial surface are input. A mesh generator is used to construct a tetrahedral model for FEM. • Mapping. This part includes two kinds of meshmapping methods, a surface-mapping method and a bone-mapping method. The deformation result is directly expressed via the skin surface by the surface-mapping method. The bone-mapping method allows simulation of interaction between bone and soft tissue.
• Simulation. The obtained forces are applied to the skin surface, according to the force-transfer principle, to simulate the deformation caused by bone movement. Using suitable boundary conditions for the biomechanical soft tissue model, the new facial appearance is computed using a linear elasticity FEM approach. • Control and display. The overall control and display structure are based on the model-viewcontroller pattern.

Data acquisition
Prior to craniofacial surgery, a pre-operative 3D CT scan and 3D photographs are acquired. The commercial product F aceSCAN 3D is used to obtain a 3D facial mesh, which describes the structure of the surface of the face as well as giving texture information. 3D CT scans are widely used as a tool to observe diseased human tissue, and provide excellent visualization of craniofacial tissue. The marching cube algorithm [8] is employed to reconstruct dense, highly detailed bone surface mesh models. The high quality skin surface acquired from the 3D photography (see Fig. 2(a)) is registered together with the bone structure, ready for simulation.

Dividing the intermediate mesh
To simulate skin tissue deformation, we need to build a volume mesh data model containing all deforming skin tissues. As volume elements, we use tetrahedra; these permit a good expression of skin tissue properties. To save processing time and memory usage, we restrict the tetrahedral mesh to the surface region that will deform during surgery, near shifting bones. Our system allows the user to confirm the deformation region near such bones via an intermediate mesh. As shown in Fig. 2

Generating tetrahedra
Our "ST Generator" algorithm is shown in Algorithm 1; it generates a high quality tetrahedral mesh from the intermediate mesh. The algorithm establishes tetrahedral elements from the surface mesh based on the negative face normal and point normal directions. The algorithm creates new points in certain offset positions along the negative normal directions of the grid points. In the same way, new points can be generated at certain offset positions from the center of each triangle surface along the negative normal direction of the triangle. Then, the algorithm creates tetrahedral elements using the following principles: • Each triangle determines a tetrahedral element using the offset point of its center, as shown in Fig. 3(a). • Given two adjacent triangles, their common edges  determine a tetrahedral element using the offset points of their centers, as shown in Fig. 3(b).
• Given two adjacent triangles, the two offset points of their centers separately determine a tetrahedral element using one of their shared points and its offset point, as shown in Fig. 3(c).

Surface-mapping method
For computational efficiency, the tetrahedral mesh used for simulation is restricted to the surface region near shifting bones. Thus, the tetrahedral mesh containing all facial tissues is unsuitable for clinical validation. Moreover, since the perception of the human face is mainly determined by the relationships between different parts of the face, it is important for surgeons to see the whole face when inspecting the effect of a maxillofacial procedure, not just a part of it. To combine the advantages of a smooth tetrahedral mesh for fast simulation and a realistic, highly detailed mesh with texture information for a good recognition, a surface-mapping method is used.

Bone-mapping method
As bone moves, facial deformation is inevitable in craniofacial surgery. Surgeons want to predict appearance changes caused by the movement of bone. Initially, we tried to build the relationship between the face mesh and bone mesh. However, we found that building a bone-related model directly is too time consuming if all the points are calculated in a time step. Our system restricts the tetrahedral mesh to the surface region near the shifting bones. For example, when the jawbone moves, tissue deformation near the upper skull bone can be ignored. We designed the relationship between the bone mesh and the intermediate mesh to take advantage of this approach. The relationship between a bone and soft tissue is determined by point normals on the bone. As the ray along the point normal hits the soft tissue, the hit point is calculated, and the relationship is built. For efficiency, an object oriented binary space partitioning (BSP) tree is used to find the hit point quickly. Obviously, most points of the bone mesh will be connected to only a limited number of face meshes. A distance threshold value must be set to prevent hyperteloric points from being linked together. The effect is shown in Fig. 4(c), where the threshold value is set to 2 cm.

Force-transfer principle
The displacement vector of the bone is regarded as the source of motion. We implement force-transfer method on the following principle: (1 i N ). We use the force vector F to represent the stress vector produced by the movement of the bone, which is a constant value for all points of the bone mesh. F i denotes the force vector at the corresponding ith point P i in the face mesh structure; θ i denotes the angle between the normal at P i and the direction of F.

Boundary conditions
In the bone-mapping method, three types of points can be distinguished in the intermediate mesh: forced points, fixed points, and free points: During simulation their movement is completely determined by the elastic force at these points. The boundary points restrict the deformation region and guarantee the stability of the model.

Simulation principle
Before analyzing the finite element model based on the tetrahedral mesh, material property values which represent attributes of particular tissue types have to be assigned; these are Young's modulus and Poisson's ratio. Measuring precise human tissue material parameters is a difficult task, and many researchers are still working on this topic [9]. Here, we assumed the tissue to be isotropic, homogeneous, and linearly elastic. Based on research by Shlivko et al. [10], we set the value of Young's modulus to 1 MPa and the value of Poisson's ratio to 0.46 during all simulations.
The generated tetrahedral mesh and boundary conditions serve as inputs for our soft tissue simulator. When moving the lower jaw, the forced points are assigned according to the bone-mapping method and the force vectors at these points are calculated based on the force-transfer principle. The simulator calculates the displacements of these points in the tetrahedral mesh at each time step. At the same time, the steady-state situation of the free points near the forced points is broken. At each time step, the coordinate states of all points in the intermediate mesh are transferred back to the skin surface mesh and shown in real time.

Control and display
Based on the model-view-controller pattern, our software framework is composed of several modules: the data model module, the view module, the simulation controller module, and the tetrahedron generation module. The data model module manages the pipeline of data and interaction between interfaces. The simulation controller applies physical rules to elements to modify their positions, velocities, accelerations, and other physical attributes. The tetrahedron generation module establishes the volume mesh composed of tetrahedral elements. The view module of our framework uses OpenGL to render meshes with appropriate structure, point coordinates, textures, and lighting information.

Results and discussion
We have tested the feasibility of the tetrahedral mesh. We also measured the accuracy of our simulation in terms of the differences between the effects of deformation, and data after operations.

Validation of ST Generator algorithm
To verify the operation of the ST Generator algorithm, we extracted the skin from CT image sequences using the marching cubes algorithm [11] and generated a tetrahedral mesh using TetGen [12].
To compare the quality of the two kinds of tetrahedral mesh, we used the aspect ratio, face angle, and dihedral angle of the data as evaluation criteria. As shown in Fig. 5, the tetrahedral meshes have similar peak values for aspect ratio, face angle, and dihedral angle. It indicates that they have the similar structural attributes and the mesh produced by the surface-based tetrahedral generator is suitable for building the tetrahedral volumetric FEM model.

Prediction of deformation
Proplanner is a commercial software product for operation design. It implements prediction of soft tissue deformation based on statistical strategies. In our experiment, we used our prediction system and this commercial system in turn to predict the shape of surface tissue when the mandible is moved to a fixed position. The difference between the results of these two systems is shown in Fig. 6.
In addition, we employed the color map method to represent the distance change between the two prediction results. As shown in Fig. 7(a), each color represents a difference value in millimetres as indicated by the color scale; hot colors such as red, orange, and yellow denote bulges in our simulation results, while cool colors such as indigo, blue, and violet denote concave parts relative to the simulation results of Proplanner.
In another experiment, a set data was collected from a patient undergoing jawbone plastic surgery. Figure 8 shows the shape of the face before and after osteotomy and advancement of the mandible. The surgical procedure greatly affected the facial    appearance. Part of the mandible was moved into a position where the upper and lower teeth can clench normally. A further small part of the mandible moved to provide a good appearance. We imitated the surgical procedure after capturing these parts of the mandible and measuring the displacement of these parts when moved into position after osteotomy as shown in Fig. 9. Figure 7(b) uses a color-map to indicate the distances between corresponding points of the predicted and actual postoperative facial skin surface. Keeping in mind the stability and predictability of maxillofacial surgery, these values show an acceptable accuracy. See Table 1: the first three columns show the number of vertices, faces, and elements used in the facial soft tissue model. The time needed to calculate the new facial appearance and to generate the tetrahedral mesh are summarized in the next two columns. In the last two columns, the mean and variance of errors are shown after simulation. The experiments ran on a standard personal computer with an Intel Core i7 CPU at 2.50 GHz and 4 GB memory.

Discussion
The test results of these experiments demonstrated that our system can provide a good prediction of the post-operative facial appearance with texture. However, the prediction results inevitably have small differences from the post-operative data. Analysis reveals some drawbacks that lead to the inaccuracy of our simulated postoperative facial appearance: • Due to the limitations of numerical calculations, our data model is based on skin surfaces and the mandible. The remaining materials such as fat tissue and muscular tissue are omitted as their presence would produce too great a density of elements. Better results would be obtained if these materials were incorporated. • The model parameters such as Young's modulus  and Poisson's ratio are recommended values; measuring precise human tissue material parameters is a difficult task. • Some materials are inevitably needed to reconnect partially cut bone or to fix bone cuts in actual surgery, but these are not modelled in the simulated surgery.

Conclusions
The results of the first experiment show that our system can achieve similar effects to a commercial system called Proplanner. In addition, our system is able to show an appearance with texture and avoid an obvious fault at the edge of the deformation area. In the second experiment, we succeeded in imitating a surgical procedure with a complete set of actual patient data. In these experiments we achieved satisfactory results of CMF surgery simulation. The test results demonstrate that fast and good prediction of post-operative facial appearance with texture is feasible. Our surgical simulation system is applicable to computer-assisted CMF surgery.
In future, we plan to implement GPU support into the simulation system to improve performance. We will also consider how to include further soft tissue information such as muscle anatomy and fat tissue in our data models. Open Access The articles published in this journal are distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
Other papers from this open access journal are available free of charge from http://www.springer.com/journal/41095. To submit a manuscript, please go to https://www. editorialmanager.com/cvmj.