CBCT-based virtual articulator system and construction method

Abstract

The present invention belongs to the field of articulator technology and specifically discloses a CBCT-based virtual articulator system and construction method. The construction method includes the following steps: acquiring CBCT image projection information from the human oral cavity and reconstructing CT slice data using the FDK algorithm; reconstructing maxillary and mandibular mesh models based on the CT slice data using the MC algorithm; performing noise reduction, coarse registration, and fine registration on the maxillary and mandibular mesh models to optimize the maxillary and mandibular mesh models; extracting the hinge axis from the optimized maxillary and mandibular mesh models; and using the hinge axis to assist in mandibular motion simulation based on the rigid body transformation of the mandibular motion. This technical solution accurately simulates the mandibular motion trajectory, extracts the hinge axis, and simplifies operation.

Description

Technical Field

[0001] The present invention belongs to the technical field of articulators and relates to a CBCT-based virtual articulator system and a construction method. Background Art

[0002] The temporomandibular joint (TMJ) is an indirect connection between the mandible and the skull. It is the only movable joint in the maxillofacial region and one of the most complex joints in the human face. Together with the masticatory muscles and tooth occlusion, it affects the trajectory of mandibular movement. Only through the coordination of the three can basic oral functions such as chewing, swallowing, speech, and facial expressions be accomplished. Oral treatments (such as denture restorations, orthodontic treatments, and bite splint treatments) often cause changes in tooth occlusal contact. For example, if mandibular movement and tooth occlusion become misaligned during treatment, this can lead to tooth wear, damage to restorations, and even temporomandibular joint disorders. Therefore, accurate preoperative recording and simulation of mandibular movement trajectory is a topic of ongoing research for dentists and is of great significance to clinical diagnosis and treatment.

[0003] Mechanical articulators have been developed for over 200 years. Equipped with joints and components similar to the human masticatory apparatus, they are commonly used to simulate mandibular movement. The hinge axis plays a crucial role in mechanical articulators. The hinge axis is the rotational axis that occurs during mandibular mouth opening and closing movements, located near the line connecting the centers of gravity of the condyles. The hinge axis effectively assists the articulator in simulating mandibular movements, such as mouth opening and closing, and protrusion. Current methods for extracting the hinge axis typically involve a doctor installing a mechanical device on the patient's face, and then drawing a dot near the ear while the patient makes small mouth opening and closing movements. With the continuous advancement of articulator technology and techniques, traditional mechanical articulators can simulate mandibular movement relatively accurately, meeting the needs of routine laboratory and clinical work. However, traditional mechanical articulators still suffer from the drawbacks of complex and cumbersome clinical operation and difficulty in intuitive and three-dimensionally simulating mandibular movement. Furthermore, patents for mechanical articulators are held in Europe and the United States, making them expensive to purchase domestically. Summary of the Invention

[0004] The purpose of the present invention is to provide a CBCT-based virtual articulator system and construction method to solve the problems that traditional mechanical articulators are complex and tedious to operate clinically and difficult to simulate mandibular movements intuitively and three-dimensionally.

[0005] In order to achieve the above-mentioned object, the basic scheme of the present invention is: a method for constructing a virtual articulator system based on CBCT, comprising the following steps:

[0006] Collect CBCT image projection information inside the human oral cavity and reconstruct CT tomographic data using the FDK algorithm;

[0007] Based on CT tomographic data, the maxillary and mandibular mesh models were reconstructed using the MC algorithm.

[0008] Perform noise removal, coarse registration, and fine registration on the maxillary and mandibular mesh models to optimize them.

[0009] Extract the hinge axis using the optimized maxillary and mandibular mesh models;

[0010] In view of the expression of mandibular movement, the hinge axis is used to assist in simulating mandibular movement from the perspective of rigid body transformation.

[0011] The working principle and beneficial effects of this basic solution are as follows: It scans the oral cavity and reconstructs a realistic 3D model of the mandible. Based on this model, the hinge axis is extracted based on the mandibular motion characteristics. Finally, by controlling the hinge axis, mandibular motion is simulated, avoiding the cumbersome operation of traditional mechanical articulators.

[0012] Furthermore, the method of reconstructing CT slice data using the FDK algorithm is as follows:

[0013] According to the three-dimensional reconstructed geometric coordinate relationship of the scanning trajectory of the CBCT image projection information, three different coordinate systems are obtained, namely the scanning object coordinate system, the detector projection data coordinate system, and the ray source coordinate system, which correspond to o-xyz, o-tz, and o-st respectively, where z is the rotation center axis, s passes through the center of the ray source and is perpendicular to the detector plane;

[0014] The formula of the FDK algorithm in the center plane is:

[0015]

[0016] u=ycosβ-xsinβ, v=xcosβ+ysinβ,

[0017] Here, p(β,a,b) represents the projection data of point X0(x0,y0,z0) in the projection direction β.

[0018]

[0019] β is the projection angle, (a, b) is the detector array coordinate, R(β) is the orbit radius, g(τ) is the filter function, c is the complex unit, ρ is the polar diameter; u, v, τ are all intermediate codes; y, x, z are the coordinates of the corresponding coordinate axes.

[0020] The FDK reconstruction algorithm is a reconstruction algorithm specially designed for the circular focus orbit of the cone beam. Its biggest advantage is that the reconstructed image effect is better when the cone angle is less than 10°.

[0021] Furthermore, the method of reconstructing the maxillary and mandibular mesh models using the MC algorithm is as follows:

[0022] Divide the upper and lower jaw area data by using a cuboid to enclose the upper and lower jaw areas, and divide the cuboid into m×n×k cubes, with the upper and lower layers of each cube covered by different faults;

[0023] Number the vertices of each cube as v = (v1, v2, v3, v4, v5, v6, v7, v8) and the edges of the cube as e = (e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12);

[0024] Based on the tissue density on the section, we calculate whether each vertex of the cube is inside or outside the upper and lower jaw areas. The cube has 8 vertices in total, which requires 8 bits to represent. There are a total of 256 possible situations.

[0025] According to the internal and external situations of the 8 vertices of the cube on the upper and lower jaws, the interpolation formula is used to calculate the positions on the edges that need to be triangulated:

[0026]

[0027] Among them, isovalue is the upper and lower tissue density thresholds, P1 and P2 represent the coordinates of the two endpoints on the edge, and V1 and V2 represent the values at these two endpoints;

[0028] Calculate the normal vector of the triangle, the three points A, B, C of the triangle, the normal vector Traverse all cubes to obtain the Mesh models of the upper and lower jaws.

[0029] The FDK reconstruction algorithm was used to obtain a 3D image of the brain's internal structure. To obtain the upper and lower jaw models, the Marching Cubes algorithm was used. The Marching Cubes algorithm extracts isosurfaces from a 3D discrete data field, primarily used to extract mesh information from objects.

[0030] Furthermore, the method for denoising the upper and lower jaw Mesh models is as follows:

[0031] The point cloud data of the upper and lower jaw Mesh models are input into the VoxelGrid voxel filter;

[0032] Set up a large cube that encloses the point cloud data. The length, width, and height of the large cube are Lx, Ly, and Lz respectively. Divide the large cube into several small cubes. Let the total number of points in the point cloud be N. Then the side length L of the small cube is:

[0033]

[0034] Among them, α is the scaling factor used to adjust the side length of the small cube, and s is the scaling coefficient;

[0035] By adjusting the parameters α and s, the number of small cubes is obtained using the rounding function, and the point cloud data is divided into m*n*k small cubes;

[0036] Traverse the point cloud position information in each small cube and calculate the centroid coordinates of the point cloud in the small cube:

[0037]

[0038] Among them, X0, Y0, Z0 are the coordinates of the center of mass, x i ,y i , z i are the coordinates of the point cloud in the small cube, is the number of point clouds falling in the small cube, i is the sequence number, i = 1, 2, ..., represents the i-th point cloud in the small cube;

[0039] By downsampling the voxel grid, the number of points in the point cloud data is reduced, and the shape characteristics of the point cloud do not change.

[0040] Use the VoxelGrid voxel filter to downsample the model, reduce the number of point clouds, remove noise, and optimize the point cloud accuracy.

[0041] Furthermore, the method for coarse registration of the maxillary and mandibular mesh models is as follows:

[0042] Select two groups of corresponding points P in the maxillary x , P y , construct the objective function and solve the rotation matrix R:

[0043] P y =RP x +t

[0044] Where t is the translation matrix;

[0045] Ask for P x , P y The centroid coordinates of :

[0046]

[0047] It's P x The coordinates of the midpoint, It's P y Coordinates of the midpoint;

[0048] Move the centroid of the two groups of points to the center of the circle. The coordinates of the center of the circle are:

[0049]

[0050]

[0051] Solve the rotation matrix:

[0052] p′ y =Rp′ x

[0053]

[0054] R is the rotation matrix, or it is obtained by scaling transformation. Scaling transformation will cause shape deformation when applied to the maxillary point cloud. Singular value decomposition is used to solve it and construct the covariance matrix H:

[0055]

[0056] Perform singular value decomposition on H:

[0057] H=U∑V T

[0058] all:

[0059] R=UV T

[0060] Where U is an m-order orthogonal matrix, V is an n-order orthogonal matrix, ∑ is an mxn rectangular diagonal matrix, and T represents a transpose operation.

[0061] The ICP algorithm is used to perform point cloud registration on the upper and lower jaws with different openings. Registration is divided into coarse registration and fine registration. Coarse registration is mainly used to provide a good initial position for fine registration.

[0062] Furthermore, the method for fine registration of the maxillary and mandibular mesh models is as follows:

[0063] Using the initial R0, t0 of the coarse registration or the R obtained in the previous iteration k-1 , t k-1 Transform the point cloud data of the mandibular Mesh model to obtain a temporary transformed point cloud. Let the temporary point cloud be P s , the target point cloud is P t ; Where R is the rotation matrix, t is the translation matrix, subscript 0 indicates the first rotation or translation operation, and subscript k-1 indicates the k-1th rotation or translation operation;

[0064] P s and P t Compare and find P s Every point in P t The nearest neighbor points in are aligned;

[0065] Determine the objective function:

[0066]

[0067] Among them, R * represents the rotation matrix, T * represents the translation matrix, argmin represents the variable value when the objective function f(x) takes the minimum value, P s For temporary point cloud, It's P t and P s The corresponding points in , t is the translation matrix, R is the rotation matrix obtained by coarse registration; i is the sequence number, i = 1, 2, ..., representing the i-th point cloud in the small cube;

[0068] Use SVD decomposition method to solve R k and t k ; R k and t k Where R is the rotation matrix, t is the translation matrix, and the subscript k represents the kth rotation or translation operation;

[0069] Calculate the optimal translation:

[0070] Let N = |P s |, N is the total number of points in the point cloud, and the loss function is:

[0071]

[0072] Taking its derivative, we get:

[0073]

[0074] Let the derivative be 0, then:

[0075]

[0076] make We can get:

[0077]

[0078] When the value of R is determined, the optimal t is solved according to the formula to minimize the loss; P t and P s The center of mass;

[0079] Calculate the optimal rotation:

[0080] Convert the temporary point cloud and target point cloud to centroid coordinates:

[0081]

[0082]

[0083] Only the rotation problem needs to be considered, and the loss can be written as:

[0084]

[0085] in: is the temporary point cloud in the centroid coordinates, is the target point cloud in the centroid coordinates,

[0086] T is the transpose operation,

[0087] Simplified to:

[0088]

[0089] Loss is rewritten as:

[0090]

[0091] From the definition of matrix multiplication and matrix trace, we can get:

[0092]

[0093] The problem is transformed into:

[0094]

[0095] make H is the covariance matrix,

[0096] Then we have:

[0097]

[0098] Perform singular value decomposition on H = U∑V T

[0099]

[0100] Among them, trace means finding the trace of the matrix, V and U are orthogonal matrices;

[0101] Since V, U, and R are all orthogonal matrices, let the transformation matrix M be:

[0102] M=V T RU

[0103] Then we have:

[0104] trace(∑V T RU)=trace(∑M)

[0105] According to the non-negative nature of singular values and the properties of orthogonal matrices: the absolute value of the elements in the orthogonal matrix is not greater than 1, and trace(∑M) is maximum only when M is the unit matrix, that is:

[0106] V T RU=E

[0107] R=V T U

[0108] The upper jaws with different openings are registered to obtain the registration matrix M, which is then applied to the mandible to unify the mandibles with different openings into one coordinate system.

[0109] The upper jaw is aligned by taking advantage of the characteristic that the upper jaw remains stationary during the opening and closing movement. Then the alignment matrix M is applied to the lower jaw, and the lower jaw in different opening states is unified into one coordinate.

[0110] Furthermore, the KD tree method is used to find the nearest neighbor points of the temporary point cloud and the target point cloud. The steps to construct the KD tree are as follows:

[0111] S1, select the x dimension, and put P t The points in the tree are selected using a quick selection algorithm in the x-dimension to find the median, which is then divided into two parts: the left and right parts. The median point is used as the root node, and the points on the left and right sides are divided into two subtrees.

[0112] S2: For the points in these two parts, change the dimension to the y-axis, use the quick selection algorithm again to select the median, and use the median point as the child node of the previous layer node. At the same time, each node is divided into two parts;

[0113] S3, change the dimension to the z-axis, repeat S2, again construct the child node according to the median of the z-dimension and divide the point into left and right parts;

[0114] Repeat S1-S3 until there are no nodes in the left and right parts.

[0115] In the process of fine registration, KD tree is used for optimization to find the nearest neighbor points of the source point cloud and the target point cloud.

[0116] Furthermore, the method for extracting the hinge axis is as follows:

[0117] Use the SIFT operator to extract key points to replace the original point cloud of the upper and lower jaw Mesh model:

[0118] The coordinates of each point in the original mandibular point cloud are expressed as p(x, y, z), and the scale space of the 3D point cloud is defined as L(x, y, z, σ), which is obtained by convolving a Gaussian kernel function G(x, y, z, σ) with the original mandibular point cloud p(x, y, z):

[0119] L(x,y,z,σ)=G(x,y,z,σ)*p(x,y,z)

[0120] in, σ is a constant;

[0121] Use the multiplication factor k to obtain different scales and construct the point cloud Gaussian pyramid. If the number of layers in each pyramid group is s, then set k s =2, use the more efficient DoG operator instead of the scale-normalized Gaussian Laplace operator to perform extreme value detection. The DoG operator calculation formula is:

[0122] D(x,y,z,k e σ)=L(x,y,z,k e+1 σ)-L(x,y,z,k e σ)

[0123] Where, e∈[0,s+2]; obtain extreme points, set a threshold, and use extreme points that are not less than the threshold as candidate key points;

[0124] By calculating the points (x j ,y j ,z j ) to the neighborhood center point (x0, y0, z0) vector modulus m (x, y, z), azimuth angle θ and pitch angle μ to confirm the key points:

[0125]

[0126] Wherein, j is the serial number, j = 1, 2, ..., arctan is the inverse tangent function, and arcsin is the inverse sine function;

[0127] The amplitude m is used as the weight, and the histogram is used to count the θ and μ in the neighborhood of the feature point. The histogram is Gaussian weighted to determine the main orientation of the key point.

[0128] Use ICP to register the temporary point cloud P of the mandible before and after rotation s and target point cloud P t , get the rotation matrix M, for the rotation matrix M, let the direction vector of the straight line l be n(n x ,n y ,n z ), passing through the point (x0,y0,z0), there is

[0129]

[0130] Where K = 1-cosθ, m = (n x x0+n y y0+n zz0),t1=(x0-n x m)K+(n z y0-n y z0)sinθ,t2=(y0-n y m)K+(n x y0-n z x0)sinθ,t3=(z0-n z m)K+(n y x0-n x y0)sinθ, we get:

[0131]

[0132]

[0133]

[0134]

[0135] Find the rotation vector n(n x ,n y ,n z ) and the rotation angle θ; n x , n y , n z is the component of the rotation vector n in each direction;

[0136] Set the corresponding points P1 and P2 to rotate around point P by angle θ in plane Plane2, and solve the rotation reference point P:

[0137] Make the direction vectors of P1 and P2:

[0138] n1=P1-P2

[0139] Solve for the midpoint of P1 and P2:

[0140]

[0141] Use n and the midpoint P0 of points P1 and P2 to make plane Plane1:

[0142] n x (xP 0x )+n y (yP 0y )+n z (zP 0z )=0

[0143] Among them, P 0x is the x-axis coordinate of point P0, P 0y is the y-axis coordinate of point P0, P 0z is the z-axis coordinate of point P0; nx , n y , n z is the coordinate of the rotation vector n; n 1x , n 1y , n 1z is the coordinate of the rotation vector n1; the subscripts x, y, and z represent the corresponding axis coordinates;

[0144] Use direction vector n1 and point P0 to make plane Plane2:

[0145] n 1x (xP 0x )+n 1y (yP 0y )+n 1z (zP 0z )=0

[0146] L is the intersection of planes Plane1 and Plane2:

[0147]

[0148] Find a point P on L such that ∠P1PP2 is θ:

[0149]

[0150] By traversing all corresponding points, a number of rotation reference point scatter points are obtained;

[0151] Fit a rough hinge axis using the least squares method:

[0152] Set the equation of the line in space:

[0153]

[0154] Perform equivalent conversion and get:

[0155]

[0156]

[0157] in: The straight line passes through the fixed point (b1, b2), and the direction vector of the straight line is (k1, k2);

[0158] Since a straight line in space is obtained by the intersection of two planes, a straight line fitting is performed on the two plane equations, and the scattered points are fitted into a straight hinge axis l:

[0159]

[0160]

[0161]

[0162]

[0163] Among them, x i , z i ,y i Represents the coordinates of the scattered points;

[0164] Use the gradient descent algorithm to solve the refined hinge axis:

[0165] The direction vector of the hinge axis l is (n x ,n y ,n z ), through the point (x0, y0, z0), the rotation angle is θ, the rotation matrix is M, and the objective function is set to:

[0166]

[0167] Among them, X i , Y i is the corresponding point;

[0168] The parameters to be solved are:

[0169] θ=(x0,y0,z0,n x ,n y ,n z ) T

[0170] so:

[0171]

[0172]

[0173]

[0174]

[0175] Among them, M ij is the data of the i-th row and j-th column in the rotation matrix M;

[0176] Iterative process:

[0177] The rotation direction and angle of the coarse hinge axis are used as the initial value of θ, and the gradient descent iteration formula is:

[0178]

[0179] Among them, θ k ,θ k-1 In the figure, θ is the rotation angle, the upper right subscripts k and k-1 are the number of iterations, and λ is the learning rate between the interval (0,1).

[0180] Corresponding points are extracted using point cloud registration and nearest neighbor search. To reduce the number of point clouds, the SIFT operator is used to extract key points, replacing the original point cloud. Leveraging the characteristics of rotational motion, the coarse hinge axis is extracted using rotational geometric relationships. This coarse hinge axis provides an excellent initial position for extracting the fine hinge axis. The extracted corresponding points and the coarse hinge axis are then passed as parameters to the gradient descent algorithm to determine the fine hinge axis. This method is simple and easy to use.

[0181] Furthermore, the method of using the hinge axis to assist in simulating mandibular movement is as follows:

[0182] Mandibular movements can be manifested in the form of mouth opening and closing, protrusion, and lateral movement;

[0183] The direction of the hinge axis is n=(n x ,n y ,n z ) passes through the condyle point Q(x0,y0,z0), and the hinge axis line is expressed as:

[0184]

[0185] According to the size of the opening amplitude, the mouth opening and closing movement is divided into small jaw opening movement and large jaw opening movement;

[0186] The trajectory of small jaw opening movement is:

[0187] P t =RP s

[0188] Among them, P s is a temporary point cloud, P t is the target point cloud, R is the point cloud (n x ,n y ,n z ,θ,x0,y0,z0) to represent the direction (n x ,n y ,n z ), the rotation matrix of the straight line through the fixed point (x0, y0, z0) with an angle θ;

[0189] When performing large jaw opening motion, in order to make the simulation smoother, the motion trajectory is further decomposed. Point Q is first translated to Q″(x1, y1, z1), and rotated around the hinge axis by an angle of θ″. Then, it is translated from Q″ to point Q′(x2, y2, z2), and rotated around the hinge axis of Q′ by an angle of θ′.

[0190] The translation matrix from point Q to point Q″ is:

[0191]

[0192] At this time, the hinge axis passes through the point Q″(x1,y1,z1), and the rotation matrix R1 is:

[0193] R1=R(n x ,n y ,n z ,θ″,x1,y1,z1)

[0194] The translation matrix from point Q″ to point Q′ is:

[0195]

[0196] At this time, the hinge axis passes through the point Q′(x2,y2,z2), so the rotation matrix R2 is:

[0197] R2=R(n x ,n y ,n z ,θ′,x2,y2,z2)

[0198] The motion trajectory of the wide jaw opening movement is:

[0199] P t =R2T2R1T1P s ;

[0200] Protrusive movement means that the mandible tilts forward, and point Q translates to point Q″(x1, y1, z1). Protrusive movement only involves translation, not rotation, so the translation matrix is constructed as follows:

[0201]

[0202] The trajectory of the forward movement is:

[0203] P t =T3P s ;

[0204] The manifestation of lateral movement is the lateral movement of the mandible, in which the left condyle S does not move, and the right condyle Q rotates around the line l perpendicular to the hinge axis by an angle θ to Q′;

[0205] Assume that the direction of line l is n l =(n lx ,n ly ,n lz ), through the condyle point S(S x ,S y ,S z ), the rotation matrix is:

[0206] R3=R(n lx ,n ly ,n lz ,θ,S x ,S y ,Sz )

[0207] The trajectory of lateral movement is:

[0208] P t =R3P s .

[0209] For the three forms of mandibular movement: opening and closing, protrusion and lateral movement, the hinge axis is used to assist in simulating mandibular movement from the perspective of rigid body transformation, which facilitates the visualization of mandibular movement.

[0210] The present invention also provides a CBCT-based virtual articulator system, comprising a data acquisition module and a construction module. The data acquisition module is used to acquire CBCT image projection information inside the human oral cavity. The output end of the data acquisition module is connected to the input end of the construction module. The construction module executes the method of the present invention to construct a virtual articulator system.

[0211] The system simulates mandibular movement by controlling the hinge axis, thus avoiding the cumbersome operation of traditional mechanical articulators. BRIEF DESCRIPTION OF THE DRAWINGS

[0212] Figure 1 1 is a flow chart of a method for constructing a virtual articulator system based on CBCT according to the present invention;

[0213] Figure 2 Schematic diagram of the geometric coordinate relationship of the three-dimensional reconstruction of the cone-beam circular scanning trajectory in the method for constructing a virtual articulator system based on CBCT of the present invention;

[0214] Figure 3 3 is a cube diagram of the MC algorithm of the method for constructing a virtual articulator system based on CBCT of the present invention;

[0215] Figure 4 2 is a cube numbering diagram of the MC algorithm of the method for constructing a virtual articulator system based on CBCT of the present invention;

[0216] Figure 5 15 basic scenarios of the MarchingCubes algorithm of the method for constructing a virtual articulator system based on CBCT according to the present invention;

[0217] Figure 6 Schematic diagram of point cloud in VoxelGrid voxel filter downsampling of the method for constructing a virtual articulator system based on CBCT of the present invention;

[0218] Figure 7 is the corresponding point P in the coarse registration process of the virtual articulator system construction method based on CBCT of the present invention. x , P y Coordinate diagram of

[0219] Figure 8 is the corresponding point P in the coarse registration process of the virtual articulator system construction method based on CBCT of the present invention. x , P y The center of mass of the coordinate diagram is moved to the center of the circle;

[0220] Figure 9 Schematic diagram of the structure of the rotation reference point P in the rough hinge axis extraction of the CBCT-based virtual articulator system construction method of the present invention. DETAILED DESCRIPTION

[0221] The following describes embodiments of the present invention in detail. Examples of the embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended only to explain the present invention and are not to be construed as limiting the present invention.

[0222] In the description of the present invention, it should be understood that the terms "longitudinal", "transverse", "up", "down", "front", "back", "left", "right", "vertical", "horizontal", "top", "bottom", "inside", "outside", etc., indicating the orientation or position relationship, are based on the orientation or position relationship shown in the accompanying drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore cannot be understood as limiting the present invention.

[0223] In the description of the present invention, unless otherwise specified and limited, it should be noted that the terms "installed", "connected" and "connected" should be understood in a broad sense. For example, it can be a mechanical connection or an electrical connection, or it can be the internal communication between two components. It can be a direct connection or an indirect connection through an intermediate medium. For ordinary technicians in this field, the specific meanings of the above terms can be understood according to the specific circumstances.

[0224] The present invention discloses a method for constructing a virtual articulator system based on CBCT. The method uses CBCT (Cone beam CT, cone beam CT, cone beam projection computerized tomography) to scan the oral cavity and reconstruct a realistic three-dimensional model of the oral mandible. On this basis, the hinge axis is extracted based on the movement characteristics of the oral mandible. Finally, by controlling the hinge axis, the mandibular movement is simulated, avoiding the cumbersome operation of traditional mechanical articulators. Figure 1 As shown, the construction method includes the following steps:

[0225] The system collects CBCT image projection information from inside the human oral cavity and reconstructs CT slice data using the FDK algorithm. To collect X-ray projections, the patient only needs to fix their head on the supporting device that comes with the CBCT device. The CBCT device then starts to rotate, the X-ray tube emits X-rays, and the flat-panel detector on the other side collects the X-ray projections.

[0226] Based on CT tomographic data, the maxillary and mandibular mesh models were reconstructed using the MC algorithm.

[0227] Preprocess the maxillary and mandibular mesh models, including noise removal, coarse registration, and fine registration, to optimize the maxillary and mandibular mesh models.

[0228] Extract the hinge axis using the optimized maxillary and mandibular mesh models;

[0229] In view of the expression of mandibular movement, the hinge axis is used to assist in simulating mandibular movement from the perspective of rigid body transformation.

[0230] In a preferred embodiment of the present invention, the method for reconstructing CT slice data using the FDK algorithm (the FDK algorithm is the most commonly used algorithm in three-dimensional cone-beam reconstruction) is as follows:

[0231] The FDK reconstruction algorithm is a reconstruction algorithm designed specifically for the circular focus orbit of cone beams. Its greatest advantage is that the reconstructed image is better when the cone angle is less than 10°. Figure 2 As shown in Figure 3, the three-dimensional reconstructed geometric coordinate relationship of the scanning trajectory of the CBCT image projection information results in three different coordinate systems, namely the scanning object coordinate system, the detector projection data coordinate system, and the ray source coordinate system, which correspond to o-xyz, o-tz, and o-st, respectively, where z is the rotation center axis, s passes through the center of the ray source, and is perpendicular to the detector plane.

[0232] The formula of the FDK algorithm in the center plane is:

[0233]

[0234] u=ycosβ-xsinβ, v=xcosβ+ysinβ,

[0235] Here, p(β,a,b) represents the projection data of point X0(x0,y0,z0) in the projection direction β.

[0236]

[0237] β is the projection angle, (a, b) is the detector array coordinate, R(β) is the orbit radius, g(τ) is the filter function, c is the complex unit, ρ is the polar radius; u, v, and τ are intermediate codes; y, x, and z are the coordinates of the corresponding axes. f(x, y, z) is the polar form of the double integral, with x, y, and z projected onto a, b, and β.

[0238] In a preferred embodiment of the present invention, the method for reconstructing the maxillary and mandibular mesh models using the MC (Marching Cubes) algorithm (also known as Isosurface Extraction, a classic surface rendering algorithm whose main essence is to approximate isosurfaces through linear interpolation in a three-dimensional discrete data field) is as follows:

[0239] The FDK reconstruction algorithm was used to obtain a 3D image of the brain's internal structure. To obtain the upper and lower jaw models, the Marching Cubes algorithm was used. The Marching Cubes algorithm extracts isosurfaces from a 3D discrete data field, primarily used to extract mesh information from objects.

[0240] Divide the upper and lower jaw area data, use a cuboid to surround the upper and lower jaw areas, and divide the cuboid into m×n×k cubes. The upper and lower layers of each cube are covered by different faults, such as Figure 3 As shown;

[0241] Number the vertices of each cube as v = (v1, v2, v3, v4, v5, v6, v7, v8) and the edges of the cube as e = (e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12), as follows: Figure 4 As shown;

[0242] Determine whether the small cube is inside or outside the upper and lower jaw areas: Based on the tissue density on the tomography, calculate whether each vertex of the cube is inside or outside the upper and lower jaw areas. The cube has 8 vertices in total, which requires 8 bits to represent. There are a total of 256 situations; Figure 5 As shown in Figure 1, there are 15 basic cases of the MarchingCubes algorithm. The other 241 cases can be achieved by rotating, mapping, and other methods of these 15 basic cases.

[0243] Interpolation calculation: Based on the internal and external positions of the eight vertices of the cube on the upper and lower jaws, the interpolation formula is used to calculate the positions on the edges that need to be triangulated:

[0244]

[0245] Among them, isovalue is the upper and lower tissue density thresholds, P1 and P2 represent the coordinates of the two endpoints on the edge, and V1 and V2 represent the values at these two endpoints;

[0246] Calculate normal direction: After the previous four steps, the point and surface information of the isosurface can be obtained. In order to further improve the display effect, the normal vector of the triangle needs to be solved. Calculate the normal vector of the triangle. The three points A, B, C of the triangle have the normal vector Traverse all cubes to obtain the Mesh models of the upper and lower jaws.

[0247] In a preferred embodiment of the present invention, the method for denoising the upper and lower jaw Mesh models is as follows:

[0248] The point cloud data of the upper and lower jaw mesh models are input into the VoxelGrid voxel filter (a downsampling filter that uses voxelization to reduce the number of point clouds);

[0249] Set up a large cube that encloses the point cloud data. The length, width, and height of the large cube are Lx, Ly, and Lz respectively. Divide the large cube into several small cubes. Let the total number of points in the point cloud be N. Then the side length L of the small cube is:

[0250]

[0251] Among them, α is the scaling factor used to adjust the side length of the small cube, and s is the scaling coefficient;

[0252] By adjusting the parameters α and s, the number of small cubes is obtained using the rounding function, and the point cloud data is divided into m*n*k small cubes;

[0253] Traverse the point cloud position information in each small cube and calculate the centroid coordinates of the point cloud in the small cube:

[0254]

[0255] Among them, X0, Y0, Z0 are the coordinates of the center of mass, x i ,y i , z i are the coordinates of the point cloud in the small cube, is the number of point clouds falling in the small cube, i is the sequence number, i = 1, 2, ..., represents the i-th point cloud in the small cube;

[0256] By downsampling the voxel grid, the number of points in the point cloud data is reduced, and the shape characteristics of the point cloud do not change, such as Figure 6 shown.

[0257] In a preferred embodiment of the present invention, the method for coarse registration of the maxillary and mandibular mesh models is as follows:

[0258] After the maxillary point cloud is downsampled by VoxelGrid, the number of point clouds has been reduced. However, in order to provide a better initial position for fine registration, coarse registration is still required. x , P y ,like Figure 7 As shown, construct the objective function and solve the rotation matrix R:

[0259] P y =RP x +t

[0260] Where t is the translation matrix;

[0261] Ask for P x , P y The center of mass coordinates, P x , P y The big one is a set of point clouds, and here is the centroid of the point cloud:

[0262]

[0263]

[0264] It's P x The coordinates of the midpoint, It's P y The coordinates of the midpoints, these points can be randomly selected;

[0265] Move the centroid of the two groups of points to the center of the circle, such as Figure 8 As shown, the coordinates of the circle center are:

[0266]

[0267]

[0268] Solve the rotation matrix:

[0269] p′ y =Rp′ x

[0270]

[0271] However, there is no guarantee that R is a rotation matrix. It may also be obtained through scaling transformation. In this case, the scaling transformation will cause shape deformation when applied to the maxillary point cloud. Singular value decomposition is used to solve it and construct the covariance matrix H:

[0272]

[0273] Perform singular value decomposition on H:

[0274] H=U∑VT

[0275] all:

[0276] R=UV T

[0277] Where U is an m-order orthogonal matrix, V is an n-order orthogonal matrix, ∑ is an mxn rectangular diagonal matrix, and T represents a transpose operation.

[0278] In a preferred embodiment of the present invention, the method for performing precise registration processing on the upper and lower jaw mesh models is as follows:

[0279] Using the initial R0, t0 of the coarse registration or the R obtained in the previous iteration k-1 , t k-1 Transform the point cloud data of the mandibular Mesh model to obtain a temporary transformed point cloud. Let the temporary point cloud be P s , the target point cloud is P t ; The target point cloud is P t ; Where R is the rotation matrix, t is the translation matrix, subscript 0 indicates the first rotation or translation operation, subscript k-1 indicates the k-1th rotation or translation operation; alignment requires rotating and translating Ps to Pt

[0280] P s and P t Compare and find P s Every point in P t The nearest neighbor point in the , is used for registration; if the nearest neighbor point search is performed directly, two cycles are required, and the computational complexity is O(|P s |*|P t |), which is time-consuming. The KD tree (KDimensionTree) method is used to find the nearest neighbor points of the temporary point cloud and the target point cloud to speed up the search.

[0281] Determine the objective function:

[0282]

[0283] Among them, R * represents the rotation matrix, T * represents the translation matrix, argmin represents the variable value when the objective function f(x) takes the minimum value, P s For temporary point cloud, It's P t and P s The corresponding points in , t is the translation matrix, R is the rotation matrix obtained by coarse registration; i is the sequence number, i = 1, 2, ..., representing the i-th point cloud in the small cube;

[0284] Use SVD (singular value decomposition) decomposition method to solve the R of this round k and t k ; R k and t k Where R is the rotation matrix, t is the translation matrix, and the subscript k represents the kth rotation or translation operation;

[0285] Calculate the optimal translation:

[0286] Let N = |P s |, N is the total number of points in the point cloud, and the loss function is:

[0287]

[0288] Taking its derivative, we get:

[0289]

[0290] Let the derivative be 0, then:

[0291]

[0292] make We can get:

[0293]

[0294] When the value of R is determined, the optimal t is solved according to the formula to minimize the loss; P t and P s The center of mass;

[0295] Calculate the optimal rotation:

[0296] Convert the temporary point cloud and target point cloud to centroid coordinates:

[0297]

[0298]

[0299] Only the rotation problem needs to be considered, and the loss can be written as:

[0300]

[0301] in: is the temporary point cloud in the centroid coordinates, is the target point cloud in the centroid coordinates,

[0302] T is the transpose operation,

[0303] Simplified to:

[0304]

[0305] Loss is rewritten as:

[0306]

[0307] From the definition of matrix multiplication and matrix trace, we can get:

[0308]

[0309] The problem turns into:

[0310]

[0311] make H is the covariance matrix,

[0312] Then we have:

[0313]

[0314] Perform singular value decomposition on H = U∑V T

[0315]

[0316] Since V, U, and R are all orthogonal matrices, let the transformation matrix M be:

[0317] M=V T RU

[0318] Then we have:

[0319] trace(∑V T RU) = trace(∑M) where trace represents the trace of the matrix, and V and U are orthogonal matrices;

[0320] According to the non-negative nature of singular values and the properties of orthogonal matrices: the absolute value of the elements in the orthogonal matrix is not greater than 1, and trace(∑M) is maximum only when M is the unit matrix, that is:

[0321] V T RU=E

[0322] R=V T U

[0323] The upper jaws with different openings are registered to obtain the registration matrix M, which is then applied to the mandible to unify the mandibles with different openings into one coordinate system.

[0324] In a preferred embodiment of the present invention, a KD tree method is used to find the nearest neighbor points of the temporary point cloud and the target point cloud. The steps of constructing the KD tree are as follows:

[0325] S1, select the x dimension, and put P t The points in the tree are selected using a quick selection algorithm in the x-dimension to find the median, which is then divided into two parts: the left and right parts. The median point is used as the root node, and the points on the left and right sides are divided into two subtrees.

[0326] S2: For the points in these two parts, change the dimension to the y-axis, use the quick selection algorithm again to select the median, and use the median point as the child node of the previous layer node. At the same time, each node is divided into two parts;

[0327] S3, change the dimension to the z-axis, repeat S2, again construct the child node according to the median of the z-dimension and divide the point into left and right parts;

[0328] Repeat S1-S3 until there are no nodes in the left and right parts.

[0329] In a preferred embodiment of the present invention, the method for extracting the hinge axis is as follows:

[0330] Point cloud registration and nearest neighbor search are used to extract corresponding points. To reduce the number of point clouds, the SI FT (Scale Invariant Feature Transform) operator is used to extract key points to replace the original point clouds of the upper and lower jaw mesh models:

[0331] The coordinates of each point in the original mandibular point cloud are expressed as p(x, y, z), and the scale space of the 3D point cloud is defined as L(x, y, z, σ), which is obtained by convolving a Gaussian kernel function G(x, y, z, σ) with the original mandibular point cloud p(x, y, z):

[0332] L(x,y,z,σ)=G(x,y,z,σ)*p(x,y,z)

[0333] in, σ is a constant;

[0334] Use the multiplication factor k to obtain different scales and construct the point cloud Gaussian pyramid. If the number of layers in each pyramid group is s, then set k s =2, use the more efficient DoG operator instead of the scale-normalized Gaussian Laplace operator to perform extreme value detection. The DoG operator calculation formula is:

[0335] D(x,y,z,k e σ)=L(x,y,z,k e+1 σ)-L(x,y,z,k e σ)

[0336] Among them, e∈[0,s+2]; the local extreme points of the DoG space constitute the key points. The method of finding the extreme points of the DoG function is to compare each point with all its adjacent points, and then determine whether the point is the maximum or minimum value in the neighborhood. At the same scale, the middle detection point needs to be compared with 26 points. At the upper and lower adjacent scales, the middle detection point needs to be compared with 54 points (27×2=54). Only when it is greater than or less than all 80 neighboring points can the middle detection point be considered an extreme point. By the same token, not all extreme points are key points. There is a concept called candidate key points. Specifically, if the threshold t=0.5 is set, then all extreme points less than the threshold are low-contrast key points that need to be eliminated, and the extreme points that pass the screening are called candidate key points. Get the extreme points, set the threshold, and use the extreme points that are not less than the threshold as candidate key points;

[0337] By calculating the points (x j ,y j ,z j ) to the neighborhood center point (x0, y0, z0) vector modulus m(x, y, z) to the neighborhood center point vector modulus, azimuth angle θ and pitch angle μ for key point confirmation:

[0338]

[0339] Where j is the sequence number, j = 1, 2, ...; arctan is the inverse tangent function, and arcsin is the inverse sine function; the amplitude m is used as the weight, and the histogram is used to count the θ and μ in the neighborhood of the feature point, and the histogram is Gaussian weighted to determine the main orientation of the key point;

[0340] Let the intermediate point cloud P mid =M*P s , intermediate point cloud P mid It is a special intermediate state. First, P mid The point cloud P before rotation is s .Through matrix transformation, so P mid With P s . The formation of a one-to-one correspondence, that is, through P mid The index of P can be found s The corresponding points in

[0341] In space, P mid With the rotated point cloud P t If the location is close, you can find P through the nearest neighbor mid With P t The corresponding point of P mid With P sThe corresponding relationship between the point cloud P before and after rotation can be obtained. s and P t Make connections.

[0342] Nearest neighbor search for corresponding points before and after rotation

[0343] By constructing the intermediate state point cloud P mid , so that the point cloud P before and after rotation s and P t Establish a connection. t and P mid The corresponding points can be found using the nearest neighbor algorithm. In order to make the corresponding points extracted more accurate, the nearest neighbor algorithm needs to be used twice to find P mid and P s The points that are closest to each other.

[0344] Use ICP to register the temporary point cloud P of the mandible before and after rotation s and target point cloud P t , get the rotation matrix M, for the rotation matrix M, let the direction vector of the straight line l be n(n x ,n y ,n z ), passing through the point (x0,y0,z0), there is

[0345]

[0346] Among them, K = 1-cosθ, K is the middle value, m = (n x x0+n y y0+n z z0),t1=(x0-n x m)K+(n z y0-n y z0)sinθ,t2=(y0-n y m)K+(n x y0-n z x0)sinθ,t3=(z0-n z m)K+(n y x0-n x y0)sinθ, we get:

[0347]

[0348]

[0349]

[0350]

[0351] Find the rotation vector n(nx ,n y ,n z ) and the rotation angle θ; n x , n y , n z is the component of the rotation vector n in each direction;

[0352] With the corresponding points, rotation direction and rotation angle, we can use geometric relationships to solve the reference point during rotation. Figure 9 As shown, set the corresponding points P1 and P2 to rotate around point P by angle θ in plane Plane2, and solve the rotation reference point P:

[0353] Make the direction vectors of P1 and P2:

[0354] n1=P1-P2

[0355] Solve for the midpoint of P1 and P2:

[0356]

[0357] Use n and the midpoint P0 of points P1 and P2 to make plane Plane1:

[0358] n x (xP 0x )+n y (yP 0y )+n z (zP 0z )=0

[0359] Among them, P 0x is the x-axis coordinate of point P0, P 0y is the y-axis coordinate of point P0, P 0z is the z-axis coordinate of point P0; n x , n y , n z is the coordinate of the rotation vector n; n 1x , n 1y , n 1z is the coordinate of the rotation vector n1; the subscripts x, y, and z represent the corresponding axis coordinates;

[0360] Use direction vector n1 and point P0 to make plane Plane2:

[0361] n 1x (xP 0x )+n 1y (yP 0y )+n 1z (zP 0z )=0

[0362] L is the intersection of planes Plane1 and Plane2:

[0363]

[0364] Find a point P on L such that ∠P1PP2 is θ:

[0365]

[0366] By traversing all corresponding points, a number of rotation reference point scatter points are obtained;

[0367] Using the characteristics of rotational motion and the geometric relationship of rotation, the rough hinge axis is fitted using the least squares method:

[0368] Set the equation of the line in space:

[0369]

[0370] Perform equivalent conversion and get:

[0371]

[0372]

[0373] in: The straight line passes through the fixed point (b1, b2), and the direction vector of the straight line is (k1, k2);

[0374] Since a straight line in space is obtained by the intersection of two planes, a straight line fitting is performed on the two plane equations, and the scattered points are fitted into a straight hinge axis l:

[0375]

[0376]

[0377]

[0378]

[0379] Among them, x i , z i ,y i Represents the coordinates of the scattered points;

[0380] The corresponding points and the rough hinge extracted above are passed as parameters to the gradient descent algorithm to solve the fine hinge axis:

[0381] Gradient descent is an optimization algorithm generally used to find the local minimum of a function. The direction of the gradient is the direction in which the objective function changes most rapidly. Therefore, during gradient descent, the optimal value can be obtained by simply iterating in the opposite direction of the gradient at a certain learning rate.

[0382] The direction vector of the hinge axis l is (n x ,n y ,n z ), through the point (x0, y0, z0), the rotation angle is θ, the rotation matrix is M, and the objective function is set to:

[0383]

[0384] Among them, X i , Y i are corresponding points in the point cloud;

[0385] The parameters to be solved are:

[0386] θ=(x0,y0,z0,n x ,n y ,n z ) T

[0387] so:

[0388]

[0389]

[0390]

[0391]

[0392] Among them, M ij is the data of the i-th row and j-th column in the rotation matrix M;

[0393] Iterative process:

[0394] The rotation direction and angle of the coarse hinge axis are used as the initial value of θ, and the gradient descent iteration formula is:

[0395]

[0396] Among them, θ k ,θ k-1 In the figure, θ is the rotation angle, the upper right subscripts k and k-1 are the number of iterations, and λ is the learning rate between the interval (0,1).

[0397] In a preferred embodiment of the present invention, the method of using the hinge axis to assist in simulating mandibular movement is as follows:

[0398] Mandibular movements can be manifested in the form of mouth opening and closing, protrusion, and lateral movement;

[0399] The direction of the hinge axis is n=(n x ,n y ,n z) passes through the condyle point Q(x0,y0,z0), and the hinge axis line is expressed as:

[0400]

[0401] According to the size of the opening amplitude, the mouth opening and closing movement is divided into small jaw opening movement and large jaw opening movement;

[0402] Assume that the mandibular joint rotates around the hinge axis by an angle of θ, so according to the hinge axis direction n=(n x ,n y ,n z ) and the condylar point Q(x0,y0,z0) to construct the rotation matrix, and the rotation matrix R can be constructed. The trajectory of the small jaw opening movement is:

[0403] P t =RP s

[0404] Among them, P s is a temporary point cloud, P t is the target point cloud, R is the point cloud (n x ,n y ,n z ,θ,x0,y0,z0) to represent the direction (n x ,n y ,n z ), the rotation matrix of the straight line through the fixed point (x0, y0, z0) with an angle θ;

[0405] When performing large jaw opening motion, in order to make the simulation smoother, the motion trajectory is further decomposed. Point Q is first translated to Q″(x1, y1, z1), and rotated around the hinge axis by an angle of θ″. Then, it is translated from Q″ to point Q′(x2, y2, z2), and rotated around the hinge axis of Q′ by an angle of θ′.

[0406] The translation matrix from point Q to point Q″ is:

[0407]

[0408] At this time, the hinge axis passes through the point Q″(x1,y1,z1), and the rotation matrix R1 is:

[0409] R1=R(n x ,n y ,n z ,θ″,x1,y1,z1)

[0410] The translation matrix from point Q″ to point Q′ is:

[0411]

[0412] At this time, the hinge axis passes through the point Q′(x2,y2,z2), so the rotation matrix R2 is:

[0413] R2=R(n x ,n y ,n z ,θ′,x2,y2,z2)

[0414] The motion trajectory of the wide jaw opening movement is:

[0415] P t =R2T2R1T1P s ;

[0416] Protrusive movement means that the mandible tilts forward, and point Q translates to point Q″(x1, y1, z1). Protrusive movement only involves translation, not rotation, so the translation matrix is constructed as follows:

[0417]

[0418] The trajectory of the forward movement is:

[0419] P t =T3P s ;

[0420] The manifestation of lateral movement is the lateral movement of the mandible, in which the left condyle S does not move, and the right condyle Q rotates around the line l perpendicular to the hinge axis by an angle θ to Q′;

[0421] Assume that the direction of line l is n l =(n lx ,n ly ,n lz ), through the condyle point S(S x ,S y ,S z ), the rotation matrix is:

[0422] R3=R(n lx ,n ly ,n lz ,θ,S x ,S y ,S z )

[0423] The trajectory of lateral movement is:

[0424] P t =R3P s .

[0425] For the three forms of mandibular movement: opening and closing, protrusion and lateral movement, the hinge axis is used to assist in simulating mandibular movement from the perspective of rigid body transformation, which facilitates the visualization of mandibular movement.

[0426] The present invention also provides a CBCT-based virtual articulator system, comprising a data acquisition module and a construction module. The data acquisition module is used to acquire CBCT image projection information inside the human oral cavity. The output end of the data acquisition module is connected to the input end of the construction module. The construction module executes the method of the present invention to construct a virtual articulator system.

[0427] This approach first uses CBCT to capture projections of the oral cavity, then uses the FDK algorithm to reconstruct 3D slices of the brain. Based on the varying tissue densities between the mandible and surrounding tissues in the CT slice data, mesh models of the mandible and upper jaw are reconstructed using the MC algorithm. Finally, the above process is repeated, with multiple scans and reconstructions of the mandibular models in different opening positions. This provides data for subsequent hinge axis extraction, and the mesh model also facilitates visualization of mandibular motion.

[0428] Throughout this specification, reference to terms such as "one embodiment," "some embodiments," "examples," "specific examples," or "some examples" means that a specific feature, structure, material, or characteristic described in conjunction with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, schematic representations of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples.

[0429] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to the embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the claims and their equivalents.

Claims

1. A method for constructing a virtual articulator system based on CBCT, characterized in that: The steps include: Collect CBCT image projection information inside the human oral cavity and reconstruct CT tomographic data using the FDK algorithm; Based on CT tomographic data, the maxillary and mandibular mesh models were reconstructed using the MC algorithm. Perform noise removal, coarse registration, and fine registration on the maxillary and mandibular mesh models to optimize them. Extract the hinge axis using the optimized maxillary and mandibular mesh models; In terms of the expression of mandibular movement, the hinge axis is used to assist in simulating mandibular movement from the perspective of rigid body transformation; The method to extract the hinge axis is as follows: Use the SIFT operator to extract key points to replace the original point cloud of the upper and lower jaw Mesh model: The coordinates of each point in the original mandibular point cloud are expressed as p(x, y, z), and the scale space of the 3D point cloud is defined as L(x, y, z, σ), which is obtained by convolving a Gaussian kernel function G(x, y, z, σ) with the original mandibular point cloud p(x, y, z): L(x,y,z,σ)=G(x,y,z,σ)*p(x,y,z) Where σ is a constant, Use the multiplication factor k to obtain different scales and construct the point cloud Gaussian pyramid. If the number of layers in each pyramid group is s, then set k s =2, use the more efficient DoG operator instead of the scale-normalized Gaussian Laplace operator to perform extreme value detection. The DoG operator calculation formula is: D(x,y,z,k e σ)=L(x,y,z,k e+1 σ)-L(x,y,z,k e σ) Where, e∈[0,s+2]; obtain extreme points, set a threshold, and use extreme points that are not less than the threshold as candidate key points; By calculating the points (x j ,y j ,z j ) to the neighborhood center point (x0, y0, z0) vector modulus m (x, y, z), azimuth angle θ and pitch angle μ to confirm the key points: Wherein, j is a serial number, j = 1, 2, ...; arctan is the inverse tangent function, and arcsin is the inverse sine function; The amplitude m is used as the weight, and the histogram is used to count the θ and μ in the neighborhood of the feature point. The histogram is Gaussian weighted to determine the main orientation of the key point. Use ICP to register the temporary point cloud P of the mandible before and after rotation s and target point cloud P t , get the rotation matrix M, for the rotation matrix M, let the direction vector of the straight line l be n(n x ,n y ,n z ), passing through the point (x0,y0,z0), there is where \(K = 1-\cos\theta\), \(m=(n x x0 + n y y0 + n z z0)\), \(t1=(x0 - n x m)K+(n z y0 - n y z0)\sin\theta\), \(t2=(y0 - n y m)K+(n x y0 - n z x0)\sin\theta\), \(t3=(z0 - n z m)K+(n y x0 - n x y0)\sin\theta\), we get: Find the rotation vector n(n x ,n y ,n z ) and the rotation angle θ; n x , n y , n z is the component of the rotation vector n in each direction; Set the corresponding points P1 and P2 to rotate around point P by angle θ in plane Plane2, and solve the rotation reference point P: Make the direction vectors of P1 and P2: n1=P1-P2 Solve for the midpoint of P1 and P2: Use n and the midpoint P0 of points P1 and P2 to make plane Plane1: n x (xP 0x )+n y (yP 0y )+n z (zP 0z )=0 Among them, P 0x is the x-axis coordinate of point P0, P 0y is the y-axis coordinate of point P0, P 0z is the z-axis coordinate of point P0; n x , n y , n z is the coordinate of the rotation vector n; n 1x , n 1y , n 1z is the coordinate of the rotation vector n1; the subscripts x, y, and z represent the corresponding axis coordinates; Use direction vector n1 and point P0 to make plane Plane2: n 1x (xP 0x )+n 1y (yP 0y )+n 1z (zP 0z )=0 L is the intersection of planes Plane1 and Plane2: Find a point P on L such that ∠P1PP2 is θ: By traversing all corresponding points, a number of rotation reference point scatter points are obtained; Fit a rough hinge axis using the least squares method: Set the equation of the line in space: Perform equivalent conversion and get: in: The straight line passes through the fixed point (b1, b2), and the direction vector of the straight line is (k1, k2); Since a straight line in space is obtained by the intersection of two planes, a straight line fitting is performed on the two plane equations, and the scattered points are fitted into a straight hinge axis l: Among them, x i , z i ,y i Represents the coordinates of the scattered points; Use the gradient descent algorithm to solve the refined hinge axis: The direction vector of the hinge axis l is (n x ,n y ,n z ), through the point (x0, y0, z0), the rotation angle is θ, the rotation matrix is M, and the objective function is set to: Among them, X i , Y i are corresponding points in the point cloud; The parameters to be solved are: θ=(x0,y0,z0,n x ,n y ,n z ) T so: Among them, M ij is the data of row i and column j in the rotation matrix M; Iterative process: The rotation direction and angle of the coarse hinge axis are used as the initial value of θ, and the gradient descent iteration formula is: Among them, θ k ,θ k-1 In the figure, θ is the rotation angle, the upper right subscripts k and k-1 are the number of iterations, and λ is the learning rate between the interval (0,1).

2. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method of reconstructing CT tomographic data using the FDK algorithm is as follows: According to the three-dimensional reconstructed geometric coordinate relationship of the scanning trajectory of the CBCT image projection information, three different coordinate systems are obtained, namely the scanning object coordinate system, the detector projection data coordinate system, and the ray source coordinate system, which correspond to o-xyz, o-tz, and o-st respectively, where z is the rotation center axis, s passes through the center of the ray source and is perpendicular to the detector plane; The formula of the FDK algorithm in the center plane is: Here, p(β,a,b) represents the projection data of point X0(x0,y0,z0) in the projection direction β. β is the projection angle, (a, b) is the detector array coordinate, R(β) is the orbit radius, g(τ) is the filter function, c is the complex unit, ρ is the polar diameter; u, v, τ are all intermediate codes; y, x, z are the coordinates of the corresponding coordinate axes.

3. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method of reconstructing the maxillary and mandibular mesh models using the MC algorithm is as follows: Divide the upper and lower jaw area data by using a cuboid to enclose the upper and lower jaw areas, and divide the cuboid into m×n×k cubes, with the upper and lower layers of each cube covered by different faults; Number the vertices of each cube as v = (v1, v2, v3, v4, v5, v6, v7, v8) and the edges of the cube as e = (e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, e11, e12); By using the tissue density on the tomogram, we calculated whether each vertex of the cube was inside or outside the upper and lower jaw areas; According to the internal and external conditions of the eight vertices of the cube on the upper and lower jaws, the interpolation formula is used to calculate the positions on the edges that need to be triangulated: Among them, isovalue is the upper and lower tissue density thresholds, P1 and P2 represent the coordinates of the two endpoints on the edge, and V1 and V2 represent the values at these two endpoints; Calculate the normal vector of the triangle, the three points A, B, C of the triangle, the normal vector Traverse all cubes to obtain the Mesh models of the upper and lower jaws.

4. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method for denoising the upper and lower jaw Mesh models is as follows: The point cloud data of the upper and lower jaw mesh models are input into the VoxelGrid voxel filter; Set up a large cube that encloses the point cloud data. The length, width, and height of the large cube are Lx, Ly, and Lz respectively. Divide the large cube into several small cubes. Let the total number of points in the point cloud be N. Then the side length L of the small cube is: Among them, α is the scaling factor used to adjust the side length of the small cube, and s is the scaling coefficient; By adjusting the parameters α and s, the number of small cubes is obtained using the rounding function, and the point cloud data is divided into m*n*k small cubes; Traverse the point cloud position information in each small cube and calculate the centroid coordinates of the point cloud in the small cube: Among them, X0, Y0, Z0 are the coordinates of the center of mass, x i ,y i , z i are the coordinates of the point cloud in the small cube, is the number of point clouds falling in the small cube, i is the sequence number, i = 1, 2, ..., represents the i-th point cloud in the small cube; By downsampling the voxel grid, the number of points in the point cloud data is reduced, and the shape characteristics of the point cloud do not change.

5. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method for coarse registration of the maxillary and mandibular mesh models is as follows: Select two groups of corresponding points P in the maxillary x , P y , construct the objective function and solve the rotation matrix R: P y =RP x +t Where t is the translation matrix; Ask for P x , P y The centroid coordinates of: It's P x The coordinates of the midpoint, It's P y Coordinates of the midpoint; Move the centroid of the two groups of points to the center of the circle. The coordinates of the center of the circle are: Solve the rotation matrix: p′ y =Rp′ x R=p′ y p′ x -1 R is the rotation matrix or obtained by scaling transformation. Scaling transformation will cause shape deformation when applied to the maxillary point cloud. Singular value decomposition is used to solve it and construct the covariance matrix H: H=p′ x p′ y T Perform singular value decomposition on H: H=U∑V T all: R=UV T Where U is an m-order orthogonal matrix, V is an n-order orthogonal matrix, ∑ is an mxn rectangular diagonal matrix, and T represents a transpose operation.

6. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method for precise registration of the upper and lower jaw Mesh models is as follows: Using the initial R0, t0 of the coarse registration or the R obtained in the previous iteration k-1 , t k-1 Transform the point cloud data of the mandibular Mesh model to obtain a temporary transformed point cloud. Let the temporary point cloud be P s , the target point cloud is P t ; Where R is the rotation matrix, t is the translation matrix, subscript 0 indicates the first rotation or translation operation, and subscript k-1 indicates the k-1th rotation or translation operation; P s and P t Compare and find P s Every point in P t The nearest neighbor points in are aligned; Determine the objective function: Among them, R * represents the rotation matrix, T * represents the translation matrix, argmin represents the variable value when the objective function f(x) takes the minimum value, P s For temporary point cloud, It's P t and P s The corresponding points in , t is the translation matrix, R is the rotation matrix obtained by coarse registration; i is the sequence number, i = 1, 2, ..., representing the i-th point cloud in the small cube; Use SVD decomposition method to solve R k and t k ; R k and t k Where R is the rotation matrix, t is the translation matrix, and the subscript k represents the kth rotation or translation operation; Calculate the optimal translation: Let N = |P s |, N is the total number of points in the point cloud, and the loss function is: Taking its derivative, we get: Let the derivative be 0, then: make We can get: When the value of R is determined, the optimal t is solved according to the formula to minimize the loss; P t and P s The center of mass; Calculate the optimal rotation: Convert the temporary point cloud and target point cloud to centroid coordinates: Only the rotation problem needs to be considered, and the loss can be written as: in: is the temporary point cloud in the centroid coordinates, is the target point cloud in the centroid coordinates, T is the transpose operation, which can be simplified to: Loss is rewritten as: From the definition of matrix multiplication and matrix trace, we can get: The problem is transformed into: make H is the covariance matrix, Then we have: Perform singular value decomposition on H = U∑V T Among them, trace means finding the trace of the matrix, V and U are orthogonal matrices; Since V, U, and R are all orthogonal matrices, let the transformation matrix M be: M=V T RU Then we have: trace(∑V T RU)=trace(∑M) According to the non-negative nature of singular values and the properties of orthogonal matrices: the absolute value of the elements in the orthogonal matrix is not greater than 1, and trace(∑M) is maximum only when M is the unit matrix, that is: V T RU=E R=V T U The upper jaws with different openings are registered to obtain the registration matrix M, which is then applied to the mandible to unify the mandibles with different openings into one coordinate system.

7. The method for constructing a virtual articulator system based on CBCT according to claim 6, wherein: The KD tree method is used to find the nearest neighbor points of the temporary point cloud and the target point cloud. The steps to construct the KD tree are as follows: S1, select the x dimension, and put P t The points in the tree are selected using a quick selection algorithm in the x-dimension to find the median, which is then divided into two parts: the left and right parts. The median point is used as the root node, and the points on the left and right sides are divided into two subtrees. S2: For the points in these two parts, change the dimension to the y-axis, use the quick selection algorithm again to select the median, and use the median point as the child node of the previous layer node. At the same time, each node is divided into two parts; S3, change the dimension to the z-axis, repeat S2, again construct the child node according to the median of the z-dimension and divide the point into left and right parts; Repeat S1-S3 until there are no nodes in the left and right parts.

8. The method for constructing a virtual articulator system based on CBCT according to claim 1, wherein: The method of using the hinge axis to assist in simulating mandibular movement is as follows: Mandibular movements can be manifested in the form of mouth opening and closing, protrusion, and lateral movement; The direction of the hinge axis is n=(n x ,n y ,n z ) passes through the condyle point Q(x0,y0,z0), and the hinge axis line is expressed as: According to the size of the opening, the mouth opening and closing movement can be divided into small jaw opening movement and large jaw opening movement. The trajectory of small jaw opening movement is: P t =RP s Among them, P s is a temporary point cloud, P t is the target point cloud, R is the point cloud (n x ,n y ,n z ,θ,x0,y0,z0) to represent the direction (n x ,n y ,n z ), the rotation matrix of the straight line through the fixed point (x0, y0, z0) with an angle θ; When performing large jaw opening motion, in order to make the simulation smoother, the motion trajectory is further decomposed, and point Q is first translated to Q″(x 1, y1, z1), rotated around the hinge axis by an angle of θ″, then translated from Q″ to point Q′(x2, y2, z2), and rotated around the hinge axis of Q′ by an angle of θ′; The translation matrix from point Q to point Q″ is: At this time, the hinge axis passes through the point Q″(x1,y1,z1), and the rotation matrix R1 is: R1=R(n x ,n y ,n z ,θ″,x1,y1,z1) The translation matrix from point Q″ to point Q′ is: At this time, the hinge axis passes through the point Q′(x2,y2,z2), so the rotation matrix R2 is: R2=R(n x ,n y ,n z ,θ′,x2,y2,z2) The motion trajectory of the wide jaw opening movement is: P t =R2T2R1T1P s 4 Protrusive movement means that the mandible tilts forward, and point Q translates to point Q″(x1, y1, z1). Protrusive movement only involves translation, not rotation, so the translation matrix is constructed as follows: The trajectory of the forward movement is: P t =T3P s ; The manifestation of lateral movement is the lateral movement of the mandible, in which the left condyle S does not move, and the right condyle Q rotates around the line l perpendicular to the hinge axis by an angle θ to Q′; Assume that the direction of line l is n l =(n lx ,n ly ,n lz ), through the condyle point S(S x ,S y ,S z ), the rotation matrix is: R3=R(n lx ,n ly ,n lz ,θ,S x ,S y ,S z ) The trajectory of lateral movement is: P t =R3P s 。 9. A CBCT-based virtual articulator system, characterized in that: The invention comprises a data acquisition module and a construction module, wherein the data acquisition module is used to acquire CBCT image projection information inside the human oral cavity, the output end of the data acquisition module is connected to the input end of the construction module, and the construction module executes the method according to any one of claims 1 to 8 to construct a virtual articulator system.

Images