Abstract
In order to extract the pixels of teeth from 3D Cone Beam Computed Tomography (CBCT) image, in this paper, a novel 3D segmentation approach based on deformable surface mode is developed for 3D tooth model reconstruction. Different forces are formulated to handle the segmentation problem by using different strategies. First, the proposed method estimates the deformation force of vertex model by simulating the deformation process of a bubble under the action of internal pressure and external force field. To handle the blurry boundary, a “braking force” is proposed deriving from the 3D gradient information calculated by transforming the Sobel operator into three-dimension representation. In addition, a “border reinforcement” strategy is developed for handling the cases with complicate structures. Moreover, the proposed method combines affine cell image decomposition (ACID) grid reparameterization technique to handle the unstable changes of topological structure and deformability during the deformation process. The proposed method was performed on 510 CBCT images. To validate the performance, the results were compared with those of two other well-studied methods. Experimental results show that the proposed approach had a good performance in handling the cases with complicate structures and blurry boundaries well, is effective to converge, and can successfully achieve the reconstruction task of various types of teeth in oral cavity.
1. Introduction
Recently, more and more people pay attention to dental healthy [1]. Cone Beam Computed Tomography (CBCT) scan technology is recommended by physicians to help diagnose dental disease because of being noninvasive, painless, cost-effective, and portable [2]. Further, it is more suitable to present the high details of skeleton even inside of teeth [3]. However, it is not intuitive and detailed enough due to the limitation of its tomography imaging characteristics.
3D scanning technology can obtain a more accurate description of exacted shape of tooth crown surface [4]. The accurate 3D modeling of mandible and the simulation of tooth movement play an important role in preoperative planning for the dental and maxillofacial surgery. To reconstruct the 3D tooth model from CBCT image slices, it is necessary to segment the boundary of each tooth from background accurately. Many segmentation methods, which are usually based on the methods described later or their combinations, thresholding-based methods, boundary-based methods, and region-based methods, have been proposed and proved to be efficient for the 2D image segmentation [5–10]. Existing segmentation methods need to utilize additional preprocessing technologies to construct the 3D model from a series of image slices. However, these methods were easily affected by the changes of curve shape and topology structure, and only few methods are capable of extending to three-dimensional models directly. Several works have been done for the 3D reconstruction of human organ from CT slices and some commercial systems are available in practical applications [11]. However, no existing system is capable of generating such a tooth model where each tooth is modeled and manipulated independently from other teeth and jaw. The global thresholding-based approach, which is widely used in existing 3D reconstruction systems, is known to be effective for automatic hard tissue segmentation. However, in most of tooth CBCT images, the tooth regions are those where some teeth touch other teeth and some are located inside of alveolar bone with similar or same intensity profile. These factors led to blurry boundaries and complicate structures in tooth CBCT image. In addition, the teeth in CBCT have great variation in shape and size, and the structures of some teeth are quite complicate. Therefore, the tooth CBCT segmentation is a difficult task.
Deformable models, which include the popular deformable contours, or snakes, and deformable surfaces, are a powerful segmentation technique designed to meet the task of medical image segmentation and have proved to be successful in boundary integration and feature extraction [12]. Based on deformable surface model, [13] developed a Brain Extraction Tool (BET) that evolves to fit the brain’s surface by the application of a set of locally adaptive model forces. The method is very fast and requires no preregistration or other preprocessing before being applied. However, its two dimensional formulation limits its ability to process 3D image data. Reference [14] developed a topology adaptive snakes (T-snakes) model that overcomes the limitations of standard deformable models. T-snakes exploited affine cell image decomposition (ACID), a theoretically sound framework that significantly extends the abilities of standard parameters snakes. ACID induces an efficient reparameterization mechanism that enables deformable contours to flow or grow into complex geometries even modifying their topology as necessary. Therefore, immersing discrete parametric snakes in ACID enables them to handle the complex structures, with a high degree of automation, efficiency, and reproducibility in many medical image analysis scenarios. In addition, an important feature of ACID is that it does not interfere with the physics-based formulation of standard snakes, preserving the intuitive user interaction mechanism and the ability of incorporate constraints through energy or force functions.
In this paper, an effective segmentation approach based on deformable surface model is proposed which is able to conduct the segmentation and reconstruction process simultaneously while working on 3D CBCT volume data directly instead of individual slices. Based on deformable triangle mesh, the proposed method estimates the deformation force of vertex model by simulating the deformation process of a bubble under the action of internal pressure and external force field. To handle the blurry boundary, a “braking force” is proposed deriving from the 3D gradient information calculated by using 3D Sobel operator. In addition, a “border reinforcement” strategy is developed for handling the cases with complicate structures by weighting the points neighboring to each vertex according to its corresponding movement. Moreover, the proposed method combines ACID grid reparameterization technique to handle the unstable changes of topological structure and deformability during the deformation process. By combining CBCT and 3D scanning data, both entire topology structure information and occlusal surface details can be integrated conveniently. Experimental results with the implemented 3D reconstruction algorithm with CBCT image demonstrate that the proposed approach is efficient to segment the 3D images with complicate structure and blurry boundaries well and can successfully achieve the reconstruction task of various types of teeth in oral cavity.
Section 1 reviews the related background knowledge. In Section 2 we briefly introduce the proposed 3D segmentation approach based on deformable triangular mesh model. The reconstruction results are given in Section 3 and Section 4 concludes the paper.
2. The Proposed Method
2.1. Overview of the Proposed Method
The purpose of the proposed method is to partition the image into object and background regions to finish the 3D model reconstruction. First, a triangular mesh model is initialized and placed into the volume data space. For each vertex, a deformation force is defined according to the image information and triangular mesh information that drives the initial triangular mesh model to slowly deform, one per a time, following forces that keep the surface well-spaced and smooth, while attempting to move toward the target’s edge. The deformation process is iterated until the deformation model fits to the target shape in volume data. Finally, a reparameterization is employed to keep the uniformity of deformed mesh model. The block diagram of the proposed method is shown in Figure 1.
2.2. Surface Model Initialization
The segmentation of the proposed method is started by specifying initial surface model. For producing the initial surface model, a tessellated sphere is generated by starting with an icosahedron and iteratively subdividing each triangle into several smaller triangles. The spherical tessellated surface is initially centered on a point selected from any slice of CBCT images produced by user’s interaction, with its radius set to half of the estimated target radius. An example of initial deformation model is shown in Figure 2.
(a)
(b)
(c)
A deformable surface model is defined by a mesh within an image domain that allows an initialized model inside the target to deform by the forces coming from the surface itself and external forces derived from image data; this model is then “fitted” to the target surface in the image when it reaches an equilibrium in which the total force acting on each vertex is zero:where is the set of vertices of the deformable surface. The external forces are based on the image data and computed by some combinations of pressure forces, or potential force, or gradient vector flow. They pull the surface toward expected target boundaries and perform the shape initialization. The internal force is obtained by some combinations of mean curve motion, elastic forces, and Laplacian smoothing. They are designed to hold the surface together, keep it from bending too much, make the surface smooth, and perform the shape regularization.
2.3. Internal Force Estimation
The proposed approach simulates the deformation process of a bubble under the action of internal pressure and external force in the gradient field of CBCT volume data. Assuming that there is a vacuum (zero pressure) outside of the bubble; there has been a quality of ideal gas inside the bubble; and the bubble surface has a tension, thus, the relationship between the pressure difference between outside and inside of the bubble and the surface curvature can be defined by Young-Laplace [15]:where is the pressure difference between outside and inside of the bubble, , are the curvature radius of two main directions, respectively. is the surface tension coefficient. It is noted that the term is two times of mean curvature , and (2) can be simplified as
At the initial time, the bubble is in a balance state. Then we fill the bubble with predetermined pressure within the time length of 0, while applying an external pressure filed to the bubble which can keep the bubble in a balance state. At this time, the difference between inside and outside of the bubble for each point can be obtained by (4), and then the pressure of additional force filed at each point is
The external pressure of vertex iswhere is the normal vector and is the unit dimension of vertices. After that, the external force filed is removed and the internal force of each vertex is
2.4. Driving Force Estimation
The driving force that drives the model to deform is provided by two aspects: internal pressure and intensity information of the image. The intensity-based thresholding is one of the most widely used image preprocessing methods because of its simplicity and ease to implement it. In this paper, we model a basic driving force based on local intensity thresholding to distinguish foreground and background firstly [13]. The external driving force iswherewhere is the weight of external driving force that controls the contribution of to the resultant force. In this paper, is selected to 0.1 by the experiment. is the local threshold defined aswhere , are the max gray value and the min gray value of pixels along the direction of vertex normal vector with distance of the image, respectively. To suppress noise, we utilize the equations as follows:where is the gray value in 2% of histogram and is the median and is the threshold that partitions the image into foreground and background derived by using Otsu thresholding method [16].
2.5. External Force Estimation
The basic driving force derived from local intensity thresholding method is based on an assumption of single target in the image. However, it may lead to incorrect segmentation due to the complex scene, and it is not suitable for separating individual tooth regions from the CBCT image where some teeth touch others and some are located inside of alveolar bone with similar intensity. To overcome the above, we propose to add an external force (braking force) to stop the deformation when the deforming model reached to the real tooth boundaries in complex scene and blurry boundaries. In this study, we use the gradient information of 3D space to model the braking force (Figure 3). The gradient at point () in 3D space is expressed as
Sobel operator is a gradient operator that is commonly used in edge detection in 2D image [17]. In the proposed method, we transform the Sobel operator into three-dimensional forms to compute the gradient information in 3D model. Sobel operator consists of two independent operators: smoothing effect perpendicular to the derivative direction: , , ; central differences along the derivative direction: , , . In 3D, the Sobel kernel can be expressed by Tensor integral [18]:where is the one-dimensional filtering core for different directions. Then the Sobel kernel along direction in 3D is given as
Similarly, the kernels of the rest of two directions can be acquired by rotating the kernel of direction. For each point in the image, three components of gradient are computed by filtering the volume data with three Sobel kernels to construct the 3D gradient field. For simplifying the computation, the gradient modules are bound to by scaling the gradient field. Figure 4(a) is a cube CT image, Figure 4(b) is the corresponding gradient field. We note that the probability of model points that just coincides with corresponding pixels is quite small. In order to obtain the gradient at any point of space, the gradient of pixel with the closest distance to the current point found by applying the nearest point interpolation method is used as the gradient of current point. Based on 3D gradient field, the braking force at point () is defined aswhere is the weight of gradient field.
(a)
(b)
2.6. Border Reinforcement
In CBCT images, the shape of some teeth is quite complicate and irregular (e.g., blind tube shape). To the initial deformation model with regular shape (Figure 5(a)), it will gradually fit to some parts of tooth surface during the deformation process (Figure 5(b)). In this situation, some parts of bubble which arrived at the target surface edge (marked with ) stopped to deform but other parts of bubble (marked with and ) will attempt to move towards the rest of target surface edge unfilled with the bubble as the pressure increasing. As the deformation continues, oversegmentation may happen.
(a)
(b)
(c)
To the oversegmentation problem caused by the complicate target structures during the deformation process, we propose a “border reinforcement” strategy by adding a weight to each pixel that increases or decreases the weights of pixels neighboring to each vertex with vertices movement situation:
The proposed border reinforcement strategy is based on the following idea: if a vertex stops to deform in an iteration, the corresponding gradient filed will be increased. If a vertex does not move during several iterations, the corresponding gradient filed will be much stronger. Based on border reinforcement, the braking force is redefined as
2.7. Deformation Force and Vertex Motion
The deformation force that drives the deformation is made up by three components: internal force , originating from the surface tension and internal pressure; external force , originating from the intensity information of image; braking force , originating from the gradient information of image. Based on three kinds of forces, the resultant force iswhere the braking force is a passive force that enables the model to stop deformation when it reaches to the target edge.
Assume that the quality of each vertex model is and the initial speed is 0. During time, the displacement of each vertex under the action of resultant force is simplified aswhere is a constraint parameter to overcome oversegmentation selected according to . Since the internal pressure in the bubble plays a leading role in resultant force, in this paper, is computed as and is selected as .
2.8. Model Reparameterization
The proposed method utilizes triangular mesh model to initialize the deformable surface model. During deformation process, the position of each vertex is changed under the action of resultant force except topological structure. After several iterations, some vertices deviate far from the initial position but some others nearly stay in the original position, which may lead the deformed meshes to be nonuniform. Nonuniform meshes increase the errors of calculating discrete attribute information (e.g., curvature) and affect the deformation of surface model seriously. Therefore, it is necessary to reparameterize the deformed model.
To overcome the limitations of standard deformable models, in our study, affine cell image decomposition (ACID) is utilized to reparameterize the deformed triangle mesh model that significantly extends the abilities of standard parameters snakes. Instead of rectangular tessellation, we partition the space into tetrahedral cells by using the Coxeter-Freudenthal triangulation in 3D simplicial decomposition. Each grid is constructed by dividing the image volume using a uniform cubic grid and subdividing each cube into six tetrahedra (Figure 6(a)). The aim of subdividing the cube grid into tetrahedra is to provide an unambiguous framework for the creation of robust consistent local polygonal approximations of the surface of anatomic structures. The polygonal approximation is constructed from the intersection of the object surface with the edges of each boundary cell. The intersection points result in one triangle or one quadrilateral (which can be subdivided into two triangles) approximating the object surface inside each boundary (Figure 6(b)), where each triangle or quadrilateral intersects a tetrahedral cell on 3 or 4 distinct edges, respectively. The triangle (or quadrilateral) separates the positive vertices of the tetrahedral cell from the negative vertices. The set of all these triangles constitutes the polygonal approximation of the object surface. With both nonsimplicial and simplicial frameworks, an approximation to any desired degree of accuracy can be obtained by decreasing the size of the grid cells.
(a)
(b)
The reparameterization process is as follows: calculating the intersection of mesh model and tetrahedron; locating the position (inside or outside of the model) of grid vertices; removing invalid intersections and reconstructing a new triangle mesh model with valid intersections; traveling all tetrahedrons and utilizing the triangles intersecting the model as the triangles of reconstructed model. Finally, we assemble the coordinates of vertices corresponding to each reconstructed triangle mesh to obtain the reparameterized model. Figure 7 is the reparameterized result of a triangular mesh model by employing ACID reparameterization. Figure 7(a) is the original triangular mesh model; it is seen that the mesh distribution of model is nonuniform. Comparing to Figure 7(a), it is seen that the mesh distribution of reparameterized model (Figure 7(b)) is improved.
(a) Original mesh model
(b) Reparameterized result
2.9. Model Smoothing
To produce a smooth and good quality mesh surface of deformed model, in this paper, a global Laplacian smoothing approach is utilized [19]. The problem of mesh smoothing is translated into finding an approximating surface with a global minimization of a surface fairing energy. By adopting the Laplacian operator in a global way instead of a local way, the smooth mesh is constructed by solving a spare linear system, which is fast and efficient, in a noniterative way with feature constraints.
2.10. Parameters Initialization
The initial parameters of the proposed method are set as follows. Assuming that the deformation is caused by the internal pressure only without any other forces, then the largest deformation force is .
To avoid the oversegmentation caused by rapid deformation, is set to , where is the pixel pitch. To simplify computing, the proposed method is operated on image coordinates (). Accordingly, .
In addition, the minimum curve radius of target deformation is also computed according to the pixel pitch (). To ensure that the deformation model can keep the expansion trend at the maximum curve point but not shrinks under the force action of itself, it is set to . Similarly, to enable the deformation to stop at the boundary region, it is set to . Accordingly, , .
3. Results
A database including 510 CBCT images is utilized as the test images for evaluating validation of the proposed method. The test images were obtained from Dental Hospital of Wu Han University and captured directly from CBCT image slices. The resolution of each image is with pixel pitch and 0.3 mm vertical spacing of cross section. Special attention is paid to protect the privacy of the patients.
In the experiment, two other well-studied approaches, BET [13] and reparameterized BET, are applied to the same images for comparison. The initial condition is set as follows: First, we choose a frame with distinct tooth boundaries from the CBCT image slices. Then, we initialized the deformation models at four seed points selected from four kinds of tooth regions: incisor, canine, the first premolar, and the first molar in the left lower jaw, respectively. With the same initial condition, we deform the initial surface model with the proposed approach, BET and reparameterized BET, and then utilize the deformed model volume after each iteration as the similarity criteria.
Figure 8 is the reconstruction results by segmenting the images with BET, reparameterized BET, and the proposed method. Comparing to the reconstruction results of two other methods, the hollow structure of pulp cavity is recognized correctly and the tubular tooth model is reconstructed more accurately by the proposed method. To the regions between roots and alveolar with the similar intensity distribution, it is also successfully distinguished by the proposed method.
Another statistical evaluation of the efficiency of the proposed method is carried out with model volume variation via iteration during 500 iterations of BET method, reparameterized BET (R-BET) method, and the proposed method in four types of tooth derived from the center incisor of left mandibular, canine, the first premolar, and the first molar, respectively. It is seen from Figure 9 that the volume variation ratio of the proposed method in different types of tooth is much higher than that of two other methods. The experiments were implemented on a Core 4 PC with 8 GB RAM, and the process time of BET was about 5–10 sec in average. The average process speed of the proposed method takes about 3–6 sec. This demonstrates that the proposed method can converge more quickly and needs less iteration to achieve a better performance.
Figure 10 is the reconstruction results of several individual teeth including one molar, one lateral incisor, and one central incisor, and then we give out the reconstruction results of all teeth in our database to evaluate the overall performance of the proposed method. It is seen that each individual tooth is effectively modeled and extracted, and successful segmentation is achieved. This demonstrates that the reconstruction result of the proposed method is accurate and reliable.
(a)
(b)
(c)
(d)
4. Conclusions
In this study, an effective 3D CBCT image segmentation approach is proposed based on deformation surface model for reconstructing the precious tooth model from the CBCT image slices. We propose a braking force that enables the segmentation process to stop the deformation when the contour of a tooth expands out to other teeth boundaries in deformation process. We also develop a “border reinforcement” strategy by weighting the points according to motion of vertices to overcome the problems that the deformed model exceeds the desired boundaries in the regions with complicate structures during the deformation process. Finally, we reparameterize the deformed model by employing ACID. In the experiment, a series CBCT image slices were employed to validate the reconstruction performance of the proposed method, and the reconstruction results of the proposed method were compared with the results of BET and reparameterized BET (R-BET) method. In addition, four types of tooth were utilized to evaluate the efficiency by evaluating model volume variation via iteration during 500 iterations. The process time of the proposed method is also more quickly than that of two other methods. Experimental results demonstrate that the proposed method can extract the 3D models from CBCT images with blurry boundaries and complicate structure well and reconstruct the 3D tooth model accurately and effectively.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
Financial support from the National Nature Science Foundation of China (NSFC) is greatly appreciated (Grant no. 61271093).