Finite Element Modeling for Orthodontic Biomechanical Simulation Based on Reverse Engineering: a Case Study

In order to improve the validity and feasibility of the solid model of oral tissue, a new method is provided based on reverse engineering. Biomechanical simulation with FEM is an important technique for orthodontic force analysis and evaluation, as well as treatment design. As the base of FEM simulation, the solid geometrical models of oral tissue including tooth, Periodontal Ligament (PDL) and alveolar bone are difficult to construct through conventional solid modeling methods because the oral models are very complicated in geometry and topology which are generally represented as triangular meshes. But in many cases, solid model is necessary for FEM performing. So how to construct the solid model of oral tissue with good quality and efficiency is a problem should be faced and solved in orthodontic biomechanical simulation. Aiming at this problem, reverse engineering modeling is adopt to transfer triangular meshes to four-sides surface model and together with techniques of medical image processing, three dimensional triangular meshes calculating, surface fitting and solid modeling, the solid geometrical model of oral cavity for FEM analysis is constructed from CT (computerized tomography) images. With a simulation case of rat molar movement test, the whole procedure of solid modeling based on reverse engineering and some main techniques are presented and the validity and feasibility are proved also by the simulation results.


INTRODUCTION
Force system supplied by orthodontic appliances is composed by complex three dimensional forces and moments in Cartesian coordinate system.Clearly known these data of forces and moments can guide to design more perfect treatment plan and make accurate tooth movements (Badawi et al., 2009).Experimental test is a valid method for acquiring the macro mechanical data applied on teeth by archwires, but special device generally is needed with force transducers (Chen et al., 2010).With these tested data, accurate control of orthodontic force and personalized treatment plan can be realized to some degree.
In the area of orthodontic biomechanics, besides the macro force system, micro mechanical data such as stress distribution on tooth, PDL and alveolar bone are important too because they are critical to know the biology properties of tooth movement, root resorption and bone remodeling (Reimann et al., 2009).Because the structure and material properties of oral cavity are very complicated, analytic methods are hard to be used to solve the problems of micro mechanics.Computer simulation with FEM can adapt complex structure and has been widely used in biomechanics, including orthodontic (Bourauel et al., 1999).In recent years, FEM has become the significant method to quantify the micro biomechanical system supplied by orthodontic appliances and to study the relationships among the tooth movements, root resorption and bone remodeling (Rudolph et al., 2001;Cattaneo et al., 2005;Penedo et al., 2010).
Geometrical model is the base of FEM, but conventional methods of solid modeling for geometry model construction in FEM analysis are not appropriate for tooth, PDL and bone.This is due to the inherent complexity of their structures.In general, there are two modeling methods existing used for FEM.One is that directly reconstruct the model from CT images via medical image processing software, such as MIMICS (Materialise NV, Belgium) (Cattaneo et al., 2005).This method uses triangular meshes as the geometrical model imported into FEM software, so only surface data can be represented, which will result difficult on some simulations between contacted solid models, such as sliding mechanics analysis (Ammar et al., 2011).Another method is transfer the triangular mesh model to solid model using combined method of reverse engineering and conventional CAD (computer aided design) method of solid modeling, such as using CAD software Unigraphics (Siements, German).For example, Li et al. (2005) proposed a routine for constructing solid model of tooth from CT scan, which uses point set on cross section of tooth surface to create spline curve and then fit the group of section curves to create NURBS (Non-uniform rational b-spline) surface and finally construct the solid model of the tooth.With this method, the solid model of tooth can be created and imported into FEM software such as ANSYS (ANSYS Inc., USA) with minor errors.But some small detail features on the tooth are ignored because the number of sections is finite, especially on the tip of crown and end of root.On the other hand, it is difficult to create PDL model based on tooth root model via offset algorithm of root surface because of the complication and small curvature existing in tooth.Moreover, because solid modeling based on section curves are incapable to the model with bifurcation, so another problem is that multi-roots tooth such as premolar and molar is very difficult to be constructed with this method because of their complicated topological structures (Ziegler et al., 2005).With section curves, the roots and crown generally are modeled separately and then try to combine them to form a whole model, which is very difficult to realize because the linkage areas are not easy to form and the errors cannot be avoid in that areas.
Reverse engineering, as a typical method for product development in industrial area (Várady et al., 1997), has been used in medical applications in recent decades (Gibson, 2005), such as customized implant development, customized surgical plan design and surgical guide development (Maravelakis et al., 2008;Hieu et al., 2010).In fact, reverse engineering can be more useful rather than supply section data of tooth for orthodontic simulation modeling.Another strategy of transferring triangular mesh model to CAD solid model via reverse engineering is a combination of region segmentation on mesh model and NURBS surface fitting based on whole model regions (Ma and Kruth, 1998).This strategy can supply a whole surface model with G 1 continuity between NURBS surface patches to CAD software, so no combination work between different parts on one model are needed and the modeling efficiency and model accuracy can be improved much (Piegl and Tiller, 2001).Taking an experimental rat molar as an example, this study will introduce the modeling procedure of multi-roots tooth, PDL and bone for FEM analysis and the work of each phase in reverse engineering, CAD modeling and FEM analysis are discussed in detail.With this case study, the validity and feasibility of proposed strategy will be proved.

METHODOLOGY
In order to construct precise solid model of tooth and PDL, reverse engineering method is adopted to create NURBS surface model from mesh model and then solid model for FEM simulation can be created easily based on the surface model.The procedure of solid model reconstruction from section data is not controllable in accuracy because the numbers of sections and discrete point sets on each section are determined by designer's experience rather than rational standard.But with the technique of regional surface fitting, the errors in the procedure of transferring triangular meshes to NURBS surface is under control since in each step the error can be pre-set.So combining the techniques of image processing, reverse engineering reconstruction based on surface fitting and solid modeling, precise solid model of the tooth with multi-roots for FEM simulation can be created.
The procedure of this solution is shown in Fig. 1, four modules are included: image processing and meshes model reconstruction with MIMICS, NURBS surface model fitting with Geomagic (Geomagic corp., USA) and solid model creating via Unigraphics and simulation in ANSYS.Between two adjacent modules, the exchanging data format is specified with different type file: first, the CT image scanned from spatial CT, Cone beam CT or micro-CT is imported into MIMICS via DICOM (Digital Imaging and Communications in Medicine) format file or BMP image file; then the mesh model is saved as STL (Stereo lithography) file; the surface model is saved as IGES (Initial Graphics Exchange Specification) file; and the solid model is exchanged by Unigraphics part file.
The procedure can be described as following: First, the DICOM file of images is imported into MIMICS and then the triangular mesh model is reconstructed through image processing including HU threshold value setting, region growing and 3D mesh calculation.Besides these auto operation tools, some manual operations such as mask editing can be adopted to create required tissue with more accuracy.Then, the mesh model is imported into Geomagic and with Fig. 1: Procedure of solid modeling from images based on reverse engineering everse engineering tools including noise remove, mesh smooth, patches constructing and surface fitting, the NURBS surface can be reconstructed.And then, the surface model is imported into Unigraphics and sewed within small error limitation to create solid model.Finally, the solid model is imported into ANSYS to preprocess and simulation.
Obviously, the primary work of this strategy lies in the creation of mesh model and surface model and the phase of solid modeling is simplified too much, which guarantees the accuracy of the reconstructed model since the forward design procedure of solid model is more appropriate for free concept design rather than complied to existing data.The BOOL operations are performed in ANSYS rather than Unigraphics to avoid bringing more errors into simulation because ANSYS has high requirements for solid model.The solid models with low quality often induce more difficult in mesh creating and even the solution cannot proceed in ANSYS.

MESH MODEL RECONSTRUCTION
In industrial applications, conventional reverse engineering starts from physical part with surface data collection via laser scanner or contact scanner.After the mass point cloud of surface is acquired, techniques of reverse modeling including noise remove, region segmentation and surface fitting are used to construct the NURBS surface model and finally the solid model of the original part is constructed (Várady et al., 1997).But in this case of medical application, the original data of tissue is scanned by CT machine and the data type is medical image, so the data processing techniques for medical model reconstruction are different from industrial applications.
Case report: background: This case is to investigate the relationship between the orthodontic stress distribution on the tooth root surface and root resorption volume during the rat orthodontic tooth movement.After acclimation for 1 week, 25 5-week-old male Sprague Dawley rats (Specific Pathogen Free level 3) fed with a powder diet and water ad libitum until sacrificed.In each animal, the maxillary left first molar was extracted under anesthesia using intraperitoneally injected chloral hydrate (2 ml/kg).A Nickel-Titanium (NiTi) coil spring 0.2 mm in diameter (IMD, Shanghai, China) as shown in Fig. 2, was ligated between the maxillary left second molar and both incisors in each rat using a 0.08-inch ligature wire.The spring was activated for approximately 1 mm to produce a continuous force of 10 g to move the maxillary left second molar forward.The ligature wire was then secured with bond adhesive (Transbond, 3 M Unitek, Monrovia, Calif.) on the incisors.Spring retention was checked daily to ensure the stability of the applied force.The incisors, molar and spring were cleaned and irrigated with tap water as needed to prevent potential trauma and irritation to the gingival and periodontal tissues.

Data collecting of micro-CT scan:
The head of each animal was scanned using an in vivo micro-CT system (SkyScan 1076, Kontich, Belgium) after spring installation.The rat was anesthetized and secured in a carbon composite animal holder and the palate was parallel to the stage.The spring was removed and reattached after scanning.The head was scanned through 180º of rotation at 0.5º step increment with 2 sec/degree with the resolution of 18 µm/pixel.The raw data were further reconstructed to provide axial cross-section images using NRecon software, version 1.4.4 (Skyscan).Approximately 1000 axial cross-section images were collected per time point for each animal and images were converted into 16-bit bitmapped TIFF images with a resolution of 1024×1024 pixels and recorded on a disc for transferring and using.

Mesh model reconstruction from images:
The scanned images are imported and processed in MIMICS v11.11.With tools of image processing, including setting threshold of HU value, region growing and 3D calculation, the mesh model of whole data is simple to be created.But if a tissue such as the single root is required to be reconstructed separately, more manual work are needed, because in the images, the border between two different tissues is not clear often.In Fig. 3, a compare of different types of CT images is given.Figure 3b is a human oral image scanned by spiral CT machine with 0.1 mm scan space and Fig. 3c is a human oral image scanned by CBCT machine, both have not clear borders between tooth, PDL and bone, consequently, many manual operations of mask editing to remove the adjacent pixels on each image are needed to do if sole tooth is required to be separated from whole model.Benefitting from high resolution, the micro-CT image in Fig. 3d shows clear boundaries between different tissues and the region of PDL is clearly emerged between the roots and bone.So constructing a single tooth is easy from micro-CT images, only the pixels of adjacent crown marked in the ellipse of Fig. 3d are needed to be erased.
The experimental molar with four roots is reconstructed with 92208 triangular meshes, in Fig. 4a and b, which includes some internal structures.Because triangular mesh has good appropriate property to complicated structure, the shape of the molar including complicated occlusal surface on crown is represented adequately (Fig. 4c).The mesh model is saved as STL file to be imported into Geomagic V10.0 for NURBS surface construction.

SURFACE MODELING VIA REVERSE ENGINEERING
In reverse engineering, in order to get 4-sides region surface patches of tooth from scanned point cloud or processed triangular mesh, different strategies are adopted in different commercial software.One of the popular strategies is a routine from point to curve and more to surface, which is adopted by some typical reverse engineering commercial platforms such as Image Ware (Siemens, German) and the reverse module in CATIA (Dassault System, France).These platforms are mainly used in industrial applications because of the relative rule structures of industrial parts.Geomagic, as a popular reverse engineering platform, adopts different strategy of utilizing surface fitting on whole parts based on mesh model, which is more appropriate for medical applications because of the complicated structures of tissue.

Mesh preprocessing:
In order to get good surface model with high quality and limited errors, the mesh model is preprocessed in Geomagic firstly (Fig. 5).From transparent display of the model in Fig. 4a, singular structure is existing in this molar model, which are some internal structures connected with the surface meshes through some blending areas.These internal meshes should be removed otherwise the fitted surface will be affected on quality and accuracy.To remove the internal structure, the connecting mesh on blending areas are deleted and the mesh model can separated to outer model as shown in Fig. 5b   Because the CT scan is not a continuous procedure, the reconstructed model from CT images is very coarse on the surface when interpolation is needed for compensate the scan space between adjacent layers.Another reason for surface coarse meshes is that most operations on image processing are discrete procedures, which makes pixels discontinuous on one layer.Geomagic supplies a reliable tool for mesh smooth with error control, which can smooth the meshes with good quality within required error limitations.The smoothed molar shown in Fig. 5d is with following error limitations: average distance between smoothed and original meshes is 0.009 mm, maximum distance is 0.049 mm and standard deviation is 0.005 mm, which are sufficient for simulation.
Region segmentation: Because parametric surface used for solid modeling such as NURBS surface is constructed on 4-sides rectangular region, the triangular meshes of the whole model should be segmented to adjacent regions with 4-sides topology.In order to guarantee any adjacent surface patches G 1 continuity fitted from under meshes, bi-cubic B-spline surface is adopted: where, ܰ ,ଷ (u) and ܰ ,ଷ (v) are the B spline base functions, ܸ , (I = 0, 1,…,m: j = 0, 1, …,n) is the control points of fitted surface through r times iterations.The u and v are knot vectors along two parametric directions: Based on parametric surface expression, the whole surface model is solved with G 1 continuity constrains (Kruth and Kerstens, 1998).Geomagic supplies automatic solution for surface fitting and also supplies some auto tools for region segmentation and mistake check.With the help of manual operations, the region segmentation can be done quickly.Figure 6 shows the segmentation procedure of this case: first, recognize and extract some feature curves on the molar by auto curvature estimation and then segment the model to 294 regions according to surface curvature automatically, as shown in Fig. 6a; but there are some errors existing which can be found by checking tool and needed to be modified manually.A typical error on region segmentation is triangular shape existing on some regions, which should be adjusted to rectangular region by moving some linking points.In Fig. 6b and c, the crosses points 1, 2 and 3 are moved to new places and the errors then are corrected.Other problems such as more than 5 sides linked with one point and too small angle between two sides also can be found automatically and then adjusted manually.The corrected regions are shown in Fig. 6d, compare to original regions in Fig. 6a, the region sides and points are arranged more regular.
Free form surface fitting: Based on region segmentation with good quality, the surface fitting can be done quickly.According to the curvature of the molar, discrete points on each side is specified to 10, which determines the numbers of interpolation points and control points of each NURBS surface patch, as shown in Fig. 7a.The fitting error is specified 1e-5, which is the error of G 1 continuity between any adjacent two surface patches and consistent with the requirements of solid modeling in Unigraphics and BOOL operations in ANSYS.The fitted surface model is shown in Fig. 7b and the error analysis between the NURBS surface model and original mesh model are displayed in Fig. 7c, which shows the good quality of the whole model with average error 0.001 mm.The surface model is saved as IGES file for importing into Unigraphics for solid modeling.

SOLID MODEL CONSTRUCTING AND PRE-PROCESS FOR FEM SIMULATION
Solid geometrical model is a universal expression in CAD, CAM (compute aided manufacturing) and FEM simulation.CAD software such as Unigraphics supplies several solid constructing methods including analytic modeling (cylinder, cone, block, etc.), feature modeling (extruding, swiping, revolving, etc.) and stitching enclosure free form surfaces patches.In simulation platforms, solid model can be used for BOOL calculations, meshing and material applied.

Solid model constructing via unigraphics:
The IGES file of surface model of molar is imported into Unigraphics NX 6.0 and then these surface patches are stitched together with 1e-5 errors to construct the solid model.If the continuous constraints between any two adjacent surface patches cannot satisfy the requirement of stitching, the solid model cannot be constructed,  In surface model, only surface is existing, inner is empty (Fig. 8a).After all of the surface patches stitched together, the inner of the molar are filled, as shown in Fig. 8b and it can be used for BOOL operation and simulation.Figure 8c is an image of Gaussian curvature distribution map of the molar surface, which shows the surface is smooth and good quality with even curvature variation.PDL modeling: It is not easy for surface model of roots to offset to create PDL, because strictly constraints on curvature for NURBS surfaces are needed for offsetting.On the other hand, mask tools in MIMCS cannot help to create ideal PDL tissue because its gray value is too low.Mesh model of molar is adopted to create PDL because triangular meshes own free property and can be offset conveniently.The distance of offsetting is 0.2 mm, which is the thickness of PDL, measured from micro-CT images.The offset mesh model is shown in Fig. 9a and with the same procedure as tooth modeling, the surface model of PDL can be created (Fig. 9b).The solid model is created based on the surface model and redundant data of crown part is cut.In simulation, alveolar bone is simplified as homogenous elastic material, so a regular block is design as substitute.The solid models of tooth, PDL and bone are shown in Fig. 9c.

Preprocessing for simulation:
The preprocessing is operated in simulation platform ANSYS after the Unigraphics part file of solid model is imported.
Utilizing overlapping of BOOL operation, the holes in the alveolar bone for the PDL and holes in the PDL for the tooth roots are created and the surplus part of PDL outside the bone is deleted.After overlapping, the unit should be uniformed to international metric m/N/Pa, so the models are scaled with factor of 1e-3 from millimeter to meter.Using mesh type solid 187, the solid models of tooth and PDL are meshed with mesh size 0.1mm and bone is meshed with size 0.5 mm.Totally 138742 meshes are created (Fig. 10).

RESULTS AND DISCUSSION
According to experiments condition, the force applied to the molar in simulation is a single direction fore with 10 g.The bone is constrained on five faces with zero displacement.The material properties of tooth, PDL and alveolar bone are referenced from the   1.
The case is a static mechanical system, with a computer workstation (intel Core 5, 8 G memory), about 40 minutes are consumed to solve.The Von Mises stress distribution of PDL and tooth roots are shown in Fig. 11a and b, the largest stress is 0.0516 Pa and 1.556 Pa.The areas with large stress exist on the distal-top side of roots, which is coincident to experimental observation of root resorption.
In this case, the solid modeling of molar is based on revere engineering with NURBS surface fitting, which conquers the challenge on modeling for multiroots with complicated topology.The modeling method based on sectional curves is capable on single root tooth (Li et al., 2005).But for multi-roots tooth with complicated topology, if each part is modeling separately, the link between the crown and roots are difficult to create.
The reverse engineering software, Geomagic, supplies powerful tools on 4-sides region segmentation and surface fitting, so the procedure of NURBS surface construction is not time consuming.In this case, the time for creating NURBS surface of whole molar is less than one hour.But in modeling method of section curves, some time-consuming steps including discrete points on original section curves, spline curves reconstruction and surface lofting are needed and each steps are complicated and needed be done carefully.

CONCLUSION
Till now, several cases including rat and human oral cavity scanned from spiral CT, CBCT and micro-CT have been completed with proposed surface fitting method.The validity of FEM simulation with this kind of model is proved by experimental results.From the applications of the proposed method, the following conclusions can be drawn: • Constructing NURBS surface directly from triangular mesh can be more efficient and feasible than conventional method.• Surface fitting from triangular mesh for whole model is more appropriate for medical application, because this method can create surface with complicated topological and geometrical structures directly.• The errors in each step are controllable, which guarantees the performability and reliability of FEM simulation.

Fig. 3 :
Fig. 2: Rat molar movement experiment Fig. 4: Mesh model (a) Rendering display, (b) Rendering and edge display, (c) Meshes on model and internal model.Then filling the holes on the outer model gets the final mesh model without singular structure, as shown in Fig. 5c.

Fig. 6 :Fig. 7 :
Fig. 6: Region segmentation for surface fitting (a) Feature curves extraction by curvature estimation and auto region plotting, (b) Error of triangular region, (c) Error correction, (d) Plotted region without errors

Fig. 9 :
Fig. 9: PDL construction (a) Offset mesh model, (b) Surface model, (c) Solid models which is the reason of high requirements needed in region segmentation for surface fitting.Figure 8 shows the difference between surface model and solid model.In surface model, only surface is existing, inner is empty (Fig.8a).After all of the surface patches stitched together, the inner of the molar are filled, as shown in Fig.8band it can be used for BOOL operation and simulation.Figure8cis an image of Gaussian curvature distribution map of the molar surface, which shows the surface is smooth and good quality with even curvature variation.
Figure8shows the difference between surface model and solid model.