3D Intelligent Scissors for Dental Mesh Segmentation

Teeth segmentation is a crucial technologic component of the digital dentistry system. The limitations of the live-wire segmentation include two aspects: (1) computing the wire as the segmentation boundary is time-consuming and (2) a great deal of interactions for dental mesh is inevitable. For overcoming these disadvantages, 3D intelligent scissors for dental mesh segmentation based on live-wire is presented. Two tensor-based anisotropic metrics for making wire lie at valleys and ridges are defined, and a timesaving anisotropic Dijkstra is adopted. Besides, to improve with the smoothness of the path tracking back by the traditional Dijkstra, a 3D midpoint smoothing algorithm is proposed. Experiments show that the method is effective for dental mesh segmentation and the proposed tool outperforms in time complexity and interactivity.


Introduction
In the context of the emerging computer-aided medical, it is possible to acquire the 3D digital model of dental and to design and manufacture dental implant guides. Dental implant guide design systems based on digital technology have developed rapidly over the years. Typical processes of dental implants design include the following: (1) obtaining the digital model of the tooth through a traditional impression and 3D scanning and then obtaining the patient's cranial CT; (2) tooth segmentation from dental mesh, tissue segmentation, and reconstruction based on CT; (3) personalized dental implant guide design; (4) engineering analysis and manufacture; (5) implement treatment. e teeth segmentation plays a crucial role in digital dental systems (steps 2, 3, and 4). e quality of the digital dental mesh is mainly dependent on the digital dental mesh acquisition method (three-dimensional scanning technology). erefore, it is crucial to design a user-friendly specialised tooth segmentation tool, whose quality of the segmentation is controllable.
ere are three types of segmentation methods based on dental mesh: automated methods, semiautomated methods, and manual methods. Automatic methods do not require user interaction and are very convenient. However, due to the lack of controlling the quality of segmentation, such methods do not meet the accuracy requirement. Although manual methods can obtain accurate results through users' intersection, they also have many shortcomings, such as tedious and time-consuming. Semiautomated methods can keep the balance between the accuracy of segmentation and the user effort for interaction. However, the general semiautomated methods are not directly suitable for dental mesh segmentation because of the unique geometry of teeth and multiple teeth arrangements on dental mode. Moreover, existing semiautomated methods proposed to handle dental mesh have the shortcoming that the positions of interactions are not suitable for improving the segmentation accuracy.
is paper aims to develop an interactive tool for dental mesh segmentation, called 3D intelligent scissors, motivated by a user-friendly segmentation tool [1,2]. e wire between two points inputted by user interaction is regarded as a part of the segmentation boundary. e main contribution of this work is threefold: (1) An improved intelligent scissors tool for triangle mesh segmentation has lower time complexity, which meets the requirements of real-time interaction and is specially optimised for dental mesh segmentation (2) e tool requires less interaction and can acquire better segmentation (3) An algorithm of 3D discrete curve smoothing on triangle mesh is proposed, which subtly transforms the complex problem of smoothing the curve on the surface into a simple problem of smoothing the curve on the plane

Related Work
e dental model is a triangle mesh and can be obtained by 3D scanning. Previous work about dental mesh segmentation and general mesh segmentation is reviewed. e proposed method depends on the geodesic line, and the review of the geodesic line is given at the end of the section.

Mesh Segmentation Approaches.
ere are some specific approaches [3][4][5][6][7] that fully exploit the dental characteristics. ese methods are either interaction-intensive or resultinaccurate. Most approaches about the dental mesh segmentation reference the methods of general mesh segmentation.
Mesh segmentation is a fundamental research topic in geometry processing and computer graphics. Numerous general mesh segmentation approaches have been proposed. e manual segmentation tool is user-unfriendly, usually as a tool for obtaining benchmark segmentation [8]. Most of the mesh segmentation algorithms are automated or semiautomated and can be briefly classified as region-based and boundary-based [9].
Boundary-based methods aim to find the boundary or contour of each segmentation part. In order to select the best boundary, Golovinskiy and Funkhouser [18] proposed a method called randomized cuts. However, due to the lack of user control, the segmentation quality is poor. Zheng and Tai [19] proposed an interactive method, according to which the user interacts with the segmentation process by drawing one or more strokes across the desired boundary. Another interactive boundary-based method, called dot scissor, is discussed in [20]. With the tool, the user's effort is reduced to placing only a single click where a cut is desired. Although the interaction needed by the methods [19,20] is little, the segmentation boundary does not fit well with the user's intentions. Zhuang et al. [1] proposed a live-wire mesh segmentation tool, with defining a wire through two endpoints entered by the user as the segmentation boundary which is the shortest path-based anisotropic metric over the surface. e method balances the ability to control results and the user effort for interaction, but the method is timeconsuming, which is intolerable as it is a real-time interactive tool. In detail, the method needs mesh embedding and local subdivision at initialisation, and computing wires is time complex because it uses the exact method [21]. Moreover, the tool fails in teeth segmentation when the seeds are far apart ( Figure 1).

Computing Geodesics Approaches.
e idea of the livewire based algorithm for mesh segmentation is to define a wire between two points on the surface as a segmentation boundary, and the critical problem is to compute the wire quickly. e geodesic-based suitable anisotropic tensor can be used as the wire, and it is a special kind of geodesic on the surface. Now, we review the work on geodesic based on triangle mesh. ere are two types of methods: exact method and estimated method [21].
Chen and Han [22] proposed an exact algorithm based on the idea of one angle one split. is algorithm complexity is O(n 2 ). MMP [23] algorithm is another method for computing exact geodesic-based continuous Dijkstra, and its average time complexity is O(n 1.5 log n) but in the worst case is O(n 2 log n). Ying et al. [24] proposed the SVG (saddle vertex graph) algorithm, which prebuilds the saddle vertex graph and can compute the shortest path through the Dijkstra algorithm. e time complexity of the algorithm is O(D n log n), D << n. In short, the exact method is time complex, although the SVG algorithm is better in time complexity. Moreover, all the above algorithms need to fulfil the triangle inequality everywhere, and violation of the triangle inequality is the typical behaviour especially for high degrees of anisotropy [25].
Estimated methods generally obtain the approximate solution of the discrete geodesic by solving the partial differential equation. Generally, it is faster than the exact methods. Among estimated methods, geodesic in heat [26,27] performs very well, but it is Euler metric based. Yang and Cohen [28] extended the method by adding the variable heat transfer coefficient, and it can lead the geodesic along the feature area. Combining the anisotropic Laplace-Beltrami operators proposed by Andreux et al. [29] with the heat method, the anisotropic heat equation can be solved and then the features sensitive geodesic can be acquired. In short, all the above algorithms are faster, but they have poor robustness because solving PDE is unstable in the condition of high degrees of anisotropy.
Another estimated method is named short-term vector Dijkstra (STVD) [25]. e method uses a new version of update function by exploiting a vector-valued short-term memory that aims to improve accuracy. However, the results are unsmooth. e conclusion of the literature review is that live-wirebased segmentation tool is user-friendly, but existing methods are time-consuming and do not exploit dental characteristics.

Our Method
Inspired by the intelligent scissors [2] in image segmentation, Zhuang et al. [1] proposed a live-wire method for mesh segmentation. In the method, the geodesic between two points on the surface is defined as a part of the segmentation boundary. However, it consumes time in computing exact geodesic and is meaningless in the segmentation application.
Similar to the live-wire segmentation, in this paper, the proposed teeth segmentation tool starts with a single click on the triangle mesh to select a seed point. en, the program traces back a smooth path between the seed point and the current position of the mouse in real time.
For real-time interaction, a method to fast compute the geodesic is presented. Besides, two anisotropic metrics are defined, which can force the geodesic line lying at the feature of the surface. Unfortunately, because of using an inexact method based graph search [30], the path along the edge is often unsmooth, so a 3D curve smoothing algorithm is proposed. e steps of obtaining a complete segmentation boundary are shown in Figure 2, which mainly include the following ones: (1) Initialization: computing the minimum and maximum curvature and the normal and the tensor-based anisotropic metric vertex (Section 3.1) 3.1. Anisotropic Metrics. In this paper, the "shortest" path between two points (geodesic) as a part of the segmentation boundary is used. Obviously, the "shortest" is not the shortest in the usual sense (Euler metric). In order to define the concept of the "shortest", a metric tensor g is introduced.
In a Riemannian manifold M, tangent plane T a M is local Euclidean representation of manifold M around point a. Moreover, in M, b is another point around a. en, the length between a and b in M based the metric tensor g a is approximately defined as follows: (1) e metric tensor g a is regarded as a local orthogonal coordinate system in the tangent plane T a M (Figure 3), where e 1 and e 2 are the eigenvectors of g a , λ 1 and λ 2 are the eigenvalues of g a , and θ is the angle formed between x and e 1 , then formula (1) can be rewritten as follows:

Computational and Mathematical Methods in Medicine
If the eigenvectors are fixed, the length is determined by λ 1 and λ 2 . erefore, controlling the value of λ 1 and λ 2 can diminish the distance between two points a and b in the feature area. Previous researches [1,31] suggested that the eigenvectors should be aligned with the local curvature directions (minimum curvature and maximum curvature).
In this paper, e 1 is aligned with the direction of minimum curvature and e 2 is aligned with the direction of maximum curvature. k min and k max denote the minimum and maximum curvature at a point of the triangle mesh. |k max | − |k min | is smaller in ridges and bigger in valleys, and |k max | + |k min | is smaller in valleys and planar area (Figure 4). e tangent direction of wire along valleys is similar to the direction of the minimum curvature ( Figure 5(a)). So for valleys, λ 1 and λ 2 are set in formula (3), which is named as valley metric. Dividing by k max makes λ 1 and λ 2 scale-independent. It is obvious λ 1 (Figure 4(b)) and sin(θ) (Figure 5(a)) are small in valleys, and according to formula (2), the above eigenvalues make the wire along the valleys shorter than that along others: For ridges, λ 1 and λ 2 are set in formula (4), which is named as ridge metric. In ridges, k max is similar to k min (Figure 4(a)), which means λ 2 is close to 1. Furthermore, in other regions, λ 2 is much more than 1. Besides, because the direction of wire along ridges is similar to the direction of maximum curvature ( Figure 5(b)), the smaller λ 2 makes the wire along the ridges shorter. ε (10 − 4 in the paper) prevents division by zero. e definitions of λ 1 and λ 2 reveal their scale-independency:

Anisotropic Dijkstra.
In order to fast compute the geodesic from all vertices to the source, an anisotropic Dijkstra algorithm is proposed. Its pseudocode is shown in Algorithm 1. e data structure of the triangle mesh of the teeth as one input of the algorithm is half-edge [32]. Every vertex of the triangle mesh has a set of triples expressed as (dist, final, and pred), where dist is the distance between the vertex and the source, and final marks whether the vertex is visited and pred represents the previous vertex of the vertex in a path.
In the initialisation, dist of source point is zero and all other vertices are infinity. A priority queue Q is employed, which is ordered by dist. In the beginning, there is only source vertex in Q. en, until Q is empty and in the process dist and pred of all vertices are updated, the main loop is running. In the main loop, the distance between two adjacent vertices-based anisotropic metric is calculated. e  (2) source Q, Q is a priority queue ordered by dist; (3) while not Q empty do (4) υ Q; (5) υ.final = true; (6) for w adjancen υ and w.final ≠ true do (7) if w.dist > υ.dist + l g (υ, w) then (8) w.dist = υ.dist + l g (υ, w); (9) w.pred = υ; (10) if not w in Q then (11) w Q; ALGORITHM 1: Anisotropic Dijkstra.
Computational and Mathematical Methods in Medicine metric tensors of the two vertices v and w are g v and g w , and distance between v and w is computed as follows: e geodesic between any vertices and the source vertex is tracked back by pred as shown in Figure 6(a). Figure 6(a), the wire we trace back is jagged, because it consists of the edges of the triangle mesh determined by the Dijkstra algorithm. However, for teeth segmentation, the geodesic should be smooth as the boundary is smooth. To achieve this goal, Polthier and Schmies [33] proposed the "Geodesic Euler Method" and "Geodesic Runge-Kutta Method" based on the gradient of the geodesic field. e methods regard the geodesic tracking back problem as an ordinary differential equation (ODE) solving problem with the initial value, which requires the geodesic field is smooth. However, the geodesic field computed by "Anisotropic Dijkstra" does not satisfy the smooth condition because of its strong anisotropy and weak accuracy. erefore, Polthier's methods cannot get the correct geodesic. Inspired by common sense that replacing the two endpoints with the midpoint makes polyline smooth on the plane, a 3D midpoint algorithm to smooth curvature on triangle mesh is proposed and the pseudocode is shown in Algorithm 2.

Smooth 3D Curve on Dental Mesh. As shown in
"Midpoint method" is not suitable for smoothing the curve on the 3D triangle mesh because the midpoint of two points in the curve may not locate on the triangle mesh in 3D space. erefore, in the proposed algorithm, k (usually 3 to 5) neighbourhoods of the curve C on triangle_mesh are extracted (Figure 6(a)), named k_ring_mesh. Using the method of least squares conformal maps (LSCM) [34], the expanded triangle mesh on the plan is named by k_ring_-mesh_2D ( Figure 6(b)). C is constructed by the vertices of triangle_mesh. e curve C is mapped to a 2D curve C 2D that is constructed by the vertices of k_ring_mesh_2D. C 2D is smoothed with the "midpoint method" to S_C 2D (Figure 6(c)). S_C 2D is mapped back to a smooth curve S_C (Figures 6(d) and 6(e)). e method mainly includes the following steps.
Firstly, the interactions of each edge of k_ring_mesh_2D and S_C 2D are computed.
Secondly, one edge of k_ring_mesh_2D intersects S_C 2D at pt_int 2D . Two vertices of this edge are defined as index 1 and index 2 , respectively. en, the information of intersection is recorded in a set of triples (index 1 , index 1 , and α), where α is calculated by the following formula: Finally, S_C 2D is mapped back to S_C by the following formula: Additionally, pt 2D and pt 3D are the vertices of k_ring_mesh_2D and k_ring_mesh_3D, respectively. And, all the pt_int 3D form S_C. s is a parameter of the algorithm to control the smooth strength.

User Interface.
Our teeth segmentation tool is developed in C++, depending on libigl [35] and OpenMesh [32]. OpenMesh and libigl mainly provide functions such as halfedge data structure, reading and writing the file, rendering, UI components, and computing some discrete geometry quantities. To begin with, the user holds the Ctrl key and presses the left button of the mouse on the triangle mesh of the teeth model and then a seed is selected. When the Ctrl key is held and the mouse is moving on the triangle mesh, the "shortest" path between the seed and current position of the mouse is acquired and rendered in real time. If the path is realised as "good," it is confirmed with the Ctrl key and the left mouse button click. When a previous path is regarded as not "good," it is cancelled by the right mouse button click with holding the Ctrl key, and the previous seed becomes the current one. It is ended by the right button click with Ctrl + Shift holding. Our program stores all the seeds. When the mouse hovers over the seed and the seed is highlighted, then the seed can be dragged to modify the path. During this process, the path is rendered in real time. When the Ctrl key is not held, the mouse is responsible for common 3D intersections including translating, rotating, and scaling.
All the tasks of teeth segmentation, including segmentation between teeth and gums, teeth and teeth, the occlusal surface of the teeth, need the path along valleys or ridges. A radio button is employed to switch the mode (along valleys or ridges). A radio button is to change the metric (formula (3) and (4)) defined in Section 3.1.

Setting.
We conducted experiments on clinical dental meshes which are acquired by traditional impression and then 3D scanning using AutoScan-DS100+. All the programs are implemented in C++ and compiled in Microsoft Visual Studio 2017 and run on a computer with 8 GB RAM, Intel ® Core ™ i7-4790 3.6 GHz CPU, and Windows 10 64 bits system. Real-time interactivity and easy-to-use are significant for the segmentation tool, so they are evaluated in our experiments by comparing with Zhuang's method. In order to evaluate our method, seven dental meshes are used and their information is listed in Table 1.

Results and Analysis.
e results of our method are shown in Figure 7. In this experiment, the valley metric is used to segment the teeth and gums from dental mesh with parameters k � 3 and s � 5. e results of Zhuang's method with his Max metric are given in Figure 7. Furthermore, the parameter settings follow the author's suggestion. Both methods contain three operations: initialisation, seedling, and tracking back path. e running time for every operation is listed in Table 2. In the experiment, the principle of seeds selection is that the wire between two seeds with the furthest distance is as close as possible to the valleys.
As shown in Figure 7, the qualities of the segmentation boundaries acquired by our and Zhuang's method are similar. However, the number of seeds acquired by our method is less than that of Zhuang's method (shown in Table 3). It is evident that, in most cases, the seeding time is much less than one second. e tracking back time is slightly longer than that of Zhuang's method, but a few milliseconds of tracking time does not affect the real-time interaction. So (4) for m = 1 : s do (5) initialize pts as a point array; (6) foreach pt i 2D in C 2D do (7) (8) update C 2D using temp_pt; (9) (10) k_ring_mesh_2D; (11) compute α based Equation (6); (12) (index 1 , index 2 , α) info_int 2D ; (13) initialize S_C as a point array; (14) foreach (index 1 , index 2 , α) in info_int 2D do (15) compute pt_int 3D based Equation (7); (16) pt_int 3D S_C foreach pt i 2D and pt i+1 in C 2D do 2D compute intersection pt_int i 2D of line pt i 2D pt i+1 and edges of 2D ALGORITHM 2: 3D midpoint smoothing algorithm. Zhuang's method using his Min and Max metric, respectively. e red points are the seeds, and there are only two seeds at the ends of the wire. e blue and dark red wires overlap each other in Figure 1(a). e wire obtained by the valley metric is closer to the real parting of teeth and gums, especially when the two seeds are far apart, as shown in Figures 1(b)-1(d). It means that the valley metric is more effective than Zhuang's method. Although the performance of ridge metric is close to Zhuang's Min metric shown in Figure 8, our method is faster than Zhuang's method. ere are two parameters in the proposed method for smoothing 3D curve. e parameter k controls the magnitude of k_ring_mesh, which is usually set between 3 and 5.
A smaller value of k may cause that the smoothing curve goes beyond the boundary of the mesh. Moreover, a larger value of k inevitably leads to reduction in efficiency because intersections between the curve and the mesh need to be computed. In our tests, when k is set to 3, the algorithm performs better. e parameter s is the smoothing strength, which affects the times of the taken midpoints. On the plane, the polyline will become a straight line by taking the midpoint enough times. It is also suitable for our 3D curve smoothing method on a triangle mesh. An appropriate value of s makes the wire smooth and not far from the original polyline ( Figure 9). In our tests, the parameter performs well from 2 to 10. Experimental results show that our method is insensitive to parameters.  Our method + ridge metric Zhuang's method + min metric (c) Figure 8: Occlusal surface segmentation of using our method + ridge metric and Zhuang's method + Min metric. e dark blue wire is obtained by our method + ridge metric, and the red one is obtained by Zhuang's method + Min metric.
Our method is proposed for dental mesh segmentation. Meanwhile, it can work well on other triangle meshes with visible valleys or ridges as shown in Figure 10.

Conclusion
An intelligent scissors tool for teeth segmentation in the dental mesh is proposed, which is inspired by Zhuang's live-wire method. Valley metric and ridge metric are defined to lead the wires along the valleys and the ridges. To quickly compute the geodesic, we adapt the Dijkstra algorithm for the anisotropic metric. erefore, in order to solve the problem that the path tracked back by Dijkstra is unsmooth, a 3D midpoint smoothing algorithm is proposed. e experiments show that the tool is effective for tasks of teeth segmentation in the dental mesh. Compared to Zhuang's  Figure 10: e blue wire is obtained by our tool, and the red points are the seeds. method, our method is better in time complexity and interactivity. However, the tool is poor in the segmentation of the mesh without prominent valleys or ridges. In future work, it is attractive to design a more versatile metric and to extend the tool for interactive texture mapping.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this article.