CN105488849A - Hybrid level set based three-dimensional tooth modeling method - Google Patents

Abstract

The invention provides a hybrid level set based three-dimensional tooth modeling method. The method comprises the following multiple steps: 1) selecting an initial slice and initializing a level set; 2) segmenting a tooth root slice image by utilizing a single-phase hybrid level set model constructed by combining priori shape constraint energy, Flux model based edge energy and priori gray based local region energy; 3) segmenting a tooth crown layer slice image by a dual-phase hybrid level set model in combination with regional competition constraint; and 4) reconstructing a contour segmentation based three-dimensional tooth model. According to the method, the hybrid level set model can effectively solve the problems of inaccuracy for edge positioning, incapability of processing an image, uneven gray level and the like in a conventional segmentation model, and an independent three-dimensional model of each tooth can be quickly and accurately constructed, so that a solid foundation is laid for making oral repair plans, biomechanical analysis and the like.

Description

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

Technical field

The invention belongs to image processing field, relate to the application of a kind of Level Set Method in three-dimensional tooth modeling.

Background technology

In recent years, along with three-dimensional digital imaging technique is in the develop rapidly in dentistry field, Computer aided decision reparation is more and more applied in the middle of oral restoration, and becomes the development trend in this field gradually.In the repair system of computing machine oral cavity, first important will obtain digitized three-dimensional tooth model exactly, and the precision of model and integrality are directly connected to the result of follow-up row's tooth, plantation, correction and biomechanical analysis.

The most frequently used method of current tooth modeling method utilizes image processing techniques to split tooth CT image sequence exactly, namely first from every one deck CT cuts into slices, is partitioned into tooth profile, then utilizes these interlayer contour reconstructions to go out tooth three-dimensional model.Because these class methods can obtain whole teeth patterning structure, for patient's pathological change of oral cavity provides complete diagnosis basis, therefore, the tooth modeling method based on CT image is more and more subject to the extensive concern of researcher.

Due to the density of oral cavity CT image Tooth and jawbone and distance all comparatively close, adopt traditional image partition method to be difficult to accurately extract the tissue contours of every tooth.The segmentation of CT image tooth is one always and is full of challenging problem.Wang Li etc. (Wang Li, Cui Jin, Han Qingkai etc., tie up solid model and set up by the tooth 3 based on CT image.Journal of Image and Graphics, 2005,10 (10): 1289-1292.) utilize binaryzation and Boundary Extraction to filter out every layer of section tooth profile key point, then utilize 3D-Delaunay tetrahedralization algorithm to obtain the solid model of whole tooth; But binaryzation processing ease produces the problem of over-segmentation or less divided.(the X.Wu such as Wu, H.Gao, H.Heo, etal.ImprovedB-splinecontourfittingusinggeneticalgorithm forthesegmentationofdentalcomputerizedtomographyimageseq uences.TheJournalofimagingscienceandtechnology, 2007,51 (4): 328-336.) utilize the tooth profile extracting every layer of section based on the B-spline curves matching of genetic algorithm, but B-spline curves cannot the problem of dental treatments topologies change.(the H.Gao such as Gao, O.Chae.IndividualtoothsegmentationfromCTimagesusinglevel setmethodwithshapeandintensityprior.PatternRecognition, 2010,43:2406-2417) movable contour model based on level set is adopted to carry out tooth segmentation, and different parted patterns is adopted to different dental layers section, the corona of every tooth and root of the tooth all can be split; This model mainly relies on the probability distribution of Image edge gradient and prior shape to carry out the evolution of Guidance Levels collection, but due to tooth Density inhomogeneity, and be easily subject to the isostructural interference of alveolar bone around, therefore rely on the gray probability of prior shape peripheral region to come contraction and the easy problem producing boundary leakage of expansion of level of control collection profile.

Summary of the invention

For overcoming above-mentioned the deficiencies in the prior art, improve efficiency and the precision of tooth segmentation, the present invention considers oral cavity CT image teeth patterning Variation Features, proposes a kind of three-dimensional tooth modeling method based on mixed-level collection.

The invention provides: a kind of three-dimensional tooth modeling method based on mixed-level collection, comprises the steps:

(1) from the CT image sequence of oral cavity, choose one as initial section, and sketch out every tooth profile on this image with initialization level set function;

(2) section of initial section is following root of the tooth layer utilizes prior shape bound energy, single-phase mixed-level collection model segmentation tooth profile based on the edge energy of Flux model, the structure that to combine based on the energy of local area of priori gray scale;

(3) the corona layer section that initial section is above utilizes calmodulin binding domain CaM to compete the bipolar mixture Level Set Models segmentation tooth profile of constraint;

(4) the tooth profile pixel after all section segmentations is converted into three-dimensional coordinate, utilizes Delaunay triangulation methodology to carry out rebuilding to obtain the triangle grid model of every tooth.

Step (1) is carried out as follows: choose in the sectioning image of tooth neck position all tooth all occur and the section of the less appearance of alveolar bone as initial section, and on this sectioning image, sketch out the general profile C of every tooth _i(i=1,2 ... n, n are tooth number) as the initial profile of level set, then utilize C _iinitialization n level set function Φ _i(i=1,2 ... n), Φ _iinitialization by point each in computed image to C _isigned distance come, that is:

${\mathrm{\Φ}}_i\left(x\right)=\left\{\begin{array}{cc}d\[\left(x\right),C_i\],& x\∈inside\left(C_i\right)\\ 0,& x\∈C_i\\ -d\[\left(x\right),C_i\],& x\∈outside\left(C_i\right)\end{array}\right.$

Wherein, d [(x), C _i] represent pixel x and curve C _ibetween Euclidean distance.

The energy functional of the single-phase mixed-level collection model of the segmentation root of the tooth layer section described in step (2) is defined as the weighted sum of prior shape bound energy, the edge energy based on Flux model, the energy of local area based on priori gray scale etc.:

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

Wherein, μ, γ, α, ν, λ are the weight coefficient of each energy term;

(3a) symbolic distance keeps ENERGY E _int(Φ), be used for ensureing the stability in level set movements process, be defined as:

$E_{\mathrm{int}}\left(\mathrm{\Φ}\right)={\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\Φ}|-1)}^{2}dx$

(3b) curve arc long smoothed energy E _length(Φ), be used for smooth level collection profile, be defined as:

$E_{\mathrm{length}}={\∫}_{\mathrm{\Ω}}|\▿H\left(\mathrm{\Φ}\right)|\mathrm{dx}$

(3c) prior shape bound energy E _prior(Φ), be used for the shape of level of control collection, the tooth profile after segmentation is at every turn mapped to contiguous slices image, and the prior shape developed as present level set function uses restraint, and its energy functional is defined as:

$E_{prior}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫}{\left(H\right(\mathrm{\Φ})-H({\mathrm{\Φ}}_p\left)\right)}^{2}dx$

Wherein Φ is the level set function of current slice, Φ _pfor level set function corresponding to rear prior shape has been split in a upper section, H (x) is Heaviside function;

(3d) based on the edge energy E of Flux model _edge(Φ), be used for detecting the outer boundary profile of tooth, be embedded into by the angle information between image gradient direction and level function gradient direction in the middle of traditional F lux model, its energy functional is defined as:

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

Wherein Δ is Laplacian operator, I _σrepresent the image after Gaussian smoothing,

$\mathrm{\ξ}\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}>0\\ 0,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}\≤0\end{array}\right.$

For gradient direction detection function, for gradient operator, represent dot product;

(3e) based on the energy of local area E of priori gray scale _region(Φ), be used for overcoming gradation of image problem of non-uniform, be embedded in the middle of regional model by priori half-tone information, its energy functional is defined as:

$E_{region}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_in}\right)}^{2}H\left(\mathrm{\Φ}\right)dx+\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_out}\right)}^2\left(1-H\left(\mathrm{\Φ}\right)\right)dx$

Wherein f _{ref_in}(x) and f _{ref_out}x r neighborhood that () is defined as reference point in prior shape is respectively at the inside and outside gray average of priori pattern curve, and reference point is the point that in prior shape, distance present image target pixel points x is nearest;

(3f) comprehensive above each energy term, is expressed as the minimization of energy functional of root of the tooth layer section mixed-level collection model:

Above-mentioned energy functional is minimized, obtains the EVOLUTION EQUATION of level set curve:

$\begin{array}{c}\frac{\∂\mathrm{\Φ}}{\∂t}=\mathrm{\μ}\[\mathrm{\Δ}\mathrm{\Φ}-div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\]\\ +\mathrm{\δ}\left(\mathrm{\Φ}\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\\ -2\mathrm{\α}\left(H\left(\mathrm{\Φ}\right)-H\left({\mathrm{\Φ}}_p\right)\right)\\ -v\mathrm{\ξ}\left(x\right)\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\end{array}\right]\end{array}$

Wherein, δ (Φ) is Dirac function.

The bipolar mixture Level Set Models of the segmentation corona layer section during step (3) is described is set up as follows: the single-phase mixed-level collection model described in utilizing initial section is split, all level set functions after segmentation are combined into the bipolar mixture level set function of two couplings according to the rule every a tooth, and add region-competitive bound energy and overlap to overcome two level sets, its energy functional is defined as:

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

Wherein, E _repulse=∫ _Ωh (Φ ₁) H (Φ ₂) dx

For region-competitive bound energy, β is used for the degree of control two level set function region overlap;

Above-mentioned energy functional is minimized, obtains two-phase level set function Φ ₁, Φ ₂eVOLUTION EQUATION be respectively: $\frac{\∂{\mathrm{\Φ}}_1}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_1-div\left(\frac{\▿{\mathrm{\Φ}}_1}{|\▿{\mathrm{\Φ}}_1|}\right)\]$

With $+\mathrm{\δ}\left({\mathrm{\Φ}}_1\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_1}{\left|\▿{\mathrm{\Φ}}_1\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_1\right)-H\left({\mathrm{\Φ}}_{p1}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_1\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_2\right)\end{array}\right]$

With $\begin{array}{c}\frac{\∂{\mathrm{\Φ}}_2}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_2-div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\]\\ +\mathrm{\δ}\left({\mathrm{\Φ}}_2\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_2\right)-H\left({\mathrm{\Φ}}_{p2}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_2\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_1\right)\end{array}\right]\end{array}$

Wherein Φ _p1, Φ _p2represent Φ respectively ₁, Φ ₂priori value; ξ ₁(x), ξ ₂x () represents Φ respectively ₁, Φ ₂gradient direction detection function during evolution:

${\mathrm{\ξ}}_1\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}\≤0\end{array}\right.$

${\mathrm{\ξ}}_2\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}\≤0\end{array}\right..$

Three-dimensional tooth described in step (4) is rebuild and is carried out as follows: the EVOLUTION EQUATION utilizing arrowband method and the above-mentioned level set function of half implicit difference Scheme Solving, extract the pixel of level set function in Φ=0 after having upgraded, obtain the tooth profile of current slice image, tooth profile pixel after all sectioning image segmentations is converted into three-dimensional coordinate, the Delaunay triangulation methodology increased based on region is utilized to rebuild, to generate the triangle grid model of every tooth.

Beneficial effect of the present invention: mixed-level collection models coupling many-sided information such as image border, regional area, priori that the present invention proposes, effectively to overcome conventional model edge local inaccurate and cannot process the problems such as gradation of image is uneven.The inventive method has less manual intervention, and there is good segmentation effect and higher accuracy rate, the three-dimensional tooth model reconstructed correctly can reflect the anatomic form of tooth, thus formulate oral restoration planning, biomechanical analysis etc. for next step and establish solid foundation, for the precision of raising oral restoration and efficiency significant.

Accompanying drawing explanation

Fig. 1 is tooth modeling technique process flow diagram of the present invention.

Fig. 2 is the level set function initialization schematic diagram of initial section.

The gradient direction schematic diagram that Fig. 3 (a) is tooth and surrounding interference structure.

The gradient direction schematic diagram that Fig. 3 (b) is level set function.

Fig. 4 is the two-phase level set function schematic diagram of corona layer sectioning image.

Fig. 5 is two-phase level-set segmentation area schematic.

Fig. 6 is that the present invention and existing method are to the segmentation effect figure of initial sectioning image.

To be the present invention and existing method to cut into slices the segmentation effect figure at position of grinding one's teeth in sleep to root of the tooth Fig. 7.

Fig. 8 is that the present invention grinds one's teeth in sleep the some section segmentation contours in position and three-dimensional reconstruction design sketch.

Embodiment

Below in conjunction with accompanying drawing, embodiments of the present invention is further illustrated:

As shown in Figure 1, specific implementation step of the present invention is as follows:

Step 1, the choosing and level set initialization of initial section: first read tooth CT image sequence, in the sectioning image of tooth neck position, then select all tooth all to occur and the section of the less appearance of alveolar bone as initial section.On initial sectioning image, sketch out the general profile C of every tooth _i(i=1,2 ... n, n are tooth number) as the initial profile of level set, as shown in Figure 2, then utilize C _i(i=1,2 ... n) initialization n level set function Φ _i(i=1,2 ... n).Φ _iinitialization by point each in computed image to C _isigned distance come, that is:

${\mathrm{\Φ}}_i\left(x\right)=\left\{\begin{array}{cc}d\[\left(x\right),C_i\],& x\∈inside\left(C_i\right)\\ 0,& x\∈C_i\\ -d\[\left(x\right),C_i\],& x\∈outside\left(C_i\right)\end{array}\right.$

Wherein, d [(x), C _i] represent some x and curve C _ibetween Euclidean distance.

Step 2, root of the tooth layer section segmentation: in root of the tooth tomographic image, tooth is vulnerable to alveolar bone and the isostructural more interference of pulp cavity.Therefore, a kind of new mixed-level collection tooth parted pattern is proposed.The energy functional of this model is primarily of part compositions such as prior shape bound energy, the edge energy based on Flux model, the energy of local areas based on priori gray scale.

(1) in priori shape constraining energy term, between contiguous slices, all relatively, therefore the tooth profile that a upper section has completed segmentation is mapped to current slice image, the prior shape as level set movements uses restraint for the shape of corresponding tooth and position.Φ is made to be the level set function of current slice, Φ _pfor level set function corresponding to rear prior shape has been split in a upper section, then prior shape bound term is defined as:

$E_{prior}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫}{\left(H\right(\mathrm{\Φ})-H({\mathrm{\Φ}}_p\left)\right)}^{2}dx$

In formula, H (x) is Heaviside function.

(2) in the edge energy based on Flux model, utilize the Flux model based on gradient vector field to come localizing objects border, utilize the angle between image gradient direction and level function gradient direction to overcome the interference problem of surrounding socket bone and pulp cavity simultaneously.

Due to image at boundary gradient and edge method vector in the same way, therefore Flux model comes positioning image edge by the edge integration of computed image gradient and curve method inner product of vectors, and the method effectively can solve the segmentation problem of weak boundary target, and Evolution Rates is fast.The maximization energy functional of Flux model is expressed as:

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

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

But, because the interference structure of Its pulp and outside is more, directly utilize above-mentioned model to split, level set can capture on the pulp cavity of Its pulp and the border such as the adjacent teeth of outside or alveolar bone simultaneously, and therefore gradient direction constraint joins in the middle of above-mentioned Flux model by the present invention.

In tooth CT image, tooth is high-brightness region, and background is low brightness area, therefore the image gradient direction of tooth and surrounding interference structure can be expressed as shown in Fig. 3 a.Due to level set function Φ defined herein just getting in level set inside, outside gets negative, then the gradient direction of Φ can be expressed as shown in Fig. 3 b.Can find out, if evolution curve will be made only to capture the outer boundary profile of tooth, then the gradient direction of image should be consistent with the gradient direction of level set function.Therefore, the functional that minimizes of edge energy of the present invention is defined as:

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

Wherein $\mathrm{\ξ}\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}>0\\ 0,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}\≤0\end{array}\right.$

For gradient direction detection function, for gradient operator, represent dot product.

(3) in the energy of local area in conjunction with priori gray scale, the local neighborhood gray average of prior shape is replaced overall average, priori half-tone information cannot be utilized to overcome conventional model and the uneven problem of gradation of image cannot be processed.

Because the intensity profile of image between contiguous slices is very identical, and tooth outline position is also close, and therefore priori half-tone information is embedded in energy model by the present invention, compares precision to improve, and its energy functional is defined as:

$E_{region}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_in}\right)}^{2}H\left(\mathrm{\Φ}\right)dx+\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_out}\right)}^2\left(1-H\left(\mathrm{\Φ}\right)\right)dx$

Wherein f _{ref_in}(x) and f _{ref_out}x r neighborhood that () is defined as reference point in prior shape is respectively at the inside and outside gray average of priori pattern curve.Reference point is determined as follows: after a upper section segmentation completes, the r neighborhood calculating each point in prior shape is at the inside and outside gray average of priori pattern curve, then in current slice image, by as a reference point for point nearest for distance present image target pixel points x in prior shape.Choosing of the r radius of neighbourhood need be determined according to the resolution of image, generally should be too not large, in order to avoid comprise too much adjacent teeth region.

(4) for ensureing the shape of level set in whole evolutionary process and stability, then add symbolic distance and keep energy:

$E_{\mathrm{int}}\left(\mathrm{\Φ}\right)={\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\Φ}|-1)}^{2}dx$

Item level and smooth with curve arc long:

$E_{\mathrm{length}}={\∫}_{\mathrm{\Ω}}|\▿H\left(\mathrm{\Φ}\right)|\mathrm{dx}.$

Like this, comprehensive each energy information above, can be expressed as the minimization of energy functional of the single-phase mixed-level collection model of segmentation root of the tooth layer section:

Minimize energy equation, the EVOLUTION EQUATION that can obtain the level set function Φ of root of the tooth layer section is:

$\begin{array}{c}\frac{\∂\mathrm{\Φ}}{\∂t}=\mathrm{\μ}\[\mathrm{\Δ}\mathrm{\Φ}-div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\]\\ +\mathrm{\δ}\left(\mathrm{\Φ}\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\\ -2\mathrm{\α}\left(H\left(\mathrm{\Φ}\right)-H\left({\mathrm{\Φ}}_p\right)\right)\\ -v\mathrm{\ξ}\left(x\right)\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\end{array}\right]\end{array}$

Wherein δ (Φ) is Dirac function.

Step 3, corona layer section parted pattern: in corona layer sectioning image, tooth is close to together gradually, gap is less, for effectively extracting the border between adjacent two teeth, the single-phase mixed-level collection model described in step 2 is utilized to carry out tooth segmentation to initial section, level set function after segmentation is combined into the bipolar mixture level set function of two couplings according to the rule every a tooth, as shown in Figure 4.

Principle is split, two level set function Φ from multi-phase horizontal set ₁, Φ ₂image can be divided into four regions, as shown in Figure 5.For { Φ ₁> 0, Φ ₂this phase of > 0}, represents the overlapping region of adjacent teeth.For overcoming Φ ₁, Φ ₂overlap in evolutionary process, in energy functional, introduce region-competitive criterion.If A ₁, A ₂be respectively level set function Φ ₁, Φ ₂the interior zone area that corresponding level set curve surrounds.According to the character of Heaviside function H (x), can by A ₁, A ₂be expressed as:

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

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

These two regions are not overlapped, then should meet:

Therefore region penalty term can be defined as:

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

Above formula is minimized and is equivalent to avoid two teeth to overlap.

Joined by region penalty term in the middle of level set energy functional, the energy equation that can obtain the bipolar mixture Level Set Models of corona layer is:

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

Wherein β is used for the degree of control two level set function region overlap, minimizes, can obtain two-phase level set function Φ to above-mentioned energy functional ₁, Φ ₂eVOLUTION EQUATION be respectively:

$\begin{array}{c}\frac{\∂{\mathrm{\Φ}}_1}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_1-div\left(\frac{\▿{\mathrm{\Φ}}_1}{|\▿{\mathrm{\Φ}}_1|}\right)\]\\ +\mathrm{\δ}\left({\mathrm{\Φ}}_1\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_1}{\left|\▿{\mathrm{\Φ}}_1\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_1\right)-H\left({\mathrm{\Φ}}_{p1}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_1\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_2\right)\end{array}\right]\end{array}$

With

$\begin{array}{c}\frac{\∂{\mathrm{\Φ}}_2}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_2-div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\]\\ +\mathrm{\δ}\left({\mathrm{\Φ}}_2\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_2\right)-H\left({\mathrm{\Φ}}_{p2}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_2\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_1\right)\end{array}\right]\end{array}$

Wherein Φ _p1, Φ _p2represent Φ respectively ₁, Φ ₂priori value; ξ ₁(x), ξ ₂x () represents Φ respectively ₁, Φ ₂gradient direction detection function during evolution, that is:

${\mathrm{\ξ}}_1\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}\≤0\end{array}\right.$

${\mathrm{\ξ}}_2\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}\≤0\end{array}\right.$

Step 4, tooth three-dimensional model are rebuild: the EVOLUTION EQUATION utilizing arrowband method and the above-mentioned level set function of half implicit difference Scheme Solving, extract the pixel of level set function in Φ=0 after having upgraded, these tooth profile pixels after all sectioning image segmentations are converted into three-dimensional coordinate, utilize Delaunay triangulation methodology to rebuild, the triangle grid model of every tooth can be obtained.

Validity of the present invention can be further illustrated by experiment below:

Using the complete patient oral cavity CT image sequence of mandibular teeth as example.Fig. 6 utilizes the present invention and Li model, the initial sectioning image of CV model to this CT sequence to carry out tooth segmentation result comparison diagram.Fig. 6 (a) each tooth initial profile polygon for drawing in initial section.Can find out, based on the Li model at edge owing to too relying on the image gradient information near initial curve, but lack gradient direction constraint, so there is the inaccurate phenomenon of a lot of tooth edge local in Fig. 6 (b).C-V model based on region relies on overall half-tone information, and the weak boundary place therefore between the nearer tooth of Fig. 6 (c) middle distance and between tooth and jawbone produces merges.Fig. 6 (d) is segmentation result of the present invention, owing to utilizing the Flux model based on gradient vector field, effectively can detect weak boundary, and by adding gradient direction and prior shape constraint, thus realize accurate positioning teeth boundary profile.

The method of people such as Fig. 7 (a) and Fig. 7 (b) difference the present invention and Gao etc. is to the segmentation result of the root of the tooth position sectioning image of this patient's second molar.Can find out, the present invention has segmentation result more accurately, but the method for the people such as Gao relies on inside or outside of curve gray probability to detect profile, because the gray scale of pulp cavity is comparatively dark, can make the curve deflection pulp cavity of detection at tooth weak boundary place.

Fig. 8 shows the segmentation result of some sections at first molar position, right side.Can find out, the segmentation result of all sectioning images can be fitted tooth objective contour well substantially.But in the section of the 160th layer, level set does not capture Its pulp groove contour, and this is mainly subject to the impact of arrowband size in level set movements equation.Due to the existence of dental surface groove, there is depression extremely in the profile that some can be caused to cut into slices.If now narrowband width obtains less, may will develop on objective contour by limit levels collection.But this does not affect last segmentation result, because target wheel profile is included in the inside of level set curve, extract final outline line by gray threshold.Show every layer of section tooth profile after to segmentation on the right of Fig. 8 and carry out the result after three-dimensional reconstruction, can find out that three-dimensional model that the present invention rebuilds correctly can reflect the anatomic form of tooth, thus provide valid data for next step computer-aided diagnosis, enable doctor formulate best operation plan for the actual conditions of patient.

The foregoing is only preferred embodiments of the present invention, not in order to limit the present invention, all any amendments done within the spirit and principles in the present invention, equivalent replacement and improvement etc., all should be included within protection scope of the present invention.

Claims (5)

1., based on a three-dimensional tooth modeling method for mixed-level collection, it is characterized in that, comprise the steps:

(1) from the CT image sequence of oral cavity, choose one as initial section, and sketch out every tooth profile on this image with initialization level set function;

(2) section of initial section is following root of the tooth layer utilizes prior shape bound energy, single-phase mixed-level collection model segmentation tooth profile based on the edge energy of Flux model, the structure that to combine based on the energy of local area of priori gray scale;

(3) the corona layer section that initial section is above utilizes calmodulin binding domain CaM to compete the bipolar mixture Level Set Models segmentation tooth profile of constraint;

(4) the tooth profile pixel after all section segmentations is converted into three-dimensional coordinate, utilizes Delaunay triangulation methodology to carry out rebuilding to obtain the triangle grid model of every tooth.

2. a kind of three-dimensional tooth modeling method based on mixed-level collection according to claim 1, it is characterized in that, step (1) is carried out as follows: choose in the sectioning image of tooth neck position all tooth all occur and the section of the less appearance of alveolar bone as initial section, and on this sectioning image, sketch out the general profile C of every tooth _i(i=1,2 ... n, n are tooth number) as the initial profile of level set, then utilize C _iinitialization n level set function Ф _i(i=1,2 ... n), Ф _iinitialization by point each in computed image to C _isigned distance come, that is:

${\mathrm{\Φ}}_i\left(x\right)=\left\{\begin{array}{cc}d\[\left(x\right),C_i\],& x\∈inside\left(C_i\right)\\ 0,& x\∈C_i\\ -d\[\left(x\right),C_i\],& x\∈outside\left(C_i\right)\end{array}\right.$

Wherein, d [(x), C _i] represent pixel x and curve C _ibetween Euclidean distance.

3. a kind of three-dimensional tooth modeling method based on mixed-level collection according to claim 1, it is characterized in that, the energy functional of the single-phase mixed-level collection model of the segmentation root of the tooth layer section described in step (2) is defined as the weighted sum of prior shape bound energy, the edge energy based on Flux model, the energy of local area based on priori gray scale etc.:

E _{root of the tooth}(Ф)=μ E _int(Ф)+γ E _length(Ф)+α E _prior(Ф)+ν E _edge(Ф)+λ E _region(Ф) wherein, μ, γ, α, ν, λ are the weight coefficient of each energy term;

(3a) symbolic distance keeps ENERGY E _int(Ф), be used for ensureing the stability in level set movements process, be defined as:

$E_{\mathrm{int}}\left(\mathrm{\Φ}\right)={\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿\mathrm{\Φ}|-1)}^{2}dx$

(3b) curve arc long smoothed energy E _length(Ф), be used for smooth level collection profile, be defined as:

$E_{length}={\∫}_{\mathrm{\Ω}}\left|\▿H\left(\mathrm{\Φ}\right)\right|dx$

(3c) prior shape bound energy E _prior(Ф), be used for the shape of level of control collection, the tooth profile after segmentation is at every turn mapped to contiguous slices image, and the prior shape developed as present level set function uses restraint, and its energy functional is defined as:

$E_{prior}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫}{\left(H\right(\mathrm{\Φ})-H({\mathrm{\Φ}}_p\left)\right)}^{2}dx$

Wherein Ф is the level set function of current slice, Ф _pfor level set function corresponding to rear prior shape has been split in a upper section, H (x) is Heaviside function;

(3d) based on the edge energy E of Flux model _edge(Ф), be used for detecting the outer boundary profile of tooth, be embedded into by the angle information between image gradient direction and level function gradient direction in the middle of traditional F lux model, its energy functional is defined as:

E _edge(Ф)＝-∫ _Ωξ△I _σ(1-H(Ф)dx

Wherein △ is Laplacian operator, I _σrepresent the image after Gaussian smoothing,

$\mathrm{\ξ}\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}>0\\ 0,& if& \▿\mathrm{\Φ}\·\▿I_{\mathrm{\σ}}\≤0\end{array}\right.$

For gradient direction detection function, for gradient operator, represent dot product;

(3e) based on the energy of local area E of priori gray scale _region(Ф), be used for overcoming gradation of image problem of non-uniform, be embedded in the middle of regional model by priori half-tone information, its energy functional is defined as:

$E_{region}\left(\mathrm{\Φ}\right)=\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_in}\right)}^{2}H\left(\mathrm{\Φ}\right)dx+\underset{\mathrm{\Ω}}{\∫\∫}{\left(I-f_{ref\_out}\right)}^2\left(1-H\left(\mathrm{\Φ}\right)\right)dx$

Wherein f _{ref_in}(x) and f _{ref_out}x r neighborhood that () is defined as reference point in prior shape is respectively at the inside and outside gray average of priori pattern curve, and reference point is the point that in prior shape, distance present image target pixel points x is nearest;

(3f) comprehensive above each energy term, is expressed as the minimization of energy functional of root of the tooth layer section mixed-level collection model:

Above-mentioned energy functional is minimized, obtains the EVOLUTION EQUATION of level set curve:

$\begin{array}{c}\frac{\∂\mathrm{\Φ}}{\∂t}=\mathrm{\μ}\[\mathrm{\Δ}\mathrm{\Φ}-div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\]\\ +\mathrm{\δ}\left(\mathrm{\Φ}\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿\mathrm{\Φ}}{\left|\▿\mathrm{\Φ}\right|}\right)\\ -2\mathrm{\α}\left(H\left(\mathrm{\Φ}\right)-H\left({\mathrm{\Φ}}_p\right)\right)\\ -v\mathrm{\ξ}\left(x\right)\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\end{array}\right]\end{array}$

Wherein, δ (Ф) is Dirac function.

4. a kind of three-dimensional tooth modeling method based on mixed-level collection according to claim 3, it is characterized in that, the bipolar mixture Level Set Models of the segmentation corona layer section during step (3) is described is set up as follows: the single-phase mixed-level collection model described in utilizing initial section is split, all level set functions after segmentation are combined into the bipolar mixture level set function of two couplings according to the rule every a tooth, and add region-competitive bound energy and overlap to overcome two level sets, its energy functional is defined as:

E _corona(Ф ₁, Ф ₂)=E _{root of the tooth}(Ф ₁)+E _{root of the tooth}(Ф ₂)+β E _repulse

Wherein,

E _repulse＝∫ _ΩE(Ф ₁)H(Ф ₂)dx

For region-competitive bound energy, β is used for the degree of control two level set function region overlap;

Above-mentioned energy functional is minimized, obtains two-phase level set function Ф ₁, Ф ₂eVOLUTION EQUATION be respectively:

$\frac{\∂{\mathrm{\Φ}}_1}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_1-div\left(\frac{\▿{\mathrm{\Φ}}_1}{|\▿{\mathrm{\Φ}}_1|}\right)\]$

With $+\mathrm{\δ}\left({\mathrm{\Φ}}_1\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_1}{\left|\▿{\mathrm{\Φ}}_1\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_1\right)-H\left({\mathrm{\Φ}}_{p1}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_1\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_2\right)\end{array}\right]$

With

$\begin{array}{c}\frac{\∂{\mathrm{\Φ}}_2}{\∂t}=\mathrm{\μ}\[{\mathrm{\Δ\Φ}}_2-div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\]\\ +\mathrm{\δ}\left({\mathrm{\Φ}}_2\right)\left[\begin{array}{c}+\mathrm{\γ}div\left(\frac{\▿{\mathrm{\Φ}}_2}{\left|\▿{\mathrm{\Φ}}_2\right|}\right)\\ -2{\mathrm{\α}}^{*}\left(H\left({\mathrm{\Φ}}_2\right)-H\left({\mathrm{\Φ}}_{p2}\right)\right)\\ -v^{*}{\mathrm{\ξ}}_2\mathrm{\Δ}I\\ -\mathrm{\λ}\[{\left(I-f_{ref\_in}\right)}^2-{\left(I-f_{ref\_iout}\right)}^2\]\\ -\mathrm{\β}H\left({\mathrm{\Φ}}_1\right)\end{array}\right]\end{array}$

Wherein Φ _p1, Φ _p2represent Ф respectively ₁, Ф ₂priori value; ξ ₁(x), ξ ₂x () represents Ф respectively ₁, Ф ₂gradient direction detection function during evolution:

${\mathrm{\ξ}}_1\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_1\·\▿I_{1\mathrm{\σ}}\≤0\end{array}\right.$

${\mathrm{\ξ}}_2\left(x\right)=\left\{\begin{array}{ccc}1,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}>0\\ 0,& if& \▿{\mathrm{\Φ}}_2\·\▿I_{2\mathrm{\σ}}\≤0\end{array}\right..$

5. a kind of three-dimensional tooth modeling method based on mixed-level collection according to claim 1, it is characterized in that, three-dimensional tooth described in step (4) is rebuild and is carried out as follows: the EVOLUTION EQUATION utilizing arrowband method and the above-mentioned level set function of half implicit difference Scheme Solving, extract the pixel of level set function in Ф=0 after having upgraded, obtain the tooth profile of current slice image, tooth profile pixel after all sectioning image segmentations is converted into three-dimensional coordinate, the Delaunay triangulation methodology increased based on region is utilized to rebuild, to generate the triangle grid model of every tooth.