Method for extracting curve skeleton line of three-dimensional tubular body model based on hierarchical graph

Abstract

This invention proposes a method for extracting the curve skeleton line of a three-dimensional tubular body model based on a hierarchical graph. The method includes: inputting three-dimensional tubular body model data; preprocessing and extracting endpoints; calculating the shortest path, calculating the geometric features between different paths, and constructing initial discrete skeleton points; calculating the distance and angle of the paths to determine the positions of candidate bifurcation points; segmenting the original data; applying spherical B-splines for fitting and optimization to obtain skeleton line data; outputting the skeleton line data, estimating the position of bifurcation points by tracking the shortest path, and using an elastic ball collision movement algorithm to obtain the precise position of the bifurcation points; and using a dynamic ball expansion algorithm with the bifurcation points as the sphere centers to divide the data model into multiple branches, laying the foundation for subsequent hierarchical graph construction and parallel computation.

Description

Technical Field

[0001] This invention belongs to the technical field of computer graphics, specifically relating to a method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph. Background Technology

[0002] Skeleton line extraction of 3D volumetric data models is an important research topic in computer graphics. The 3D curve skeleton is a one-dimensional representation of the 3D model and has topological isomorphism. Using skeleton lines can achieve dimensionality reduction of the model and reduce complexity and redundancy.

[0003] Skeletonization provides an effective and compact method to obtain the skeleton of an object by reducing its dimensionality, while preserving the topological structure and geometric properties of the data. Since Blum proposed the concept of the central axis in the 1960s, the academic community has made many expansions in skeleton line extraction. The idea of ​​skeletonization is to reduce the size of objects, simplify two-dimensional objects into one-dimensional curves, and three-dimensional objects into two-dimensional or one-dimensional curves, and make the skeleton have properties such as centering, smoothness, and topological isomorphism.

[0004] Skeleton algorithms vary considerably, but from a computational and strategic perspective, they can be broadly categorized into five types: The first type is morphological thinning algorithms. This method is based on iterative erosion propagation using a Blum-like burning model simulated in digital voxels until a one-dimensional skeleton line is obtained. This type of method is simple to implement and performs well with smooth data, but it cannot guarantee the single-pixel integrity, smoothness, and centering of the skeleton, and is sensitive to noise, easily generating irrelevant branches. The second type is distance field methods. This method calculates the distance from all points to the boundary, then finds local ridges, connects all ridges, and prunes to obtain the skeleton line. This method is mainly suitable for 3D graphics of tubular objects. For complex graphics, the obtained skeleton has poor centering and smoothness, but this method is computationally very efficient, can reconstruct the original object, and is simple to implement, possesses translation and rotation invariance, and can guarantee topological connectivity. Therefore, it is often used in the forward work of some algorithms. The third type is geometric approximation methods, which mainly include… Voronoi diagrams and Reeb diagrams are used to calculate the central plane of the model and obtain the skeleton lines. This type of method can directly generate one-dimensional skeleton lines without relying on the central plane. However, due to the approximation by geometric methods, the computational complexity is high, and the time consumed in extracting the skeleton of complex models is long, so it has not been widely used. The fourth type is the field function method, such as forming an electric field or potential energy field at the boundary of the object and extracting the skeleton based on the distribution of the electric field inside the object. This method averages the influence of noise points, is not sensitive to boundary noise, and the obtained skeleton points are relatively continuous. Another drawback of this method is that the calculation process involves first and second derivatives, which makes the calculation unstable, so it cannot be widely used. The fifth type is to approximate the skeleton lines by connecting the ridge points of the data model using the shortest path algorithm. The main advantage of this type of method is that the shortest cost path has scale independence and has good performance in noise detection. At the same time, the algorithm is easy to implement. However, this method has the characteristics of poor centering and poor smoothness.

[0005] Since there is no unified and universally accepted mathematical definition for skeleton lines, and there are differences between different definitions, a large number of algorithms have been developed to optimize and iterate on different skeleton line definitions. Choosing a suitable skeletonization method for a specific application is a huge challenge. However, in most applications, people want the skeleton of the model to be robust, smooth, and centered under different digitization and imaging conditions, and to be able to achieve skeleton path tracking, allowing for acceptable reconstruction of the original object. However, previous algorithms have certain limitations and do not fully possess all of the above properties.

[0006] To address the shortcomings of the above algorithms, this invention proposes a method for extracting the curve skeleton line of a three-dimensional tubular body model based on a hierarchical graph. Summary of the Invention

[0007] This invention proposes a method for extracting the curve skeleton line of a three-dimensional tubular body model based on a hierarchical graph. This algorithm can maintain the skeleton line's smoothness, centering, and single-pixel characteristics.

[0008] The present invention adopts the following technical solution:

[0009] A method for extracting the curve skeleton line of a 3D tubular model based on a hierarchical graph, comprising:

[0010] Step 1, input the 3D tubular model data O;

[0011] Step 2: Preprocess the input model data O, and then extract the endpoints VE of the three-dimensional tubular body;

[0012] Step 3: Calculate the shortest path between any endpoints p, q∈VE on the three-dimensional tubular body, calculate the geometric features between different paths, and construct the initial discrete skeleton points;

[0013] Step 4: Calculate the distance and angle between two paths with the same starting point but different ending points on the three-dimensional tubular body to determine the location of candidate bifurcation points;

[0014] Step 5: Optimize and adjust the position of the bifurcation point using the dynamic ball expansion algorithm to obtain the final bifurcation point; establish a hierarchical topology, and apply a shrinkage algorithm to segment the original data based on the topology;

[0015] Step 6: For each segment of data on the three-dimensional tubular body, obtain the initial skeleton line through the distance transformation field; apply spherical B-splines for fitting and optimization to obtain the skeleton line data;

[0016] Step 7: Output skeleton line data.

[0017] Furthermore, in step 2, a neighborhood exploration algorithm is used to extract the endpoints of the three-dimensional tubular structure:

[0018] Calculate the number of voxels contained in the foreground points of a voxel's 6-neighborhood and 26-neighborhood. Let m be the number of foreground points in the 6-neighborhood of the object voxel, and n be the number of foreground points in the 26-neighborhood of the object voxel. The candidate point set of endpoints is represented as:

[0019] {x|N6(x)≤m}∩{x|N 26 (x)≤n}......(1),

[0020] Where x represents a data point in model data 0, and N6(x) and N 26 (x) are the 6-neighborhood and 26-neighborhood of data point x, respectively. The candidate endpoint set is obtained through the expression in equation (1).

[0021] The K-means method is used to cluster candidate points at the endpoints, defined as follows:

[0022]

[0023]

[0024] Among them, E b μ is the inter-class distance. i μ j Let E be the mean of the i-th and j-th clusters, respectively, and k be the number of clusters. ω x is the intra-class distance. i μ represents each data point. x The cluster mean is used, and the loss function is:

[0025] E = E b +E ω ......(4);

[0026] Among them, the inter-class distance and intra-class distance determine the average profile coefficient. The clustering effect is evaluated by the average profile coefficient. The number of clusters k in equation (2) is calculated based on the average profile coefficient.

[0027] The K-means algorithm is used for iterative calculation until the termination condition is met.

[0028] Using the center of each cluster as an endpoint voxel, we form the endpoint set V. E .

[0029] Furthermore, in step 4, for any two shortest paths with the same starting point but different ending points: π1 = <p0,p1,p2,...p l-1 , >π2= <q0,q1,q2,...q l-1 The two discrete paths are fitted using B-spline curves, denoted as C1(t) and C2(t). In C2(t), find a point q′ that minimizes the distance of the path p′q′, where D represents the minimum distance between two points, and θ is the angle between the tangent vectors of the corresponding points of the projected point and the original point. When D and θ are greater than a given threshold and When the time is reached, stop the iteration, determine the candidate points of the current p′q′ point as the bifurcation point, and successively obtain the candidate points of any two paths with the same starting point but different ending points to obtain the candidate point data of the bifurcation point;

[0030] K-means clustering is used to cluster candidate path points, and each cluster center is used as a candidate point for branching points. iTo expand the ignition point range, if the curves of different ignition points collide, a neighborhood relationship is established between the two ignition points. All candidate points with a degree less than 3 are deleted, and the remaining points constitute the candidate points in the current level. All bifurcation points in the outermost layer are found, and the bifurcation points in the next layer are iteratively obtained using the outermost bifurcation points as new endpoints until all bifurcation points are obtained, forming a set V of candidate bifurcation points. C .

[0031] Furthermore, in step 5, for each bifurcation point v i ∈V C The electric force at that point is defined as:

[0032]

[0033] Where n represents all the pairs of v i Distance ||p i -v i The number of points whose ||≤σ, where σ is a given threshold, s i Represents each point p on the boundary j Point to v i unit vector, p j and v i The potential energy between the two points is the Euclidean distance between them, where m represents the weight parameter. The ideal position of the bifurcation point is defined as the point where the internal field strength is 0. For each change of the candidate bifurcation point, the direction of movement of vi is calculated according to the electrostatic field model, and the distance of movement is calculated by the method of elastic ball, thus completing one movement.

[0034] Define the center of the elastic sphere as A sphere of radius r is:

[0035] (x-x0) 2 +(y-y0) 2 +(z-z0) 2 =r 2 ,......(6);

[0036] For each v i ∈V C Define σ(r,M(v) i ) represents the distance the ball's center moves after each collision. Where r is the current radius of the sphere, M(v i ) indicates that v i Let α and β be the number of voxels contained in the sphere at the center, and α and β be coefficients. As the radius increases, after the current sphere collides with the boundary, the electric potential field value at that point and the direction of movement are recalculated, and the new position of the sphere's center is determined. Defined as x1 = x0 + σs x y1=y0+σs y z1=z0+σs z Through continuous iteration until a stable state is reached, the ideal location of the bifurcation point is finally obtained. The center of each dynamic ball constitutes the final set V of bifurcation points. B .

[0037] Furthermore, in step 6, obtaining the initial skeleton line through the distance transformation field includes:

[0038] Calculate the minimum distance from each voxel to the boundary. Define d(p) as the distance value between each voxel and the boundary. The distance transformation of any voxel p is equivalent to a central sphere with radius d(p). A maximum central sphere centered on voxel p is defined as a sphere that can contain point p within the model. After calculating the distance transformation field, extract the center of the maximum central sphere to obtain the initial skeleton point.

[0039] The burning grass model is used to find the branch endpoints. The two sets covered_points and fire_front are calculated iteratively to simulate the fire transmission process. Here, fire_front represents the set of voxel points that the fire is currently burning, and covered_points represents the set of voxel points that the fire has already burned. All voxel points are included in the covered_points set as a termination condition. The time t when the outermost voxel point of covered_points is transmitted is recorded.

[0040] The initial skeleton points obtained by the distance field method and the endpoints obtained by the burning grass model method are used as the initial points of the skeleton lines of the final branches.

[0041] Furthermore, step 6 includes fitting the skeleton line using a B-spline curve:

[0042] Using cubic spline basis functions, discrete points are fitted. For m+1 data points q0, q1, ..., q... m Find a cubic B-spline curve:

[0043]

[0044] Satisfying p0 = q0, p1 = q m The remaining data points q i (i = 1, 2, ..., m-1) are calculated using the least squares method, and the optimization function is:

[0045]

[0046] The optimization function is about n-1 control vertices dj A minimum value of (j = 1, 2, ..., n-1), where Let be the parameter value for the data point, which is determined by accumulating the chord length parameter.

[0047]

[0048] Will Substituting the values ​​into the calculation, we get:

[0049]

[0050] Solving equation (10) is a least squares problem. To minimize the objective function f, we need to minimize the n+1 control vertices d. j The derivative of (j = 1, 2, ..., n-1) is 0, where d j The l-th derivative is:

[0051]

[0052] A system of n-1 equations with n-1 unknowns was obtained. The solution was obtained by Gaussian elimination, resulting in a cubic B-spline curve. This curve was expressed by parameters, and the discrete skeleton points were fitted into a skeleton line.

[0053] Furthermore, in step 6, spherical B-splines are applied for fitting optimization, and Hausdorff distance is introduced to optimize the original splines. The objective function for optimization is:

[0054]

[0055] d h Let represent the Hausdorff distance from all scattered points to the B-spline curve. The final optimization function is:

[0056]

[0057] Wherein, V(X) i (t, θ), Y i (t, θ), Z i (t, θ) represents the projection point of v(x, y, z) onto the surface of the sphere B spline;

[0058] The quasi-Newton algorithm is used to iteratively optimize and adjust the coordinates, number, and radius of the spheres controlling the sphere B-spline curve. Based on the error accuracy, the optimized sphere B-spline curve is obtained.

[0059] Furthermore, in step 1, the model data is a sequence of discrete points described using three-dimensional coordinates.

[0060] Furthermore, in step 2, the input model data O is preprocessed, including denoising and downsampling.

[0061] Furthermore, step 6 includes optimizing the shortest total time for all voxel points to be covered, thereby obtaining the outermost voxel point of this combustion as the endpoint. The optimized formula is as follows:

[0062]

[0063] Among them, L(p i ) indicates the distance from the ignition point to p i The time elapsed since the point was reached.

[0064] Compared with the prior art, the superior effects of the present invention are as follows:

[0065] 1. The method for extracting the curve skeleton line of a three-dimensional tubular body model based on a hierarchical graph, as described in this invention, transforms the skeleton extraction problem into a volume model equivalence class problem by leveraging the compactness characteristics of the volume model, thereby achieving topological isomorphism between the skeleton and the volume model;

[0066] 2. The curve skeleton line extraction method of the three-dimensional tubular body model based on the hierarchical graph described in this invention estimates the position of the bifurcation point by tracking the shortest path and obtains the precise position of the bifurcation point by using an elastic ball collision and movement algorithm; the data model is divided into multiple branches by a dynamic ball expansion algorithm with the bifurcation point as the ball center, laying the foundation for the subsequent construction of the hierarchical graph and parallel computing;

[0067] 3. The method for extracting the skeleton line of a three-dimensional tubular body model based on a hierarchical graph described in this invention first uses a B-spline approximation algorithm to fit discrete skeleton points to obtain an initial continuous skeleton line, and then applies a spherical B-spline algorithm to optimize the position of the fitted curve to make it have better centering, while obtaining the corresponding radius of the skeleton line. Attached Figure Description

[0068] Figure 1 This is a flowchart of the method for extracting the curve skeleton line of the three-dimensional tubular body model of the hierarchical diagram described in the embodiments of the present invention. Detailed Implementation

[0069] To better understand the above-mentioned objectives, features and advantages of the present invention, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be noted that, unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other.

[0070] Example

[0071] like Figure 1As shown, the method for extracting the skeleton lines of a 3D tubular body model based on a hierarchical graph extracts skeleton lines from the original data. In this embodiment, a divide-and-conquer approach is adopted to find the open cover of the original body model O, dividing the original data model into multiple branches K. i and satisfy and The collection of all branches constitutes the open cover of the original model. Skeleton lines are extracted in each branch, and finally the skeleton lines of the original data are formed.

[0072] For topology construction, this embodiment adopts the hierarchical graph method to ensure the topological isomorphism between the skeleton lines and the original data, and gives the corresponding definition, which leads to the definition of hierarchical graph;

[0073] Definition 1: G V (VE, W) represents the graph topology of the volume data, where V represents each voxel, E represents the connection between each voxel and its 26 neighbors, and W is a function defined on each edge, the value of which is the distance between the two points corresponding to that edge.

[0074] Definition 2: Define V S V is the set of all points on the skeleton lines, consisting of three types of points. B V E and V R Composition, V B A branching point is a node with three or more neighboring nodes, V E For an endpoint, the endpoint has only one neighbor node, V R For a normal node, there are two adjacent nodes. By definition, |V S |=|V R |+|V B |+|V B |;

[0075] Definition 3: Define an edge BE as the shortest path directly connecting a branch point and its endpoint, and an edge BB as the shortest path directly connecting any two branch points. Each branch is defined as the path corresponding to either edge BE or edge BB. V The volume data in the graph, therefore the branch graph G B It is defined as the set of all branches, that is

[0076] Definition 4: Define the graph topology of the skeleton lines S of the solid model O as G. C G C The vertex is v∈V S , side e i，j ∈{BB,BE}, representing v i v j They are adjacent, lv Defined as skeleton line graph topology G C The level of the midpoint is set to one, and the level of other points is the minimum level of all adjacent points plus one. e Defined as skeleton line graph topology G C The level of an edge is defined as the smaller of the levels of the two endpoints corresponding to that edge. For branching graph The level is the level of the skeleton line corresponding to the branch diagram;

[0077] Definition 5: Subgraph at level n It is the global graph G V a subset of in Define a hierarchy diagram A hierarchical diagram can be considered a partition of the original model, with the i-th layer being... This represents the set of all branches in the graph that are at level i.

[0078] In terms of topology, the topological structure constructed based on the above methods ensures the hierarchical graph. It is topologically isomorphic to the original data, thus ensuring the skeleton lines corresponding to the hierarchical graph. It is also isomorphic to the original data topology. In order to construct a hierarchical graph, an edge shrinking algorithm is used to recursively extract the branch graph by layering and splitting the data, defining... For skeleton topology Delete e ij This edge, simultaneously in the volume model Delete e ij The corresponding subgraph, after one contraction, will have all edges l at level 1. e Delete, and simultaneously split the corresponding branch. The resulting branches form a branch graph with a level of 1. Then, using the bifurcation point as the new endpoint, the skeleton lines are iteratively deleted. The corresponding edges and their corresponding branch graphs are obtained by deleting all branches, thus obtaining all the branch graphs of the original data, and these branch graphs preserve the topological structure of the original data.

[0079] Specifically, the method for extracting the curve skeleton line of the three-dimensional tubular body model based on the hierarchical graph includes topology construction and geometry optimization;

[0080] In topology construction, this embodiment preserves the topological structure of the original data and segments the original data by constructing a hierarchical graph, constructing an open cover of the original data. During the construction of the hierarchical graph, the endpoints and bifurcation points of the data model are extracted, specifically including:

[0081] For endpoint determination, this embodiment proposes two methods: neighborhood exploration and local segmentation-diffusion. The neighborhood exploration method uses N6, N 26 The neighborhood exploration algorithm extracts endpoints. For a voxel, it calculates the number of voxels contained in its 6-neighborhood and 26-neighborhood foreground points. Compared to ordinary points, the foreground points of voxels in the 26-neighborhood of an endpoint are necessarily incomplete and sparse. Let m be the number of foreground points in the 6-neighborhood of the object voxel, and n be the number of foreground points in the 26-neighborhood of the object voxel. The candidate points of the endpoints satisfy {x|N6(x)≤m}∩{x|N 26 (x)≤n}; In the experiment, the best results were obtained when the thresholds m=3 and n=9 were set. Through the neighborhood exploration method, many voxel points that meet the conditions will be selected in the endpoint part. In order to improve the accuracy and correctness of endpoint identification, the K-means method is used to cluster the candidate endpoint points. It is hoped that the intra-class distance of the sample is minimized and the distance between different clusters is maximized. Therefore, the following is defined: Where E b μ is the inter-class distance. i μ j This represents the mean among different clusters, where k is the number of clusters. Among them, E ω x is the intra-class distance. i μ represents each data point. x The cluster mean is given by the loss function E = Ec b +E ω Because the distance between clusters is relatively large, the parameter k can be easily and automatically calculated based on the average silhouette coefficient. Then, the K-means algorithm is used for iterative calculation until the movement of the cluster centers is less than the required error, reaching the termination condition. Finally, the center of each cluster is used as the endpoint voxel to form the endpoint set V. E The local segmentation diffusion method involves establishing local separators and determining endpoints based on the topological connections of these separators. The separator mentioned in this embodiment is a minimal separator; that is, if some data points are deleted from the current separator, its properties will no longer be satisfied. For the volume data of a tubular structure, the shape of the established local separator ∑ is similar to the inscribed disk of branched data. To establish a separator ∑, the boundary points of the disk are first sampled on the input voxel model. Then, the model expands in all directions centered on each sampled point. Multiple separators collide during expansion, eventually stopping. The boundaries of these collisions form the separators. However, the current separator is not the smallest separator. This embodiment explores ∑-v for each element in ∑ through connectivity tests between multiple separators. iIf the properties of the separator can still be guaranteed, iterative deletion is performed; otherwise, it is retained until the smallest separator is found. By setting different initial values ​​for the separators, an adjacency matrix for each individual separator is established, preserving the topological relationships between the separators. If the separator ∑ i , ∑ j By data grouping F k Connected, then we can in ∑ i and ∑ j An edge is established between the two separators, assuming they are adjacent, and their connection relationship is recorded in the adjacency matrix. Clearly, in the tubular model, the degree of the separator at the endpoint is 1, while the degree of separators at other locations will be greater than 1. Therefore, the endpoint separator is found by determining the degree of each separator in the adjacency matrix. This ensures that each segmented part F... i To make the sample size sufficiently small, the number of separators needs to be increased. Encryption is achieved by increasing the number of sampling points. After obtaining the separators and segmentation results at the endpoints, the Manhattan distance dist(x) between any two points in the separators is calculated. i -x j )=||x i -x j || 2 To obtain a new endpoint;

[0082] While both methods mentioned in this embodiment can acquire endpoints, their applicable conditions differ. The neighborhood exploration method is based on digital topology principles, traversing all voxel points according to neighborhood rules. This method requires high data integrity and smoothness, and parameter adjustments are needed for different types of data to obtain better exploration results. Its advantage is its computational simplicity and ease of implementation. The local segmentation diffusion algorithm, on the other hand, is based on a local segmenter, which is based on random sampling, generating a separator for each sampling point. This method has good robustness to different smoothness levels and different types of data, especially for tubular objects with large variations in thickness. Because this method is based on data point sampling, there may be no corresponding sampling points on branches with a small number of voxels, making it impossible to acquire the endpoints of the corresponding branches. To ensure that endpoints can be acquired on thin branches, the sampling density needs to be increased, which in turn increases the computational load.

[0083] The bifurcation point is crucial for preserving the topology. In this embodiment, the shortest path between any two endpoints is found, and the bifurcation point is estimated based on the shortest path to obtain the bifurcation point.

[0084] For any two shortest paths with the same starting point but different ending points, π1 = <p0,p1,p2,...p l-1, >π2= <q0,q1,q2,...q l-1 To facilitate the calculation of the projection and tangent vector between the two curves and reduce the calculation error between discrete points, B-spline curves are used to fit the two discrete paths, denoted as C1(t) and C2(t) respectively. In C2(t), find a point q′ that minimizes the distance of path p′q′. Then D represents the minimum distance between the two points. For finding the bifurcation point, consider the distance D between the two projected points, the angle cosθ between the tangent vectors of the corresponding points of the projected and original points, and the positional relationship between the line connecting the original and projected points and the current voxel background point. When D and θ are greater than a given threshold and When the time is reached, stop the iteration, determine the current p′q′ point as a candidate point for the bifurcation point, and then sequentially obtain the candidate points for any two paths with the same starting point but different ending points.

[0085] K-means clustering was applied to cluster all candidate path points, and each cluster center was used as a candidate for bifurcation points. Because the candidate bifurcation points were relatively close together, some erroneous clustering results existed. To remove these erroneous clustering structures, a simulated burning method was used, with each cluster center C... i To expand the ignition point range, if the curves of different ignition points collide, a neighborhood relationship is established between the two ignition points. All candidate points with a degree less than 3 are deleted, and the remaining points constitute the candidate points in the current level. Following this method, all bifurcation points in the outermost layer can be found. Similarly, using the outermost bifurcation points as new endpoints, the bifurcation points in the next layer are iteratively obtained until all bifurcation points are obtained, forming a set V of candidate bifurcation points. C ;

[0086] Based on the optimization of bifurcation points using dynamic spheres, it's important to note that bifurcation points obtained through K-means may not be centered within the voxel model, potentially resulting in edge-grabbing. Therefore, to ensure the centrality of bifurcation points, a dynamic collision sphere and electrostatic field method are employed for further optimization. The electrostatic field model assumes that the object's surface is filled with the same type of charge, ultimately forming a comprehensive electric field within the model. For each bifurcation point v... i ∈V C The electric force at that point is defined as: Where n represents all the pairs of v i Distance ||p i -v i The number of points whose ||≤σ, where σ is a given threshold, and s i Represents each point p on the boundary j Point to v i unit vector, p jand v i The potential energy between the two points is the Euclidean distance between them, where m represents the weighting parameter. Therefore, the ideal position of the final bifurcation point is defined as the point where the internal field strength is 0. For each change of the candidate bifurcation point, the direction of movement of vi is calculated according to the electrostatic field model, and then the movement distance is calculated by the method of elastic ball, thus completing one movement.

[0087] In the first case, simply increasing the radius of the sphere will transform the first case into the second, where the sphere is tangent to the boundary of the data. At this point, the position of the sphere's center needs to be adjusted, defined as... The sphere of radius r is (x-x0). 2 +(y-y0) 2 +(z-z0) 2 =r 2 For each v i ∈V C Define σ(r), M(v) i ) represents the distance the ball's center moves after each collision. Where r is the radius of the current sphere, M(v i ) indicates that v i Let α and β be the number of voxels contained in the sphere at the center, and α and β be coefficients. As the radius increases, after the current sphere collides with the boundary, the electric potential field value at that point and the direction of movement are recalculated, thus determining the new position of the sphere's center. It is defined as: x1 = x0 + σs x y1=y0+σs y z1=z0+σs z At this point, scenario two transforms into scenario one. Through continuous iteration, a stable state is reached, ultimately yielding the ideal location of the bifurcation point. The center of each dynamic sphere constitutes the final set V of bifurcation points. B ;

[0088] The hierarchical graph is constructed and the skeleton lines are initialized. Since the skeleton lines are roads connecting the points, there must be a shortest path between any two endpoints. Furthermore, paths with the same starting point but different ending points will always pass through the first branching point and then separate. Therefore, by finding the shortest path between any two endpoints, the outermost branching point can be obtained. Then, using the outermost branching point as the endpoint, the internal branching points are iteratively calculated. Based on the hierarchy of the branching points, the corresponding subgraphs and the hierarchical graph are constructed. Finally, the data is partitioned according to the topological relationships of the hierarchical graph, using a divide-and-conquer approach, and the skeleton lines are extracted in each branch.

[0089] By establishing a sphere separator with the bifurcation point as the sphere center, the initial radius of the separator is set to 1. The radius of the sphere is continuously increased until multiple branches are no longer connected. For volume data in the separator, it is simultaneously assigned to all adjacent branches to ensure the integrity of the branches. For each data branch, this paper uses the distance transformation field method to initialize the skeleton line. The distance transformation field is to calculate the minimum distance from each voxel point to the boundary. d(p) is used to define the distance value between each voxel point and the boundary. Therefore, the distance transformation of any voxel p can be interpreted as a central sphere with a radius of d(p). A maximum central sphere centered on voxel point p is defined as a sphere within the model that can contain point p. After calculating the distance transformation field, the center of the maximum central sphere is extracted to initialize the skeleton points.

[0090] For each segment after splitting, the branch k i It is necessary to calculate k i For the two endpoints, since the data after segmentation is no longer smooth (due to loss in the segmentation and uncertainty in the shape at the segmentation point), the neighborhood method cannot be used to determine the endpoints. For branches with relatively small curvature, the Euclidean distance between any two voxel points can be directly maximized, and then the two voxels with the farthest distance are used to approximate the two endpoints. To solve the problem of finding the endpoints of branches with large curvature, a method for finding branch endpoints is proposed based on the burning grass model proposed by Blum. The algorithm calculates the time for the entire branch to burn completely by taking each voxel point as the ignition point, and then optimizes the shortest time to determine the endpoints. The algorithm is mainly divided into two parts: one is to simulate the process of randomly selecting an ignition point on the model M and spreading it to the surrounding area; the other is the cutoff case, that is, to verify that all points have been transmitted to the end point. Part of the algorithm is based on breadth-first search, iteratively calculating the two sets, covered_points and fire_front, to simulate the fire transmission process. Here, fire_front represents the set of voxels currently burning, and covered_points represents the set of voxels that have already been burned. The termination condition is that all voxels are included in the covered_points set, meaning all points have been burned. At this point, the time t at which the outermost voxel in covered_points is transmitted is recorded. The purpose of the method described in this embodiment is to optimize the shortest total time for all voxels to be covered, thereby obtaining the outermost voxel as the endpoint of this burning event. The optimization formula is as follows: Among them, L(p i ) indicates the distance from the ignition point to p i The time elapsed at each point is used to determine the initial skeleton points obtained by the distance field method and the endpoints obtained by the fire method, which are then used as the initial points of the final branch skeleton lines.

[0091] In geometric optimization, after obtaining the branch graph and corresponding skeleton line of the original data, the initial skeleton line needs to be optimized to ensure its smoothness and centering. This process is divided into two parts. First, the B-spline fitting algorithm is applied to fit the original discrete skeleton line into a cubic spline curve to ensure the smoothness of the skeleton line. Then, the spherical B-spline optimization algorithm is applied to optimize the position of the curve to ensure the centering of the skeleton line.

[0092] B-spline curve fitting of the skeleton line is achieved through a least-squares approximation algorithm, smoothing the skeleton line and converting the discrete skeleton line into a parametric spline curve. A cubic spline basis is used as the basis function to fit discrete points. For m+1 data points q0, q1, q2...q... m Find a cubic B-spline curve Satisfying p0=q0,p m =q m The remaining data points q i (i = 1, 2, 3, ..., m-1) are calculated using the least squares method, and the optimization function is: This function is about n-1 control vertices d j A minimum value of (j = 1, 2, 3, ..., n-1) is obtained, thus yielding a system of n-1 equations with n-1 unknowns. This system is solved using Gaussian elimination, resulting in a cubic B-spline curve. This curve is expressed parametrically, fitting the discrete skeleton points into a skeleton line, thus ensuring the smoothness of the skeleton line.

[0093] Using cubic spline basis functions, we fit discrete points. For m+1 data points q0, q1, ..., q m Find a cubic B-spline curve:

[0094]

[0095] Satisfying p0 = q0, p1 = q m The remaining data points q i (i = 1, 2, ..., m-1) are calculated using the least squares method, and the optimization function is:

[0096]

[0097] This function is about n-1 control vertices d j A minimum value of (j = 1, 2, ..., n-1), where i = 0, 1, ..., m are the parameter values ​​for the data points. These parameter values ​​are determined by accumulating the chord length parameter. Let:

[0098]

[0099] Will Substituting the values ​​into the calculation, we get:

[0100]

[0101] To optimize this cost function using the least squares method, to minimize the objective function f, we need to minimize the n+1 control vertices d. j The derivative of (j = 1, 2, ..., n-1) is 0, and its first derivative is:

[0102]

[0103] Thus, a system of n-1 equations containing n-1 unknowns was obtained. The solution was obtained by Gaussian elimination. After the solution, a cubic B-spline curve was obtained. This curve was expressed by parameters and the discrete skeleton points were fitted into a skeleton line, which ensured the smoothness of the skeleton line.

[0104] To optimize the skeleton line of a sphere B-spline curve, Hausdorff distance is introduced to optimize the original spline in order to ensure good centering of the skeleton line in the data. The sphere B-spline algorithm consists of two parts: the first part calculates the initial radius of the sphere B-spline and fits the branch to the initial sphere B-spline curve; the second part calculates the average Hausdorff distance to optimize the skeleton line centering. Where d h (Q, (t) is the Hausdorff distance of the discrete point set, and the spherical B-spline curve can be expressed by the following two parts:

[0105]

[0106] The initial radius of the spherical B-spline curve is initialized using the elastic sphere method. First, the model boundary is determined. When calculating the boundary points of the voxel model, the 26-neighbor detection method is still used. For any branch volume data k obtained through segmentation... i Iterate through all voxel points v∈k i Then, calculate the 26-neighborhood of v. If all points in the 26-neighborhood are foreground points, then the point must be surrounded by other foreground points, meaning it is not a boundary point. Otherwise, if the 26-neighborhood contains background points, then the point is considered a boundary point, marked, and added to the boundary point set. After traversing all the data, k is obtained. i For all boundary data, with the points on the B-spline curve as the centers of the elastic sphere, the radius of the sphere is continuously expanded. By determining the positional relationship between the sphere and the boundary, it is known that the elastic sphere is tangent to the boundary, and thus the initial radius is determined, and the initial sphere B-spline curve is generated.

[0107] To ensure good centering of the skeleton lines in the data, Hausdorff distance is introduced to optimize the original splines. The objective function for optimization is:

[0108] Let represent the Hausdorff distance from all scattered points to the B-spline curve. The final optimization function is:

[0109]

[0110] Where V(X) i (t, θ), Y i (t, θ), Z i (t, θ) represents the projection point of v(x, y, z) onto the surface of the spherical B-spline. For the initial spherical B-spline curve, a quasi-Newton algorithm is used to iteratively optimize and adjust the coordinates, number, and radius of the control spheres of the spherical B-spline curve. Then, based on the error accuracy, the optimized spherical B-spline curve is finally obtained. Through the optimization of the spherical B-spline, the B-spline curve at this time has the characteristics of centering, smoothness, and single pixel. Finally, based on the topological structure at the time of segmentation, the skeleton lines of all branches are merged again to form the skeleton line S of the final original data model.

[0111] This invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the invention as claimed. The scope of protection of this invention is defined by the appended claims.

Claims

1. A method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph, characterized in that, include: Step 1, input the 3D tubular model data O; Step 2: Preprocess the input model data O, and then extract the endpoints VE of the three-dimensional tubular body; Step 3: Calculate the shortest path between any endpoints p, q ∈ VE on the three-dimensional tubular body, calculate the geometric features between different paths, and construct the initial discrete skeleton points; Step 4: Calculate the distance and angle between two paths with the same starting point but different ending points on the 3D tubular body to determine the location of candidate bifurcation points: Fitting two discrete paths using B-spline curves is represented as follows: ,for ,exist Find a point in Satisfy path The distance is minimized by , where D represents the minimum distance between the two points. Let D be the angle between the tangent vectors of the corresponding points of the projected point and the original point. Greater than a given threshold and When the iteration stops, the current iteration is set to... Candidate points for determining the bifurcation point are obtained by sequentially finding candidate points for any two paths with the same starting point but different ending points, thus obtaining candidate point data for the bifurcation point. The K-means clustering method is used to cluster the candidate path points, and each cluster center is used as a candidate point for branching points. To expand the ignition point range, if the curves of different ignition points collide, a neighborhood relationship is established between the two ignition points. All candidate points with a degree less than 3 are deleted, and the remaining points constitute the candidate points in the current level. All bifurcation points in the outermost layer are found, and the bifurcation points in the next layer are iteratively obtained using the outermost bifurcation points as new endpoints until all bifurcation points are obtained, forming a set of candidate bifurcation points. ; Step 5: Optimize and adjust the position of the bifurcation point using the dynamic ball expansion algorithm to obtain the final bifurcation point; establish a hierarchical topology, and apply a shrinkage algorithm to segment the original data based on the topology: For each branching point The electric force at that point is defined as: ……（5）； Where n represents all the pairs of... Distance || - || ≤ The number of points, For a given threshold, Represents each point on the boundary point to unit vector, and The potential energy between the two points is the Euclidean distance between them, where m represents the weight parameter. The ideal position of the bifurcation point is defined as the point where the internal field strength is 0. For each change of the candidate bifurcation point, the direction of movement of vi is calculated according to the electrostatic field model, and the distance of movement is calculated by the method of elastic ball, thus completing one movement. Define the center of the elastic sphere as A sphere with radius r is: ,……（6）； For each ,definition The distance the ball's center moves after each collision. Where r is the radius of the current sphere. Indicated by The number of voxels contained in the sphere at the center. , As a coefficient, after the current sphere collides with the boundary, the electric potential field value at that point and the direction of movement are recalculated, and the new position of the sphere's center is determined. Defined as , , Through continuous iteration until a stable state is reached, the ideal location of the bifurcation point is finally obtained. The center of each dynamic ball constitutes the final set of bifurcation points. ; Step 6: For each segment of data on the three-dimensional tubular body, obtain the initial skeleton line through the distance transformation field; The skeleton line data were obtained by fitting and optimizing the B-spline. Step 7: Output skeleton line data.

2. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 1, characterized in that, In step 2, the endpoints of the 3D tubular structure are extracted using a neighborhood exploration algorithm: Calculate the number of voxels contained in the foreground points of a voxel's 6-neighborhood and 26-neighborhood. Let m be the number of foreground points in the 6-neighborhood of the object voxel, and n be the number of foreground points in the 26-neighborhood of the object voxel. The candidate point set of endpoints is represented as: ……（1）， in, For the data points in model data O, and Data points The 6-neighborhood and 26-neighborhood are used to obtain the candidate endpoint set through the expression (1); The K-means method is used to cluster candidate points at the endpoints, defined as follows: ……（2）， ……（3）， in, For inter-class distance, , ... Intra-class distance, Represents each data point, The cluster mean is used, and the loss function is: ……（4）； Among them, the inter-class distance and intra-class distance determine the average profile coefficient. The quality of clustering is evaluated by the average profile coefficient. The number of clusters k in equation (2) is calculated based on the average profile coefficient. The K-means algorithm is used for iterative calculation until the termination condition is met. Using the center of each cluster as an endpoint voxel, an endpoint set is formed. .

3. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 1, characterized in that, In step 6, obtaining the initial skeleton line through the distance transformation field includes: Calculate the minimum distance from each voxel to the boundary, and define... For each voxel point, the distance to the boundary is the value of any voxel. The distance transformation is equivalent to a radius of The central sphere, a voxel-based sphere The largest central sphere centered on the model is defined as the sphere that cannot be contained within the model. After calculating the distance transformation field, the center of the largest central sphere is extracted to obtain the initial skeleton point; The burning grass model is used to find the branch endpoints. The two sets covered_points and fire_front are calculated iteratively to simulate the fire transmission process. Here, fire_front represents the set of voxel points that the fire is currently burning, and covered_points represents the set of voxel points that the fire has already burned. All voxel points are included in the covered_points set as a termination condition. The time t when the outermost voxel point of covered_points is transmitted is recorded. The initial skeleton points obtained by the distance field method and the endpoints obtained by the burning grass model method are used as the initial points of the skeleton lines of the final branches.

4. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 3, characterized in that, Step 6 includes fitting the skeleton line using a B-spline curve: Using cubic spline basis as the basis function, the discrete points are fitted. For m+1 data points... Find a cubic B-spline curve: ……（7）； satisfy , The remaining data points The optimization function is calculated using the least squares method: ……（8）； The optimization function is about n-1 control vertices. A minimum value of , where Let be the parameter value for the data point, which is determined by accumulating the chord length parameter. ……（9）； Will Substituting the values ​​into the calculation, we get: ……（10）； Solving equation (10) is a least squares problem. To minimize the objective function f, we need to minimize the number of control vertices (n-1). The derivative is 0, where The l-th derivative is: ……（11）； A system of n-1 equations with n-1 unknowns was obtained. The solution was obtained by Gaussian elimination, resulting in a cubic B-spline curve. This curve was expressed by parameters, and the discrete skeleton points were fitted into a skeleton line.

5. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 4, characterized in that, In step 6, spherical B-splines are applied for fitting optimization. Hausdorff distance is introduced to optimize the original splines. The objective function for optimization is: ……（12）； Let represent the Hausdorff distance from all scattered points to the B-spline curve. The final optimization function is: ……（13）； in, express The projection point on the surface of the sphere B spline; The quasi-Newton algorithm is used to iteratively optimize and adjust the coordinates, number, and radius of the spheres controlling the sphere B-spline curve. Based on the error accuracy, the optimized sphere B-spline curve is obtained.

6. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 1, characterized in that, In step 2, the input model data O is preprocessed, including denoising and downsampling.

7. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 1, characterized in that, In step 1, the model data is a sequence of discrete points described using three-dimensional coordinates.

8. The method for extracting the curve skeleton line of a three-dimensional tubular model based on a hierarchical graph according to claim 5, characterized in that, Step 6 includes optimizing the shortest total time for all voxel points to be covered, and then obtaining the outermost voxel point of this combustion as the endpoint. The optimization formula is as follows: ……（14）； in, Indicates from the source ignition point to The time elapsed since the point was reached.

Images