Three-dimensional tooth modeling method based on mixed level set
Technical Field
The invention belongs to the field of image processing, and relates to application of a level set method in three-dimensional tooth modeling.
Background
In recent years, with the rapid development of three-dimensional digital imaging technology in the field of oral medicine, computer-aided diagnosis and treatment repair is increasingly applied to oral repair and gradually becomes the development trend in the field. In a computer oral restoration system, it is important to acquire a digitized three-dimensional tooth model, and the accuracy and integrity of the model are directly related to the results of subsequent tooth arrangement, implantation, orthodontics, and biomechanical analysis.
At present, the most common method for tooth modeling is to segment the tooth CT image sequence by using an image processing technique, i.e., first segment the tooth contour from each layer of CT slices, and then reconstruct the tooth three-dimensional model by using these inter-layer contours. Because the method can obtain the whole tooth shape structure and provide complete diagnosis basis for the oral lesion of a patient, the tooth modeling method based on the CT image is more and more concerned by researchers.
Because the density and distance of teeth and jaw bone in the oral cavity CT image are relatively close, it is difficult to accurately extract the tissue contour of each tooth by adopting the traditional image segmentation method. CT image tooth segmentation has been a challenging problem. Wangli et al (Wangli, Chinay, Konjac et al, tooth 3-dimensional solid model establishment based on CT images, Chinese image graphic report, 2005, 10 (10): 1289-; however, the binarization operation is prone to generate over-segmentation or under-segmentation. Wu et al (X.Wu, H.Gao, H.Heo, et al. improved B-spline conforming using genetic algorithm for the segmentation of dental computer science image sequences. the Journal of imaging science and technology, 2007, 51 (4): 328-336.) utilize a B-spline curve fitting based on a genetic algorithm to extract the tooth profile of each slice, but the B-spline curve cannot deal with the problem of tooth topology changes. Gao et al (H.Gao, O.Chae.Indai. dental promotion from CT images using level set method with shape and sensitivity. Pattern Recognition, 2010, 43: 2406-; the model mainly depends on the probability distribution of image edge gradient and prior shape to guide the evolution of the level set, but because the tooth density is not uniform and the surrounding is easily interfered by structures such as alveolar bone, the problem of boundary leakage is easily generated by controlling the contraction and expansion of the level set outline by the gray probability of the surrounding area of the prior shape.
Disclosure of Invention
In order to overcome the defects of the prior art and improve the efficiency and the precision of tooth segmentation, the invention provides a three-dimensional tooth modeling method based on a mixed level set by comprehensively considering the tooth shape change characteristics of an oral CT image.
The present invention provides: a three-dimensional tooth modeling method based on a mixed level set comprises the following steps:
(1) selecting one of the oral cavity CT image sequences as an initial slice, and drawing a contour of each tooth on the image to initialize a level set function;
(2) cutting the root layer below the initial section by utilizing a single-phase mixed level set model constructed by combining prior shape constraint energy, edge energy based on a Flux model and local region energy based on prior gray level;
(3) cutting the tooth profile of the dental crown layer slice above the initial slice by using a biphase mixed level set model in competitive constraint of a binding area;
(4) and converting all the tooth profile pixel points after the section segmentation into three-dimensional coordinates, and reconstructing by using a Delaunay triangulation method to obtain a triangular mesh model of each tooth.
The step (1) is carried out according to the following steps: selecting a section with all teeth and less alveolar bone from the section image of the tooth neck part as a starting section, and drawing an approximate outline C of each tooth on the section imageiN (i ═ 1,2,. n, n is the number of teeth) as the initial contour of the level set, and then C is usediInitializing n level set functions Φi(i=1,2,...n),ΦiBy calculating each point on the image to CiIs completed by the signed distance of:
wherein d [ (x), Ci]Representing pixel point x and curve CiThe euclidean distance between them.
Defining an energy functional of the single-phase mixed level set model for segmenting the dental root layer slice in the step (2) as a weighted sum of prior shape constraint energy, edge energy based on a Flux model, local region energy based on prior gray level and the like:
Etooth root(Φ)=μEint(Φ)+γElength(Φ)+αEprior(Φ)+vEedge(Φ)+λEregion(Φ)
Wherein mu, gamma, α, nu and lambda are weight coefficients of each energy item;
(3a) symbol distance maintenance energy Eint(Φ) to ensure stability during level set evolution, defined as:
(3b) curve arc length smoothing energy Elength(Φ) for smoothing the level set contour, defined as:
(3c) prior shape constrained energy Eprior(Φ) for controlling the shape of the level set, mapping the tooth profile after each segmentation to the adjacent slice image, and constraining the tooth profile as the prior shape of the current level set function evolution, wherein the energy functional is defined as:
where Φ is the level set function of the current slice, ΦpA level set function corresponding to the prior shape after the last slice is segmented, wherein H (x) is a Heaviside function;
(3d) flux model-based edge energy Eedge(Φ) for detecting the contour of the outer boundary of the tooth, embedding angle information between the gradient direction of the image and the gradient direction of the horizontal function in the imageIn the unified Flux model, the energy functional is defined as:
Eedge(Φ)=-∫ΩξΔIσ(1-H(Φ)dx
wherein Δ is Laplacian operator, IσRepresenting the image after the gaussian smoothing and,
for the purpose of the gradient direction detection function,as gradient operator,. represents the dot product;
(3e) local region energy E based on prior gray scaleregion(Φ) for overcoming the problem of uneven image gray scale, embedding the prior gray scale information into a region model, wherein an energy functional is defined as:
wherein f isref_in(x) And fref_out(x) Respectively defining r neighborhood of a reference point on the prior shape as the gray average values inside and outside the prior shape curve, wherein the reference point is the point closest to a current image target pixel point x on the prior shape;
(3f) and (3) combining the energy terms, and expressing the minimized energy functional of the root layer slice mixed level set model as follows:
minimizing the energy functional to obtain an evolution equation of a level set curve:
where δ (Φ) is the Dirac function.
The biphase mixed level set model for segmenting the dental crown layer slices in the step (3) is established according to the following steps: the initial slice is segmented by using the single-phase mixed level set model, all segmented level set functions are combined into two coupled biphase mixed level set functions according to the rule of every other tooth, and region competition constraint energy is added to overcome the two level sets to generate overlap, wherein the energy functional is defined as:
Edental crown(Φ1,Φ2)=ETooth root(Φ1)+ETooth root(Φ2)+βErepulse
Wherein E isrepulse=∫ΩH(Φ1)H(Φ2)dx
To constrain the energy for region competition, β is used to control the degree of overlap of the two level set function regions;
minimizing the energy functional to obtain a biphase level set function phi1,Φ2Respectively as follows:
and
and
wherein phip1,Φp2Each represents phi1,Φ2ξ1(x)、ξ2(x) Each represents phi1,Φ2Gradient direction detection function during evolution:
the three-dimensional tooth reconstruction in the step (4) is carried out according to the following steps: solving an evolution equation of the level set function by using a narrow-band method and a semi-implicit difference scheme, extracting pixel points of the updated level set function at phi of 0 to obtain a tooth profile of the current slice image, converting the pixel points of the tooth profile after all the slice images are segmented into three-dimensional coordinates, and reconstructing by using a Delaunay triangulation method based on region growth to generate a triangular mesh model of each tooth.
The invention has the beneficial effects that: the mixed level set model provided by the invention combines information in various aspects such as image edges, local areas, priori knowledge and the like, and can effectively solve the problems that the traditional model is inaccurate in edge positioning and cannot process uneven image gray scale and the like. The method has less manual intervention, better segmentation effect and higher accuracy, and the reconstructed three-dimensional tooth model can correctly reflect the anatomical form of the tooth, thereby laying a solid foundation for further making oral restoration planning, biomechanical analysis and the like and having important significance for improving the accuracy and efficiency of oral restoration.
Drawings
FIG. 1 is a flow chart of a dental modeling technique of the present invention.
Fig. 2 is a diagram illustrating the initialization of the level set function for a start slice.
Fig. 3(a) is a schematic diagram showing the gradient direction of the tooth and the surrounding interfering structures.
FIG. 3(b) is a schematic diagram showing the gradient direction of the level set function.
FIG. 4 is a graphical representation of a biphasic level set function of a slice image of a crown layer.
Fig. 5 is a schematic diagram of a bi-phase level set partition.
FIG. 6 is a diagram illustrating the segmentation effect of the present invention and the prior art method on the initial slice image.
FIG. 7 is a graph showing the effect of the present invention on the segmentation of the molar region of a root segment.
FIG. 8 is a diagram of the segmentation contour and three-dimensional reconstruction effect of a plurality of slices at the molar part of the present invention.
Detailed Description
The embodiments of the invention will be further described with reference to the accompanying drawings in which:
as shown in fig. 1, the specific implementation steps of the present invention are as follows:
step 1, selection of initial slices and initialization of a level set: first, a tooth CT image sequence is read, and then a slice in which all teeth appear and alveolar bones rarely appear is selected as a starting slice from the slice images of the cervical region. On the initial slice image, the approximate contour C of each tooth is outlinedi(i 1, 2.. n, n is the number of teeth) as the initial contour of the level set, as shown in fig. 2, and then using Ci(i ═ 1,2,. n) initializing n level set functions Φi(i=1,2,...n)。ΦiBy calculating each point on the image to CiIs completed by the signed distance of:
wherein d [ (x), Ci]Represents point x and curve CiThe euclidean distance between them.
Step 2, cutting root laminas: in the root canal image, the tooth is susceptible to more interference from structures such as alveolar bone and pulp cavity. Therefore, a new mixed level set tooth segmentation model is proposed. The energy functional of the model mainly comprises prior shape constraint energy, Flux model-based edge energy, prior gray level-based local region energy and the like.
(1) In the prior shape constraint energy term, because the shapes and the positions of corresponding teeth are relatively close to each other between adjacent slices, the tooth profile of the segmented previous slice is mapped to the current slice image and is constrained as the prior shape of level set evolution. Let Φ be the level set function of the current slice, ΦpAnd defining the prior shape constraint term as the level set function corresponding to the prior shape after the previous slice is segmented:
wherein H (x) is a Heaviside function.
(2) In the Flux-based edge energy, a gradient vector field-based Flux model is used to locate the target boundary, while the angle between the image gradient direction and the horizontal function gradient direction is used to overcome the interference problem of the surrounding alveolar bone and pulp cavity.
Because the gradient of the image at the boundary is in the same direction as the normal vector of the edge, the Flux model positions the edge of the image by calculating the edge integral of the inner product of the gradient of the image and the normal vector of the curve, and the method can effectively solve the problem of segmentation of the weak boundary target and has high evolution speed. The maximum energy functional of the Flux model is expressed as:
EFlux(Φ)=∫ΩΔI(1-H(Φ)dx
in the formula, Delta is Laplacian operator, IσRepresenting a gaussian smoothed image.
However, since there are many interfering structures inside and outside the tooth, the division is performed directly by using the above model, and the level set is captured at the same time on the boundaries of the pulp cavity inside the tooth and the adjacent teeth or alveolar bone outside the tooth, so the present invention adds the gradient direction constraint to the above Flux model.
In the tooth CT image, the tooth is a high-brightness region, and the background is a low-brightness region, so the image gradient direction of the tooth and the surrounding interfering structure can be represented as shown in fig. 3 a. Since the level set function Φ as defined herein takes positive inside and negative outside the level set, the gradient direction of Φ can be represented as shown in fig. 3 b. It can be seen that if the evolution curve is to capture only the outer boundary contours of the tooth, the gradient direction of the image should coincide with the gradient direction of the level set function. Therefore, the minimization functional of the edge energy of the present invention is defined as:
Eedge(Φ)=∫ΩξΔI(1-H(Φ)dx
wherein
For the purpose of the gradient direction detection function,for the gradient operator,. represents the dot product.
(3) In the local area energy combined with the prior gray scale, the global mean value is replaced by the local neighborhood gray scale mean value of the prior shape, so that the problems that the prior gray scale information cannot be utilized and the image gray scale unevenness cannot be processed in the traditional model are solved.
Because the gray distribution of the images between adjacent slices is very same and the tooth profile positions are also close, the prior gray information is embedded into an energy model to improve the comparison accuracy, and the energy functional is defined as:
wherein f isref_in(x) And fref_out(x) Defined as the mean gray values of the r neighborhoods of the reference points on the prior shape inside and outside the prior shape curve respectively. The reference point is determined as follows: after the previous slice is divided, calculating the gray average value of the r neighborhood of each point in the prior shape inside and outside the prior shape curve, and then taking the point which is closest to the target pixel point x of the current image in the prior shape as a reference point in the current slice image. The radius of the r neighborhood is chosen according to the resolution of the image and should not be too large in order to avoid containing too many neighboring tooth regions.
(4) To ensure the shape and stability of the level set in the whole evolution process, the symbol distance is added to maintain energy:
and a curve arc length smoothing term:
thus, by integrating the above energy information, the energy minimization functional of the single-phase mixed level set model for segmenting the root layer slices can be expressed as:
minimizing the energy equation to obtain an evolution equation of the level set function phi of the root layer slice as follows:
where δ (Φ) is the Dirac function.
Step 3, cutting the dental crown layer into slices and cutting the model: in the dental crown layer slice image, teeth are gradually connected together, the gap is small, in order to effectively extract the boundary between two adjacent teeth, the single-phase mixed level set model described in step 2 is used for tooth segmentation on the initial slice, and the segmented level set functions are combined into two coupled biphase mixed level set functions according to the rule of every other tooth, as shown in fig. 4.
Known from the multi-phase level set partitioning principle, two level set functions phi1、Φ2The image may be divided into four regions as shown in fig. 5. For [ phi ]1>0,Φ2> 0, representing the overlapping area of adjacent teeth. To overcome phi1、Φ2Overlap is generated in the evolution process, and a regional competition criterion is introduced into the energy functional. Let A1、A2Respectively as a level set function phi1、Φ2The corresponding horizontal set curve encloses an inner area. A may be determined by the nature of the Heaviside function H (x)1,A2Respectively expressed as:
A1=area(Φ1≥0)=∫ΩH(Φ1)dx
A2=area(Φ2≥0)=∫ΩH(Φ2)dx
to make these two regions not overlap, it should satisfy:
the region penalty term can therefore be defined as:
Erepulse=∫ΩH(Φl)H(Φ2)dx
minimizing the above equation amounts to avoiding two teeth from overlapping.
Adding the region penalty term into the level set energy functional to obtain an energy equation of the biphase mixed level set model of the dental crown layer, wherein the energy equation comprises the following steps:
Edental crown(Φ1,Φ2)=ETooth root(Φ1)+ETooth root(Φ2)+βErepulse
β is used to control the overlapping degree of the two level set function regions, and the energy functional is minimized to obtain a biphase level set function phi1,Φ2Respectively as follows:
and
wherein phip1,Φp2Each represents phi1,Φ2ξ1(x)、ξ2(x) Each represents phi1,Φ2Gradient direction detection function in evolution, namely:
step 4, reconstructing a tooth three-dimensional model: solving an evolution equation of the level set function by using a narrow-band method and a semi-implicit difference scheme, extracting pixel points of the updated level set function at phi 0, converting the tooth profile pixel points after all the slice images are segmented into three-dimensional coordinates, and reconstructing by using a Delaunay triangulation method to obtain a triangular mesh model of each tooth.
The effectiveness of the present invention can be further illustrated by the following experiments:
take a CT image sequence of a patient's mouth with intact mandibular teeth as an example. FIG. 6 is a graph comparing the tooth segmentation results of the initial slice images of the CT sequence using the present invention and Li and CV models. Fig. 6(a) is an initial contour polygon of each tooth drawn on the starting slice. It can be seen that the edge-based Li model has many phenomena of inaccurate tooth edge positioning in fig. 6(b) due to the fact that it is too dependent on image gradient information near the initial curve, but lacks the gradient direction constraint. The region-based C-V model relies on global gray scale information, thus creating fusion between closely spaced teeth and at weak boundaries between teeth and jaws in fig. 6 (C). Fig. 6(d) shows the segmentation result of the present invention, which can effectively detect weak boundaries by using a Flux model based on a gradient vector field, and accurately position the tooth boundary contour by adding the gradient direction and the prior shape constraint.
FIGS. 7(a) and 7(b) are the results of the segmentation of the root site slice images of the second molar of the patient using the method of the present invention and Gao et al, respectively. It can be seen that the present invention has a more accurate segmentation result, but the method of Gao et al detects the profile by means of the probability of the inner and outer gray levels of the curve, and the detected curve is biased to the pulp chamber at the weak boundary of the tooth because the gray level of the pulp chamber is darker.
Fig. 8 shows the result of segmentation of some slices at the right first molar site. It can be seen that the segmentation results for substantially all slice images fit well to the contour of the dental target. However, in slice 160, the level set does not capture the tooth's internal groove contour, which is mainly affected by the size of the narrow band in the level set evolution equation. The presence of the grooves in the tooth surface can cause extreme dishing of the contours of some of the slices. If the narrow bandwidth is now made smaller, it may be possible to limit the evolution of the level set onto the target profile. However, this does not affect the final segmentation result, and since the target contour is included in the level set curve, the final contour can be extracted by the grayscale threshold. The right side of fig. 8 shows the result of three-dimensional reconstruction of the profile of each layer of sliced teeth after segmentation, and it can be seen that the three-dimensional model reconstructed by the invention can correctly reflect the anatomical form of the teeth, thereby providing effective data for the next computer-aided diagnosis and enabling doctors to make the optimal surgical plan according to the actual situation of patients.
The above description is only a preferred embodiment of the present invention, and should not be taken as limiting the invention, and any modifications, equivalents, improvements, etc. made within the spirit and principle of the present invention should be included in the scope of the present invention.