CN105279762A - An oral cavity soft and hard tissue CT sequence and three-dimensional grid model registration method - Google Patents

Abstract

The invention provides an oral cavity soft and hard tissue CT sequence and three-dimensional grid model registration method comprising the following steps of: firstly manually dividing a biting impression model to obtain the part without a pedestal in the biting impression model; secondly, according to rough tooth outlines of user editing and calibrating, dividing input CT sequences to obtain the coordinates of the tooth outlines in the CT sequences, and constructing a hard tissue point set with the point coordinates of tooth hard tissue outlines of input layers of CT sequences; thirdly using iterative closest point algorithm for registration of the point set coordinates of the biting impression model and the point set coordinates of the hard tissue outlines to obtain the rough registration position in three-dimensional space of the biting impression model; fourthly, using the rough registration position obtained through the ICP algorithm as the initial result of accurate registration, constructing initial energy function values for the given CT sequences and the position of a three-dimensional grid model of initial registration, and optimizing the energy function through the quasi-Newton method to obtain an accurate registration result. The method enables superposition display of registration of CT sequences and biting impression models and can measure the thickness of soft tissue in CT sequences.

Description

Oral cavity soft and hard tissue CT sequence and three-dimensional grid model registration method

Technical Field

The invention relates to the field of image processing, computer vision and augmented reality, in particular to a registration method of an oral soft and hard tissue CT sequence and a three-dimensional grid model.

Background

Currently, some researchers have conducted research on related technologies. In image segmentation, in 2005, Chun-MingLi, Connecticut, proposed an improved Level-set (Level-set) method based on variational approach [5], which does not require re-initialization of the user-specified initial boundaries. This method has the following advantages over the conventional level set method: a larger timestamp value can be adopted to improve the convergence rate of the algorithm and reduce the calculation time; a generalized calibration area boundary can be used as an initial boundary without being limited by the fact that the boundary is necessarily a distance function; can be efficiently implemented using a finite difference method. Experimental results show that the method has a good segmentation effect even on relatively fuzzy object boundaries. However, the method has a low fault tolerance rate for the rough boundary initially calibrated manually by the user, and is prone to be trapped in the situation that the rough boundary cannot be converged to the correct boundary. In 2010, Chun-MingLi, university of Connecticut, improved [5] using a distance-uniformizing method [6 ]. In the improved level set method, the consistency of the level set function is maintained during its evolution, and the level set boundary with a value of zero is driven towards the correct object contour by the gradient flow direction of the distance-uniformizing term that minimizes the energy function. The distance uniformization term is defined as a potential function, so that the evolution direction of the level set has a unique FAB (ForwardBackWard) fusion effect, and the correctness of the level set boundary with the value of zero is guaranteed.

In 2010, OksamChae at the university of celebration, korea proposed a contour segmentation method for tooth CT sequences based on a level set method [7 ]. The method adopts single-level set segmentation and multi-level set segmentation to different parts (crown and root) of the tooth respectively, and overcomes the defect that the common horizontal method has poor segmentation effect on different objects with too close outline boundaries. Meanwhile, the method utilizes the front background gray and the back background gray near the outline boundary in the tooth CT slice as prior information, adds a region gray consistency item in an energy function of the traditional level set method, and realizes the evolution of the level set boundary with a zero driving value to a correct tooth boundary when the energy function is minimized. Compared with the traditional level set method, the method has certain improvement on the segmentation of the tooth CT sequence, but the calculation speed is too slow.

In terms of point set alignment, the iterative closest points (iterative closest points) algorithm is a registration algorithm for minimizing the difference between two point sets. The ICP algorithm continuously optimizes the similarity between the two point sets by iteratively calculating the matching relationship between the two point sets and minimizing the distance, thereby realizing the effect of point set registration.

2010, the German Kasikong institute of technology proposes an expert knowledge-based method for analyzing the augmented reality position in the oral implant operation. The method constructs a context-sensitive virtual reality system by acquiring local context information of the surgical instrument, and identifies the stage of the operation by tracking the positions of the human and the surgical instrument, thereby effectively visualizing the result of the virtual-real fusion. A schematic of this process is shown in figure 1.

An operation assisting system using augmented reality and virtual reality technology is introduced at the international society of ISMAR in 2011, munich, germany, and the application of the system in 100 clinical operation operations is analyzed in detail. The composition of the system is shown in figure 2.

In 2006, p.ljung, the university of snowman, sweden, proposed a virtual dissection method based on interactive large-scale, high-resolution CT scan data. The method adopts a world-leading volume rendering method to realize rendering of model data of upper G level on a GPU, and utilizes a data dimension reduction method based on a transformation function which adopts a multi-resolution rendering technology and a hierarchical detail selection method, as shown in figure 3.

Disclosure of Invention

In light of the above-mentioned practical needs and key issues, the present invention is directed to: the method for registering the oral soft and hard tissue CT sequence and the three-dimensional grid model can be used for registering the three-dimensional grid model into a local coordinate system of the CT sequence by utilizing the head CT sequence and the three-dimensional grid model of the human teeth.

The technical scheme adopted by the invention is as follows: editing the marked rough outline of the tooth by a user, and performing pixel segmentation on an input CT sequence by adopting a level set-based image segmentation method to obtain coordinates of the outline of the tooth in the CT sequence; registering point set coordinates of the input tooth three-dimensional grid model and the tooth outline point set coordinates obtained in the step by an ICP (inductively coupled plasma) algorithm in a local coordinate space of the CT sequence to obtain a rough position of the tooth three-dimensional grid model in the local coordinate space of the CT sequence; and performing iterative optimization on the tooth three-dimensional grid model by using the gray information of the tooth part in the CT sequence by using an energy function optimization method based on a quasi-Newton method to obtain the optimal position of the three-dimensional grid model in the CT sequence coordinates.

The rough outline of the calibrated tooth is edited by a user, and the polygon vertexes of the outline of the tooth are selected and determined mainly through visual recognition of the user for a plurality of pictures in a CT sequence. And (3) segmenting in the plane coordinate of each picture calibrated by the user to obtain the fine contour of the tooth by adopting an image segmentation method based on a level set method.

And (3) CT sequence point set registration, which is to obtain the three-dimensional coordinates of points on the tooth profile in a CT coordinate space by mainly utilizing the pixel interval information of the CT sequence and the fine contour of the tooth profile obtained in the image segmentation step, and perform point set registration with the three-dimensional mesh model of the tooth so as to transform the three-dimensional mesh model of the tooth to the approximately correct position in the CT coordinate space.

The fine registration of the three-dimensional grid model mainly utilizes the gray information of the CT sequence to construct a gray gradient matrix of soft and hard tissues of the CT sequence. According to the density difference of soft and hard tissues in the oral cavity and the characteristics of CT sequence scanning, the corresponding gray scale value near the surface of the hard tissue dental crown has a large change rate, namely the corresponding gray scale gradient value is the local maximum value. And continuously iterating and optimizing the position of the three-dimensional grid model by adopting a quasi-Newton method through the constraint that the gray gradient of the corresponding CT pixel coordinate in the CT sequence coordinate system where each vertex of the three-dimensional grid model is located is the maximum, thereby finally obtaining the accurate coordinate of the three-dimensional grid model in the CT sequence coordinate system.

Compared with the prior art, the invention has the beneficial effects that: 1, the invention completely utilizes the CT sequence data of the head and the corresponding bite model data without any other data; 2, the invention fully utilizes the gray information of the head CT sequence, and improves the registration precision on the basis of using the iterative closest point to register the algorithm; 3, iterative optimization is carried out by adopting a quasi-Newton method, and the optimal registration position can be rapidly calculated;

drawings

FIG. 1 is a method for analyzing augmented reality position in oral implant surgery based on expert knowledge, proposed by the acarbo institute of technology, Germany.

Fig. 2 shows a surgical assistance system using augmented reality and virtual reality technology, proposed by the university of munich, germany.

Fig. 3 a virtual dissection method of CT scan data proposed by the university of snowy forest, sweden.

FIG. 4 is a diagram illustrating a process of the present invention.

Fig. 5 is a cross-sectional view of the registration result of the present invention.

FIG. 6 is a flow chart of the architecture of the present invention;

fig. 7 is a flow chart of the registration algorithm of the present invention.

Detailed Description

The present invention will be described in detail with reference to the accompanying drawings.

See the energy function definition proposed by the present invention. The registration energy function of the oral soft and hard tissue CT sequence and the three-dimensional grid model is an energy function which is based on the density distribution of the CT sequence and the rotation and translation changes of the three-dimensional grid model in the space and is used for describing the registration degree of the CT sequence and the three-dimensional grid model.

The CT sequence A to be registered and three-dimensional grid data B are set, wherein the origin point coordinate of A in three-dimensional space is { x₀,y₀,z₀Pixel pitch of three directions is { sx }₀,sy₀,sz₀The coordinates { x, y, z } in A located in the volume will be mapped to the coordinates { x, y, z } in three-dimensional space₀+sx₀x,y₀+sy₀y,z₀+sz₀z }. Similarly, for any vertex in B, assuming its coordinates in three-dimensional space are { px, py, pz }, its volume coordinates mapped into the CT sequence are $\{\frac{\mathrm{px}-x_0}{{\mathrm{sx}}_0},\frac{\mathrm{py}-y_0}{{\mathrm{sy}}_0},\frac{\mathrm{pz}-z_0}{{\mathrm{sz}}_0}\}.$

For CT sequence, the density distribution function is set as D (x, y, z), where x, y and z are the volume coordinates of CT sequence, and the distribution function of density gradient can be obtainedAccording to the density distribution of the actual CT sequence, the part with lower density corresponds to the soft tissue, air and the like of the patient, and the part with higher density corresponds to the hard tissue of the teeth of the patient. The density of the CT sequence at the boundary between the soft and hard tissues of the patient will vary significantly, i.e., for points x, y, z on the contour of the hard tissue,die ofLocal maxima are taken.

provided with rotational-translational transformation { α, beta, gamma, t₁,t₂,t₃wherein α, β, γ represent rotation angles in three directions of x, y, z, respectively, and t₁,t₂,t₃Representing the translation in x, y, z directions, respectively, the corresponding rotation matrix is R ═ Rx × Ry × Rz, wherein $Rx=\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& -\sin \mathrm{\α}\\ 0& \sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right],Ry=\left[\begin{array}{ccc}\cos \mathrm{\β}& 0& -\sin \mathrm{\β}\\ 0& 1& 0\\ \sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right],Rz=\left[\begin{array}{ccc}\cos \mathrm{\γ}& -\sin \mathrm{\γ}& 0\\ \sin \mathrm{\γ}& \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right],$ Can obtain $R=\left[\begin{array}{ccc}\cos \mathrm{\β}\cos \mathrm{\γ}& -\cos \mathrm{\β}\sin \mathrm{\γ}& -\sin \mathrm{\β}\\ -\cos \mathrm{\γ}\sin \mathrm{\α}\sin \mathrm{\β}+\cos \mathrm{\α}\cos \mathrm{\γ}& \cos \mathrm{\α}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}& -\cos \mathrm{\β}\sin \mathrm{\α}\\ \cos \mathrm{\α}\cos \mathrm{\γ}\sin \mathrm{\β}+\sin \mathrm{\α}\sin \mathrm{\γ}& \cos \mathrm{\γ}\sin \mathrm{\α}-\cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}& \cos \mathrm{\α}\cos \mathrm{\β}\end{array}\right].$ The coordinates of the new point obtained by the rotation translation transformation of the point { x, y, z } on the three-dimensional mesh model are $\{x^{\′},y^{\′},z^{\′}\}=transpose\left(R*\left[\begin{array}{c}x\\ y\\ z\end{array}\right]+\left[\begin{array}{c}{t}_1\\ t_2\\ t_3\end{array}\right]\right),$ Calculating volume coordinates for new point mapping into CT sequenceThen its corresponding CT sequence density gradient amplitude isDeriving a definition of a registration energy functionthis function reflects the variables { α, β, γ, t in the rotational translation₁,t₂,t₃the registration degree of the three-dimensional grid model and the CT sequence is lower, the lower the value of E is, the better the registration result is, therefore, the parameters { α, beta, gamma, t of E are optimized by iteration₁,t₂,t₃can obtain { alpha, β, gamma, t }₁,t₂,t₃I.e. the best value for registration.

Referring to fig. 7, a flow chart of the registration method of the present invention is shown. According to definition of energy functionthe energy function is obtained with respect to { α, β, γ, t₁,t₂,t₃Expression of

$\begin{array}{c}|\∂D(t_1+\frac{\cos \mathrm{\β}\cos \mathrm{\γ}\left(-x_0+{\mathrm{vx}}^{\′}\right)}{{\mathrm{sx}}_0}-\frac{\cos \mathrm{\β}\sin \mathrm{\γ}\left(-y_0+{\mathrm{vy}}^{\′}\right)}{{\mathrm{sy}}_0}-\frac{\sin \mathrm{\β}\left(-z_0+{\mathrm{vz}}^{\′}\right)}{{\mathrm{sz}}_0},\\ t2+\frac{\left(-\cos \mathrm{\γ}\sin \mathrm{\α}\sin \mathrm{\β}+\cos \mathrm{\α}\sin \mathrm{\γ}\right)\left(-x_0+{\mathrm{vx}}^{\′}\right)}{{\mathrm{sx}}_0}+\frac{\left(\cos \mathrm{\α}\cos \mathrm{\γ}+\sin \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}\right)\left(-y_0+{\mathrm{vy}}^{\′}\right)}{{\mathrm{sy}}_0}-\frac{\cos \mathrm{\β}\sin \mathrm{\α}\left(-z_0+{\mathrm{vz}}^{\′}\right)}{{\mathrm{sz}}_0},\end{array}$ $t3+\frac{\left(\cos \mathrm{\α}\cos \mathrm{\γ}\sin \mathrm{\β}+\sin \mathrm{\α}\sin \mathrm{\γ}\right)\left(-x_0+{\mathrm{vx}}^{\′}\right)}{{\mathrm{sx}}_0}+\frac{\left(\cos \mathrm{\γ}\sin \mathrm{\α}-\cos \mathrm{\α}\sin \mathrm{\β}\sin \mathrm{\γ}\right)\left(-y_0+{\mathrm{vy}}^{\′}\right)}{{\mathrm{sy}}_0}+\frac{\cos \mathrm{\α}\cos \mathrm{\β}\left(-z_0+{\mathrm{vz}}^{\′}\right)}{{\mathrm{sz}}_0}\left)\right|$ . Let the current search position be x_n＝{α_n,β_n,γ_n,t_1n,t_2n,t_3nGet x by definition of E_nPartial derivatives of the components ofAnd E is with respect to x_nJacobian of_nWhile E can be determined with respect to x_nHessian matrix H_n. The updated formula of the line search is x_n+1＝x_n-a_nd_nWherein d is_nIs the direction of the line search, anda_nis the step size of the search.

When the BFGS quasi-Newton method is adopted to carry out line search, an iterative method is adopted to approximate the valueApproximate substitutionHas a recurrence formula ofWherein y is_n+1＝x_n+1-x_n，s_n+1＝J_n+1-J_n. To this end, for x_nThe optimized line search direction of (2) has already been determined, the remaining problem is to determine the step size of the line search. An appropriate search step size may be determined according to the StrongWolfe-Powell criterion such that the step size is neither too large to result in a deviation from the optimal value, nor too small to trap in a local extremum and reduce the search efficiency.

Referring to fig. 6, a step diagram of the present invention is shown.

Step one, displaying the pretreatment operation to be carried out on the mold-biting model. In practical application, the bite model is segmented by adopting a manual model segmentation method to obtain a model without a base, so that the interference of noise data of a non-human tissue structure is eliminated by separating from original data, and the registration accuracy is improved.

And step two, calibrating the tooth profile of each layer of the CT sequence by a user. And drawing the boundary by reading the polygon vertex of the tooth profile in the CT sequence selected point by the user.

And thirdly, accurately segmenting the CT sequence by adopting a level set method and utilizing the rough calibration result of the user on the CT sequence in the last step to obtain the fine contour of the tooth in each CT sequence. According to the predecessor method, the level set function should be kept as a directed distance function, i.e. within a given scale space (Ω, d), given a function f and a curved surface Ω in space_cThe following expression is satisfied: $f\left(x\right)=\left\{\begin{array}{cc}d(x,{\mathrm{\Ω}}^c)& x\∈\mathrm{\Ω}\\ -d(x,\mathrm{\Ω})& x\∈{\mathrm{\Ω}}^c\end{array}\right\},$ wherein, $d(x,\mathrm{\Ω})=\underset{y\∈\mathrm{\Ω}}{inf}d(x,y),$ d (x, Ω) is the distance from x to the curved surface Ω.

To ensure that f is a directed distance function, defineMeasure the degree of difference between f and the real directional distance function, whereinIs the laplacian operator. For the image segmentation problem, the following boundary detection function is defined:wherein G is_σRepresenting a gaussian kernel with sigma as the standard deviation, and I representing the CT image to be processed. Then defining an energy function $E\left(f\right)=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(\right|\▿f|-1)}^{2}dxdy+\mathrm{\λ}{\∫}_{\mathrm{\Ω}}g\mathrm{\σ}\left(f\right)|\▿f|dxdy+\mathrm{\ν}{\∫}_{\mathrm{\Ω}}gH\left(-f\right)dxdy,$ Wherein σ is a dirac function, H is a Heaviside function, and μ, λ, and ν are parameters of each item. By aiming at thisThe gradient of f is calculated from the energy function and takenAnd continuously iterating and updating until f converges as the updating direction of f, wherein the zero horizontal line of the obtained function f is the boundary obtained by segmentation. And combining the pixel spacing information of the CT sequence to obtain the three-dimensional coordinates of the tooth profile in the CT sequence coordinate space.

And step four, adopting an iterative closest point algorithm to finely divide the tooth profile of each layer of the CT sequence. The iterative closest point algorithm, ICP, is an abbreviation for iterative closestpoint, used for registration of a set of points with a set of points. And aligning the three-dimensional grid model to the CT sequence coordinate space by using the known point set coordinates of the tooth three-dimensional grid model and adopting an iterative closest point algorithm. Since the point set contained in the tooth profile of each layer of the CT sequence obtained by the level set segmentation method in the previous step is not the same as the point set coordinates of the three-dimensional mesh model, only the rough position of the three-dimensional mesh model can be obtained in this step.

Step five, constructing an energy equation of the three-dimensional grid model point set corresponding to the sum of the gradient values of the CT sequence by utilizing the gray information of the CT sequence, and obtaining a result R by iterating the nearest point algorithm in the step₀,t₀And as an initial result of registration, performing iterative optimization on the values of R and t by using a quasi-Newton method to obtain a final registration position.

The above description is only a few basic descriptions of the present invention, and any equivalent changes made according to the technical solutions of the present invention should fall within the protection scope of the present invention.

Claims (5)

1. A registration method of an oral soft and hard tissue CT sequence and a three-dimensional grid model is characterized by comprising the following steps:

1) manually dividing the mold biting model to obtain a part without a base in the mold biting model;

2) the method comprises the steps of utilizing a user to edit a calibrated tooth rough contour, segmenting an input CT sequence, obtaining coordinates of the tooth contour in the CT sequence, and constructing a hard tissue point set by using point coordinates of the tooth hard tissue contour in each layer of the input CT sequence;

3) registering the point set coordinates of the bite model and the hard tissue contour point set coordinates by adopting an iterative closest point algorithm to obtain a rough registration position of the bite model in a three-dimensional space;

4) and adopting a rough registration position obtained by an ICP (inductively coupled plasma) algorithm as an initial result of accurate registration, constructing an initial energy function value aiming at a given CT (computed tomography) sequence and the position of the initially registered three-dimensional grid model, and optimizing the energy function by adopting a quasi-Newton method to obtain an accurate registration result.

2. the method of claim 1, wherein the energy function is { α, β, γ, t₁,t₂,t₃describing the rotation translation transformation of the three-dimensional grid model to be registered to the accurate registration position by six variables, and adopting the rotation translation transformation of the three-dimensional grid model { α, beta, gamma, t on the basis of the prior condition that the gray gradient amplitude of the CT sequence is maximum at the hard tissue boundary₁,t₂,t₃The opposite number of the sum of the gray gradient amplitudes of the CT sequence corresponding to the coordinates of all the points under the points is defined as an energy function.

3. the method for registering an oral soft and hard tissue CT sequence with a three-dimensional grid model according to claim 2, wherein the parameters { α, β, γ, t) are optimized by using a quasi-Newton method in the step 4)₁,t₂,t₃the rotation angle of the bite-mold three-dimensional model in the three directions of x, y and z is { alpha, β gamma }, and the rotation angle of the bite-mold three-dimensional model in the three directions of x, y and z is { t }₁,t₂,t₃And determining the search direction by using an approximate value of a hessian matrix of the variables according to a BFGS method, and determining the value of the search step length based on a strong Wolfe-Powell condition.

4. The method for registering an oral soft and hard tissue CT sequence and a three-dimensional mesh model according to claim 1, characterized in that: further comprising the steps of: the step 2 further comprises:

1) by utilizing an interactive oral soft and hard tissue segmentation tool, a user can sequentially select rough boundary vertexes of hard tissue outlines in a current CT sequence slice;

2) selecting a slice with an obvious hard tissue contour in a CT sequence as an initial frame of marking, and marking a plurality of manually selected rough hard tissue contours of the obvious teeth by using the segmentation tool;

3) accurately dividing the rough hard tissue contour of the CT serial slice by adopting a level set method to obtain the precise contour of the rough hard tissue contour, outwards expanding a new contour obtained by a pixel by the precise contour to be used as the rough contour of the next adjacent CT serial slice, also accurately dividing by adopting the level set method, and dividing by adopting the same method to obtain the hard tissue contours of a plurality of CT serial slices;

4) and (3) constructing the hard tissue contour data of a plurality of slices of the CT sequence obtained in the previous step to obtain a point set of the corresponding hard tissue contour in a three-dimensional space.

5. The method for registering an oral soft and hard tissue CT sequence with a three-dimensional mesh model according to claim 3, wherein the method comprises the following steps: and 4) adopting an ICP (inductively coupled plasma) algorithm, taking the point set of the bite model as an input point set to be registered, taking the three-dimensional point set of the oral cavity hard tissue outline obtained in the step 3 as a reference registration point set, calculating the distance of each point of the target point set from each point on the source point set, matching each point with the nearest point of the target point set, and repeating the steps until the mean square error is smaller than a specified threshold value, so as to obtain the approximate registration position of the bite model in the space.