# Single-Tooth Modeling  for 3D Dental Model

Authors: Tianran Yuan; Wenhe Liao; Ning Dai; Xiaosheng Cheng; Qing Yu
Journal: International Journal of Biomedical Imaging
 (2010)

---

## Abstract

An integrated single-tooth modeling scheme is proposed for the 3D dental model acquired by optical digitizers. The cores of the modeling scheme are fusion regions extraction, single tooth shape restoration, and single tooth separation. According to the “valley” shape-like characters of the fusion regions between two adjoining teeth, the regions of the 3D dental model are analyzed and classified based on the minimum curvatures of the surface. The single tooth shape is restored according to the bioinformation along the hole boundary, which is generated after the fusion region being removed. By using the extracted boundary from the blending regions between the teeth and soft tissues as reference, the teeth can be separated from the 3D dental model one by one correctly. Experimental results show that the proposed method can achieve satisfying modeling results with high-degree approximation of the real tooth and meet the requirements of clinical oral medicine.


---

## Body

## 1. Introduction

As the 3D (three dimensional) dental model can be acquired easily through different kinds of intra or extra oral measurement methods including optical digitizers [1–3], CT (CBCT) [4, 5] and MRI [6], CAD (Computer-Aided Design)/CAM (Computer-Aided Manufacturing) has been introduced to dentistry and achieved great success in clinical applications [7–10] such as orthodontics, oral and maxillofacial surgery. Dental restorations can be designed and manufactured much more easily compared with traditional complex and labor-intensive process. Pre- or postsurgery simulation can be used to achieve assessment of dental skeletal relationships and facial aesthetics, audit orthodontic outcomes with regard to soft and hard tissues, and direct 3D treatment planning.In general, the 3D dental models (including 3D single-tooth) used in CAD/CAM dentistry system are mostly obtained by optical digitizers, which is typically represented by using a watertight triangular mesh. The 3D dental model is an integral model without obvious blending boundary between the single-tooth and the soft tissues. Two adjoining teeth sometimes are fused together and without obvious tooth gap, due to teeth overlapping, lower measurement precision, and limited resolution triangulating methods during digitizing step. In order to satisfy the prerequisites of manufacturing the dental restorations and assessing the virtual dental behaviors, the teeth have to be independent of each other and keep the original shape of the real tooth. The accurate single-tooth shape restoration and extraction techniques for the 3D dental model play a vital role in CAD/CAM dentistry system.Although the surface of the 3D dental model is extremely irregular and complex, the fusion regions between adjoining teeth and the blending regions between teeth and soft tissues are distributed like “valleys” on the 3D dental model. So, the regions of the 3D dental model can be analyzed quantitatively by applying the corresponding geometric differential component [11, 12]—minimum curvature. The regions identified based on the geometric differential component are the clustering of vertices with similar curvature behavior, which may also include the nontarget regions. In graphics field, the target regions are usually selected by using the window polygon selection mapping method [13], which is difficult to deal with the feature regions of the 3D dental models with complex surface. In this paper, we propose a spatial polygon selection method, of which the edge is straight “line” on the 3D dental model surface.After the fusion regions are selected and removed, the corresponding holes will be generated on the 3D dental model. Researchers have done lots of work on shape restoration for the triangle mesh models. The existing approaches can be classified into two main categories: the nongeometric [14–16] and geometric [17–21]. (1) Non-geometric methods are mainly based on the attributes of the boundary and its n-ring neighbor vertices, in order to reconstruct the field function [14, 15] or implicit surface [16] which can describe the missing part approximately. The corresponding restoration surface patch is generated by using the isosurface extraction method [22]. The restoration result of the nongeometric methods is unique, which cannot achieve restoration with given continuity, and the overall efficiency of these kinds of algorithm is low. (2) In geometric methods, the hole boundary is triangulated based on the mapping plane [20] or spatial triangulation method [18] to get an initial surface patch firstly, and then the initial surface patch is refined and reshaped to obtain the restoration surface patch. The key of the geometric methods is the triangulation of the hole boundary and the following reshaping adjustment.The blending region between two adjacent teeth with obvious tooth gap is similar to a flipped “saddle” shape surface, of which the left and right sides reflect the local shape of the corresponding single-tooth, respectively. Because the surface patch reconstructed by using the existing shape restoration method represents the “whole” instead of the “partial” nature of the model, if the holes formed after the fusion regions being removed are directly filled without being further processed, we will get the restoration results similar to the original model which fails to satisfy the biocharacteristics of the single-tooth (see Figure1(c)). In this paper, we propose a single-tooth shape restoration approach: the hole is firstly divided into two subholes and triangulated separately by using the occlusal plane as the reference; secondly, the triangulation result corresponding to each subhole is subdivided and reshaped as a whole according to the biocharacteristics of the single-tooth.A general scheme of the single-tooth modeling.(a) Zoom view of the local region corresponding to the left first and second molars. (b) Fusion region identified and removed. (c) Restoration result with the hole filled directly (d) Restoration result that satisfies the biocharacteristics of the single-tooth.
(a)(b)(c)(d)After the 3D dental model has been shape restored, the single-tooth can be extracted from the 3D dental model. Various techniques [23–25] have been proposed to segment 3D dental model, which are based on the plane view image information of the 3D dental model. The above methods are limited to segment dental models with mild malocclusion, and missed interstices or wrong cuts will be introduced when dealing with models with severe malocclusion. In this paper, we propose a segmentation boundary extraction method, which is applied directly on the 3D dental model and can separate the single-tooth from the 3D dental model correctly.The single-tooth modeling techniques of the 3D dental model are very important and nontrivial (see Figure1). In this paper, we present an integrated modeling scheme, which mainly includes the following steps.(1)
Digitize the 3D dental model through extra or intra oral measurement method.(2)
Analyze, select, and remove the fusion regions between the adjacent teeth.(3)
Restore the shape of the single-tooth.(4)
Analyze and select the blending region between the teeth and soft tissues.(5)
Extract the segmentation boundary and separate the tooth from the 3D dental model.


## 2. Digital Dental Model Acquisition

Traditional measuring devices used to measure dental casts including dividers, calipers have provided the standard of plaster model analysis [26, 27], but the manual measurement techniques have disadvantages of being time consuming, inaccuracte, and capable of making linear measurements only in a few locations. With advances in computer and optical technology, the dental cast can be digitized through various scanning techniques [1–6]. The 3D dental model can benefit CAD/CAM dentistry in accuracy, efficiency, and ease of measurement of tooth size, arch form, and its dimensions.In this paper, the 3D dental model is scanned from plaster models with a commercially available 3D scanner MCS-30 [28] depending on the structured light technique. A video camera records the structured light distortions after it has been projected over the study models, and then the computer processes the recorded images and merges them together to create a complete 3D dental model. The precision of the 3D scanner MCS-30 with 1280*1024 image resolution can reach 10 μm. The average triangle numbers of the mesh that can meet the clinical precision requirement are usually no less than 20 thousands. The 3D dental model is represented by using a watertight or 2-manifold triangular mesh and usually stored as (Stereo-lithographic) STL or (Virtual Reality Modeling Language) VRML format.


## 3. Feature Regions Analysis and Extraction

### 3.1. Notation

LetMbe the 2-manifold triangular mesh corresponding to surfaceSembedded inR3, V={v1,v2,…,vn} denote the set of vertices in M·n⃗vi represents the unit normal vector of vertex vi. We define NeiV1(i)as 1-ring neighbors of vertexvi, and get n-ring neighbors NeiVn(i) through recursively enlarging the radius of the current neighborhood:(1)NeiV1(i)={i}+{j∣∃edge(vi,vj)},NeiVn(i)=NeiV1(NeiVn-1(i))(n>1)NeiT1(i)is defined as the 1-ring neighboring triangles that share vertexvi·|NeiV1(i)|, |NeiT1(i)|denotes the set size ofNeiV1(i)andNeiT1(i), respectively.


### 3.2. Differential Characteristics Analysis of the 3D Dental Model

Letpbe a point on surfaceS. Consider all curves Ci onSpassing through the pointp. Every such Ci has an associated curvature κigiven atp. Of those curvaturesκi, at least one is characterized as maximal κmaxand one as minimalκmin, and these two curvatures are known as the principal curvatures ofS. In mathematics, the minimum curvature κmin is used to describe the hills (κmin>0) and valleys (κmin<0) of the 3D models, while the maximum curvature is used to describe ridges (κmax>0) and depressions (κmax<0).After a detailed analysis of the 3D dental model’ bioshape characteristics, we find that the blending region between the teeth and soft tissues, the fusion regions between adjoining teeth, and regions including alveolar bone ridges are distributed like “valleys” on the 3D dental model, while the regions corresponding to the incisal edges, cusp tips are “ridges” like. So, the feature regions of the 3D dental model can be classified quantitatively by using the principal curvature information.For the smooth triangular mesh model with uniform triangles, the second-order differential components can be solved by using the corresponding discrete differential geometry operators with guarantee accuracy, which are constructed based on the Laplace-Beltrami operator and spherical mapping methods [11]. But when the triangle shape is irregular and the mesh model is noisy, the calculation results will have much deviation compared with the real value. In this paper, we propose a local surface fitting based method used to estimate the second-order differential properties, which is proved to be robust and accurate in the following experiments. The local shape of any arbitrary complex surface can be described approximately by anm(m≥2)-order polynomial surface:(2)S(u,v)=(u,v,ϕ(u,v)),ϕ(u,v)=∑ak,sf(uk,vs),0≤k+s≤m,k≥0,s≥0,
where ak,sis the polynomial coefficients,f(uk,vs)=ukvs, (2) is the parameter representation of the m(m≥2)-order polynomial surface in the local coordinate system. For vertex viof the mesh model, the corresponding local coordinate systemPuvϕ is determined as follows:Letvibe the origin point of the local coordinate system. ϕaxis coincides with the normal n⃗viof vertexvi.u,vare orthogonal to each other in the tangents plane Tiof vertex vi. Whenϕaxis is paralleled to z-axis of the absolute coordinate systemOxyzafter being applied for by a series of rotation and translation operation, u,vare also paralleled to x,y, respectively.LetKNb(vi)={p1,…,pk}denote the knearest neighbors of viin its n-ring neighborsNeiVn(i). We apply the method proposed by Meyer et al. [11] to calculate the discrete normal vector of the triangular mesh:(3)n⃗vi=14Amixed∑j∈NeiV1(i)n(cotαij+cotβij)(vi-vj),
where αijand βijare two angles  opposite to the edge eijby which viand vjare connected. Amixed is the weighted summation of triangle areas from NeiT1(i) of vertexvi. After n⃗vi is obtained, we can get the mapped vertices KNb(vi)Puvϕ={q1,…,qk} ofKNb(vi)={p1,…,pk} in the local coordinate system Puvϕ,qk=(uk,vk,ϕk). After Puvϕ and the vertices KNb(vi)Puvϕ={q1,…,qk} are determined, the local surfaceS(u,v)fitted to KNb(vi)Puvϕ can be obtained by using the weighted least square method. In this paper, the corresponding least square error is
(4)δ=∑qj∈KNb(vi)Puvϕ|ϕ(uj,vj)-ϕj|2·e-di,j/max(di,j),
where di,j=∥qi-qj∥,1≤i,j≤k. In order to make the local surface be solved with high efficiency, m=2 is applied in this paper. When k= 16~20, the local surface can achieve a better approximation of the real shape. According to the first and second fundamental forms ofS(u,v), the Gaussian curvatureκG and mean curvature κH at u=0,v=0can be solved. Because the differential characters of the mesh model at vertex vi can be  substituted by the differential characters of the local surface S(u,v) at u=0,v=0, the minimum curvature κmin(vi) of vertex vi can be calculated by the following equation:(5)κmin(vi)=kH-kH2-kG.In order to make the solved minimum curvature by (5) reflect the regional characters much more accurately, the curvature values have to be smoothed and denoised further:
(6)κmin(vi)=∑j∈nei1(i)θi,jκmin(vj),
where θi,j=1/∥vi-vj∥.We draw a color map of the minimum curvature values as shown in Figure2 to visualize where the high- and low-curvature areas locate. The highest and lowest curvatures are corresponding to the red and blue color, respectively, the remains are assigned colors between red and green according to the curvature values. As can be seen from Figure 2, the region marked with blue can include the fusion regions and blending regions clearly.Figure 2
Minimum curvature color map of the 3D dental modelWe compare the curvature evaluation method proposed in this paper with that of Meyer et al. [11] by using the torus model:(7)r(u,v)=r(x(u,v),y(u,v),z(u,v)),x(u,v)=(R+rcosv)cosu,y(u,v)=(R+rcosv)sinu,0≤u≤2π,0≤v≤2π,z(u,v)=rsinv,
where R is the wheel radius and r is the tube radius. In this paper, we choose R=2andr=1 as the torus parameters. We obtain the corresponding noisy torus model by adding Gaussian noise with noise level h=0.1,0.3,0.5,0.7,0.9,1.0, respectively. Figure 3 shows that the mean and Gaussian curvatures are robust to noise, and the estimation results (see Figure 4) are much more stable and reliable than those of Meyer et al. [11] when the differential components of noise model are calculated by the method proposed above.(a) Torus model withR=2andr=1. (b) The accurate mean curvature plot of (a). (c) Gaussian noisy model withh=0.5; (d) Mean curvature plot of (c) by Meyer et al. [11]. (e) Mean curvature plot of (c) by the method in this paper.
(a)(b)(c)(d)(e)Average absolute error comparison. (a) Mean curvature. (b) Gaussian curvature.
(a)(b)


### 3.3. Post Processing of the Feature Regions

After the 3D dental model has been analyzed based on the minimum curvature information, the regions of the 3D dental model can be classified and extracted according to the given curvature threshold. However, the feature regions extracted usually contain small pieces and holes (see Figure5(a)). The small pieces contained in the feature regions can be identified efficiently according to the vertex neighboring relationship, and can be removed from the feature regions automatically when the vertex number of the small pieces is less than the given value. We use the mathematical morphology operation extended from the image field to fill the small holes and smooth the feature regions boundaries. There are also four main operators such asdilation,erosion,openandclose included in the 3D mathematical morphology operation [29].(a) Feature regions extracted corresponding to vertices marked with blue in the minimum curvature color map. (b) Feature region before and after being applied for opening and closing operation; (c) Feature regions of (a) after being processed further.
(a)(b)(c)LetF′denote the index set of the vertices in the feature regions. ∀j∈F'⇒vj∈M=(v1,…,vi,…,vn), (1≤j≤n), thedilationanderosion morphology operators corresponding to the 3D models are defined as follows: (8)dilationn(F′)={k∣∃j∈F′:k∈NeiVn(j)},erosionn(F′)={k∣NeiVn(k)∈F′}.Dilation operation is used to “attract” the vertices unmarked as feature vertices but lying inside or at the boundary of the feature regions and can still keep the “shape” of the feature region during dilating. Erosion operation is used to delete undesired branches and will make the feature regions seem much more smooth and thin. We obtain the opening operation by consecutively dilating and then eroding the feature region. The closing operation is obtained by swapping the applying order:(9)openn(F′)=erosionn(F′)∘dilationn(F′),closen(F′)=dilationn(F′)∘erosionn(F′).Multiple application of opening and closing operation can filter out the noise and artifacts of the feature regions effectively. Figure5(b) shows the feature region after being applied for opening and closing operation.


### 3.4. Spatial Polygon Selection Method

As can be seen from Figure5(c), after being further processed, the feature regionsF′can include the fusion regions between the adjacent teeth and the blending regions between the teeth and the soft tissues completely. However, because the feature regions extracted according to the given threshold are the clustering of the vertices with similar curvature behavior, the nontarget regions such as areas including alveolar bone ridges are also extracted. In order to obtain the target regions accurately, the feature regions have to be divided and selected interactively. In this paper, we propose a spatial polygon selection method, of which the edge is straight “line” on the 3D dental model surface.Letpsandpedenote the starting position and the destination on the triangular meshM. The spatial “line” α=(ps,p1⋯pn,pe)between these two vertices on the triangular mesh model can be solved approximately by using the direction tracing method described as follows:Letni⃗, ni+1⃗be the normal ofpi,pi+1. AssumeQ with unit length is the object being moved on the surface of the triangular mesh. If pi is the current position, andpi+1 is the next position Qgoing along the direction pipe⃗ as shown in Figure 6(a), we do not change the direction pipe⃗and make sure that the pose of Q parallels to the normal ni⃗ when Qis moving in the interior of a triangle or along an edge until it get topi+1. When Qis atpi+1, we change the moving direction pipe⃗intopi+1pe⃗, and the pose ni⃗ into ni+1⃗. Because (p0,…,pi-1,pi,pi+1,…,pn) is piecewise linear continuous and lies over the triangular mesh, piandpi+1must belong to the same triangle. The line segmentpipi+1⃗ is also the intersection between the triangle which includes pi,pi+1and the normal section πat pithrough pipe⃗. We can obtain pi+1 according to pi, pe and ni⃗as shown in Figure 6(a). The incident triangles of pi intersecting with the normal section πmay be more than one sometime. Let pi+11, pi+12⋯ denote the intersection points. We choose pi+1* as pi+1 when the angle between pipi+1*⃗ and pipe⃗ is the smallest of all thepi+1*∈(pi+11,pi+12⋯). Beginning with ps, p1,p2⋯pi⋯is solved in turn. When pi and pe are in the same triangle, we get the 3D “line”α=(ps,p1⋯pn,pe)(see Figure 6).(a) Concept of the direction tracing method (b) Examples of “line” solved on the planar and spatial surface.
(a)(b)As shown in Figure7(a), the edges of the spatial polygon can be determined by a series of vertices on the 3D model, which are selected interactively according to the profile of the target region. The triangles are marked with “select”, of which the three vertices fall into the inner of the spatial polygon together.(a) Fusion region extraction example by using the spatial polygon selection method; (b) Results of the fusion regions being extracted and removed.
(a)(b)


## 3.1. Notation

LetMbe the 2-manifold triangular mesh corresponding to surfaceSembedded inR3, V={v1,v2,…,vn} denote the set of vertices in M·n⃗vi represents the unit normal vector of vertex vi. We define NeiV1(i)as 1-ring neighbors of vertexvi, and get n-ring neighbors NeiVn(i) through recursively enlarging the radius of the current neighborhood:(1)NeiV1(i)={i}+{j∣∃edge(vi,vj)},NeiVn(i)=NeiV1(NeiVn-1(i))(n>1)NeiT1(i)is defined as the 1-ring neighboring triangles that share vertexvi·|NeiV1(i)|, |NeiT1(i)|denotes the set size ofNeiV1(i)andNeiT1(i), respectively.


## 3.2. Differential Characteristics Analysis of the 3D Dental Model

Letpbe a point on surfaceS. Consider all curves Ci onSpassing through the pointp. Every such Ci has an associated curvature κigiven atp. Of those curvaturesκi, at least one is characterized as maximal κmaxand one as minimalκmin, and these two curvatures are known as the principal curvatures ofS. In mathematics, the minimum curvature κmin is used to describe the hills (κmin>0) and valleys (κmin<0) of the 3D models, while the maximum curvature is used to describe ridges (κmax>0) and depressions (κmax<0).After a detailed analysis of the 3D dental model’ bioshape characteristics, we find that the blending region between the teeth and soft tissues, the fusion regions between adjoining teeth, and regions including alveolar bone ridges are distributed like “valleys” on the 3D dental model, while the regions corresponding to the incisal edges, cusp tips are “ridges” like. So, the feature regions of the 3D dental model can be classified quantitatively by using the principal curvature information.For the smooth triangular mesh model with uniform triangles, the second-order differential components can be solved by using the corresponding discrete differential geometry operators with guarantee accuracy, which are constructed based on the Laplace-Beltrami operator and spherical mapping methods [11]. But when the triangle shape is irregular and the mesh model is noisy, the calculation results will have much deviation compared with the real value. In this paper, we propose a local surface fitting based method used to estimate the second-order differential properties, which is proved to be robust and accurate in the following experiments. The local shape of any arbitrary complex surface can be described approximately by anm(m≥2)-order polynomial surface:(2)S(u,v)=(u,v,ϕ(u,v)),ϕ(u,v)=∑ak,sf(uk,vs),0≤k+s≤m,k≥0,s≥0,
where ak,sis the polynomial coefficients,f(uk,vs)=ukvs, (2) is the parameter representation of the m(m≥2)-order polynomial surface in the local coordinate system. For vertex viof the mesh model, the corresponding local coordinate systemPuvϕ is determined as follows:Letvibe the origin point of the local coordinate system. ϕaxis coincides with the normal n⃗viof vertexvi.u,vare orthogonal to each other in the tangents plane Tiof vertex vi. Whenϕaxis is paralleled to z-axis of the absolute coordinate systemOxyzafter being applied for by a series of rotation and translation operation, u,vare also paralleled to x,y, respectively.LetKNb(vi)={p1,…,pk}denote the knearest neighbors of viin its n-ring neighborsNeiVn(i). We apply the method proposed by Meyer et al. [11] to calculate the discrete normal vector of the triangular mesh:(3)n⃗vi=14Amixed∑j∈NeiV1(i)n(cotαij+cotβij)(vi-vj),
where αijand βijare two angles  opposite to the edge eijby which viand vjare connected. Amixed is the weighted summation of triangle areas from NeiT1(i) of vertexvi. After n⃗vi is obtained, we can get the mapped vertices KNb(vi)Puvϕ={q1,…,qk} ofKNb(vi)={p1,…,pk} in the local coordinate system Puvϕ,qk=(uk,vk,ϕk). After Puvϕ and the vertices KNb(vi)Puvϕ={q1,…,qk} are determined, the local surfaceS(u,v)fitted to KNb(vi)Puvϕ can be obtained by using the weighted least square method. In this paper, the corresponding least square error is
(4)δ=∑qj∈KNb(vi)Puvϕ|ϕ(uj,vj)-ϕj|2·e-di,j/max(di,j),
where di,j=∥qi-qj∥,1≤i,j≤k. In order to make the local surface be solved with high efficiency, m=2 is applied in this paper. When k= 16~20, the local surface can achieve a better approximation of the real shape. According to the first and second fundamental forms ofS(u,v), the Gaussian curvatureκG and mean curvature κH at u=0,v=0can be solved. Because the differential characters of the mesh model at vertex vi can be  substituted by the differential characters of the local surface S(u,v) at u=0,v=0, the minimum curvature κmin(vi) of vertex vi can be calculated by the following equation:(5)κmin(vi)=kH-kH2-kG.In order to make the solved minimum curvature by (5) reflect the regional characters much more accurately, the curvature values have to be smoothed and denoised further:
(6)κmin(vi)=∑j∈nei1(i)θi,jκmin(vj),
where θi,j=1/∥vi-vj∥.We draw a color map of the minimum curvature values as shown in Figure2 to visualize where the high- and low-curvature areas locate. The highest and lowest curvatures are corresponding to the red and blue color, respectively, the remains are assigned colors between red and green according to the curvature values. As can be seen from Figure 2, the region marked with blue can include the fusion regions and blending regions clearly.Figure 2
Minimum curvature color map of the 3D dental modelWe compare the curvature evaluation method proposed in this paper with that of Meyer et al. [11] by using the torus model:(7)r(u,v)=r(x(u,v),y(u,v),z(u,v)),x(u,v)=(R+rcosv)cosu,y(u,v)=(R+rcosv)sinu,0≤u≤2π,0≤v≤2π,z(u,v)=rsinv,
where R is the wheel radius and r is the tube radius. In this paper, we choose R=2andr=1 as the torus parameters. We obtain the corresponding noisy torus model by adding Gaussian noise with noise level h=0.1,0.3,0.5,0.7,0.9,1.0, respectively. Figure 3 shows that the mean and Gaussian curvatures are robust to noise, and the estimation results (see Figure 4) are much more stable and reliable than those of Meyer et al. [11] when the differential components of noise model are calculated by the method proposed above.(a) Torus model withR=2andr=1. (b) The accurate mean curvature plot of (a). (c) Gaussian noisy model withh=0.5; (d) Mean curvature plot of (c) by Meyer et al. [11]. (e) Mean curvature plot of (c) by the method in this paper.
(a)(b)(c)(d)(e)Average absolute error comparison. (a) Mean curvature. (b) Gaussian curvature.
(a)(b)


## 3.3. Post Processing of the Feature Regions

After the 3D dental model has been analyzed based on the minimum curvature information, the regions of the 3D dental model can be classified and extracted according to the given curvature threshold. However, the feature regions extracted usually contain small pieces and holes (see Figure5(a)). The small pieces contained in the feature regions can be identified efficiently according to the vertex neighboring relationship, and can be removed from the feature regions automatically when the vertex number of the small pieces is less than the given value. We use the mathematical morphology operation extended from the image field to fill the small holes and smooth the feature regions boundaries. There are also four main operators such asdilation,erosion,openandclose included in the 3D mathematical morphology operation [29].(a) Feature regions extracted corresponding to vertices marked with blue in the minimum curvature color map. (b) Feature region before and after being applied for opening and closing operation; (c) Feature regions of (a) after being processed further.
(a)(b)(c)LetF′denote the index set of the vertices in the feature regions. ∀j∈F'⇒vj∈M=(v1,…,vi,…,vn), (1≤j≤n), thedilationanderosion morphology operators corresponding to the 3D models are defined as follows: (8)dilationn(F′)={k∣∃j∈F′:k∈NeiVn(j)},erosionn(F′)={k∣NeiVn(k)∈F′}.Dilation operation is used to “attract” the vertices unmarked as feature vertices but lying inside or at the boundary of the feature regions and can still keep the “shape” of the feature region during dilating. Erosion operation is used to delete undesired branches and will make the feature regions seem much more smooth and thin. We obtain the opening operation by consecutively dilating and then eroding the feature region. The closing operation is obtained by swapping the applying order:(9)openn(F′)=erosionn(F′)∘dilationn(F′),closen(F′)=dilationn(F′)∘erosionn(F′).Multiple application of opening and closing operation can filter out the noise and artifacts of the feature regions effectively. Figure5(b) shows the feature region after being applied for opening and closing operation.


## 3.4. Spatial Polygon Selection Method

As can be seen from Figure5(c), after being further processed, the feature regionsF′can include the fusion regions between the adjacent teeth and the blending regions between the teeth and the soft tissues completely. However, because the feature regions extracted according to the given threshold are the clustering of the vertices with similar curvature behavior, the nontarget regions such as areas including alveolar bone ridges are also extracted. In order to obtain the target regions accurately, the feature regions have to be divided and selected interactively. In this paper, we propose a spatial polygon selection method, of which the edge is straight “line” on the 3D dental model surface.Letpsandpedenote the starting position and the destination on the triangular meshM. The spatial “line” α=(ps,p1⋯pn,pe)between these two vertices on the triangular mesh model can be solved approximately by using the direction tracing method described as follows:Letni⃗, ni+1⃗be the normal ofpi,pi+1. AssumeQ with unit length is the object being moved on the surface of the triangular mesh. If pi is the current position, andpi+1 is the next position Qgoing along the direction pipe⃗ as shown in Figure 6(a), we do not change the direction pipe⃗and make sure that the pose of Q parallels to the normal ni⃗ when Qis moving in the interior of a triangle or along an edge until it get topi+1. When Qis atpi+1, we change the moving direction pipe⃗intopi+1pe⃗, and the pose ni⃗ into ni+1⃗. Because (p0,…,pi-1,pi,pi+1,…,pn) is piecewise linear continuous and lies over the triangular mesh, piandpi+1must belong to the same triangle. The line segmentpipi+1⃗ is also the intersection between the triangle which includes pi,pi+1and the normal section πat pithrough pipe⃗. We can obtain pi+1 according to pi, pe and ni⃗as shown in Figure 6(a). The incident triangles of pi intersecting with the normal section πmay be more than one sometime. Let pi+11, pi+12⋯ denote the intersection points. We choose pi+1* as pi+1 when the angle between pipi+1*⃗ and pipe⃗ is the smallest of all thepi+1*∈(pi+11,pi+12⋯). Beginning with ps, p1,p2⋯pi⋯is solved in turn. When pi and pe are in the same triangle, we get the 3D “line”α=(ps,p1⋯pn,pe)(see Figure 6).(a) Concept of the direction tracing method (b) Examples of “line” solved on the planar and spatial surface.
(a)(b)As shown in Figure7(a), the edges of the spatial polygon can be determined by a series of vertices on the 3D model, which are selected interactively according to the profile of the target region. The triangles are marked with “select”, of which the three vertices fall into the inner of the spatial polygon together.(a) Fusion region extraction example by using the spatial polygon selection method; (b) Results of the fusion regions being extracted and removed.
(a)(b)


## 4. Single Tooth Shape Restoration

After the fusion regions have been removed, the corresponding holes are generated on the 3D dental model. The holes are typical “saddle shape” and each one is shared by two adjoining single-teeth (see Figure1(b)). If the holes are directly filled without being further processed, we will get the restored model similar to the original (see Figure 1(c)). The failure reason is that the hole belongs to two teeth which are adjoining but independent from each other. If the hole is filled as a whole, the boundary information of the two independent teeth will be diffused into the same restoration surface patch averagely, and cannot reflect the biocharacteristics of the single-tooth. In this paper, we present a single-tooth shape restoration approach. In order to achieve a best approximation of the original tooth, the restoration process has to satisfy the following prerequisites.(a)
The restoration surface patch corresponding to the missing part should be reconstructed in a way that is minimally distinguishable from the surrounding regions, and should also preserve the sampling density of the original 3D dental model.(b)
There is no interference between the restored teeth according to the independence properties of the teeth. The blending region between the adjacent teeth has to be natural and continuous.LetBdenote the hole boundary of the 3D dental model. P represents the filling patch forB. Pmin,PrefineandPdeformare the surface patches corresponding to different filling stages: spanning triangulation, refinement, deformation.Pfinalrepresents the final restoration result, which meets the restoration prerequisites.


### 4.1. Hole Boundary Triangulation

In geometric methods, the hole boundary is mostly triangulated based on the mapping plane [20] or spatial triangulation method [18] to get an initial surface patch. The mapping plane triangulation methods convert the 3D hole boundary into the 2D polygon by projecting it onto the mapping plane, which is fitted to the boundary vertices by the least square method. The mapping plane triangulation methods can achieve satisfying results in dealing with the simple regular hole, which is homeomorphic to a disc after projection. But for the complex hole with sharp curvature changes along the boundary, there will appear self-intersection in the projected 2D polygon. Barequet and Sharir [18] give an interesting solution of the 3D polygon triangulating problem. The spatial triangulation method [18] has order of O(N3)time complexity (Nis number of the boundary vertices), which is adaptable to deal with the boundary with small vertices number, but difficult to triangulate the big hole. In this paper, we proposed a spatial triangulation method based on local optimized weight rule, in which various influencing factors that may affect the triangulation are considered completely.We defineΩ:B3→Las the weight function, where B={v1b,v2b,…,vnb},Lis the weight set and Ωassigns a weight for each triangle with three consecutive vertices ofB. Let A(vi)denote the sum of the adjacent angles of the current vertexvi(10)A(vi)=∑j∈NeiT1(i)Aj.
During the triangulation process, when0<A(vib)<απ, after the new triangle (vi-1b,vib,vi+1b)is added as shown in Figure 8, sharp corners or triangle with interior angle close to πwill be formed. In order to avoid such situations to appear, the candidate triangle(vi-1b,vib,vi+1b) should be assigned a weight llesswith lower choice priority when 0<A(vib)<απ. we found empirically that α=1.2can yield good results.Abnormal results of the triangle being added whenA(vib)<απ. (a) Form sharp corner; (b) Generate triangle with interior angle close to π.
(a)(b)In the 2-manifold triangular mesh model, the better number of neighboring triangles is usually 5 to 8. In order to avoid too many new generated triangles converge at the same boundary vertex, the number of the vertex 1-ring neighboring triangles has to be limited during triangulation. So, when|NeiT1(i)|>8, the vertexvib should be removed firstly, which means that the candidate triangle (vi-1b,vib,vi+1b) should be assigned a weightlbigger with higher choice priority. At the same time, when the current vertex’ 1-ring neighboring triangles are projected onto their own tangent plane, there should be no intersection between the projected edges except at the current vertex itself. So, the new added triangle (vi-1b,vib,vi+1b) has to satisfy the nonintersection projection condition at vi-1b,vib,vi+1b simultaneously.When|NeiT1(i)|≤8,A(vib)>απ, the weight of the candidate triangle (vi-1b,vib,vi+1b) should be determined by its own geometric attributes such as edge length, area, and interior angle. In order to obtain a triangulation surface patch with moderate internal changes, the edges should be distributed along the boundary averagely similar to a curtain covering at a window, and the vertices of the edges should be the pairs relative nearest to each other in space. So, the candidate triangle (vi-1b,vib,vi+1b) should be weighted according to its corresponding edge length. The smaller the perimeter of the triangle is, the higher choice priority it will have.Based on the analysis of the influencing factor which will affect the triangulation results, weight functionsΩ,lless, and lbigger are described as follows:(11)Ω(vjb,vib,vkb)={-(∥eij∥+∥eik∥+∥ejk∥)when|NeiT1(i)|≤8,A(vib)>απ,llesswhen0<A(vib)<απ,|NeiT1(i)|≤8,lbiggerwhen|NeiT1(i)|>8,-∞whenintersectionafterprojection,
wherelless(vjb,vib,vkb)=-(π/Asum)*RC, lbigger(vjb,vib,vkb)=(|NeiT1(i)|/8)*RC, andRCis the radius of the model’ bounding sphere. We apply the following procedure to implement the triangulation process:(12)vn+1b=v1b,vn+2b=v2bStep 1.
Compute all the weightsΩ(vi-1b,vib,vi+1b) according to the weight function given above for each triangle(vi-1b,vib,vi+1b) with three consecutive vertices ofB, and insert the weights intoL in which the weight is sorted using an AVL tree.Step 2.
Select the maximumlmax from the weight setL, and insert its corresponding triangle (vi-1b,vib,vi+1b) intoMC. Remove the weights of the trianglesΩ(vi-2b,vi-1b,vib), Ω(vi-1b,vib,vi+1b) and Ω(vib,vi+1b,vi+2b) fromL that include vertexvib. Eliminate vertex vibfromB, and then BH={v1b,v2b,…,vi-2b,vi-1b,vi+1b,vi+2b,…,vnb}. Compute the weightsΩ(vi-2b,vi-1b,vi+1b),Ω(vi-2b,vi-1b,vi+1b) of the triangle(vi-2b,vi-1b,vi+1b),(vi-1b,vi+1b,vi+2b), and insert them into L.Step 3.
Execute Step2 iteratively until the vertex number of Bis less than three, and obtain the initial surface patchPmin.


### 4.2. Subhole Division

In order to reconstruct the surface patch with the shape of the flip “saddle”, the hole boundaryBhas to be divided into two subholesB1, B2 as shown in Figure 11(a). Each subhole is corresponding to its own tooth. The end points of the bridge edge, by which the hole boundary is bridged to form two separate subholes, are two points farthest to the occlusal plane on the buccal and lingual side of the hole boundary respectively, and can be selected automatically by using the occlusal plane as the reference. In this paper, the occlusal plane is fitted with four reference points (including the buccal cusp tips of the left and right first molars, and the mesiobuccal points of the left and right first permanent molars) as shown in Figure 9.(a) B1, B2, B3, and B4 are the corresponding reference points. (b) The occlusal plane fitted to the four reference points.
(a)(b)The subholeB1(B2) is first filled in with a local optimized triangulation P1min(P2min) of its 3D contour (see Figure 11(b)). The initial subfilling surface patchesP1min, P2min combine a complete initial filling surface patchPmin for Btogether (see Figure 11(c)).


### 4.3. Refinement

Because the edges in the initial surface patchPmin are the direct connections between boundary vertices, the surface patch Pmin has to be refined according to the boundary information to obtain a further surface patchPrefine, which approximates the density of the surrounding mesh. The mesh density is usually measured based on the average length of the edges. In this paper, the bigger triangle (vi,vj,vk) is split into three smaller ones by using “1-3” face splitting method, in which the new added vertex is the centroid vc=(vi+vj+vk)/3of the triangle(vi,vj,vk), and the interior edges are relaxed while splitting to maintain a Delaunay-like triangulation (see Figure 10).(a) Refinement and optimization mechanism. (b) Surface patch before and after being refined.
(a)(b)(a) Hole bridged to divide into two separate subholes. (b) Triangulation of the subhole. The corresponding shaded and wireframe surface patches. (c) After completely filling of the entire hole. (d) After refinement. (e) After reshaping adjustment.
(a)(b)(c)(d)(e)


### 4.4. Reshaping

Prefineis still a surface patch withC0continuity both at the boundary and the internal. The surface patchPrefinehas to be reshaped in order to generate a surface patch, which can both reflect the local characteristics of the missing part and have a good degree of visual reality. In this paper, we design a reshaping adjustment scheme based on the discrete Euler-Lagrange equation. The reshaping adjustment scheme is described as follows.LetS:Ω→R3be the smooth surface corresponding toM.S*denotes the k-order partial derivatives, andδΩstands for the surface boundary. The quadratic energy function [30] for the surface is  (13)Ek(S)=∫Fk(Su⋯u,Su⋯uv,…,Sv⋯v).
In order to actually compute the solution to the above optimization problem, we apply variational calculus to derive the corresponding Euler-Lagrange equation which characterizes the minimizers of(13)(14)ΔkS(x)=0,x∈Ω∖δΩΔjS(x)=bj(x),x∈δΩ,j<k,
where Δis the Laplace Operator, bj(j<k) is the boundary constraints. In order to ensure the efficiency and stability of algorithms, the value range of kare limited to1≤k≤3. When we use a triangle mesh as the underlying surface representation, the Laplace operator is discretized as (15)Δ(vi)=2Area(vi)∑j∈NeiV1(i)(cotαij+cotβij)(vi-vj),
whereArea(vi)is the area sum of the vertex’ 1-ring neighboring triangles, and αijandβijare two angles opposite to the edgeeij. The higher order Laplace operator can be solved iteratively  (16)Δ¯k(v)=Δ(Δ¯k-1(v)).
And then, (14) becomes a linear equation with sparse matrix  (17)[Δ¯k0∣IF](PF)=(0F),
where P=(v1,…,vp) are the free vertices interior of the surface patch. F=(f1,…,fF) are the constraint vertices with Ck-1boundary continuity. For k=1,k=2, and k=3, the surface solved from (17) is corresponding to a membrane with minimization surface area, a thin plate with minimization bending, and a surface with minimization curvature variation, respectively.According to the geometric characteristics of the tooth surface, the surface patch obtained after being reshaped should be a surface with minimum bending variation. So, the constraint parameterkis assigned the value 2 in this paper. During the deformation stage, the triangles that have greater shape change should be refined again. If we apply Pdeformas the final result Pfinaldirectly, small interference will appear between the adjoining teeth sometimes (see Figure 12). We use the following equation to control the deformation degree:(18)Pfinal=Prefine+λ(Pdeform-Prefine),0<λ≤1.The values ofλ determines the deformation degree of the restoration patch. From left to right: the refined mesh patch deformed withλ=0.7,λ=0.85, andλ=1.0.
(a)(b)(c)Figure12 shows the restoration results with λassigned different values. The value of λ determines the deformation degree of the restoration patch.


## 4.1. Hole Boundary Triangulation

In geometric methods, the hole boundary is mostly triangulated based on the mapping plane [20] or spatial triangulation method [18] to get an initial surface patch. The mapping plane triangulation methods convert the 3D hole boundary into the 2D polygon by projecting it onto the mapping plane, which is fitted to the boundary vertices by the least square method. The mapping plane triangulation methods can achieve satisfying results in dealing with the simple regular hole, which is homeomorphic to a disc after projection. But for the complex hole with sharp curvature changes along the boundary, there will appear self-intersection in the projected 2D polygon. Barequet and Sharir [18] give an interesting solution of the 3D polygon triangulating problem. The spatial triangulation method [18] has order of O(N3)time complexity (Nis number of the boundary vertices), which is adaptable to deal with the boundary with small vertices number, but difficult to triangulate the big hole. In this paper, we proposed a spatial triangulation method based on local optimized weight rule, in which various influencing factors that may affect the triangulation are considered completely.We defineΩ:B3→Las the weight function, where B={v1b,v2b,…,vnb},Lis the weight set and Ωassigns a weight for each triangle with three consecutive vertices ofB. Let A(vi)denote the sum of the adjacent angles of the current vertexvi(10)A(vi)=∑j∈NeiT1(i)Aj.
During the triangulation process, when0<A(vib)<απ, after the new triangle (vi-1b,vib,vi+1b)is added as shown in Figure 8, sharp corners or triangle with interior angle close to πwill be formed. In order to avoid such situations to appear, the candidate triangle(vi-1b,vib,vi+1b) should be assigned a weight llesswith lower choice priority when 0<A(vib)<απ. we found empirically that α=1.2can yield good results.Abnormal results of the triangle being added whenA(vib)<απ. (a) Form sharp corner; (b) Generate triangle with interior angle close to π.
(a)(b)In the 2-manifold triangular mesh model, the better number of neighboring triangles is usually 5 to 8. In order to avoid too many new generated triangles converge at the same boundary vertex, the number of the vertex 1-ring neighboring triangles has to be limited during triangulation. So, when|NeiT1(i)|>8, the vertexvib should be removed firstly, which means that the candidate triangle (vi-1b,vib,vi+1b) should be assigned a weightlbigger with higher choice priority. At the same time, when the current vertex’ 1-ring neighboring triangles are projected onto their own tangent plane, there should be no intersection between the projected edges except at the current vertex itself. So, the new added triangle (vi-1b,vib,vi+1b) has to satisfy the nonintersection projection condition at vi-1b,vib,vi+1b simultaneously.When|NeiT1(i)|≤8,A(vib)>απ, the weight of the candidate triangle (vi-1b,vib,vi+1b) should be determined by its own geometric attributes such as edge length, area, and interior angle. In order to obtain a triangulation surface patch with moderate internal changes, the edges should be distributed along the boundary averagely similar to a curtain covering at a window, and the vertices of the edges should be the pairs relative nearest to each other in space. So, the candidate triangle (vi-1b,vib,vi+1b) should be weighted according to its corresponding edge length. The smaller the perimeter of the triangle is, the higher choice priority it will have.Based on the analysis of the influencing factor which will affect the triangulation results, weight functionsΩ,lless, and lbigger are described as follows:(11)Ω(vjb,vib,vkb)={-(∥eij∥+∥eik∥+∥ejk∥)when|NeiT1(i)|≤8,A(vib)>απ,llesswhen0<A(vib)<απ,|NeiT1(i)|≤8,lbiggerwhen|NeiT1(i)|>8,-∞whenintersectionafterprojection,
wherelless(vjb,vib,vkb)=-(π/Asum)*RC, lbigger(vjb,vib,vkb)=(|NeiT1(i)|/8)*RC, andRCis the radius of the model’ bounding sphere. We apply the following procedure to implement the triangulation process:(12)vn+1b=v1b,vn+2b=v2bStep 1.
Compute all the weightsΩ(vi-1b,vib,vi+1b) according to the weight function given above for each triangle(vi-1b,vib,vi+1b) with three consecutive vertices ofB, and insert the weights intoL in which the weight is sorted using an AVL tree.Step 2.
Select the maximumlmax from the weight setL, and insert its corresponding triangle (vi-1b,vib,vi+1b) intoMC. Remove the weights of the trianglesΩ(vi-2b,vi-1b,vib), Ω(vi-1b,vib,vi+1b) and Ω(vib,vi+1b,vi+2b) fromL that include vertexvib. Eliminate vertex vibfromB, and then BH={v1b,v2b,…,vi-2b,vi-1b,vi+1b,vi+2b,…,vnb}. Compute the weightsΩ(vi-2b,vi-1b,vi+1b),Ω(vi-2b,vi-1b,vi+1b) of the triangle(vi-2b,vi-1b,vi+1b),(vi-1b,vi+1b,vi+2b), and insert them into L.Step 3.
Execute Step2 iteratively until the vertex number of Bis less than three, and obtain the initial surface patchPmin.


## 4.2. Subhole Division

In order to reconstruct the surface patch with the shape of the flip “saddle”, the hole boundaryBhas to be divided into two subholesB1, B2 as shown in Figure 11(a). Each subhole is corresponding to its own tooth. The end points of the bridge edge, by which the hole boundary is bridged to form two separate subholes, are two points farthest to the occlusal plane on the buccal and lingual side of the hole boundary respectively, and can be selected automatically by using the occlusal plane as the reference. In this paper, the occlusal plane is fitted with four reference points (including the buccal cusp tips of the left and right first molars, and the mesiobuccal points of the left and right first permanent molars) as shown in Figure 9.(a) B1, B2, B3, and B4 are the corresponding reference points. (b) The occlusal plane fitted to the four reference points.
(a)(b)The subholeB1(B2) is first filled in with a local optimized triangulation P1min(P2min) of its 3D contour (see Figure 11(b)). The initial subfilling surface patchesP1min, P2min combine a complete initial filling surface patchPmin for Btogether (see Figure 11(c)).


## 4.3. Refinement

Because the edges in the initial surface patchPmin are the direct connections between boundary vertices, the surface patch Pmin has to be refined according to the boundary information to obtain a further surface patchPrefine, which approximates the density of the surrounding mesh. The mesh density is usually measured based on the average length of the edges. In this paper, the bigger triangle (vi,vj,vk) is split into three smaller ones by using “1-3” face splitting method, in which the new added vertex is the centroid vc=(vi+vj+vk)/3of the triangle(vi,vj,vk), and the interior edges are relaxed while splitting to maintain a Delaunay-like triangulation (see Figure 10).(a) Refinement and optimization mechanism. (b) Surface patch before and after being refined.
(a)(b)(a) Hole bridged to divide into two separate subholes. (b) Triangulation of the subhole. The corresponding shaded and wireframe surface patches. (c) After completely filling of the entire hole. (d) After refinement. (e) After reshaping adjustment.
(a)(b)(c)(d)(e)


## 4.4. Reshaping

Prefineis still a surface patch withC0continuity both at the boundary and the internal. The surface patchPrefinehas to be reshaped in order to generate a surface patch, which can both reflect the local characteristics of the missing part and have a good degree of visual reality. In this paper, we design a reshaping adjustment scheme based on the discrete Euler-Lagrange equation. The reshaping adjustment scheme is described as follows.LetS:Ω→R3be the smooth surface corresponding toM.S*denotes the k-order partial derivatives, andδΩstands for the surface boundary. The quadratic energy function [30] for the surface is  (13)Ek(S)=∫Fk(Su⋯u,Su⋯uv,…,Sv⋯v).
In order to actually compute the solution to the above optimization problem, we apply variational calculus to derive the corresponding Euler-Lagrange equation which characterizes the minimizers of(13)(14)ΔkS(x)=0,x∈Ω∖δΩΔjS(x)=bj(x),x∈δΩ,j<k,
where Δis the Laplace Operator, bj(j<k) is the boundary constraints. In order to ensure the efficiency and stability of algorithms, the value range of kare limited to1≤k≤3. When we use a triangle mesh as the underlying surface representation, the Laplace operator is discretized as (15)Δ(vi)=2Area(vi)∑j∈NeiV1(i)(cotαij+cotβij)(vi-vj),
whereArea(vi)is the area sum of the vertex’ 1-ring neighboring triangles, and αijandβijare two angles opposite to the edgeeij. The higher order Laplace operator can be solved iteratively  (16)Δ¯k(v)=Δ(Δ¯k-1(v)).
And then, (14) becomes a linear equation with sparse matrix  (17)[Δ¯k0∣IF](PF)=(0F),
where P=(v1,…,vp) are the free vertices interior of the surface patch. F=(f1,…,fF) are the constraint vertices with Ck-1boundary continuity. For k=1,k=2, and k=3, the surface solved from (17) is corresponding to a membrane with minimization surface area, a thin plate with minimization bending, and a surface with minimization curvature variation, respectively.According to the geometric characteristics of the tooth surface, the surface patch obtained after being reshaped should be a surface with minimum bending variation. So, the constraint parameterkis assigned the value 2 in this paper. During the deformation stage, the triangles that have greater shape change should be refined again. If we apply Pdeformas the final result Pfinaldirectly, small interference will appear between the adjoining teeth sometimes (see Figure 12). We use the following equation to control the deformation degree:(18)Pfinal=Prefine+λ(Pdeform-Prefine),0<λ≤1.The values ofλ determines the deformation degree of the restoration patch. From left to right: the refined mesh patch deformed withλ=0.7,λ=0.85, andλ=1.0.
(a)(b)(c)Figure12 shows the restoration results with λassigned different values. The value of λ determines the deformation degree of the restoration patch.


## 5. Single Tooth Extraction

After the 3D dental model has been shape restored, the single-tooth can be extracted from the 3D dental model. The differential information of the 3D dental model is re-analyzed and processed by using the methods proposed above. As can be seen from Figure14(a), the feature regions identified based on the minimum curvature value can include the blending regions completely, and has already possessed the coarse profile of the segmentation boundary. The feature regions are still too coarse to be accepted as the segmentation boundary. We have to peel the vertices of the region boundary inward until obtaining its skeleton with width of one vertex. The key step of the boundary extraction is how to judge a vertex should be peeled or not, and the skeleton must follow the original topology of the feature region. In this paper, the segmentation boundary extraction procedure is designed based on the vertex complexity [29]∀i∈F′; letCnei(i)denote the 1-ring neighborhood of vertex vi ordered counter clockwise. ∀k∈Cnei(i); if k∈F′at the same time, we recordCnei(i)k=1 orCnei(i)k=0. With the above assumption, the vertex complexity CP(i) of vi is defined as follows (see Figure 13): (19)CP(i)=∑k=1k≤m|Cnei(i)k-Cnei(i)k+1|.Figure 13
Examples of vertex’ complexity computation.Procedure of segmentation boundary extraction and single-tooth separation. (a) 3D dental model with extracted feature regions. (b) after peeling. (c) after pruning. (d) single-tooth separation.
(a)(b)(c)(d)IfCP(i)≥4, vertex viis defined to be complex. IfCnei(i)⊆F′, vertex viis defined as center vertex. The neighbor of the center vertex is called satellite vertex, when its corresponding complexity is no less than zero. During the boundary extracting, the center vertex and complex vertex are marked as feature vertices that should be preserved. If the center vertex is removed, small close ring will be formed in the inner of the feature region, and the regional connectivity will be undermined if the complex vertex is removed. The set of satellite vertices is denoted byFS, center vertices byFC, and complex vertices byFCP. Then, we obtain the set of candidate vertices FDthat will be removed as follows.(20)FD=FS∩FCP∪FC¯.We remove one vertex from the candidate setFD each time, and recalculate its neighboring vertex’ complexity simultaneously. The set FS,FC,FCP and FD are updated after each removing. The removing and updating operation is iterated until the “shape” of the feature regions does not change anymore.The skeleton obtained after being applied for the above operation also contains unnecessary open branches as shown in Figure14(b). Because the segmentation boundary used to extract the single-tooth is a set of closed rings. The branches can be identified and pruned by deleting the line segment from the skeleton iteratively, which has at least one endpoints only connected with itself. Sometimes, there will be small redundant close rings existing, which is need to be removed interactively. After pruning, we obtain the segmentation boundary as shown in Figure 14(c). Figure 14(d) shows the single-tooth extracted according to the segmentation boundary.


## 6. Experimental Results and Analysis

In order to verify the validity and adaptability of the proposed method, we have conducted a series of experiments on various types of 3D models. Figures16 and 17 show the modeling results of the two typical kinds of 3D dental models (see Table 3) including model with normal tooth arrangement and model with severe malocclusion. Table 1 shows the detailed information of the 3D dental models including bounding box size, vertex/triangle numbers before and after shape modeling.Table 1
The detailed information of the 3D dental models.
Bounding box size (unit: mm)Vertex (V)/Triangle numbers (T)Before shape modelingAfter shape modelingModel in Figure1668.3*38.1*50.05131508 (V)/262140 (T)148020 (V)/295164 (T)Model in Figure1767.4*28*49.8145548 (V)/290116(T)186849 (V)/372608 (T)As can be seen from Figure2, the minimum curvature calculation method proposed in this paper can detect the fusion regions effectively. After the 3D dental model has been analyzed quantitatively based on the minimum curvature, and processed further by applying the morphology operation, we can extract the target regions according to the corresponding regional characteristics (see Figure 7(b)).Figure6 shows that the spatial “line” solved by the direction tracing method proposed is an approximate geodesic curve, which has a linear time complexity ofO(n), where n is the vertex number of the “line”. The polygon selection is real-time, and the time consumed can be omitted. So, the target regions can be selected fast and intuitively (see Figure 7).When we use the weight rule proposed in this paper to triangulate the hole boundary, at the initial stage, because the number of the neighboring triangles is small, the boundary is triangulated primarily based on the adjacent angles and the perimeter of the candidate triangle. As shown in Figure18(b), the initial stage is also a process used to eliminate the saw tooth and smooth the boundary. As the boundary is triangulated continuously, the weight rule will select a vertex at the corner with the highest choice priority as the forwarding location. The two vertices of the new added edge usually have much higher choice priority than the rest of the boundary, which will drive the triangulation forward until a curtain like surface patch is covering at the boundary (see Figures 18(c) and 18(d)). The weight rule divides the triangulation process into boundary smoothing and boundary zipping approximately, by which a uniform and natural triangulation surface patch can be obtained.The time complexity of the proposed method is O(N log N) (N is the number of the boundary vertices).As can be seen from Figures10 and 11, the refinement surface patch can achieve a similar mesh density with the original model, which can avoid the situation of irregular triangles to appear when the surface patch is applied by the reshaping adjustment operation. During the reshaping adjustment stage, in order to control the deformation degree, the parameter λ was introduced in (18) to ensure the restored surface patch satisfies both the continuality and noninterference conditions. We apply the method proposed by Park [31] to detect the self-intersection. The parameter value ofλ is limited to the range from 0.8 to 1.0 based on a great deal of experimental analysis, and the adjustment stepτshould not be bigger than 0.01. Then, the deformation degree can be adjusted from λ=1 toλ=1-k*τ automatically until there is no intersection existing. We apply the incremental least squares method [32] to solve the reshaping matrix, which can reach the rate of 50000 vertices per seconds on the personal computer with P4, 2.4 GHz processor.We compared the triangulation quality (see Figure15) and efficiency (see Table 2) with the method proposed by Barequet and Sharir [18]. As can be seen from Figure 15, the method proposed in this paper can deal with complex holes with much more uniform triangulation result than that of Barequet and Sharir [18]. The triangulation efficiency of [18] is measured in minutes, and it takes no less than half an hour to deal with the 13–15 holes of the 3D dental model (without single-tooth missing), which cannot meet the actual efficiency needs. Statistical results of Table 2 show that the average triangulation and restoration efficiency can reach 20000 V/s and 4000 V/s, respectively by using the method proposed in this paper. The vertex number of the restoration surface patch is usually 1/30~1/40 of the original 3D dental model. So, the whole shape modeling procedure of the 3D dental model can be complete in 2~3 minutes.Table 2
Time consumption and number of vertices (V)/triangles (T) generated at different stages in shape restoration withC1 continuity.
Number of boundary verticesGeneratePminGeneratePrefine (including time consumedGeneratePfinal (including time consumedO (N3) time complexity of Barequet and Sharir [18]O(N log N) time complexity of method in this paperduring generatingPmin)during generatingPrefine)24191.40s (241V/239T)0.015s (241V/239T)0.297s (4581V/8826T)1.005s (5264V/10182T)437191.844s (437V/435T)0.016s (437V/435T)0.392s (6031V/11636T)1.232s (6844V/12962T)5871408.14s (857V/855T)0.031s (857V/855T)0.173s (2663V/4628T)0.677s (2733V/4828T)176912762.3s (1769V/1767T)0.062s (1769V/1767T)0.328s (4559V/7342T)1.104s (4743V/7986T)Table 3
Models corresponding to different boundaries in Table2 before and after shape restoration.Triangulation results comparison. (a) Bunny model with complex hole. (b) Triangulation result by Barequet and Sharir [18]. (c) Triangulation result by the method in this paper.
(a)(b)(c)Figure 16
Single-tooth modeling results of 3D dental model with normal tooth arrangement. (a) The initial 3D dental model. (b) Minimum curvature color map. (c) Separation results shown as a whole corresponding to the initial 3D dental model. (d) After the fusion regions being deleted. (e) After single-tooth shape being restored. (f) Separation results corresponding to the 3D dental model after shape restoration. (g) The separated tooth in (c) displayed, respectively. (h) The separated tooth in (f) displayed, respectively.Figure 17
Single-tooth Modeling results of 3D dental model with severe malocclusion. (a)The initial 3D dental model. (b) Minimum curvature color map. (c) Separation results shown as a whole corresponding to the initial 3D dental model. (d) After the fusion regions being deleted. (e) After single-tooth shape being restored. (f) Separation results corresponding to the 3D dental model after shape restoration. (g) The separated tooth in (c) displayed, respectively. (h) The separated tooth in (f) displayed, respectively.Triangulation mechanism. (a) Initial hole. (b) Eliminating saw tooth and smooth the boundary. (c) Closing the hole in the zipping way. (d) The triangulation results.
(a)(b)(c)(d)Figure14 shows that the blending regions between the teeth and the soft tissues can be extracted completely. Because the skeleton of the nontarget regions such as areas including alveolar bone ridges is open branches, which can be removed automatically in the pruning stage, the nontarget regions donot need to be removed interactively. The regions can be used to the extract segmentation boundary directly.The separated teeth in Figures16 and 17 show that the modeling techniques proposed in this paper can restore the shape of the single-tooth and segment the teeth correctly.In order to test the accuracy, we have done lots of experiments on comparing the restored tooth with its corresponding plaster single-tooth. Statistical results show that the radial deviation between this two models is usually ranging from 0 to 50 um, which can meet the clinical requirements.


## 7. Conclusion

In this paper, we proposed an integrated single-tooth modeling scheme, which is mainly composed of fusion regions extraction, single-tooth shape restoration and separation. As can be seen from the above examples, the modeling results can satisfy the biocharacteristics of the real tooth. Unlike the method based on plan-view range image of teeth, we directly compute bioinformation needed on the 3D dental model. We have demonstrated that the modeling scheme can achieve satisfying modeling results with high degree approximation of the original tooth and meet the requirements of clinical oral medicine.



---

*Source: 535329-2010-06-20.xml*