Surface extraction method of dual octree based on industrial CT adaptive volume data

Abstract

In order to obtain a high-precision industrial CT contour surface grid for reverse engineering, the invention provides a dual octree surface extraction method based on industrial CT adaptive volume data, which comprises the following steps of: directly constructing an adaptive octree volume data grid from an industrial X-ray computed tomography projection drawing sequence; and based on this, a high-precision and crack-free interpolation contour surface grid is obtained. In the aspect of contour surface grid extraction, the Marching Cubes method is easy to generate aliasing and oscillating ambiguous surfaces around an accurate surface, and the Marching Cubes method is based on CT discrete body data and is inevitably influenced by CT reconstruction artifacts and discretization errors. According to the method, firstly, an adaptive volume data grid is constructed recursively, voxel values at nodes are obtained by reconstructing a projection drawing through an FDK algorithm, then vertexes are calculated through the QEFs technology, dual grids are constructed, and finally, a Marching Cubes lookup table is utilized to accurately extract the surface contour.

Description

Technical Field

[0001] The invention belongs to the technical field of industrial CT digital image processing and relates to a surface extraction method based on a dual octree of industrial CT adaptive volume data. Background Art

[0002] CT volume data is a discrete data field, and uniformly spaced volume data distance fields have been widely used, but there is a major disadvantage: it must rely on high-resolution scan data to reproduce thin features and sharp edges on the model surface, but the initial volume data field resolution obtained by current CT scanning equipment is limited. If these features are captured by subdividing the grid density, since the sampling rate of the structured grid is consistent, the volume data obtained in areas where high resolution is not required, such as areas with little or no surface changes, will have too fine grid particles, wasting a lot of memory space, and since the structured grid must traverse the entire grid to find the boundary contour, the isosurface extraction process is very slow.

[0003] This situation can be reduced by using octrees to construct adaptive volume data. However, the main disadvantage of the traditional octree method is that when the surface is extracted by isosurfaces, since the isosurfaces are located between cells of different sizes, the triangle vertices on both sides are calculated with different resolutions during calculation, and cracks may appear. In order to avoid generating cracks, the strategy of this paper is to construct a dual mesh. Under this method, the vertex of each triangle in the dual surface is located inside the cell, rather than on the edge of the cell as in the classical method. Therefore, each surface cell has a unique associated vertex, so there are no cracks on the surface. The CT value of the octree node is calculated using the method of calculating the analytical differential of the CT value directly from the X-ray projection value. Compared with calculating the discrete differential of the CT volume, the CT value is less affected by CT reconstruction artifacts. Experimental analysis shows that the algorithm can extract the surface more accurately than the traditional method. Summary of the invention

[0004] The main purpose of the present invention is to provide a surface extraction method based on a dual octree of industrial CT adaptive volume data. The octree is used to replace the traditional regular grid in volume data construction, and the dual graph is combined to extract the contour, which can realize adaptive surface extraction of CT data and reproduce the sharp edges of CT scan data. Compared with the traditional MC algorithm, the accuracy is significantly improved.

[0005] The technical solution of the present invention:

[0006] The surface extraction method based on dual octree of industrial CT adaptive volume data has the following steps:

[0007] Step 1: Input CT image data. Use Shepp-Logan convolution filter and perform deconvolution filtering along the x-axis using function h(x), as shown in the following formula:

[0008]

[0009] Among them, β(q) = cos(γ), γ is the angle between the point p where the X-ray penetrates the object and the y-axis of the turntable, S θ (t) represents the projection intensity at a given angle, represents the coordinates of the projection point in the x direction, shares the y coordinate with p, h(t) is the Shepp-Logan reconstruction filter function, and sets h(t) = 2 / (π 2 δ 2 -4t 2 ), δ is the pixel size of the projected image.

[0010] Step 2: Use the octree structure to construct an adaptive volume data field. A bottom-up approach is used to construct an adaptive volume data field in order to reduce extraction time and improve surface accuracy. The CT value f(p) and gradient value of the octree node The FDK algorithm is used for calculation, as shown in formula (2) and formula (3):

[0011]

[0012] Step 3: Construct dual mesh vertices. Since each node of the adaptive octree generated in the second step has attributes such as gradient and voxel value, the QEF energy function is used to define the features to extract the accurate dual mesh vertices of the surface:

[0013]

[0014]

[0015] Among them, ω i is the weight factor, which normalizes the contribution of each tangent plane to vertex calculation so that they have equal weights, expressed as formula (5). is a node (x i ,y i ,z i ) is the gradient value at T i (p) is the implicit function f(x,y,z) at point p i The equation of the tangent plane at p is given by the function f(x,y,z) at p i It is obtained from the first-order Taylor expansion at , which is defined as formula (6):

[0016] Solving linear equations Get the dual vertex p, which generates an internal feature vertex in the dual mesh of each octree node.

[0017] Step 4: Establish and connect the dual graph.

[0018] The dual graph is built recursively. The specific steps are as follows: If the traversal reaches a leaf node, the subsequent recursive call is stopped, otherwise the octree will recurse along the parent node to the eight child nodes. Topology creation actually generates the topological structure of the dual grid, which is topologically dual to the octree.

[0019] Step 5: Surface extraction with sharp edges preserved Since linear interpolation may lead to inaccurate isosurfaces and visible reconstruction artifacts, the present invention uses the Hermite second-order interpolation technique to generate more accurate isosurfaces.

[0020] After converting the above problem, the linear system in the formula is solved by the conventional equation using the least square method, which is expressed by formula (6):

[0021]

[0022] The calculated value of x is obtained by interpolating the quadratic function to obtain the value of the point, and finally the marching cube algorithm (MarchingCubes) is used to extract the contour model.

[0023] The beneficial effects of the present invention are as follows:

[0024] The present invention proposes a method for constructing an adaptive octree volume data grid directly from an industrial X-ray computed tomography (ICT) projection image sequence using an analytical reconstruction method, and obtains a high-precision and crack-free interpolation isosurface grid based on this data. The proposed method first recursively constructs an adaptive volume data grid from top to bottom, reconstructs the projection image by the FDK algorithm to calculate the voxel values ​​at the nodes, then calculates the vertices and constructs the dual grid by the QEFs technology, and finally uses the Hermite interpolation technology to accurately extract the surface contour. The size accuracy and surface accuracy of the generated grid by this method are better than those of the traditional method. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 is the result of SL filtering of the initial image;

[0026] Figure 2 It is a schematic diagram of cone beam CT scanning projection;

[0027] Figure 3 It is a schematic diagram of the geometric relationship between the imaginary plane penetrated by the ray source and the projection plane;

[0028] Figure 4 This is a schematic diagram of dual graph construction;

[0029] Figure 5 It is a schematic diagram of QEF minimization;

[0030] Figure 6 is the final extracted geometric model. DETAILED DESCRIPTION

[0031] In order to make the purpose, technical solution and advantages of the present invention more clearly described, the technical solution of the present invention is fully described below by taking the stepped cylindrical data set obtained by cone-beam CT scanning as an example, combined with the drawings in the embodiments of the present invention.

[0032] 1) Input CT image data. Input image and filter it using Shepp-Logan convolution filter. The result is as follows: Figure 1 As shown, the deconvolution filter uses the function h(x) to perform deconvolution along the x-axis, as shown in the following formula:

[0033]

[0034] The convolution function is set to h(t) = 2 / (π 2 δ 2 -4t 2 ), δ is the pixel size of the projected image on the scaled detector.

[0035] 2) Construct an adaptive octree for volume data. A bottom-up approach is used to construct an adaptive volume data field using an octree structure, that is, the surface is sampled more densely in the area of ​​interest so that the surface exhibits rapidly changing curvature, while sparser sampling is used in other locations to reduce extraction time while improving surface accuracy.

[0036] Calculate the voxel value of the octree sampling point. The CT value can be obtained from the projection map generated by cone beam CT using filtered back projection as follows Figure 2 As shown. Project point p onto the xy plane ( Figure 3 ) is defined as:

[0037] P(p)=α(p)(ps)(1)

[0038] α(p)=|os| / (|os|+(p) z )(2)

[0039] Based on (1) and (2), the FDK algorithm is then used to calculate the node CT value f(p) and the gradient value

[0040]

[0041] 3) Construct dual mesh vertices. Use QEF energy function to define features such as Figure 5 As shown, the tangent plane of the graph f(x, y, z) is calculated on the grid of points sampled on c, and its form is:

[0042] E(p)=∑ i [ni ·(pp i )] 2 (4)

[0043] Its geometric meaning is to define the point p to the plane set {p i , n i}The sum of the bi-norms of the distances.

[0044] Rewrite (4) as follows:

[0045]

[0046] Among them, v i is the weight factor, which normalizes the contribution of each tangent plane to vertex calculation so that they have equal weights, expressed as formula (6). is a node (x i ,y i ,z i ) is the gradient value at T i (p) is the implicit function f(x,y,z) at point p i The equation of the tangent plane at p is given by the function f(x,y,z) at p i It is obtained from the first-order Taylor expansion at , which is defined as formula (7):

[0047] Solving linear equations The dual vertex p is obtained, which generates an internal feature vertex in the dual grid of each octree node. This paper solves it by the solution closest to the center of the voxel grid. The pseudo-inverse matrix involved can be obtained by correcting the singular values ​​of its SVD decomposition. This algorithm will generate an internal feature vertex in each dual grid in the octree node.

[0048] 4) Establishment and connection of the dual graph. Based on the result of step 3), we know that there must be a characteristic point in the grid where the implicit function sign changes. Topology creation actually generates the topological structure of the dual grid, which is topologically dual to the octree. Figure 4 Shows the topology of a dual mesh created from a 2D quadtree, where each vertex of the dual mesh is located at the center of its cube.

[0049] 5) Surface extraction with sharp edges preserved. The properties of the new mesh constructed in step 5) make it possible to extract contours using marching cubes. Since linear interpolation leads to inaccurate isosurfaces with visible reconstruction artifacts, this paper uses the Hermite second-order interpolation technique to generate more accurate isosurfaces.

[0050] After converting the above problem, use the conventional equation to solve the linear system in the formula using the least squares method.

[0051]

[0052] Finally, the value of x is calculated by quadratic function interpolation and the surface contour is connected according to the Marching Cubes lookup table. The result is as follows Figure 6 shown.

Claims

1. A surface extraction method based on dual octree of industrial CT adaptive volume data, characterized in that: Here are the steps: Step 1: Input CT image data; use Shepp-Logan convolution filter and perform deconvolution filtering along the x-axis using function h(x), as shown in the following formula: Among them, β(q) = cos(γ), γ is the angle between the point p where the X-ray penetrates the object and the y-axis of the turntable, S θ (t) represents the projection intensity at a given angle, represents the coordinates of the projection point in the x direction, shares the y coordinate with p, h(t) is the Shepp-Logan reconstruction filter function, and δ is the pixel size of the projection image; Step 2: Use the octree structure to construct an adaptive volume data field; Use a bottom-up approach to construct an adaptive volume data field, where the CT value f(p) and gradient value of the octree node The FDK algorithm is used for calculation, as shown in formula (2) and formula (3): Step 3: Construct dual mesh vertices; Since each node of the adaptive octree generated in the second step has gradient and voxel value attributes, the QEF energy function is used to define the features to extract the accurate dual mesh vertices of the surface: Among them, ω i is the weight factor that normalizes the contribution of each tangent plane to vertex calculation, expressed as formula (5); is a node (x i ,y i ,z i ) at the gradient value; T i (p) is the implicit function f(x,y,z) at point p i The equation of the tangent plane at p is given by the function f(x,y,z) at p i It is obtained from the first-order Taylor expansion at , which is defined as formula (6): Solving linear equations Get the dual vertex p, which generates an internal feature vertex in the dual mesh of each octree node; Step 4: Establishment and connection of dual graph; The dual graph is established recursively; Step 5: Surface extraction with sharp edges preserved.

2. The surface extraction method based on dual octree of industrial CT adaptive volume data according to claim 1, characterized in that: In the first step, h(t)=2 / (π 2 δ 2 -4t 2 ).

3. The surface extraction method based on dual octree of industrial CT adaptive volume data according to claim 1 or 2, characterized in that: The fifth step is specifically as follows: the present invention uses the Hermite second-order interpolation technology to generate a more accurate isosurface; After converting the above problem, the linear system in the formula is solved by the conventional equation using the least square method, which is expressed by formula (6): The calculated value of x is obtained by interpolating the quadratic function to obtain the value of the point, and finally the marching cubes algorithm is used to extract the contour model.

4. The surface extraction method based on dual octree of industrial CT adaptive volume data according to claim 1 or 2, characterized in that: The fourth step is as follows: if a leaf node is traversed, the subsequent recursive calls are stopped, otherwise the octree will recurse along the parent node to the eight child nodes; the topology creation actually generates a topological structure of a dual grid, which is topologically dual to the octree.

5. The surface extraction method based on dual octree of industrial CT adaptive volume data as claimed in claim 3, characterized in that: The fourth step is as follows: if a leaf node is traversed, the subsequent recursive calls are stopped, otherwise the octree will recurse along the parent node to the eight child nodes; the topology creation actually generates a topological structure of a dual grid, which is topologically dual to the octree.

Images