CN105488849B - A kind of three-dimensional tooth modeling method based on mixed-level collection - Google Patents

Abstract

The present invention provides a three-dimensional tooth modeling method based on a mixed level set, comprising the following steps: 1) selection of a starting slice and level set initialization; 2) utilizing prior shape constraint energy, edge energy based on a Flux model, The single-phase mixed level set model based on the combination of prior grayscale and local region energy to segment the root layer slice image; 3) The two-phase mixed level set model combined with regional competition constraints to segment the crown layer slice image; 4) Based on the segmentation contour 3D model reconstruction of teeth. The mixed level set model of the present invention can effectively overcome the problems of inaccurate edge positioning of traditional segmentation models and the inability to deal with uneven gray levels of images, etc., and quickly and accurately build an independent three-dimensional model of each tooth, so as to provide a basis for the formulation of oral restoration planning, biomechanics, etc. A solid foundation of analysis, etc.

Description

A 3D Teeth Modeling Method Based on Mixed Level Set

technical field

The invention belongs to the field of image processing and relates to the application of a level set method in three-dimensional tooth modeling.

Background technique

In recent years, with the rapid development of three-dimensional digital imaging technology in the field of stomatology, more and more computer-aided diagnosis and treatment have been applied in oral restoration, and it has gradually become the development trend in this field. In the computer oral restoration system, the first important thing is to obtain a digital 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 commonly used method for tooth modeling is to use image processing technology to segment tooth CT image sequences, that is, firstly segment tooth contours from each layer of CT slices, and then use these inter-layer contours to reconstruct a three-dimensional tooth model. Because this type of method can obtain the shape and structure of the entire tooth and provide a complete diagnosis basis for patients' oral lesions, the tooth modeling method based on CT images has attracted more and more attention from researchers.

Because the densities and distances of teeth and jaws in oral CT images are relatively close, it is difficult to accurately extract the tissue contour of each tooth using traditional image segmentation methods. Teeth segmentation in CT images has always been a challenging subject. Wang Li et al. (Wang Li, Cui Jin, Han Qingkai, etc. Establishment of a 3D solid model of teeth based on CT images. Chinese Journal of Image and Graphics, 2005, 10(10): 1289-1292.) used binarization and boundary extraction to screen The key points of the tooth contour of each slice are obtained, and then the solid model of the whole tooth is obtained by using the 3D-Delaunay tetrahedronization algorithm; however, the binarization operation is prone to over-segmentation or under-segmentation. Wu et al. (X. Wu, H. Gao, H. Heo, et al. Improved B-spline contourfitting using genetic algorithm for the segmentation of dental computerized tomography image sequences. The Journal of imaging science and technology, 2007, 51(4): 328-336.) B-spline curve fitting based on genetic algorithm is used to extract the tooth profile of each slice, but the B-spline curve cannot deal with the change of tooth topology. Gao et al. (H.Gao, O.Chae.Individual toothsegmentation from CT images using level set method with shape and intensityprior.Pattern Recognition, 2010, 43: 2406-2417) used a level set-based active contour model for tooth segmentation, and Different segmentation models are used for different tooth layer slices, and the crown and root of each tooth can be segmented; the model mainly relies on the image edge gradient and the probability distribution of the prior shape to guide the evolution of the level set, but due to the tooth density It is not uniform, and the surrounding area is easily disturbed by structures such as alveolar bone, so relying on the gray probability of the surrounding area of the prior shape to control the shrinkage and expansion of the level set contour is prone to the problem of boundary leakage.

Contents of the invention

In order to overcome the shortcomings of the above-mentioned prior art and improve the efficiency and accuracy of tooth segmentation, the present invention comprehensively considers the characteristics of tooth shape changes in oral CT images, and proposes a 3D tooth modeling method based on mixed level sets.

The invention provides: a method for modeling three-dimensional teeth based on mixed level sets, comprising the following steps:

(1) Select one from the oral CT image sequence as the starting slice, and draw the outline of each tooth on the image to initialize the level set function;

(2) The root layer slice below the initial slice uses the single-phase mixed level set model constructed by combining the prior shape constraint energy, the edge energy based on the Flux model, and the local area energy based on the prior gray level to segment the tooth contour;

(3) The crown layer slice above the initial slice is segmented by the two-phase mixed level set model combined with regional competition constraints;

(4) Transform the pixel points of the tooth contour after all slices into three-dimensional coordinates, and use the Delaunay triangulation method to reconstruct to obtain the triangular mesh model of each tooth.

Step (1) is carried out as follows: select a slice in which all teeth appear and less alveolar bone appears in the slice image of the dental neck as the starting slice, and outline each tooth on the slice image The approximate contour C _i (i=1, 2,...n, n is the number of teeth) is used as the initial contour of the level set, and then use C _i to initialize n level set functions Φ _i (i=1, 2,. ..n), the initialization of Φ _i is done by calculating the signed distance from each point on the image to C _i , namely:

Wherein, d[(x), C _i ] represents the Euclidean distance between the pixel point x and the curve C _i .

The energy functional of the single-phase mixed level set model for segmenting the root layer slices described in step (2) is defined as the weighted sum of the prior shape constraint energy, the edge energy based on the Flux model, and the local area energy based on the prior grayscale, etc. :

E _{tooth root} (Φ)＝μE _int (Φ)+γE _length (Φ)+αE _prior (Φ)+vE _edge (Φ)+λE _region (Φ)

Among them, μ, γ, α, ν, λ are the weight coefficients of each energy item;

(3a) The signed distance preserving energy E _int (Φ), which is used to ensure the stability in the evolution process of the level set, is defined as:

(3b) Curve arc length smoothing energy E _length (Φ), used to smooth the level set profile, defined as:

(3c) The prior shape constraint energy E _prior (Φ), used to control the shape of the level set, maps the tooth contour after each segmentation to the adjacent slice image, and constrains it as the prior shape of the evolution of the current level set function, Its energy functional is defined as:

Among them, Φ is the level set function of the current slice, Φ _p is the level set function corresponding to the prior shape after the previous slice is segmented, and H(x) is the Heaviside function;

(3d) The edge energy E _edge (Φ) based on the Flux model is used to detect the outer boundary contour of the tooth, and the angle information between the image gradient direction and the horizontal function gradient direction is embedded into the traditional Flux model, and its energy functional definition for:

E _edge (Φ)＝- _∫Ω ξΔI _σ (1-H(Φ)dx

Where Δ is the Laplacian operator, I _σ represents the image after Gaussian smoothing,

is the gradient direction detection function, is the gradient operator, and represents the dot product;

(3e) The local region energy E _region (Φ) based on the prior gray scale is used to overcome the problem of uneven gray scale of the image, and embed the prior gray scale information into the regional model, and its energy functional is defined as:

Among them, f _{ref_in} (x) and f _{ref_out} (x) are respectively defined as the gray mean value of the r neighborhood of the reference point on the prior shape inside and outside the prior shape curve, and the reference point is the distance from the current image target pixel on the prior shape point x nearest point;

(3f) Combining the above energy items, the minimized energy functional of the mixed level set model of the root slice is expressed as:

The above energy functional is minimized to obtain the evolution equation of the level set curve:

Among them, δ(Φ) is the Dirac function.

The two-phase mixed level set model of the segmented crown layer slice described in step (3) is established as follows: the initial slice is segmented using the single-phase mixed level set model, and all level set functions after segmentation According to the rule of every other tooth, it is combined into two coupled two-phase mixed level set functions, and the area competition constraint energy is added to overcome the overlapping of the two level sets, and its energy functional is defined as:

E _{tooth crown} (Φ ₁ , Φ ₂ ) = E _{tooth root} (Φ ₁ ) + E _{tooth root} (Φ ₂ ) + βE _repulse

Among them, E _repulse ＝∫ _Ω H(Φ ₁ )H(Φ ₂ )dx

Constrain energy for regional competition, and β is used to control the overlapping degree of two level set functions;

Minimize the above energy functional to obtain the evolution equations of the two-phase level set functions Φ ₁ and Φ ₂ respectively:

with

with

Among them, Φ _p1 and Φ _p2 represent the prior values of Φ ₁ and Φ ₂ respectively; ξ ₁ (x) and ξ ₂ (x) represent the gradient direction detection function when Φ ₁ and Φ ₂ evolve respectively:

The three-dimensional tooth reconstruction described in step (4) is carried out according to the following steps: use the narrow-band method and the semi-implicit difference scheme to solve the evolution equation of the above-mentioned level set function, extract the pixel point of the updated level set function at Φ=0, and obtain For the tooth contour of the current slice image, convert the tooth contour pixel points of all slice images into three-dimensional coordinates, and use the region growing-based Delaunay triangulation method for reconstruction to generate a triangular mesh model of each tooth.

Beneficial effects of the present invention: the mixed level set model proposed by the present invention combines various information such as image edges, local regions, prior knowledge, etc., and can effectively overcome the problems of inaccurate edge positioning and inability to deal with uneven gray levels of images in traditional models . The method of the present invention has less manual intervention, and has better segmentation effect and higher accuracy rate, and the reconstructed three-dimensional tooth model can correctly reflect the anatomical shape of the tooth, so as to formulate oral restoration planning and biomechanical analysis for the next step It is of great significance to improve the accuracy and efficiency of dental restoration.

Description of drawings

Fig. 1 is a flow chart of the tooth modeling technique of the present invention.

Figure 2 is a schematic diagram of the initialization of the level set function of the initial slice.

Fig. 3(a) is a schematic diagram of the gradient directions of teeth and surrounding interference structures.

Figure 3(b) is a schematic diagram of the gradient direction of the level set function.

Fig. 4 is a schematic diagram of a two-phase level set function of a slice image of a dental crown layer.

Fig. 5 is a schematic diagram of biphasic level set segmentation regions.

Fig. 6 is an effect diagram of segmentation of the initial slice image by the present invention and the existing method.

Fig. 7 is a diagram showing the segmentation effect of the molar part of the tooth root slice by the present invention and the existing method.

Fig. 8 is a diagram showing segmentation contours and three-dimensional reconstruction effects of several slices of molar parts according to the present invention.

detailed description

Embodiments of the present invention will be further described below in conjunction with accompanying drawings:

As shown in Figure 1, the specific implementation steps of the present invention are as follows:

Step 1. Selection of the starting slice and initialization of the level set: first read the tooth CT image sequence, and then select a slice in which all the teeth appear and the alveolar bone rarely appears in the slice image of the dental neck as the starting point slice. On the initial slice image, outline the approximate outline C _i (i=1, 2, ... n, n is the number of teeth) of each tooth as the initial outline of the level set, as shown in Figure 2, and then use C _i (i=1, 2,...n) initializes n level set functions Φ _i (i=1, 2,...n). The initialization of _Φi is done by computing the signed distance from each point on the image to _Ci , namely:

Among them, d[(x), C _i ] represents the Euclidean distance between the point x and the curve C _i .

Step 2. Root layer slice segmentation: In the root layer image, the teeth are easily disturbed by structures such as alveolar bone and pulp cavity. Therefore, a new hybrid level set tooth segmentation model is proposed. The energy functional of the model is mainly composed of prior shape constraint energy, edge energy based on Flux model, and local area energy based on prior gray level.

(1) In the prior shape constraint energy item, since the shape and position of the corresponding teeth are relatively close between adjacent slices, the tooth contour that has been segmented in the previous slice is mapped to the current slice image as a level set The evolution of the prior shape is constrained. Let Φ be the level set function of the current slice, and Φ _p be the level set function corresponding to the prior shape after the previous slice is segmented, then the prior shape constraint term is defined as:

In the formula, H(x) is the Heaviside function.

(2) In edge energy based on the Flux model, use the Flux model based on the gradient vector field to locate the target boundary, and at the same time use the angle between the gradient direction of the image and the gradient direction of the horizontal function to overcome the surrounding alveolar bone and pulp cavity interference problem.

Since the gradient of the image at the boundary is in the same direction as the edge normal vector, the Flux model locates the edge of the image by calculating the edge integral of the inner product of the image gradient and the curve normal vector. This method can effectively solve the segmentation problem of weak boundary targets, and the evolution speed is fast . The maximized energy functional of the Flux model is expressed as:

E _Flux (Φ)＝ _∫Ω ΔI(1-H(Φ)dx

In the formula, Δ is the Laplacian operator, and I _σ represents the image after Gaussian smoothing.

However, since there are many interfering structures inside and outside the tooth, the above-mentioned model is directly used for segmentation, and the horizontal assembly can be captured on the pulp cavity inside the tooth and the boundary of the adjacent tooth or alveolar bone outside, so the present invention divides the gradient direction Constraints are added to the Flux model above.

In the tooth CT image, the tooth is a high-brightness area, and the background is a low-brightness area, so the image gradient direction of the tooth and surrounding interference structures can be expressed as shown in Figure 3a. Since the level set function Φ defined in this paper is positive inside the level set and negative outside, the gradient direction of Φ can be expressed as shown in Figure 3b. It can be seen that if the evolution curve only captures the outer boundary contour of the tooth, the gradient direction of the image should be consistent with the gradient direction of the level set function. Therefore, the minimized functional of the edge energy of the present invention is defined as:

E _edge (Φ)＝ _∫Ω ξΔI(1-H(Φ)dx

in

is the gradient direction detection function, is the gradient operator, and represents the dot product.

(3) In the local area energy combined with the prior gray level, the local neighborhood gray level mean of the prior shape is replaced by the global mean value to overcome the inability of the traditional model to use the prior gray level information and the inhomogeneity of the gray level of the image question.

Since the grayscale distribution of images between adjacent slices is very similar, and the tooth contour positions are also close, the present invention embeds the prior grayscale information into the energy model to improve the comparison accuracy, and its energy functional is defined as:

where f _{ref_in} (x) and f _{ref_out} (x) are defined as the mean gray values of the r neighborhood of the reference point on the prior shape inside and outside the prior shape curve, respectively. The reference point is determined according to the following method: after the segmentation of the previous slice is completed, the gray mean value of the r neighborhood of each point on the prior shape inside and outside the prior shape curve is calculated, and then in the current slice image, the first The point on the test shape that is closest to the target pixel point x of the current image is used as a reference point. The selection of r neighborhood radius depends on the resolution of the image, and generally should not be too large, so as not to include too many adjacent tooth areas.

(4) In order to ensure the shape and stability of the level set throughout the evolution process, add the signed distance to preserve energy:

and a curvilinear arc-length smoothing term:

In this way, based on the above energy information, the minimized energy functional of the single-phase mixed level set model for segmenting root slices can be expressed as:

By minimizing the energy equation, the evolution equation of the level set function Φ of the root slice can be obtained as:

Among them, δ(Φ) is the Dirac function.

Step 3. Segmentation model of the crown slice: In the slice image of the crown, the teeth are gradually connected together with small gaps. In order to effectively extract the boundary between two adjacent teeth, use the method described in step 2 for the initial slice. The single-phase mixed level set model of the tooth is segmented, and the divided level set function is combined into two coupled two-phase mixed level set functions according to the rule of every other tooth, as shown in Figure 4.

According to the principle of multiphase level set segmentation, two level set functions Φ ₁ and Φ ₂ can divide the image into four regions, as shown in Fig. 5 . For the phase {Φ ₁ >0, Φ ₂ >0}, it represents the overlapping area of adjacent teeth. In order to overcome the overlapping of Φ ₁ and Φ ₂ in the evolution process, the regional competition criterion is introduced in the energy functional. Let A ₁ and A ₂ be the internal area enclosed by the level set curves corresponding to the level set functions Φ ₁ and Φ ₂ respectively. According to the properties of the Heaviside function H(x), A ₁ and A ₂ can be expressed as:

A ₁ ＝area(Φ ₁ ≥0)＝∫ _Ω H(Φ ₁ )dx

A ₂ ＝area(Φ ₂ ≥0)＝∫ _Ω H(Φ ₂ )dx

In order for these two regions not to overlap, it should satisfy:

Therefore, the area penalty term can be defined as:

E _repulse ＝∫ _Ω H(Φ _l )H(Φ ₂ )dx

Minimizing the above formula is equivalent to avoiding the overlap of two teeth.

Adding the area penalty term to the level set energy functional, the energy equation of the two-phase mixed level set model of the dental crown can be obtained as:

E _{tooth crown} (Φ ₁ , Φ ₂ ) = E _{tooth root} (Φ ₁ ) + E _{tooth root} (Φ ₂ ) + βE _repulse

Among them, β is used to control the overlapping degree of the two level set functions, and the above energy functional is minimized, and the evolution equations of the two-phase level set functions Φ ₁ and Φ ₂ can be obtained as follows:

with

Among them, Φ _p1 and Φ _p2 represent the prior values of Φ ₁ and Φ ₂ respectively; ξ ₁ (x) and ξ ₂ (x) respectively represent the gradient direction detection function when Φ ₁ and Φ ₂ evolve, namely:

Step 4. Reconstruction of the three-dimensional tooth model: use the narrow-band method and the semi-implicit difference scheme to solve the evolution equation of the above level set function, extract the updated level set function at the pixel point of Φ=0, and divide all slice images into these The pixel points of the tooth contour are transformed into three-dimensional coordinates, and the Delaunay triangulation method is used for reconstruction to obtain the triangular mesh model of each tooth.

Effectiveness of the present invention can further illustrate by following experiment:

Take the oral CT image sequence of a patient with complete mandibular teeth as an example. Fig. 6 is a comparison diagram of tooth segmentation results using the present invention, Li model, and CV model for the initial slice image of the CT sequence. Fig. 6(a) is the initial outline polygon of each tooth drawn on the initial slice. It can be seen that the edge-based Li model relies too much on the image gradient information near the initial curve, but lacks gradient direction constraints, so there are many inaccurate tooth edge positioning in Figure 6(b). The region-based C-V model relies on global grayscale information, thus producing fusions between teeth that are closer together and weak boundaries between teeth and jaws in Fig. 6(c). Fig. 6(d) is the segmentation result of the present invention. Due to the use of the Flux model based on the gradient vector field, the weak boundary can be effectively detected, and the tooth boundary contour can be accurately positioned by adding gradient direction and prior shape constraints.

Fig. 7(a) and Fig. 7(b) are the segmentation results of the root slice image of the patient's second molar by the method of the present invention and Gao et al. respectively. It can be seen that the present invention has more accurate segmentation results, but Gao et al.’s method relies on the gray-scale probability inside and outside the curve to detect the contour. Since the gray scale of the pulp cavity is dark, the detected curve will be biased at the weak boundary of the tooth pulp cavity.

Figure 8 shows the segmentation results of some slices at the right first molar site. It can be seen that the segmentation results of basically all slice images can fit the tooth target contour well. However, in the slice of the 160th layer, the level set did not capture the groove profile inside the tooth, which was mainly affected by the size of the narrow band in the level set evolution equation. Due to the presence of grooves on the tooth surface, the contours of some sections were extremely concave. If the narrow band width is made smaller at this time, it may limit the evolution of the level set to the target contour. But this does not affect the final segmentation result. Since the target contour line is included in the level set curve, the final contour line can be extracted through the gray threshold. The right side of Fig. 8 shows the results of three-dimensional reconstruction of the segmented tooth profile of each slice. It can be seen that the reconstructed three-dimensional model of the present invention can correctly reflect the anatomical shape of the tooth, thereby providing effective data for the next step of computer-aided diagnosis. To enable doctors to formulate the best surgical plan according to the actual situation of the patient.

The above description is only a preferred example of the present invention, and is not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention within.

Claims (4)

1. A three-dimensional tooth modeling method based on mixed level sets, is characterized in that, comprises the steps: (1) Select one from the oral CT image sequence as the starting slice, and draw the outline of each tooth on the image to initialize the level set function; (2) The root layer slice below the initial slice uses the single-phase mixed level set model constructed by combining the prior shape constraint energy, the edge energy based on the Flux model, and the local area energy based on the prior gray level to segment the tooth contour; (3) The crown layer slice above the initial slice is segmented by the two-phase mixed level set model combined with regional competition constraints; (4) Convert the pixel points of the tooth contour after all the slice segmentations into three-dimensional coordinates, and use the Delaunay triangulation method to reconstruct to obtain the triangular mesh model of each tooth, (5) The energy functional of the single-phase mixed level set model for segmenting root slices described in step (2) is defined as prior shape-constrained energy, edge energy based on the Flux model, local area energy based on prior grayscale, etc. The weighted sum of: E _{tooth root} (Φ)＝μE _int (Φ)+γE _length (Φ)+αE _prior (Φ)+vE _edge (Φ)+λE _region (Φ) Among them, μ, γ, α, ν, λ are the weight coefficients of each energy item, Φ is the level set function of the current slice; (3a) The signed distance preserving energy E _int (Φ), which is used to ensure the stability in the evolution process of the level set, is defined as: <mrow><msub><mi>E</mi><mi>int</mi></msub><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mo>=</mo><msub><mo>&amp;Integral;</mo><mi>&amp;Omega;</mi></msub><mfrac><mn>1</mn><mn>2</mn></mfrac><msup><mrow><mo>(</mo><mo>|</mo><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi></mrow><mo>|</mo><mo>-</mo><mn>1</mn><mo>)</mo></mrow><mn>2</mn></msup><mi>d</mi><mi>x</mi></mrow> in, is the gradient operator, (3b) Curve arc length smoothing energy E _length (Φ), used to smooth the level set profile, defined as: <mrow><msub><mi>E</mi><mrow><mi>l</mi><mi>e</mi><mi>n</mi><mi>g</mi><mi>t</mi><mi>h</mi></mrow></msub><mo>=</mo><msub><mo>&amp;Integral;</mo><mi>&amp;Omega;</mi></msub><mo>|</mo><mo>&amp;dtri;</mo><mi>H</mi><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mo>|</mo><mi>d</mi><mi>x</mi></mrow> Among them, H(x) is the Heaviside function (3c) The prior shape constraint energy E _prior (Φ), used to control the shape of the level set, maps the tooth contour after each segmentation to the adjacent slice image, and constrains it as the prior shape of the evolution of the current level set function, Its energy functional is defined as: <mrow><msub><mi>E</mi><mrow><mi>p</mi><mi>r</mi><mi>i</mi><mi>o</mi><mi>r</mi></mrow></msub><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mo>=</mo><munder><mo>&amp;Integral;</mo><mi>&amp;Omega;</mi></munder><msup><mrow><mo>(</mo><mi>H</mi><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo><mo>-</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mi>p</mi></msub><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mi>d</mi><mi>x</mi></mrow> Where Φ _p is the level set function corresponding to the prior shape after the previous slice is segmented; (3d) The edge energy E _edge (Φ) based on the Flux model is used to detect the outer boundary contour of the tooth, and the angle information between the image gradient direction and the horizontal function gradient direction is embedded into the traditional Flux model, and its energy functional definition for: E _edge (Φ)＝- _∫Ω ξ(x)ΔI _σ (1-H(Φ)dx Where Δ is the Laplacian operator, I _σ represents the image after Gaussian smoothing, <mrow><mi>&amp;xi;</mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mn>1</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mi>&amp;sigma;</mi></msub><mo>&gt;</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mi>&amp;sigma;</mi></msub><mo>&amp;le;</mo><mn>0</mn></mrow></mtd></mtr></mtable></mfenced></mrow> is the gradient direction detection function, and represents the dot product; (3e) The local region energy E _region (Φ) based on the prior gray scale is used to overcome the problem of uneven gray scale of the image, and embed the prior gray scale information into the regional model, and its energy functional is defined as: <mrow><msub><mi>E</mi><mrow><mi>r</mi><mi>e</mi><mi>g</mi><mi>i</mi><mi>o</mi><mi>n</mi></mrow></msub><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mo>=</mo><munder><mrow><mo>&amp;Integral;</mo><mo>&amp;Integral;</mo></mrow><mi>&amp;Omega;</mi></munder><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>n</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mi>H</mi><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mi>d</mi><mi>x</mi><mo>+</mo><munder><mrow><mo>&amp;Integral;</mo><mo>&amp;Integral;</mo></mrow><mi>&amp;Omega;</mi></munder><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>o</mi><mi>u</mi><mi>t</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mrow><mo>(</mo><mn>1</mn><mo>-</mo><mi>H</mi><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo><mo>)</mo></mrow><mi>d</mi><mi>x</mi></mrow> Among them, f _{ref_in} (x) and f _{ref_out} (x) are respectively defined as the gray mean value of the r neighborhood of the reference point on the prior shape inside and outside the prior shape curve, and the reference point is the distance from the current image target pixel on the prior shape point x nearest point; (3f) Combining the above energy items, the minimized energy functional of the mixed level set model of the root slice is expressed as: The above energy functional is minimized to obtain the evolution equation of the level set curve: <mfenced open = "" close = ""><mtable><mtr><mtd><mrow><mfrac><mrow><mo>&amp;part;</mo><mi>&amp;Phi;</mi></mrow><mrow><mo>&amp;part;</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mi>&amp;mu;</mi><mo>&amp;lsqb;</mo><mi>&amp;Delta;</mi><mi>&amp;Phi;</mi><mo>-</mo><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi></mrow><mrow><mo>|</mo><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><mi>&amp;delta;</mi><mrow><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo></mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><mo>+</mo><mi>&amp;gamma;</mi><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi></mrow><mrow><mo>|</mo><mrow><mo>&amp;dtri;</mo><mi>&amp;Phi;</mi></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mn>2</mn><mi>&amp;alpha;</mi><mrow><mo>(</mo><mi>H</mi><mo>(</mo><mi>&amp;Phi;</mi><mo>)</mo><mo>-</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mi>p</mi></msub><mo>)</mo><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>v</mi><mi>&amp;xi;</mrow>mi><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>&amp;Delta;I</mi><mi>&amp;sigma;</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>&amp;lambda;</mi><mo>&amp;lsqb;</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>n</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>o</mi><mi>u</mi><mi>t</mi></mrow></msub><mo>(</mo><mi>x</mi>mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>&amp;rsqb;</mo></mrow></mtd></mtr></mtable></mfenced></mrow></mtd></mtr></mtable></mfenced> Among them, δ(Φ) is the Dirac function. 2. A kind of three-dimensional tooth modeling method based on mixed level set according to claim 1, it is characterized in that, step (1) is carried out as follows: in the slice image of tooth neck position, select a piece of all teeth The slice that appears and the alveolar bone rarely appears is used as the starting slice, and the rough outline C _i of each tooth is drawn on the slice image, i=1,2,...n,n is the number of teeth, as a level set , then use C _i to initialize n level set functions Φ _i , i=1, 2,...n, the initialization of Φ _i is done by calculating the signed distance from each point on the image to C _i , namely: <mrow><msub><mi>&amp;Phi;</mi><mi>i</mi></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mi>d</mi><mo>&amp;lsqb;</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>,</mo><msub><mi>C</mi><mi>i</mi></msub><mo>&amp;rsqb;</mo><mo>,</mo></mrow></mtd><mtd><mrow><mi>x</mi><mo>&amp;Element;</mo><mi>i</mi><mi>n</mi><mi>s</mi><mi>i</mi><mi>d</mi><mi>e</mi><mrow><mo>(</mo><msub><mi>C</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>x</mi><mo>&amp;Element;</mo><msub><mi>C</mi><mi>i</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>d</mi><mo>&amp;lsqb;</mo><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>,</mo><msub><mi>C</mi><mi>i</mi></msub><mo>&amp;rsqb;</mo><mo>,</mo></mrow></mtd><mtd><mrow><mi>x</mi><mo>&amp;Element;</mo><mi>o</mi><mi>u</mi><mi>t</mi><mi>s</mi><mi>i</mi><mi>d</mi><mi>e</mi><mrow><mo>(</mo><msub><mi>C</mi><mi>i</mi></msub><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced></mrow> Wherein, d[(x),C _i ] represents the Euclidean distance between the pixel point x and the curve C _i . 3. A kind of three-dimensional tooth modeling method based on mixed level set according to claim 1, it is characterized in that, the two-phase mixed level set model of the segmented tooth crown layer slice described in step (3) is established according to the following steps : The initial slice is segmented using the single-phase mixed level set model, and all the level set functions after segmentation are combined into two coupled two-phase mixed level set functions according to the rule of every other tooth, and added The region competes to constrain the energy to overcome the overlapping of two level sets, and its energy functional is defined as: E _{tooth crown} (Φ ₁ ,Φ ₂ )=E _{tooth root} (Φ ₁ )+E _{tooth root} (Φ ₂ )+βE _repulse in, E _repulse ＝∫ _Ω H(Φ ₁ )H(Φ ₂ )dx Constrain energy for regional competition, and β is used to control the overlapping degree of two level set functions; Minimize the above energy functional to obtain the evolution equations of the two-phase level set functions Φ ₁ and Φ ₂ respectively: <mfenced open = "" close = ""><mtable><mtr><mtd><mrow><mfrac><mrow><mo>&amp;part;</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub></mrow><mrow><mo>&amp;part;</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mi>&amp;mu;</mi><mo>&amp;lsqb;</mo><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>1</mn></msub><mo>-</mo><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>1</mn></msub></mrow><mrow><mo>|</mo><mrow><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>1</mn></msub></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><mi>&amp;delta;</mi><mrow><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub><mo>)</mo></mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><mo>+</mo><mi>&amp;gamma;</mi><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub></mrow><mrow><mo>|</mo><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mn>2</mn><mi>&amp;alpha;</mi><mrow><mo>(</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub><mo>)</mo><mo>-</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mrow><mi>p</mi><mn>1</mn></mrow></msub><mo>)</mo><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><msub><mi>v&amp;xi;</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>&amp;Delta;I</mi><mi>&amp;sigma;</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>&amp;lambda;</mi><mo>&amp;lsqb;</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>n</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>o</mi><mi>u</mi><mi>t</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>&amp;rsqb;</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>&amp;beta;</mi><mi>H</mi><mrow><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced></mrow></mtd></mtr></mtable></mfenced> <mfenced open = "" close = ""><mtable><mtr><mtd><mrow><mfrac><mrow><mo>&amp;part;</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub></mrow><mrow><mo>&amp;part;</mo><mi>t</mi></mrow></mfrac><mo>=</mo><mi>&amp;mu;</mi><mo>&amp;lsqb;</mo><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>2</mn></msub><mo>-</mo><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>2</mn></msub></mrow><mrow><mo>|</mo><mrow><msub><mi>&amp;Delta;&amp;Phi;</mi><mn>2</mn></msub></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow><mo>&amp;rsqb;</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>+</mo><mi>&amp;delta;</mi><mrow><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub><mo>)</mo></mrow><mfenced open = "[" close = "]"><mtable><mtr><mtd><mrow><mo>+</mo><mi>&amp;gamma;</mi><mi>d</mi><mi>i</mi><mi>v</mi><mrow><mo>(</mo><mfrac><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub></mrow><mrow><mo>|</mo><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub></mrow><mo>|</mo></mrow></mfrac><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mn>2</mn><mi>&amp;alpha;</mi><mrow><mo>(</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub><mo>)</mo><mo>-</mo><mi>H</mi><mo>(</mo><msub><mi>&amp;Phi;</mi><mrow><mi>p</mi><mn>2</mn></mrow></msub><mo>)</mo><mo>)</mo></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><msub><mi>v&amp;xi;</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><msub><mi>&amp;Delta;I</mi><mi>&amp;sigma;</mi></msub></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>&amp;lambda;</mi><mo>&amp;lsqb;</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>n</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>-</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msub><mi>f</mi><mrow><mi>r</mi><mi>e</mi><mi>f</mi><mo>_</mo><mi>i</mi><mi>o</mi><mi>u</mi><mi>t</mi></mrow></msub><mo>(</mo><mi>x</mi><mo>)</mo><mo>)</mo></mrow><mn>2</mn></msup><mo>&amp;rsqb;</mo></mrow></mtd></mtr><mtr><mtd><mrow><mo>-</mo><mi>&amp;beta;</mi><mi>H</mi><mrow><mo>(</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub><mo>)</mo></mrow></mrow></mtd></mtr></mtable></mfenced></mrow></mtd></mtr></mtable></mfenced> Among them, Φ _p1 and Φ _p2 represent the prior values of Φ ₁ and Φ ₂ respectively; ξ ₁ (x) and ξ ₂ (x) represent the gradient direction detection function when Φ ₁ and Φ ₂ evolve respectively: <mrow><msub><mi>&amp;xi;</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mn>1</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mrow><mn>1</mn><mi>&amp;sigma;</mi></mrow></msub><mo>&gt;</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>1</mn></msub><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mrow><mn>1</mn><mi>&amp;sigma;</mi></mrow></msub><mo>&amp;le;</mo><mn>0</mn></mrow></mtd></mtr></mtable></mfenced></mrow> <mrow><msub><mi>&amp;xi;</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>x</mi><mo>)</mo></mrow><mo>=</mo><mfenced open = "{" close = ""><mtable><mtr><mtd><mrow><mn>1</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mrow><mn>2</mn><mi>&amp;sigma;</mi></mrow></msub><mo>&gt;</mo><mn>0</mn></mrow></mtd></mtr><mtr><mtd><mrow><mn>0</mn><mo>,</mo></mrow></mtd><mtd><mrow><mi>i</mi><mi>f</mi></mrow></mtd><mtd><mrow><mo>&amp;dtri;</mo><msub><mi>&amp;Phi;</mi><mn>2</mn></msub><mo>&amp;CenterDot;</mo><mo>&amp;dtri;</mo><msub><mi>I</mi><mrow><mn>2</mn><mi>&amp;sigma;</mi></mrow></msub><mo>&amp;le;</mo><mn>0</mn></mrow></mtd></mtr></mtable></mfenced><mo>.</mo></mrow> 4. A kind of three-dimensional tooth modeling method based on mixed level set according to claim 1, it is characterized in that, the three-dimensional tooth reconstruction described in step (4) is carried out as follows: Utilize narrow-band method and semi-implicit difference scheme Solve the evolution equation of the above-mentioned level set function, extract the pixel point of the updated level set function at Φ=0, obtain the tooth contour of the current slice image, convert the tooth contour pixel points of all slice images into three-dimensional coordinates, use Reconstruction was performed based on the region-growing Delaunay triangulation method to generate a triangular mesh model of each tooth.