Disclosure of Invention
According to the problems in the prior art, the invention discloses a modeling method for automatic registration of curved bone fracture and pre-bending of an internal fixation steel plate, which specifically comprises the following steps:
s1: extracting the main direction of a broken bone model point set by adopting a PCA algorithm, selecting clustering seed points along the main direction, clustering the clustering seed point set of the model by adopting a K-Means clustering method to calculate the central point of each cluster, and fitting a Bezier curve by utilizing the central points of the clusters to obtain the trend curve information of the broken bone model;
s2: constructing capsule-shaped mapping models at two ends of the fractured bone model respectively, and mapping points in the fractured bone model onto the mapping models along the normal vector direction of the points;
s3: clustering the feature point set of the mapping model by adopting an OPTIC algorithm to the mapping model constructed in S2, and performing feature matching on the feature point set of the mapping model to realize coarse registration of the two mapping models;
s4: extracting and segmenting a point set of a section part of the fractured bone model: carrying out regional self-growth on the feature points in the coarse registration to obtain an expansion point set, and carrying out matching search in the expansion point set to obtain a cross-section point set;
s5: and accurately registering the fractured bone model according to the two sectional point sets obtained by segmentation: extracting the main direction of the cross-section point set by using a PCA algorithm, performing pre-registration according to the main direction, and performing fine registration by using an ICP algorithm;
s6: and selecting a control point at the fracture part of the fractured bone model, and thickening the triangular patch within the range of the control point to obtain the steel plate model.
Further, the following method is specifically adopted in S2:
calculating tangent lines at two ends of the curve according to the trend curve of the fractured bone model obtained by fitting, respectively constructing four capsule-shaped mapping models by taking the four tangent lines as axes, and respectively mapping point sets in the range of the mapping models onto the mapping models along the normal vector direction of the point sets because the four mapping models have the same size, so as to obtain four mapping models containing the mapping point sets.
Further, the coarse registration process of the two mapping models in S3 specifically adopts the following manner:
performing mirror transformation on one of the two mapping models, matching the two mapping models in a coordinate system, performing grid division on the surface part of a cylinder in the mapping models, judging the attribution of grids of a point set of the cylinder part in a mapping point set, dividing all the point sets of the cylinder part into the grids, counting the number of mapping points in each grid, extracting all the points in the grids with the number larger than a threshold value as feature points of the mapping point set, clustering the feature point set by adopting an OPTIC algorithm to extract a clustering center, traversing the clustering center point, calculating the matching degree between clusters, extracting two clustering point sets with the maximum matching degree and obtaining a rotation angle, and rotating a broken bone model according to the rotation angle to achieve coarse registration.
Further, in S4, the following method is specifically adopted for the extraction and segmentation process of the point set of the broken bone model section part: and reversely mapping the two clustering point sets extracted in the rough registration, finding a point set corresponding to the two clustering point sets on the surface of the fractured bone model as a section characteristic point, carrying out regional self-growth on the section characteristic point set to obtain two extension point sets, and searching for a matching point according to the coordinate relation and the normal vector relation between the two points to obtain the point sets of the two sections.
By adopting the technical scheme, the automatic registration of the bent bone fracture and the pre-bending modeling method of the internal fixing steel plate can be realized, the steel plate is fitted according to the bent bone model after the registration, the algorithm provided by the invention can be used for full-automatically registering the bent bone, the robustness is high, the method can be suitable for various complicated sections, the precision is high, and the steel plate model required by the operation can be obtained through the pre-operation fitting.
Detailed Description
In order to make the technical solutions and advantages of the present invention clearer, the following describes the technical solutions in the embodiments of the present invention clearly and completely with reference to the drawings in the embodiments of the present invention:
the modeling method for automatic registration of the bent bone fracture and pre-bending of the internal fixation steel plate as shown in fig. 1 specifically comprises the following steps:
as shown in fig. 2: s1: and extracting the main direction of the broken bone model point set by utilizing a PCA algorithm, selecting clustering seed points along the main direction, clustering the clustering seed point set of the model by adopting a K-Means clustering method to calculate the central point of each cluster, and fitting a Bezier curve by utilizing the central points of the clusters to obtain the trend curve information of the bent bone model. The concrete mode is as follows: and extracting the main direction of the broken bone model point set by using a PCA algorithm. The extraction effect is shown in fig. 3, PCA is a forming algorithm for extracting the principal direction of a point set, and the method specifically comprises the following steps:
calculating the coordinates of the center point of the point set: respectively calculating the average value of X, Y, Z coordinates of all points in the fractured bone model to obtain the coordinate p of the central point of the fractured bone point set0(x0,y0,z0)。
Feature centralization: constructing the space coordinates of the broken bone point set into a 3 multiplied by n matrix
n is the total number of points, and the coordinate of each point in the matrix and the coordinate p of the central point are respectively calculated
0To obtain an updated matrix a, i.e.
Calculating a covariance matrix: multiplying the matrix A and the transposed matrix A 'thereof to obtain a covariance matrix M, namely M-AA'
Solving an eigenvalue and an eigenvector of the covariance matrix M: and the vector corresponding to the maximum characteristic value is the main direction of the model point set, namely the direction of the extracted straight line.
Equidistant 5 seed points are selected along a straight line, and then K-Means clustering is carried out by taking the 5 seed points as central points. The clustering effect is shown in fig. 4. K-Means is a classic clustering algorithm, and the specific steps are as follows:
and respectively establishing 5 clusters by taking the 5 seed points as central points, calculating all points in the point set to obtain the seed point closest to the seed point, and adding the point into the cluster where the seed point is located.
After all the points are classified, the central points of the 5 classes are recalculated and used as new seed points to perform the next iteration. And stopping iteration until the set iteration times are reached.
The center points of the clusters are used to fit a Bezier curve. Calculating a curve expression, wherein a parameter formula of the n-order Bezier curve is as follows:
wherein t is a parameter with a value range of [0,1 ], n is the order of the Bezier curve, P
iRepresenting the coordinates of the ith central point.
And (4) performing interpolation, namely performing interpolation on the curve, namely calculating to obtain the coordinates of a plurality of points on the curve by giving different t, and performing approximate fitting to obtain a complete curve. The final effect is shown in fig. 2.
S2: capsule-shaped mapping models are respectively constructed at two ends of the fractured bone model, and points in the fractured bone model are mapped onto the mapping models along the normal vector direction of the points. The concrete mode is as follows: calculating tangent lines at two ends of the curve, and moving the broken bone model when each model is constructed to facilitate construction of each mapping model, so that the end point of the curve is positioned at the origin of coordinates, and the tangential direction is the positive direction of the x axis.
The shape of the constructed mapping model is similar to that of a capsule, the middle part is cylindrical, the two ends are hemispherical, the geometric center of the model is the end point of a broken bone model curve, and the axis of the model is the tangent line obtained by the calculation. Setting the radius of a transverse circle of the cylindrical part of the mapping model as r and the length of the cylindrical part as l, and then setting the coordinate equation of the mapping model as
The set of points is mapped. And (3) mapping the point sets in the model range along the direction of the normal vector at the position of the point sets respectively, namely calculating the intersection point of the straight line where the normal vector is located and the mapping model, taking the intersection point as the mapping point of the point, repeating the steps, and constructing to obtain 4 mapping models in total, as shown in figure 5.
S3: and clustering the feature point set of the mapping model by adopting an OPTIC algorithm to the mapping model constructed in the S2, and performing feature matching on the feature point set of the mapping model to realize the coarse registration of the two mapping models. The specific process is as follows: performing mirror transformation on one of the two mapping models, matching the two mapping models in a coordinate system, performing grid division on the surface part of a cylinder in the mapping models, judging the attribution of grids of a point set of the cylinder part in a mapping point set, dividing all the point sets of the cylinder part into grids, counting the number of the mapping points in each grid, extracting all points in the grids with the number larger than a threshold value as feature points of the mapping point set, clustering the feature point set by applying an OPTIC algorithm, extracting a clustering center, traversing the clustering center points, calculating the matching degree between clusters, extracting two clusters with the maximum matching degree, acquiring a rotation angle, and rotating a broken bone model according to the rotation angle to achieve coarse registration. The specific algorithm comprises 4 steps:
s31: and (5) mirroring. Because the two matched models are mirror-symmetrical in the initial construction, mirror processing is required to be performed first for the convenience of subsequent matching. In the method, the geometric centers of the mapping models are all at the origin, so that the coordinates of the X axis of each point in the mapping point set only need to be changed into opposite numbers during the mirror image processing.
S32: and extracting the characteristic points. The surface portion of the cylinder in the mapping model is gridded and divided into n along the height (parallel to the X-axis direction) of the cylinder in the longitudinal direction
xA section divided into n along the circumference of the side surface
yIntervals, so that the total number of meshes is N
grid=n
x×n
y. Then, the grid attribution of the cylindrical part in the mapping point set is judged, and the number of the grid to which the mapping point belongs is set as (x, y, z) if the space coordinate of the mapping point is (x, y, z)
l is the length of the cylinder, and theta is a positive angle taking the positive direction of the Y axis as a starting edge and taking the connecting line of the point and the circle center of the cross section as a final edge, as shown in figure 7. Dividing the point set of the cylindrical part into grids, counting the number of mapping points in each grid, wherein the number is greater than a threshold value (a)Set as 3), all points in the grid are extracted as feature points of the mapping point set, and the extraction effect is shown in fig. 8.
S33: OPTICS clustering. And (5) clustering the feature point set by applying an OPTIC algorithm. The concrete steps
(1) Traversing the feature point set, searching the adjacent point (the point with the distance less than the threshold epsilon) of each point, counting the number of the adjacent points corresponding to each point, if the number is more than the threshold M, the point is a core point, calculating the core distance of the point, and defining the core distance of the point as p
Wherein N is
ε(p) represents a set of proximate points for point p,
represents N
ε(p) the point nearest to point i in p.
(2) Two queues are established, an ordered queue and a result queue. Traversing all points in the point set, selecting an unprocessed core point to be added into the result queue, and adding all the adjacent points of the point into the ordered queue. Calculating the reachable distance of each point in the ordered queue from the core point, from p0To p1Is defined as
Arranging the points in the ordered queue according to the reachable distance in an ascending order, traversing the points if the ordered queue is not empty, taking out one point each time, adding the point into the result queue if the point is the core point and is not in the result queue, adding the point close to the point into the ordered queue, updating the reachable distance if the point close to the point is already in the ordered queue and has a smaller reachable distance, and finally reordering the ordered queue. And repeating the process until all the feature points are traversed.
(3) And outputting a clustering result. Firstly, setting a threshold radius epsilon ', sequentially taking out a point from a result queue, and if the reachable distance of the point is not more than epsilon', the point belongs to the current category; if the reachable distance of the point is greater than epsilon' and the point is a core point, a new cluster is created; if the reachable distance of the point is greater than ε' and is not a core point, then the point is a noise point. The results of the clustering are shown in fig. 9.
S34: and extracting a clustering center. Firstly, merging the OPTIC clustering results by adopting a hierarchical clustering method, merging the categories with closer distances, and deleting the categories with less element numbers. The specific method comprises the following steps: traversing the central points of the OPTICS clusters according to the numbering sequence, traversing the central points of the subsequent classes for the central point of each class, and merging the two classes if the distance between the two central points is less than a threshold value (set as 4.0). After clustering is completed, the categories with the element number smaller than the threshold value (set as 20) are traversed and removed, and finally the updated clusters are numbered again. After hierarchical clustering is completed, the coordinate center point of each category is calculated, and the result is shown in fig. 10.
S35: and carrying out feature matching on the mapping point set. The method comprises the following specific steps:
(1) the vector from the center point of the transverse circle where the clustering center point is located to the center point is
Will be provided with
And rotating counterclockwise around the X axis to the positive direction of the Z axis by an angle alpha, calculating the corresponding alpha of the central point of each cluster, sequencing all alpha values, removing coupling, and adding the alpha values into an angle queue Q.
(2) When comparing the feature point sets of two mapping models, the cross section type is determined according to the number of the point sets of the head (end point spherical part) of the mapping model. If the number of head point sets of both mapping models is greater than the threshold (set to 300), the following method is performed: traversing the head point set in the mapping model II and for each point p
iFinding the nearest point in the head point set of the mapping model I and calculating the nearest distance s
i,s
iThe shortest distance on a sphere is not the shortest distance in space. And introducing a matching distance S to represent the matching degree of the two mapping point sets, wherein the smaller S represents the higher matching degree. In the method of studying a set of head points,
wherein n is
hThe number of points in the map point set is two. If the head point sets of the two mapping models are not both larger than the threshold value, analyzing the cylindrical part of the mapping models, wherein the method comprises the following steps: the clustering central point angle queues of the two mapping models are respectively Q
1And Q
2To Q, pair
1Q
2All angle combinations in (a) are traversed, and the currently traversed angle combination is assumed to be (alpha)
i,β
j) In which α is
i∈Q
1,β
j∈Q
2Then the mapping model is rotated counterclockwise by α
iRotate the mapping model two counterclockwise by beta
jAnd will be alpha
i,β
jAnd setting the cluster C, C' where the corresponding central point is as the current cluster. If the number of the point sets in C and C' exceeds the threshold value (set to 500), traversing the cylindrical part point set in the mapping model II, and for each point p
i' finding its neighborhood points in the cylindrical part point set of the mapping model one can improve the retrieval efficiency by using the mesh model established in the previous step. And respectively calculating the corresponding matching distance S for each angle combination. Firstly, judging the condition of high matching of C and C', the calculation formula is
Wherein a is
1Represents the number of matching points in C' in C, b
1Representing the number of all point sets with matching points in C, a
2b
2The same is true. If S'<0.1, then S ═ S'. Otherwise
m1And m2Are respectively provided withRepresents the total number of points in C and C'.
(3) Traversing all the angle combinations and calculating the matching distance S, wherein the minimum angle combination corresponding to S is the optimal registration angle of the two mapping models. And performing pairwise matching analysis on 4 mapping point sets of the two broken bone models, and performing 4 groups of matching analysis, wherein the two mapping point sets with the minimum matching distance are the mapping point sets of the two sections, and the corresponding registration angle is the angle of coarse registration.
(4) And (4) rotating the broken bone model anticlockwise around the X axis according to the angle obtained in the step (3), and then properly translating to enable the coarse registration result to be more accurate. The distance of translation is determined by the number n of head point sets of the two mapping modelsh1And nh2To calculate.
When t is less than or equal to 0, no movement is needed;
when t is>At 0, a translation vector v is calculatedt=[t*0.016 0 t*0.008]T。
Pressing the broken bone model twice according to vtThe coarse registration can be achieved by performing the translation, and the result of the coarse registration is shown in fig. 6.
S4: the point set of the section part of the broken bone model is extracted and segmented as shown in fig. 7-13: and carrying out regional self-growth on the feature points in the coarse registration to obtain an extended point set, and carrying out matching search in the extended point set to obtain a cross-section point set. In this step, on the basis of the rough registration, the point set of the cross section of the fractured bone model is accurately extracted and segmented, and the segmentation result is shown in fig. 11. The method comprises the following specific steps:
s41: the cross-sectional feature points undergo regional self-growth. In the course of coarse registration, two sets of clustering points C and C' that are optimally matched can be obtained, and in this step, these two sets of clustering points are reversely mapped back to the fractured bone model (as shown in fig. 6), so that feature point sets F of two cross sections can be obtained1And F1' then, these two feature point sets are subjected to a plurality of rounds (set as 5) of region self-growth to obtain an expanded point set F2And F2', as in FIG. 12.
S42: to F2And F2' conducting a match search. The search strategy is: traverse F2Calculating the included angle beta between the normal vector at the point and the tangent line at the section if cos beta>0.35 directly extracting the point to a section point set, otherwise setting a coordinate threshold value r120.0 and vector threshold r21.0, at F for data point p of the current study2In' go through, if F2' the point p ' in the above satisfies the following three conditions at the same time, and p ' are matching points:
the squared distance of the coordinates is less than a threshold value, i.e. (x)1-x2)2+(y1-y2)2+(z1-z2)2<r1
②(xp+xp′)2+(yp+yp′)2+(zp+zp′)2<r2,
(xp,yp,zp)(xp′,yp′,zp') is the normal vector at points p and p', respectively
And the included angle of the normal vectors at the points p and p' is gamma, and cos gamma is greater than 0.
If at least 5 p' can be found for the point p to satisfy the above condition, then p is added to the cross-section point set.
S43: respectively make F2And F2Taking the two sets of point points as a target point set to perform two traversal searches, and finally extracting to obtain two cross section point sets P1And P2The extraction results are shown in fig. 11.
S5: and accurately registering the fractured bone model according to the two sectional point sets obtained by segmentation: and extracting the main direction of the cross-section point set by using a PCA algorithm, performing pre-registration according to the main direction, and performing fine registration by using an ICP algorithm.
Pre-registration process: the pre-registration aims to enable the sections to be approximately matched, so that the subsequent fine registration is facilitated, and the method comprises the following steps: firstly, extracting the main directions of two cross-section point sets by utilizing a PCA algorithm, taking the increasing direction of an X axis as the positive direction, and taking the two main directionsThe direction vectors are all adjusted to positive directions. Then extracting the projection vectors v of the two main direction vectors on the y-o-z plane1v2Fixing the fractured bone 2, rotating the fractured bone 1 with the principal axis as the axis until v1v2Parallel, the two sections can then be made to substantially coincide in the direction of the main axis. The result of the pre-registration is shown in fig. 13.
And (3) fine registration process: on the basis of pre-registration, accurate registration is carried out by utilizing an ICP algorithm, and the specific process is that (1) two part point sets are respectively marked as U and P. (2) Calculating the closest point, that is, for each point in the set U, finding the corresponding point closest to the point in the set P, and setting the new point set consisting of the corresponding points in the set P as Q ═ QiI is 0,1, 2. (3) Using the least mean square method, the registration between the sets of points U and Q is computed such that a registration transformation matrix R, T is obtained, where R is a 3 x 3 rotation matrix and T is a 3 x 1 translation matrix. (4) Calculating coordinate transformation, i.e. for set U, using registration transformation matrix R, T to make coordinate transformation to obtain new point set U1I.e. U1RU + T. (5) Calculate U1The root mean square error between Q is ended if the root mean square error is less than a preset limit value e, otherwise, a point set U is used1And replacing U, and repeating the steps. Fig. 14 shows the effect of the point set fine registration, and fig. 15 shows the result of the final registration of the two fractured bone models.
S6: and selecting a control point at the fracture part of the fractured bone model, and thickening the triangular patch within the range of the control point to obtain the steel plate model. And (3) clicking near the spliced fractured bone section by a user as shown in fig. 16, determining the approximate shape and size of the steel plate model, recording the clicked triangular plane value, selecting all surface triangular patches in the range, calculating the normal magnitude of each triangular patch, and recording the normal magnitude. Thickening each plane to a certain extent according to the normal vector direction of each plane, and filling the gap position of each plane. The obtained thickened part is the three-dimensional data of the simulated steel plate model, and can be exported and output as a result. Fig. 17 shows the effect of steel plate fitting.
The method applies various algorithms such as curve fitting, clustering, characteristic point mapping and the like, can efficiently and accurately realize automatic splicing of the bent bones, realizes full-automatic registration and steel plate pre-bending of a bent bone model, and does not need manual operation.
The above description is only for the preferred embodiment of the present invention, but the scope of the present invention is not limited thereto, and any person skilled in the art should be considered to be within the technical scope of the present invention, and the technical solutions and the inventive concepts thereof according to the present invention should be equivalent or changed within the scope of the present invention.