Interactive Tooth Separation from Dental Model Using Segmentation Field

Tooth segmentation on dental model is an essential step of computer-aided-design systems for orthodontic virtual treatment planning. However, fast and accurate identifying cutting boundary to separate teeth from dental model still remains a challenge, due to various geometrical shapes of teeth, complex tooth arrangements, different dental model qualities, and varying degrees of crowding problems. Most segmentation approaches presented before are not able to achieve a balance between fine segmentation results and simple operating procedures with less time consumption. In this article, we present a novel, effective and efficient framework that achieves tooth segmentation based on a segmentation field, which is solved by a linear system defined by a discrete Laplace-Beltrami operator with Dirichlet boundary conditions. A set of contour lines are sampled from the smooth scalar field, and candidate cutting boundaries can be detected from concave regions with large variations of field data. The sensitivity to concave seams of the segmentation field facilitates effective tooth partition, as well as avoids obtaining appropriate curvature threshold value, which is unreliable in some case. Our tooth segmentation algorithm is robust to dental models with low quality, as well as is effective to dental models with different levels of crowding problems. The experiments, including segmentation tests of varying dental models with different complexity, experiments on dental meshes with different modeling resolutions and surface noises and comparison between our method and the morphologic skeleton segmentation method are conducted, thus demonstrating the effectiveness of our method.


Introduction
In recent years, orthodontic computer-aided-design (CAD) systems are widely used by clinicians to prepare their orthodontics surgery. These CAD systems exploit hardware-supported computer graphics technology to effectively and efficiently plan diagnosis traditionally done manually. These systems provide information precise enough to be used in diagnosis and prognosis, thus free clinical dentists from repeated works and facilitate accurate treatment planning [1][2][3][4].
As well as in traditionally case, after acquiring digital dental model which is generated from scanning teeth of patient, dentist firstly need to accurately separate teeth from the digitized them. The image based segmentation method, however, lacks three dimensional geometrical information, complex interstices or non-convex points are difficult to extract, which leads to inaccurate cutting between adjacent teeth.
In the article, we proposed a novel, effective and efficient teeth segmentation method based on segmentation field. The segmentation field is a scalar-valued field, which is solved by a linear system defined by a discrete Laplace-Beltrami operator with Dirichlet boundary conditions, which are imposed at a set of constraint points. A concavity-sensitive weighting scheme is applied to calculation of segmentation field, such scheme make the segmentation field exhibits significant greater variation at concave regions of dental mesh.
Through identifying maximum concavity variation, the final segmentation boundaries can be extracted from these contour lines that are sampled from segmentation field.
Overall, the aims and contributions of our dental segmentation method are included as follows: • Compare to other common Laplacian weighting schemes, the adjusted weighting scheme is introduced to partition teeth from dental models accurately and efficiently.
• A strategy to assign the constraint points is proposed, according to different geometrical shapes of teeth.
• An easy-to-use interactive tool is developed to segment dental models with minimal user interactions.

Materials and Methods
A series of strategies are applied by us to accomplish the purposes mentioned in previous section, the overview of our framework as illustrated in Fig 1. The scanned dental model is firstly preprocessed to identify and to group the feature points in dental mesh. These feature points, are assigned to segmentation field as boundary constraints, have the ability to indicate useful geometry information of dental model for tooth segmentation. Through user interaction, after  this, a point on mesh is clicked to indicate which tooth the user want to separate from dental mesh, thus the target boundary constraints and the background boundary constraints are assigned according to the relevant location of the clicked point. Then, the segmentation field, whose scalar values are colored from red to blue, is computed based on the concavity-aware weighting scheme. Afterwards, a set of contour lines are sampled from the segmentation field, and the cutting boundaries between tooth-tooth and tooth-gingiva are detected from these contour lines in the concave regions. Finally, all the teeth are separated and colored from the input dental mesh based on the cutting boundaries.

Definition of the segmentation field
The segmentation field is actually a harmonic field which has been exploited to mesh segmentation [19][20][21]32]. In our method, partitioning tooth from dental model relies on these desirable properties of segmentation field. The segmentation field is a smooth scalar field and is so sensitive to concave regions that tooth boundaries can easily extracted from the contour lines which lie on concave regions. Let M = (V, E) be the given dental model, where V is the set of vertex indices, E is the set of edges and v i the location of vertex i in space. A segmentation field can be denoted by u, a solution to the Poisson equation Δu = 0, which is subjected to Dirichlet boundary conditions imposed by sites S, which can be denoted by S 2 {1, 2, . . ., n}. Δ is the Laplace-Beltrami operator, which can be discretized and resulted in a symmetric Laplacian matrix. Such symmetric operator, which firstly presented by Pinkall and Polthier [33], leads to a symmetric and positive-definite linear system, allowing for fast Cholesky factorization. The site constraints are assigned in segmentation field computation, which can be described as where s i is the prescribed value of segmentation field at site i. To obtain constrained segmentation field, we minimize the membrane energy [34] u ¼ arg min leading to the linear system where L is the concavity-aware Laplacian matrix, P is the diagonal penalty matrix The penalty factor α is a large constant used to tweak the importance of constraint satisfaction (α = 10 8 in our experiments), D is a diagonal matrix of the row sums of W, W is the weight matrix which is defined by ω ij , and the selection of weighting scheme and scalar weight assigned to ω ij are discussed in the following subsection.

Selection of weighting scheme in segmentation field
Harmonic fields with different Laplacian weighting schemes have been introduced for different application purposes, such as construction of reduced deformation models [35], transformation propagation [36] and shape approximation [37]. However, for the harmonic field that is suited for dental mesh segmentation, there are two conditions should be satisfied as follows: • The segmentation field should be sufficiently smooth so that the smooth contour lines can be extracted from such field without post-processing, this means that the cutting boundaries to separate target teeth from dental model are clear and can be directly and easily identified.
• Considering the segmentation purpose for dental models, the segmentation field is able to identify the shape variation and is sensitive to concavity regions of models where the candidate cutting boundaries lie on.
To fulfill the conditions mentioned above, the small prescribed values are assigned on the edges lying on concave regions, which leads to the scalar values of segmentation field change abruptly along concave regions, through setting the concavity-aware weighting scheme as where |e ij | denotes the length of edge e ij , which is constituted by vertex i and j, G i and G j is the Gaussian curvature value of vertex i and j respectively. γ is a small constant (we choose γ = 0.0001 for all examples demonstrated in the manuscript) to prevent the divisor in Eq (8) to be zero. We make a large variation in field data in concave regions, which leads to a phenomenon that these concave regions can be distinguished from other areas on dental mesh, through setting small constant β (β = 0.01) when either vertex i or j lies on concave regions. Vertex i is treated as a concave vertex if one or more its adjacent vertexes hold the following inequality where v adj and n adj is the position and normal of the vertex adjacent to the given vertex i, whose position and normal vector are denoted by v i and n i respectively, and θ = 0.001 is the empirical threshold value. Besides to concavity-aware weighting scheme, the segmentation field is also related to selection of site constraints S, relevant contents are discussed in the following subsection.

Site constraints in segmentation field
As described in Eq (1), a prescribed value s i of the segmentation field is assigned at site constraints i. At our case, the set of constraint sites are divided into two types, target set S t and background set S b , and Eq (1) can be extended as follow: where 1 and 0 are the maximum and minimum constraint values assigned to. The paths connecting maximum constraints and minimum constraints should pass through the desirable cutting boundaries. To a incisor with relevant smooth geometrical shape, a single segmentation field with two boundary constraints is sufficiently to separate from model. However, if we need to partition a molar with complex geometrical shape, such segmentation field with only two boundary constraints cannot accurately observe the candidate cutting lines. Fig 2 shows the single segmentation field with two constraint sites, which are illustrated by the red sphere (target constraint) and blue sphere (background constraint) in Fig 2A. It is clear that the color variation in concave seam, where is located between tooth and gingiva, is small, as shown in Fig  2B, so that the final segmentation result is not accurate (as illustrated in Fig 2C). In contrast, the multiple segmentation field which is constituted by four pairs of boundary constraints is able to find the concave seam we need to extract, as shown in Fig 3. Obviously, the color variation in the seam ( Fig 3B) is so large that the final segmentation result (Fig 3C) is accurate. Such comparison also shows that the amount and locations of boundary constraints are critical for the accurate segmentation results. An efficient way to identify the locations and amount of boundary constraints is to employ the method introduced by Yokesh Kumar et al. [38], in which the feature points can be automatically detected (see figure "Model with Feature Points" in Fig 1 with feature points denoted by gray spheres). By detection of tooth interstices [30], all of feature points can be grouped into N groups, where N is the amount of teeth. Hence, the feature points in same group are represented by the same color spheres (red spheres or blue spheres), as shown in figure "Grouped Feature Points" in Fig 1. With less user interaction, these processed feature points are automatically treated as the target constraints or background constraints, relevant contents are presented in the following section. Cutting boundary selection Extracting cutting boundaries in concave regions is the main purpose to segment dental model. The segmentation field has linear variation on model mesh so that a set of contour lines can be easily sampled from dental model (Fig 4A), and a large variation of field data can be detected in concave regions (Fig 4B). These properties of segmentation field provide an efficient mean to extract cutting boundaries from the set of contour linesfor teeth separation on dental model,  Tooth Separation Using Segmentation Field through computing the field gradient magnitudes of these contour lines and identifying the one which possess the maximum gradient.
In Fig 4C, the histogram illustrates the distribution of the amount of mesh vertex and how many vertices fall into each small interval. As shown in the histogram, the segmentation field data are mapped into interval [0, 1]. The distribution of vertex are mainly located in three intervals, (0.04, 0.22), (0.78, 0.96) and (0.46, 0.64). The three scopes correspond to the target tooth, the adjacent tooth and the rest area on dental mesh, which are colored by red, blue and green rectangles respectively. The interval spacings, I a and I b in Fig 4C, between these intervals are so obviously, which means that there are fewer vertex distributing in the wider intervals and the cutting boundaries should be located in such intervals. The Bin max (denoted by red circle in Fig 4C) is the bin that contains the maximum amount of vertices, thus Bin max provides a median value V med to split the interval [0, 1] into two intervals [0, V med ] and (V med , 1]. By picking two contour lines with maximum gradient magnitudes in intervals [0, V med ] and (V med , 1] respectively, the target tooth and its adjacent tooth can be easily segmented from model.

User interface
An easy-to-use interactive tool for dental model segmentation is indispensable in an interactive segmentation approach. In our method, teeth can effectively separated from dental model by using only a single mouse click without any additional user input. In order to indicate which tooth need to be partitioned, the interactive segmentation tool (as the cross cursor showed in Fig 5) should be placed on the target tooth. Owing to the feature points are automatically detected and grouped in advance, we only need to find the two feature point groups closet to point p, which is selected by user interaction tool and is denoted by yellow sphere in Fig 5), through computing the geodesic distance between point p and feature points (denoted by blue and red spheres in Fig 5 shows). The feature points located in the closest group are assigned as the target constraints, while the feature points belong to the next-closest one are treated as the background constraints. Thus, the target tooth and its adjacent tooth can be segmented from dental model simultaneously with two cutting contour lines. It means that we only need to take only N/2 times (when N is even) or N/2 + 1 times (when N is odd) for our interactive operation to partition all teeth from dental model, if N is the amount of teeth.

Segmentation experiments
Four types of experiments were devised • Experiments on dental models with different levels of crowding problems.
• Experiments on dental meshes with different modeling resolutions and surface noises.
• Comparative experiments conducted using a classical morphologic skeleton approach and our approach.
• Comparative experiments conducted using other common Laplacian weighting scheme and our weighting scheme. Different levels of crowding problems. To testify the effectiveness and accuracy of our method, ten dental models, including five maxilla models (B, C, F, G, I) and five mandible models (A, D, E, H, J), are segmented, which are illustrated by Fig 6. In these ten dental models, case A and G, as well as case B and D are obtained from the Internet, and the rest of dental models are obtained from patients who need medical treatments with the same scan device. According to the crowding levels, these dental models are classified into three types, including "mild crowding", "moderate crowding" and "severe crowding", as Fig 6 shows. Since these dental models are collected from different sources, the scales and qualities of models are different, which can be described in Table 1.
Different levels of modeling resolutions and surface noises. The segmentation field in our method is sensitive to concave regions and is stable in the presence of different levels of tessellations and mesh noises. In Fig 7, we tested three dental models with different mesh resolutions and found that our method is largely sensitive to different tessellations on dental models,  Vertices  151022  113489  451622  226232  399264  511478  151057  537571  572587  334925   Triangles  302118  226843  903240  452460  798524  1022952  302118  1075146  1145170  669852 doi:10.1371/journal.pone.0161159.t001 The segmentation results of employing our approach on three dental models with different surface resolutions. From up to down ((A) to (C)), each dental model is with 301964, 50000 and 10000 faces, respectively. even to coarse tessellation as illustrated in Fig 7C. In addition, our dental model segmentation method is also insensitive to surface noises of dental models, and similar segmentation results are obtained under different scales of mesh noises in Fig 8. Comparison with classical morphologic skeleton segmentation approach. By running source code that is available from [28], a classical morphologic skeleton segmentation approach can be used to compare with our approach. Under the exactly same experimental conditions with same dental models, the comparison of segmentation results are shown in Fig 9. Obviously, the striking failures, including "under-segmentation" and "over-segmentation" in Fig 9A  and 9C, are avoided in our method as show in Fig 9B and 9D. On the other hand, a huge difference in time consumption between segmentation method of [28] and ours are also be detected in the comparison results. The experimental results, as shown in Fig 9A and 9C, cast 42629ms and 159208ms respectively, under framework adopted in [28], while the experimental results, as shown in Fig 9B and 9D, cast 1319ms and 3279ms respectively, under our segmentation framework. Overall, the difference between method of [28] and ours in errors and time consumption are closely related to the curvature threshold value issue. In [28], the curvature threshold value is automatically adjusted, through comparing the amount of segment parts and number of tooth that user entered for equality. Such mechanism is very time-costing due to the Tooth Separation Using Segmentation Field fact that, in order to find a satisfied threshold value, algorithm need to traverse whole the mesh vertex set multiple times with complex computation process. Conversely, corresponding iteration traversal times in our method is limited. The segmentation field computation, which is the most time consuming operation in the procedure of our segmentation framework, only need to traverse mesh vertices once. In addition to this, the appropriate curvature threshold values are difficult to obtain in some cases, the segmentation failures occurred in Fig 9A and 9C are resulted from the inappropriate threshold values. But, the threshold selection can be avoided in our method, such advantage is a significant property to reflect the superiority of our segmentation method, as the comparisons shown.
Comparison with other Laplacian weighting scheme. The other common Laplacian weighting schemes [39,40], are implemented and tested with our segmentation field framework. In Fig 10, we can see that, compare with the classical cotangent weighting scheme [39] ( Fig 10G and 10H),a large color variation can be obviously detected in the concave seams of dental model under the weighting scheme in our method (Fig 10A and 10B) and the weighting scheme adopted by [40] (Fig 10D and 10E). In addition, the concavity-sensitive weighting scheme make the field value variation mainly gathered at narrower concave areas with denser contour lines, as illustrated in Fig 10B, compare to the weighting scheme adopted by [40] (Fig  10E), leading to better segmentation result (Fig 10C and 10F).

Evaluation
The time consumption of our method on the ten dental models in Fig 6 was recorded. In Fig  11, the blue, red, green and purple lines indicate the time consumed by segmentation field computation, cutting line selection, tooth separation and coloring and total time respectively. The  Table 1. Our segmentation method is implemented as a C plus plus plug-in for Meshlab [41], a open source system for the processing and editing of three dimensional triangular meshes, with Intel Core i5-4200H CPU @2.8G Hz and 4GB RAM.

Conclusion
In this article, a novel, effective and efficient dental-target segmentation approach is presented. A segmentation field which is sensitive to concave seams of dental model is computed, through less user interaction to assign constraint sites. The cutting boundaries to separate target tooth from model mesh can be directly extracted from the smooth scalar segmentation field. As shown in the experiment results and comparisons, the properties of our approach is able to present segmentation results with less errors and less time consumption, regardless of different levels of crowding problems and different model qualities. Our future work will involve avoiding user interaction to develop a fully automatic segmentation algorithm, as well as speeding up the operation of our algorithm through optimization algorithm and applying parallel computing technique.

Author Contributions
Conceived and designed the experiments: ZL HW.
Performed the experiments: HW.