CN107689049B

CN107689049B - Tooth preparation body restoration model characteristic line extraction method

Info

Publication number: CN107689049B
Application number: CN201610634704.3A
Authority: CN
Inventors: 王煜; 沈永辉; 高键; 夏鸿建; 马杰
Original assignee: Foshan Nuowei Technology Co ltd
Current assignee: Foshan Nuowei Technology Co ltd
Priority date: 2016-08-03
Filing date: 2016-08-03
Publication date: 2020-06-16
Anticipated expiration: 2036-08-03
Also published as: CN107689049A

Abstract

The invention discloses a tooth preparation body restoration model characteristic line extraction method, which comprises the following steps: picking up key feature points in a feature area of a preparation body, then searching out all optimal candidate nodes among the key feature points by adopting a heuristic search algorithm according to a comprehensive constraint control function, and finally sequentially connecting all the optimal candidate nodes to form an initial feature path among the key feature points, wherein the comprehensive constraint control function takes the distance, the direction and the features as comprehensive constraint conditions; editing initial characteristic paths among the key characteristic points; and performing form optimization on the edited initial characteristic path by adopting an improved characteristic path fairing algorithm to obtain a characteristic line of the tooth preparation body restoration model, wherein the improved characteristic path fairing algorithm directly completes the form optimization of the initial characteristic path in a three-dimensional space according to the space topological form of the preparation body grid. The invention has the advantages of accuracy, stability, high efficiency and guaranteed extraction quality, and can be widely applied to the field of false tooth repair.

Description

Tooth preparation body restoration model characteristic line extraction method

Technical Field

The invention relates to the field of false tooth restoration, in particular to a tooth preparation body restoration model characteristic line extraction method.

Background

The noun explains:

triangular mesh: one type of polygonal mesh is a data structure used in computer graphics to model various irregular objects. The surface of an object in the real world is intuitively formed by curved surfaces; in the computer world, only discrete structures can be used to simulate real continuous things. A real world surface is actually composed of numerous small polygonal patches in a computer.

Tooth preparation body: in the digital manufacturing process of the crown and the inlay, a dentist needs to prepare an affected tooth in advance (for short, preparation), and the preparation process of the tooth is that the dentist prepares a tooth preparation body which meets the form requirement on the affected tooth by using a high-speed electric drill according to medical knowledge.

Characteristic information: the information is information that can represent a significant shape of a model, and mainly includes elements such as a ridge line, a valley line, an inflection point, and a boundary line. In the differential geometry, the feature information of the curved surface is generally calculated by a differential amount such as a principal curvature, a principal direction, and a principal curvature extremal coefficient.

Characteristic line: a line formed by the connection of feature points containing feature information.

In recent years, the application of CAD/CAM technology in the field of prosthodontia has been rapidly developed. Compared with the traditional manual restoration mode, the oral CAD/CAM technology remarkably reduces the labor intensity of dental technicians, greatly improves the restoration efficiency, has better restoration quality and meets the design requirement of personalized restoration. The characteristic lines of the tooth preparation body comprise a crown preparation body neck margin line and an inlay preparation body hole margin line, are reference lines of full crown restoration and inlay restoration, the extraction quality of the characteristic lines directly influences the modeling precision of the restoration body, and the characteristic lines have important influence on the health of gum in relation to the comfort level of a patient wearing the restoration body. Therefore, the extraction of the characteristic line of the tooth preparation body is a key link of the design of the digital dental prosthesis modeling. The preparation body generally has complex geometric appearance characteristics, the characteristic quality of the preparation body depends on the preparation level of a doctor and the quality of a grid obtained by measurement, and surface defects, noise data and the like can influence the extraction of characteristic lines. At present, the characteristic line extraction method on the mesh curved surface is mainly divided into two main categories, namely an automatic extraction method and a semi-automatic extraction method.

1. An automatic extraction method.

Et al in the name of "Extraction of feature lines on standardized surface operatorsThe paper proposes a method for extracting a characteristic line of a triangular mesh curved surface by adopting a morphological operator, and the method firstly sets a threshold value according to an average curvature value to divide a mesh model into two regions: a low curvature area and a high curvature area (namely an initial characteristic area), then, a semantic operation in image analysis is used for expansion and corrosion treatment, the initial characteristic area is enhanced, and finally, a characteristic line is extracted from the enhanced initial characteristic area. The method reduces noise data and manual intervention, but the setting of the model threshold value can be determined only by repeated tests, and is easily limited by the number and density of triangular plates, and the accuracy is not high.

The document named "triangular mesh characteristic line extraction based on Morse-Smale complex" by the Yangzhe-Smale et al proposes a characteristic line extraction algorithm based on Morse-Smale complex, which constructs an index function by calculating the curvature of the mesh vertex, establishes MS complex based on the index function, defines the significance as a control parameter for judging the importance of the characteristic line, and finally deletes the secondary characteristics in turn by the process of simplifying the complex, thereby realizing the automatic extraction of the interconnected characteristic lines. The setting of the intensity and the significance of the characteristic line in the method also needs to be determined through experiments, and the accurate characteristic line cannot be extracted for the model with characteristic interference.

2. A semi-automatic extraction method.

Liu Shengland et al put forward a characteristic line optimization extraction method based on active contour model in the thesis entitled "optimize characteristic line on mesh curved surface with active contour model": firstly, interactively selecting a plurality of initial characteristic points, and quickly determining an initial characteristic line by utilizing a tracking projection method; then calculating the active contour model energy of the feature points on the initial feature line; and finally, moving the characteristic points on the initial characteristic line to a position with extremely small energy after multiple iterations, and realizing optimization.

Kudzuvine et al proposed a similar characteristic line extraction method to Liu Shengan et al in a paper entitled "a new characteristic line extraction method on mesh curved surface", but it did not use active contour model energy to optimize, but found out the point with the maximum mean curvature in the n-ring neighborhood of the characteristic sampling point as the new characteristic point, and then utilized these new characteristic points to fit the B spline curve and projected the B spline curve on the triangular mesh curved surface.

Both the methods proposed by Liu Shenglan et al and Kudzuvine flash et al need to be mapped to a two-dimensional plane, only curvature information is used, and the method is only suitable for small-range form optimization, so that the form between two manually selected feature points cannot be changed greatly, and the extraction quality of feature lines is difficult to ensure; more feature points are required to be selected to extract an accurate feature line, and the efficiency is low.

In summary, the current method for extracting the characteristic lines on the mesh curved surface mainly has the following disadvantages:

(1) the automatic characteristic line extraction method cannot accurately extract the characteristic lines of the tooth preparation body in a local characteristic interference area, is easily influenced by the quality of a grid model, and is not stable enough.

(2) The semi-automatic characteristic line extracting method needs to manually select a plurality of characteristic points to extract an accurate characteristic line, and the shape of a grid model between two manually selected characteristic points cannot be changed violently.

Disclosure of Invention

To solve the above technical problems, the present invention aims to: the tooth preparation body restoration model characteristic line extraction method is accurate, stable, high in efficiency and guaranteed in extraction quality.

The technical scheme adopted by the invention is as follows:

a tooth preparation body restoration model characteristic line extraction method comprises the following steps:

s1, picking up key feature points in a feature area of the preparation body, then searching out all optimal candidate nodes among the key feature points by adopting a heuristic search algorithm according to a comprehensive constraint control function, and finally sequentially connecting all the optimal candidate nodes to form an initial feature path among the key feature points, wherein the comprehensive constraint control function takes distance, direction and features as comprehensive constraint conditions;

s2, editing the initial characteristic path among the key characteristic points to correct the deviation of the initial characteristic path;

and S3, performing form optimization on the edited initial characteristic path by adopting an improved characteristic path fairing algorithm to obtain a characteristic line of the tooth preparation body restoration model, wherein the improved characteristic path fairing algorithm directly completes the form optimization of the initial characteristic path in a three-dimensional space according to the space topological form of the preparation body grid.

Further, the step S1 includes:

s11, picking up the starting point V of the characteristic path in the characteristic area of the tooth preparation body_sAnd a terminal point V_eWill V_sAnd V_sCost information F of_sStoring in associated container NodeMap and storing V_sAdding into a candidate queue NodeQueue, wherein V_sThe cost information is obtained by calculation according to a comprehensive constraint control function;

s12, fetching NodeMap first element V_topAnd judging V_top＝V_eIf yes, ending the search process, and going to step S14; otherwise, go to step S13;

s13 at V_topAs the current search center point, find out V_topOne-ring neighborhood point set Neibourrhood { V₁,V₂,…,V_oAnd emptying the associated container NodeMap, traversing each point in Neiburhood, storing the point which is not in the candidate queue NodeQueue in the Neiburhood and corresponding cost information into the associated container NodeMap, and finally returning to the step S12, wherein o is V_topTotal number of neighborhood points of a ring;

s14, sequentially connecting the starting points V_sAnd a terminal point V_eAnd constructing an initial characteristic path by all path nodes in between.

Further, the step S11 includes:

marking possible starting points and end points of the characteristic path in the characteristic area of the tooth preparation body;

judging whether the starting point and the end point of the mark are the mesh vertexes, if so, taking the starting point and the end point of the mark as the starting point V of the pick-up_sAnd a terminal point V_eOtherwise, it is markedSelecting a point with the highest main curvature value from two end points of the mesh edge where the starting point or the end point is located and three vertexes of the triangular plate where the starting point or the end point is located as a picked starting point V_sOr a terminal point V_e；

Calculating a starting point V according to the comprehensive constraint control function_sCost information F of_sThe expression of the comprehensive constraint control function f (g) is f (g) ═ α f_dis+β*f_dir+γ*(f_c+f_t)，F_sF(s), where f (g) is the current candidate point V_gOf the complex constraint control function of f_disFor the current candidate point V_gAnd end point V_eHas a Euclidean geometric distance between them, f_dirFor the current candidate point V_gWith the current search center point V_cThe angle change value of the current path direction relative to the previous path direction is constructed; f. of_cFor the current candidate point V_gInverse of maximum principal curvature, f_tFor the current candidate point V_gRelative to the current search center point V_cMinimum principal direction T_cα, β and gamma are respectively a distance control factor, a direction control factor and a characteristic control factor;

will V_sAnd V_sCost information F of_sStoring in associated container NodeMap and storing V_sAnd adding the candidate queue NodeQueue.

Further, f in the comprehensive constraint control function f (g)_cAnd f_tThe calculation process of (2) is as follows:

estimating the normal vector N of the current point p_iNormal vector N of the current point p_iThe estimation formula of (c) is:

wherein the current point p is the current candidate point V_gOr current search center point V_c，l_jAnd n_jRespectively representing the opposite side and the normal vector of the jth triangular plate of the current point p in a ring neighborhood, wherein m represents the number of grid vertexes of the current point p in the ring neighborhood;

normal vector N according to current point p_iEstablishing a local coordinate system p-uvw at the mesh verticesThe local coordinate system p-uvw is obtained by the global coordinate system O-XYZ through the following operations: translating the point O to the point p, and rotating the Z axis to be equal to the normal vector N_iOverlapping, and taking the X axis and the Y axis as u and v respectively;

determining a corresponding local cubic surface equation f (x, y) according to 2-ring neighborhood points of grid vertexes in a local coordinate system p-uvw, then solving coefficient values of f (x, y) through a least square fitting method, and further determining a corresponding Weingarten curvature matrix W₁The expression of the local cubic surface equation f (x, y) is:

wherein A, B, C, D, E, F and G are both coefficient values of f (x, y),

for Weingarten curvature matrix W₁Performing singular value decomposition to obtain two eigenvalues corresponding to the maximum principal curvature and the minimum principal curvature of the vertex of the current mesh, and then obtaining two eigenvectors t corresponding to the two eigenvalues₁And t₂Determining a maximum principal direction t of a current mesh vertex_maxAnd a minimum principal direction t_minWherein, t₁、t₂、t_maxAnd t_minAre respectively: t is t₁＝(t₁₁,t₁₂)^T，t₂＝(t₂₁,t₂₂)^T，t_max＝t₁₁x+t₁₂y，t_min＝t₂₁x+t₂₂y；

According to t₁、t₂、t_maxAnd t_minDetermining f in the complex constraint control function f (g)_cAnd f_tThe value of (c).

Further, the step S13 includes:

s131 with V_topAs the current search center point, find out V_topOne-ring neighborhood point set Neibourrhood { V₁,V₂,…,V_oAnd emptying the associated container NodeMap;

s132, using NeiAny point V in bourhood_nTraversing and analyzing the current candidate point, and finally storing the points which are not in the candidate queue NodeQueue in the Neibourhood and corresponding cost information into an associated container NodeMap;

and S133, judging whether the traversal analysis of all the points in the Neibourhood is finished, if so, returning to the step S12, otherwise, returning to the step S132.

Further, the step S132 includes:

calculating the current candidate point V according to the comprehensive constraint control function f (g)_nCost information Cost of (1);

judging the current candidate point V_nWhether the current candidate point V is in the candidate queue NodeQueue or not is judged, if yes, the current candidate point V is indicated_nThe candidate points are searched as candidate points before, no processing is carried out on the candidate points at this time, and the next candidate point in Neiburhood is continuously selected to calculate the cost information of the candidate points; otherwise, then V is set_nAdding into a candidate queue NodeQueue and ordering F_nAfter Cost, V is finally added_nAnd F_nAnd storing the data into the associated container NodeMap.

Further, the step S2 specifically includes:

and judging whether the initial characteristic path between the key characteristic points has deviation with an expected characteristic line, if so, editing the initial characteristic path between the key characteristic points by adopting point editing operation, otherwise, executing the step S3, wherein the point editing operation comprises but is not limited to delete point operation, insert point operation and move point operation.

Further, the step S3 includes:

s31, traversing the edited initial characteristic path { V }₁,V₂,…,V_hAll nodes in the path, each node V in the path is calculated_iCorresponding association point set Connect (V)_i-1,V_i,V_i+1) And has a structure V_iTo the previous node V_i-1Direction vector Dir of₁And V_iTo a subsequent node V_i+1Direction vector Dir of₂Wherein h is the total number of the nodes of the path;

s32, calculating Dir₁And Dir₂Angle (V) between_i) Size;

s33, Angle (V) according to the included Angle_i) Size judgment whether node V is needed_i-1And V_i+1A new local path is constructed to carry out the fairing of the local path, if so, a node V is arranged_i-1And V_i+1Constructing a new local path; otherwise, it is not at node V_i-1And V_i+1A new local path is constructed in between.

Further, in the step S33:

s331, Angle (V) according to included Angle_i) Size calculation node V_iThe fairing weight P (V)_i) And node V_iThe fairing weight P (V)_i) Storing the node V in a fairing weight set W_iThe fairing weight P (V)_i) The calculation formula of (2) is as follows:

wherein epsilon is a set angle parameter, and Pass represents that the node V is not considered_iThe fairing weight at this time;

s332, extracting the maximum weight W in the weight set W_maxAnd obtain w_maxCorresponding node V_iAnd associated point set Connect (V)_i-1,V_i,V_i+1)；

S333 according to w_maxCorresponding node V_iAt node V_i-1And V_i+1Constructing a new local path;

s334, continuously taking the maximum weight W from the rest weight value set W_maxAnd then returning to the step S333 until the weight set W is empty, thereby completing the morphological optimization of the initial feature path.

Further, the step S333 includes:

s3331, judging node V_iWhether the vertex belongs to the mesh or the point on the mesh edge, if the vertex belongs to the mesh, executing the step S3332, otherwise, executing the step S3333;

s3332, traverse node V_iA ring neighborhood point set E_i{X₁,X₂,…,X_qAt each point in the } a vector V is calculated_iX_iAnd vector V_i-1X_iAngle and vector V between_iX_iAnd vector V_i+1X_iAngle therebetween, the sum of these two angles theta_iThen, theta is determined_iIf it is less than pi, then obtain node V first_iCrossing edge V_iX_iSet C of_iThen construct a transit node V_i-1And V_i+1And parallel to a given plane normal vector

Then, the plane Pane and the set C are calculated_iAnd finally, constructing a new local path according to the calculated intersection point: if plane Pane and set C_iIf there is an intersection point between the partially intersected edges in (1), then the set C is calculated_iMiddle and Pane non-intersection point intersecting edge non-node end point and node end point V_iThe distance from the plane Pane is determined by taking the end point nearest to the plane Pane as the intersection point of the intersection edges and the plane, and taking the node V_i-1And V_i+1And set C_iThe intersection points of all the intersected edges and the Pane are sequentially connected to serve as a new local path; if plane Pane and set C_iAll the intersecting edges in the node V have intersection points, the node V is connected_i-1And V_i+1And set C_iThe intersection points of all the intersected edges and the Pane are sequentially connected to serve as a new local path; otherwise, directly connecting the node V_i-1And V_i+1Connected as a new local path, where q is E_iTotal number of points of (2), X_iIs a V_iCorresponding one-ring neighborhood point, set C_iElement in accordance with theta_iOrdering from small to large, given the normal vector of the plane

The value of (d) is the weighted average of the trigonometric vectors of the intersecting edges;

s3333, with node V_iThe edge is used as an intersecting edge, and a passing node V is constructed_i-1And V_i+1And parallel to a given plane normal vector

Then calculating the intersection point of the plane Pane and the intersecting edge, if the intersection point of the plane Pane and the intersecting edge does not exist, calculating the node endpoint V_iAnd node V_iThe distance between two end points of the edge and the plane P is defined as the intersection point of the end points closest to the plane

The value of (d) is the weighted average of the trigonometric vectors of the intersecting edges; then, the node V is connected_iDeleting from the characteristic Path, adding the intersection points into the characteristic Path in sequence at the current position, and finally, comparing the node V with the node V_iAnd (3) carrying out corresponding updating on related nodes and weight information thereof: first, a node V is deleted from the set W of values_iAssociated point set Connect (V)_i-1,V_i,V_i+1) The fairing weight of each node in the system is updated, and then the node V is updated_iIs concentrated in two adjacent nodes V_i-1And V_i+1Corresponding associated point set and fairing weight, and finally connecting (V) the updated associated point set_i-1,V_i,V_i+1) Storing the corresponding fairing weight value into a weight value set W, wherein the updated association point set Connect (V)_i-1,V_i,V_i+1) The corresponding calculation formula of the fairing weight value is as follows:

wherein delete is the delete operation, updatde is the update operation, V_i+2Is a V_i+1The latter node of (2).

The invention has the beneficial effects that: based on a heuristic characteristic line extraction technology, extraction of characteristic lines of a preparation body model with local characteristic interference areas can be realized only by picking up key characteristic points in a characteristic area of the preparation body, designing a comprehensive constraint control function by taking distance, direction and characteristics as comprehensive constraints and extracting initial characteristic lines among the key characteristic points by combining a heuristic search algorithm, so that the preparation body model with various quality conditions can be better dealt with, and the accuracy and the stability of characteristic line extraction are improved; by adopting a heuristic search algorithm, complete extraction of the preparation body characteristic line can be realized without picking up a large number of key characteristic points, and the efficiency is higher; and performing form optimization on the edited initial characteristic path by adopting an improved characteristic path fairing algorithm, performing form optimization on the initial characteristic line in a three-dimensional space according to the topological form of the grid space, and not mapping to a two-dimensional plane, thereby ensuring the extraction quality of the characteristic line and ensuring the finally extracted characteristic line to have an accurate form and be smooth and natural.

Drawings

FIG. 1 is an overall flow chart of a method for extracting characteristic lines of a dental preparation restoration model according to the present invention;

FIG. 2 is a schematic diagram of a complex constraint control function of the present invention;

FIG. 3 is a schematic diagram of a local coordinate system established during the calculation of discrete differential by local surface fitting according to the present invention;

FIG. 4 is a schematic diagram of four cases of the feature path shape optimization of the present invention;

FIG. 5 is a morphological comparison graph before and after optimization of the feature lines extracted from the coronal preparation model according to the present invention;

FIG. 6 is a graph showing the effect of extracting the cervical ridge line of the incisor preparation using the method of the present invention;

FIG. 7 is a graph showing the effect of extracting the cervical ridge line of the canine tooth preparation using the method of the present invention;

FIG. 8 is a graph of the effect of extracting the cervical ridge line of a molar preparation using the method of the present invention;

FIG. 9 is a diagram of the effect of hole edge line extraction of a three-sided inlay preparation using the method of the present invention.

Detailed Description

Referring to fig. 1, a method for extracting a characteristic line of a tooth preparation restoration model includes the following steps:

Further preferably, the step S1 includes:

Further preferably, the step S11 includes:

judging whether the starting point and the end point of the mark are the mesh vertexes, if so, taking the starting point and the end point of the mark as the starting point V of the pick-up_sAnd a terminal point V_eAnd on the contrary, selecting a point with the highest main curvature value from the two end points of the grid edge where the start point or the end point of the mark is located and the three vertexes of the triangular plate as the picked start point V_sOr a terminal point V_e；

Further in a preferred embodiment, f of the complex constraint control function f (g)_cAnd f_tThe calculation process of (2) is as follows:

normal vector N according to current point p_iEstablishing a local coordinate system p-uvw at the grid vertex, the local coordinate system p-uvw being obtained by performing the following operations on a global coordinate system O-XYZ: translating the point O to the point p, and rotating the Z axis to be equal to the normal vector N_iOverlapping, and taking the X axis and the Y axis as u and v respectively;

wherein A, B, C, D, E, F and G are both coefficient values of f (x, y),

Further preferably, the step S13 includes:

s132, using any point V in Neibourhood_nTraversing and analyzing the current candidate point, and finally storing the points which are not in the candidate queue NodeQueue in the Neibourhood and corresponding cost information into an associated container NodeMap;

Further preferably, the step S132 includes:

Further, as a preferred embodiment, the step S2 is specifically:

Further preferably, the step S3 includes:

s32, calculating Dir₁And Dir₂Angle (V) between_i) Size;

Further preferably, in step S33:

S333 according to w_maxCorresponding node V_iAt node V_i-1And V_i+1Between them construct a newThe local path of (a);

Further, in a preferred embodiment, the step S333 includes:

Then, the plane Pane and the set C are calculated_iAnd finally, constructing a new local path according to the calculated intersection point: if plane Pane and set C_iIf there is an intersection point between the partially intersected edges in (1), then the set C is calculated_iMiddle and Pane non-intersection point intersecting edge non-node end point and node end point V_iThe distance from the plane Pane is determined by taking the end point nearest to the plane Pane as the intersection point of the intersection edges and the plane, and taking the node V_i-1And V_i+1And set C_iThe intersection points of all the intersected edges and the Pane are sequentially connected to serve as a new local path; if plane Pane and set C_iAll the intersecting edges in the node V have intersection points, the node V is connected_i-1And V_i+1And set C_iWherein all the intersection points of the intersection edges and the Pane are connected in sequence to serve as a new local partA path; otherwise, directly connecting the node V_i-1And V_i+1Connected as a new local path, where q is E_iTotal number of points of (2), X_iIs a V_iCorresponding one-ring neighborhood point, set C_iElement in accordance with theta_iOrdering from small to large, given the normal vector of the plane

The invention will be further explained and explained with reference to the drawings and the embodiments in the description.

Example one

Aiming at the defects of the existing characteristic line extraction technology, the invention provides a heuristic tooth preparation model characteristic line extraction method facing an oral doctor, which can realize the extraction of various common preparation model characteristic lines only by manually selecting 3-5 characteristic points, can still realize the extraction of the characteristic lines in a local characteristic interference area, and adopts a new form optimization algorithm to ensure that the extracted characteristic line form is more accurate and smooth and natural.

As shown in fig. 1, the method for extracting the characteristic line of the dental preparation restoration model of the present invention includes the steps of:

step one, key feature points are interactively picked up in a neck edge or hole edge feature region of the prepared body, and then an initial feature path between two adjacent key feature points is extracted.

The path nodes between any two adjacent key feature points are obtained through a heuristic search algorithm. The principle of the heuristic search algorithm is as follows: starting from the starting point, the searching process evaluates the cost of each candidate point by using a control function, determines the position of the next searching point by comparing the sizes of the costs, jumps to the position, further searches the current searching position, and gradually reaches the final target point. And the path between the starting point and the target point is formed by connecting a series of intermediate nodes in sequence, and all the intermediate nodes are optimal candidate points obtained by comparing the cost by using a control function.

The problem of feature path extraction between two adjacent key feature points can be converted into: and constructing a characteristic path between two vertexes on the surface of the triangular mesh. By V_sRepresents the search engineBeginning point, V_eIndicates a search target point, V_cRepresenting the current search center point. At the time of first search, with V_sAs a search center point. And in the searching process, evaluating the cost obtained by calculating each candidate point in a ring of neighborhood point sets of the searching central point through a control function. With the current search center point V_cSearching candidate point V with minimum cost_nAs a new search center point of the next step, a search center point V_cContinuously updating with the searching process. When searching for the center point V_cFor searching for a target point V_eAnd then, ending the whole searching process. The optimal candidate point determines the path direction, and how to accurately filter the optimal candidate point depends on the design of the control function f (g), so the core of the problem is the design of the control function.

The invention designs a distance control function f by combining the current candidate point, the exploration central point, the position information and the characteristic information of the exploration target point_disA direction control function f_dirCharacteristic control function f_cAnd f_tThe three control functions are specifically defined as follows:

f_dis＝||v_n,v_e|| (1)

f_dir＝θ_dir(2)

f_dis: representing the euclidean geometrical distance between the current candidate point and the target point. The function ensures that the search can finally reach the target point, ensures that the search rule meets the convergence, and avoids the dead cycle of oscillation divergence. f. of_disThe smaller the value of (c), the closer the search for the characteristic path is to the final target point, the better the parameter.

f_dir: representing the angle change of the current path direction constructed by the current candidate point and the search center point relative to the previous path direction, such as the angle theta in fig. 2_dirAs shown. The function can keep the smoothness degree of the search path and reduce the generation of a sawtooth path. Theta_dirThe smaller, theThe closer the current path direction is to the previous path direction, the more gradual the change among the characteristic path nodes is; further, when θ_dirWhen the value of (2) is greater than 90 degrees, the direction of the current candidate path is opposite to the overall direction of the path, and at the moment, the candidate points are directly deleted from the candidate point set, so that the search efficiency is improved, and backtracking exploration is prevented.

f_c: representing the inverse of the maximum principal curvature K of the current candidate point. f. of_t: representing the minimum principal direction T of the current candidate point_nMinimum principal direction T relative to the current search center point_cE.g. angle theta in fig. 2_tAs shown. The feature control function ensures that the searched path can describe the salient features of the feature region to the maximum extent. (f)_c+f_t) The smaller the parameter, the better the parameter.

The triangular mesh curved surface is a piece-divided linear curved surface, cannot be accurately represented by a parameter equation or an implicit equation, can be generally approximated to a mesh form by a curved surface fitting method, and approximately represents the characteristic information of the triangular mesh by using the differential geometric characteristics of the fitted curved surface. In this embodiment, a method based on local surface fitting is used to calculate the principal curvature and principal direction of discrete differential quantities of the preparation mesh model, and the specific process is as follows:

(1) a local coordinate system is established at the mesh vertices.

The local fitting method needs to establish a local coordinate system p-uvw at each grid vertex, then fit the current point p and its neighborhood points in the local coordinate system with a polynomial surface equation f (x, y), and calculate the differential quantity of the discrete grid vertex by using a polynomial surface, as shown in fig. 3. The local coordinate system of the present embodiment can be obtained by the following method:

1) estimating the normal vector N of the p point_iWhere N is_iThe calculation method of (2) adopts a simplified method of 'opposite side balance' of the following formula:

wherein l_j、n_iRespectively representing the opposite sides of the jth triangle in a current point-ring neighborhoodVector of normal; m represents the number of grid vertexes in the current point-ring neighborhood.

2) Let the global coordinate system O-XYZ perform the following operations to establish the local coordinate system p-uvw: translating the point O to the point p, and rotating the Z axis to be equal to the normal vector N_iThe X and Y axes are taken as u and v, respectively, when they are coincident.

(2) Defining local cubic surface equations at mesh vertices

In order to more accurately calculate the characteristic information of the mesh vertex, the invention adopts 2-ring neighborhood points of the current vertex to fit a cubic polynomial curved surface, and the cubic polynomial curved surface equation f (x, y) adopted by the invention can be expressed as:

f (x, y) corresponding Weingarten curvature matrix W₁Can be expressed as:

to solve for W₁The data points in the local range { V } can be added to the unknown item in (1)₁,V₂,…,V_lTransforming to the local coordinate system of the current point, and then solving the coefficient value of f (x, y) by using the least square fitting method, wherein:

X＝(A,B,C,D,E)^T

wherein, { n_x,n_y,n_zIs any point { x ] on the curved surface_i,y_i,z_iThe normal vector of.

(3) Calculating the main curvature and main direction of the grid vertex.

Substituting coefficient A, B, C in X into W₁In, to W₁Performing singular value decomposition to obtain two eigenvalues and two eigenvectors t corresponding to the two eigenvalues₁、t₂. Wherein the maximum eigenvalue corresponds to the maximum principal curvature of the current mesh vertex and the eigenvector t₁The minimum eigenvalue corresponds to the minimum principal curvature of the current mesh vertex and the eigenvector t₂. So that the maximum and minimum principal directions t of the vertices of the current mesh_max、t_minThe following method is adopted for calculation:

the complex constraint control function f (g) can be expressed as:

f(g)＝α*f_dis+β*f_dir+γ*(f_c+f_t) (10)

α, β, gamma denote the weight factors of the three control functions, since f_dis、f_dir、(f_c+f_t) The three control functions represent three physical quantities of length, angle and curvature respectively and represent different quantitative meanings, so that the determination of the weight factors among the three control functions is an important problem. In order to ensure that the normalization processing can be respectively carried out on the three control functions on the same order of magnitude, the normalized object range is a ring neighborhood point set of the current exploration center point, and the weight factors of the three control functions are all 1 at the moment. The normalized three control functions can be rewritten as follows:

the maximum and minimum values in equations (11), (12) and (13) represent the maximum and minimum values of the relevant data (including the euclidean geometric distance d, the angular variation dir, the inverse c of the maximum principal curvature, and the minimum principal direction variation t) in the set of one-ring neighborhood points of the current search center point.

Based on the introduction of the comprehensive constraint control function and combined with a heuristic search algorithm, the extraction steps of the initial characteristic path of the invention are as follows:

(1) picking up a characteristic path starting point V in a characteristic region of a dental preparation_sEnd point V_eWill V_sAnd V_sCost information F of_sStoring the node map into a related container (the interior of which is sorted according to the cost of the node), wherein the cost information F_sCalculated from the formula (10), and V is calculated_sAnd adding the candidate queue NodeQueue.

V_s、V_eAll the mark points are positioned on the triangular mesh curved surface, and in the actual extraction process, the mark points picked up may be mesh vertexes, or may be points on mesh edges or points inside triangular plates. According to the requirement of the heuristic search algorithm, the mark points belong to the grid vertexes, and when the mark points are positioned on the grid edge or inside the triangular plate, one point with the highest main curvature value can be selected from two end points of the grid edge where the mark points are positioned and three vertexes of the triangular plate where the mark points are positioned as the initial point and the final point of indirect search.

(2) Taking out the first element V of NodeMap_top. If V_top＝V_eIf yes, ending the searching process and turning to the step (4); otherwise with V_topAs the current search center point, find out V_topOne-ring neighborhood point set Neibourrhood { V₁,V₂,…,V_oAnd emptying the associated container NodeMap and ordering V_nRepresents any point in Neiburhood, and V is calculated from the formula (10)_nCost information Cost, then for each point V in neibourwood_nAll were analyzed as follows:

(2.1) if V_nAlready in the candidate queue NodeQueue, V is indicated_nThe candidate point is searched as a candidate point, no processing is carried out on the candidate point at the moment, and the next candidate point in Neibourhood is continuously evaluatedCost condition of point selection;

(2.2) if V_nIf not in the candidate queue NodeQueue, V is set_nJoin in queue NodeQueue, and F_nAfter Cost, V is finally added_nAnd V_nCost information F of_nStoring the data into a related container NodeMap;

(3) judging whether the traversal of all the points in the Neiburhood is finished or not, if so, executing the step (4), otherwise, returning to execute the step (2);

(4) and finishing searching the characteristic path, and connecting all path nodes between the target point and the tail end point in sequence to construct an initial characteristic path.

Step two, editing the characteristic path: and adding point editing operation to the deviation part of the local path in the preliminarily extracted characteristic path for correction.

And respectively generating a characteristic path between two adjacent key characteristic points of all the key characteristic points to obtain a complete characteristic path. Considering the situation that the characteristic line of a part of tooth preparations is long and the form change is large, especially in a local characteristic interference area, an accurate characteristic line may not be extracted from the preliminarily picked key characteristic points, and then the local deviation of the characteristic line needs to be corrected by adopting an interactive mode of point editing. When the extracted feature line deviates from the feature line desired by the operator, a corresponding point edit operation may be added to correct the feature line. Among them, the point editing operation has three types: a delete point operation, an insert point operation, and a move point operation.

And step three, optimizing the path form of the characteristic line.

The preliminarily extracted feature paths are composed of orderly connected grid edges, so that the form of the feature paths lacks good smooth effect and is most likely to have jagged paths. In order to meet the fairing requirement of the characteristic line of the tooth preparation body, the invention designs a new characteristic path fairing algorithm to optimize the form of the characteristic path, so that the characteristic path after fairing is completely attached to the grid curved surface.

Path for feature Path V₁,V₂,…,V_hIndicates that any node V in the Path is connected_iAnd the previous node V connected with it in sequence_i-1The latter node V_i+1Is defined as node V_iAssociated point set Connect (V)_i-1,V_i,V_i+1) Then, there is a one-to-one correspondence between each node and its associated point set: v_i→Connect(V_i-1,V_i,V_i+1). Construction node V_iTo node V_i-1、V_i+1Dir1 and Dir2, and calculating the Angle (V) between Dir1 and Dir2_i) Can be obtained by Angle (V)_i) Whether it needs to be at the node V or not is judged_i-1And V_i+1And constructing a new path to perform local path fairing. Therefore, when constructing the fairing weight function, the new feature path fairing algorithm only needs to consider the Angle (V) factor between the current node and the front and rear nodes_i) And (4) finishing.

The calculation formula for defining the fairing weight function P is as follows:

epsilon is a set angle parameter and controls the whole fairing intensity. The smaller epsilon, the better the fairing effect, but the lower the corresponding efficiency, so that epsilon is generally 175/180, and the good fairing treatment effect can be obtained.

Based on the theoretical basis, the characteristic path form optimization method comprises the following steps:

(1) traversing all nodes in the Path of the characteristic Path, and calculating each node V_iCorresponding association point set Connect (V)_i-1,V_i,V_i+1) And node V_iSubstituting formula (14), calculating weight P (V)_i) Then the weight P (V)_i) Storing the weight value into a weight value set W;

(2) taking out the maximum value W in the weight set W_maxTo obtain the corresponding node V_iAnd associated point set Connect (V)_i-1,V_i,V_i+1) Then w is required to be_maxCorresponding node V_iMoved to the appropriate position to effect smoothing of the local path, V_iFinding new positionThe method comprises the following steps:

(2.1) when node V_iWhen belonging to a mesh vertex, node V is set as shown in FIG. 4(a)_iOne ring neighborhood point set of (1) is marked as E_i{X₁,X₂,…,X_qGet through E_iFor each point in (a), calculating a vector V_iX_iAnd V_i-1X_i、V_i+1X_iSum of included angles theta of_i. If theta is greater than theta_iLess than pi, the edge V is formed_iX_iCalled the intersecting edge, to obtain a node V_iSet of intersecting edges C_iThe elements in the set of intersecting edges are according to a vector V_iX_iAnd V_i-1X_iThe included angles of (a) are ordered from small to large. Then a pass node V is constructed_i-1And V_i+1And parallel to a given plane normal vector N_iPlane Pane of (1); calculating the intersection point of the plane Pane and the intersecting edge if C_iWhere there are partially intersecting edges that do not intersect with the plane Pane, as shown in fig. 4(b), non-node end points and node end points V of these intersecting edges are calculated_iThe distance to the plane Pane is defined as the intersection point of the end points closest to the plane. There may be a special case where the set of intersecting edges C_iEmpty, as shown in fig. 4(c), there is no intersection point at this time. If plane Pane and set C_iIf the partial intersecting edges in (1) have intersection points, the node V is connected_i-1And V_i+1And set C_iThe intersections of all the intersecting edges and the Pane are sequentially connected together to serve as a new local path, and the new local path is shown by gray thick black lines in FIG. 4 (b); if plane Pane and set C_iAll the intersecting edges in the node V have intersection points, the node V is connected_i-1And V_i+1And set C_iThe intersections of all the intersecting edges and the Pane are sequentially connected together to serve as a new local path, and the new local path is shown by gray thick black lines in FIG. 4 (a); if plane Pane and set C_iAll the intersecting edges in the node V have no intersection point, the node V is directly connected_i-1And V_i+1Connected as a new partial path as shown by the gray bold black line in fig. 4 (c).

For the construction of a planar Pane,

the value of (d) is the weighted average of the trigonometric vectors at which the intersecting edges lie. N is a radical of_iThe topological form of the local characteristic path can be reflected, and the form change trend of the local characteristic path is kept; v_i-1And V_i+1Can ensure the smoothness requirement of the local characteristic path to ensure the node V_iThe included angle between the front node and the rear node is larger than the set threshold value.

(2.2) when node V_iWhen located on the edge of the mesh, as shown in FIG. 4(d), with node V_iThe edge is used as an intersecting edge, the construction method of the plane Pane is the same as that in the step (2.1), the intersection point of the plane Pane and the intersecting edge is calculated, and if the intersection point does not exist, the node endpoint V is calculated_iAnd node V_iThe distance between the two end points of the edge and the plane Pane is taken as the intersection point. Then, the node V is connected_iDeleted from the characteristic Path and added the intersection point to the characteristic Path in order at the current position, thus connecting with the node V_iThe relevant nodes and the weight information thereof need to be updated correspondingly: first, a node V is deleted from the set W of values_iThe associated point of (2) centralizes the weight of each node, and then updates the node V_iIs concentrated in two adjacent nodes V_i-1And V_i+1Corresponding associated point set and weight, storing the updated weight of the associated point set into the weight set W, and updating the associated point set Connect (V)_i-1,V_i,V_i+1) The corresponding calculation formula of the fairing weight value is as follows:

(3) and (5) iteratively executing the step (2) until the weight set W is empty, thereby finishing the fairing processing of the characteristic path.

In order to prevent endless iteration trend which may be generated by extremely distorted feature paths, an iteration mark can be added to each node on the feature paths to record the number of times of node movement, so as to avoid repeated access.

According to the above algorithm, the morphological graphs of the feature lines extracted by the coronal preparation model before and after optimization are respectively shown in fig. 5(a) and fig. 5(b), and it can be seen that the feature lines after morphological optimization are smoother and more natural than those before optimization.

The final effect of extracting characteristic lines of incisors, canine teeth, molar teeth and three-side inlay preparation models by adopting the method of the invention is shown in figures 6-9, and the characteristic lines extracted by the method of the invention are very smooth and natural and are suitable for extracting tooth models with different qualities.

The invention provides a heuristic tooth preparation model characteristic line extraction method facing an oral doctor, which has the following advantages:

(1) the method for fitting the curved surface based on the local cubic polynomial is adopted to calculate the principal curvature and the principal direction of the vertex of the preparation mesh model, and a characteristic control function is designed according to the maximum principal curvature and the minimum principal direction, so that the characteristic information of the preparation mesh model can be accurately expressed.

(2) A distance control function, a direction control function and a characteristic control function are designed according to the current candidate node and the relevant information of the key characteristic points, the optimal candidate node is calculated through the control function, all the optimal candidate nodes between the two key characteristic points are calculated by combining a heuristic search strategy, extraction of the characteristic line of the preparation body model with the local characteristic interference area can be achieved, a large number of key characteristic points do not need to be picked up, complete extraction of the characteristic line can be achieved only by 3-5 key characteristic points, and the extracted characteristic line is more accurate in shape.

(3) The method has the advantages that the form optimization of the initial characteristic path is directly realized in the three-dimensional space according to the topological form of the grid space of the preparation body, the mapping to a two-dimensional plane is not needed, compared with the existing characteristic line optimization mode of constructing the plane, projecting to the plane, smoothing in the plane and mapping back to the grid, the calculation process is direct, the error is small, the extraction quality of the characteristic line is guaranteed, and the finally extracted characteristic line is accurate in form and smooth and natural.

While the preferred embodiments of the present invention have been illustrated and described, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims

1. A tooth preparation body restoration model characteristic line extraction method is characterized by comprising the following steps: the method comprises the following steps:

s3, performing morphological optimization on the edited initial characteristic path by adopting an improved characteristic path fairing algorithm to obtain a characteristic line of the tooth preparation body restoration model, wherein the improved characteristic path fairing algorithm directly completes morphological optimization of the initial characteristic path in a three-dimensional space according to the spatial topological morphology of a preparation body grid;

the step S3 includes:

s32, calculating Dir₁And Dir₂Angle (V) between_i) Size;

s33, Angle (V) according to the included Angle_i) Size judgment whether node V is needed_i-1And V_i+1A new local path is constructed to carry out smoothing of the local path, if so, the local path is smoothedNode V_i-1And V_i+1Constructing a new local path; otherwise, it is not at node V_i-1And V_i+1Constructing a new local path;

the Angle (V) according to the included Angle_i) Size judgment whether node V is needed_i-1And V_i+1A new local path is constructed to carry out the fairing of the local path, if so, a node V is arranged_i-1And V_i+1The step of constructing a new local path "specifically includes:

according to the Angle (V)_i) Size calculation node V_iThe fairing weight P (V)_i) And node V_iThe fairing weight P (V)_i) Storing the data into a fairing weight value set W;

taking out the maximum value W in the weight set W_maxTo obtain the corresponding node V_iAnd associated point set Connect (V)_i-1,V_i,V_i+1) Then w is required to be_maxCorresponding node V_iMoved to the appropriate position to effect smoothing of the local path, V_iThe method for finding the new position is as follows:

if node V_iIf the vertex belongs to the mesh, traverse the node V_iA ring neighborhood point set E_i{X₁,X₂,…,X_qAt each point in the } a vector V is calculated_iX_iAnd vector V_i-1X_iAngle and vector V between_iX_iAnd vector V_i+1X_iAngle therebetween, the sum of these two angles theta_iThen, theta is determined_iIf the value is less than pi, obtaining a node V if the value is less than pi_iCrossing edge V_iX_iSet C of_iThen construct a transit node V_i-1And V_i+1And parallel to a given plane normal vector

Then, the plane Pane and the set C are calculated_iAnd finally, constructing a new local path according to the calculated intersection point: if plane Pane and set C_iHas an intersection at the part of the intersecting edgesPoint, then calculate set C_iMiddle and Pane non-intersection point intersecting edge non-node end point and node end point V_iThe distance from the plane Pane is determined by taking the end point nearest to the plane Pane as the intersection point of the intersection edges and the plane, and taking the node V_i-1And V_i+1And set C_iThe intersection points of all the intersected edges and the Pane are sequentially connected to serve as a new local path; if plane Pane and set C_iAll the intersecting edges in the node V have intersection points, the node V is connected_i-1And V_i+1And set C_iThe intersection points of all the intersected edges and the Pane are sequentially connected to serve as a new local path; otherwise, directly connecting the node V_i-1And V_i+1Connected as a new local path, where q is E_iTotal number of points of (2), X_iIs a V_iCorresponding one-ring neighborhood point, set C_iElement in accordance with theta_iOrdering from small to large, given the normal vector of the plane

if node V_iPoints belonging to the edges of the grid, in node V_iThe edge is used as an intersecting edge, and a passing node V is constructed_i-1And V_i+1And parallel to a given plane normal vector

The value of (d) is the weighted average of the trigonometric vectors of the intersecting edges; then, the node V is connected_iDeleting the Path from the characteristic Path, and adding the intersection points to the characteristic Path in sequence at the current positionIn Path, finally, the AND node V_iAnd (3) carrying out corresponding updating on related nodes and weight information thereof: first, a node V is deleted from the set W of values_iAssociated point set Connect (V)_i-1,V_i,V_i+1) The fairing weight of each node in the system is updated, and then the node V is updated_iIs concentrated in two adjacent nodes V_i-1And V_i+1Corresponding associated point set and fairing weight, and finally connecting (V) the updated associated point set_i-1,V_i,V_i+1) Storing the corresponding fairing weight value into a weight value set W, wherein the updated association point set Connect (V)_i-1,V_i,V_i+1) The corresponding calculation formula of the fairing weight value is as follows:

2. The method for extracting the characteristic line of the dental preparation restoration model according to claim 1, wherein: the step S1 includes:

s13 at V_topAs the current search center point, find out V_topOne-ring neighborhood point set Neibourrhood { V₁,V₂,…,V_oAnd the associated container NodeMap is emptied, then each point in neibourwood is traversed,storing points which are not in the candidate queue NodeQueue in the Neibourhood and corresponding cost information into the associated container NodeMap, and finally returning to the step S12, wherein o is V_topTotal number of neighborhood points of a ring;

3. The method for extracting the characteristic line of the dental preparation restoration model according to claim 2, wherein: the step S11 includes:

4. The method according to claim 3, wherein the method comprises: f in the comprehensive constraint control function f (g)_cAnd f_tThe calculation process of (2) is as follows:

wherein A, B, C, D, E, F and G are both coefficient values of f (x, y),

for Weingarten curvature matrix W₁Singular value decomposition is carried out to obtain the maximum principal curvature and the minimum principal curvature of the vertex of the current meshTwo eigenvalues corresponding to the principal curvatures, and then two eigenvectors t corresponding to the two eigenvalues₁And t₂Determining a maximum principal direction t of a current mesh vertex_maxAnd a minimum principal direction t_minWherein, t₁、t₂、t_maxAnd t_minAre respectively: t is t₁＝(t₁₁,t₁₂)^T，t₂＝(t₂₁,t₂₂)^T，t_max＝t₁₁x+t₁₂y，t_min＝t₂₁x+t₂₂y；

5. The method according to claim 3, wherein the method comprises: the step S13 includes:

6. The method for extracting the characteristic line of the dental preparation restoration model according to claim 5, wherein: the step S132 includes:

judging the current candidate point V_nWhether the queue is in the candidate queue NodeQueue or not is judged, if yes, the current queue NodeQueue is indicatedCandidate point V_nThe candidate points are searched as candidate points before, no processing is carried out on the candidate points at this time, and the next candidate point in Neiburhood is continuously selected to calculate the cost information of the candidate points; otherwise, then V is set_nAdding into a candidate queue NodeQueue and ordering F_nAfter Cost, V is finally added_nAnd F_nAnd storing the data into the associated container NodeMap.

7. The method for extracting the characteristic line of the dental preparation restoration model according to claim 1, wherein: the step S2 specifically includes:

and judging whether the initial characteristic path between the key characteristic points has deviation with an expected characteristic line, if so, editing the initial characteristic path between the key characteristic points by adopting point editing operation, otherwise, executing the step S3, wherein the point editing operation comprises but is not limited to point deleting operation, point inserting operation and point moving operation.