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

Abstract

In order to obtain high-precision industrial CT isosurface mesh for reverse engineering, the present application proposes a surface extraction method based on dual octree of adaptive volume data of industrial CT, directly constructs adaptive octree volume data grid from industrial X-ray computed tomography projection graph sequence, and based on this, high-precision and crack-free interpolation isosurface mesh is obtained. In terms of extracting isosurface mesh, the Marching Cubes method is easy to generate ambiguous surfaces with aliasing and oscillation around the accurate surface, and such methods based on CT discrete volume data cannot avoid the influence brought by CT reconstruction artifacts and discretization errors. The method firstly recursively constructs adaptive volume data grid, the voxel value at the node is obtained by reconstructing the projection graph by FDK algorithm, then the vertex is calculated and the dual grid is constructed by QEFs technology, finally the surface profile is accurately extracted by using Marching Cubes lookup table.

Description

TECHNICAL FIELD

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

[0002] CT volume data is a discrete data field, and a uniformly spaced volume data distance field has been widely used, but has a major drawback: it must rely on high-resolution scanning data to reproduce thin features and sharp edges of the model surface, but the resolution of the initial volume data field obtained by the current CT scanning device is limited. If the grid density is subdivided to capture these features, since the sampling rate of the structured grid is uniform, the resulting volume data has little or no change in the surface in the region that does not require high resolution, and the grid is too fine, wasting a large amount of memory space. Moreover, 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 an octree to build adaptive volume data. However, the main drawback of the traditional octree method is that when the surface is extracted by the isosurface, since the isosurface is located between different size cells, the vertices of the triangles on both sides are calculated at different resolutions, which may cause cracks. In order to avoid generating cracks, the strategy in this paper is to build a dual grid. Under this method, the vertices of each triangle in the dual surface are located inside the cell, rather than on the edge of the cell as in the classic method. Therefore, each surface cell has a unique associated vertex, so there are no cracks on the surface, and the calculation of the CT value of the octree node uses the analytical differential method to calculate the CT value directly from the X-ray projection value, which is less affected by the CT reconstruction artifacts compared to the discrete differential calculation of the CT volume. Experimental analysis shows that this algorithm can extract the surface more accurately than the traditional method. SUMMARY

[0004] The main purpose of the present application is to provide a surface extraction method of a dual octree based on adaptive volume data of industrial CT, which uses an octree instead of a traditional regular grid for volume data construction, and combines a dual graph to extract the contour, so as to realize adaptive surface extraction of CT data and reproduce sharp edges of CT scanning data, and has obvious improvement in accuracy compared with the traditional MC algorithm.

[0005] The technical scheme of the present application is as follows:

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

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

[0008]

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

[0010] Second step: Constructing the adaptive volume data field by using the octree structure. The adaptive volume data field is constructed by using the bottom-up method so as to reduce the extraction time and improve the surface precision. The CT value f(p) and the gradient value The FDK algorithm is used for calculation, as shown in equations (2) and (3):

[0011]

[0012] Third step: Constructing the dual grid vertex. Since each node of the adaptive octree generated in the second step has attributes such as the gradient and the voxel value, the QEF energy function is used to define the feature so as to extract the accurate dual grid vertex of the surface:

[0013]

[0014]

[0015] where ω i is a weight factor, which normalizes the contribution of each tangent plane to the vertex calculation, so that it has equal weight, which is represented by equation (5). is the gradient value at the node (x i ,y i ,z i ). T i (p) is the tangent plane equation of the implicit function f(x,y,z) at the point p i , which is obtained from the first-order Taylor expansion of the function f(x,y,z) at p i , and is defined as equation (6):

[0016] Solving the linear equation set to obtain the dual vertex p, which generates an internal feature vertex in the dual grid of each octree node.

[0017] Fourth step: Establishment and connection of the dual graph.

[0018] The specific steps are as follows: if a leaf node is traversed, subsequent recursive calls are stopped, otherwise the octree will recursively go along the parent node to eight child nodes. The topological structure generated by the topology creation is actually the topological structure of the dual grid, which is topologically dual to the octree.

[0019] Step 5: surface extraction of sharp edges. Since linear interpolation can cause inaccurate isosurfaces and visible reconstruction artifacts, the present application uses Hermite second-order interpolation technology to generate more accurate isosurfaces.

[0020] The above problem is converted into a linear system in the equation, which is solved by the least square method using a conventional equation, as represented by equation (6):

[0021]

[0022] The value of the point is calculated by quadratic function interpolation, and finally the marching cubes algorithm is used to extract the contour model.

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

[0024] The present application proposes a method of constructing an adaptive octree volume data grid directly from an industrial X-ray computed tomography (ICT) projection graph sequence using an analytical reconstruction method, and obtaining a high-precision and crack-free interpolated isosurface grid based on the data. The proposed method first constructs an adaptive volume data grid recursively from top to bottom, calculates the voxel value at the node by the FDK algorithm, then calculates the vertex by the QEFs technology and constructs a dual grid, and finally accurately extracts the surface contour by the Hermite interpolation technology. The size accuracy and surface accuracy of the grid generated by the method are superior to those of the conventional method. BRIEF DESCRIPTION OF DRAWINGS

[0025] Figure 1 is the result of the initial image filtered by S-L;

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

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

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

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

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

[0031] To make the description of the object, technical solutions and advantages of the present application clearer, the following will take the stepped cylinder data set obtained by cone beam CT scanning as an example, and combine the accompanying drawings in the embodiments of the present application to make a complete description of the technical solutions in the present application.

[0032] 1) Input CT image data. Input image, filter using Shepp-Logan convolution filter, the result is shown in Figure 1 , wherein the deconvolution filter uses function h(x) for deconvolution along the x axis, as shown in the following formula:

[0033]

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

[0035] 2) Construct adaptive octree of volume data. Adopt bottom-up method to construct adaptive volume data field using octree structure, that is, more dense sampling is used on the surface in the region of interest to make the surface present rapidly changing curvature, and more sparse sampling is used in other positions to reduce extraction time while improving surface precision.

[0036] Calculate voxel value of octree sampling point. CT value can be obtained from the projection map generated by cone beam CT using filter back projection based generation as shown in Figure 2 . The projection operator P of point p projected to xy plane Figure 3 ) is defined as:

[0037] P(p)=a(p)(p-s) (1)

[0038] a(p)=|o-s| / (|o-s|+(p z ) (2)

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

[0040]

[0041] 3) Construct double grid vertex. Adopt QEF energy function to define features as shown in Figure 5 , calculate tangent plane of f(x,y,z) map on the grid of points sampled on c, which is in the form of:

[0042] E(p)=∑ i [ni • (p - p i ] 2 (4)

[0043] Geometrically, it defines the sum of the squared distance of point p to the set of planes {p i , n i}.

[0044] Rewrite (4) as:

[0045]

[0046] where v i is a weight factor that normalizes the contribution of each tangent plane to the vertex so that each has equal weight, which is given by (6). is the gradient value at node (x i , y i , z i ). T i (p) is the tangent plane equation of the implicit function f(x, y, z) at point p i , which is obtained from the first order Taylor expansion of f(x, y, z) at p i , and is defined as (7).

[0047] Solving the linear system of equations yields the dual vertex p, which produces an interior feature vertex in each dual grid of the octree node. In this paper, the solution is obtained by the closest distance to the grid center. The pseudo-inverse matrix involved can be obtained by the correction of singular values of its SVD decomposition. This algorithm produces an interior feature vertex in each dual grid of 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 feature point in the grid where the sign of the implicit function changes. The topology creation actually generates the topology of the dual grid, which is topologically dual to the octree. Figure 4 The topology of the dual grid created from a two-dimensional quadtree is shown, where each vertex of the dual grid is located at the center of its cube.

[0049] 5) Surface extraction preserving sharp edges. The properties of the new grid constructed in step 5) make it possible to use marching cubes to extract the contour. Since linear interpolation can lead to inaccurate isosurfaces and visible reconstruction artifacts, this paper uses Hermite second-order interpolation techniques to generate more accurate isosurfaces.

[0050] Convert the above problem to solve the linear system in the equation using the least squares method

[0051]

[0052] The values of x are finally computed by quadratic function interpolation for the points and the surface contours are connected according to the Marching Cubes lookup table, resulting in Figure 6 as shown.

Claims

1. A surface extraction method based on dual octree of industrial CT adaptive volume data, characterized by: Here are the steps: Step 1: Input CT image data; use the Shepp-Logan convolution filter and perform deconvolution filtering along the x-axis using the function h(x), as shown below: Where, β(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 projected image; Step 2: Use the octree structure to construct the adaptive volume data field; use the bottom-up method to construct the 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, and generate an internal feature vertex in the dual mesh of each octree node; Step 4: Establish and connect the dual graph; The dual graph is built 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: using the Hermite second-order interpolation technique to generate an isosurface; After converting the above problem, we can use the conventional equation to solve the linear system in the least square method, which is expressed by equation (6): The calculated value of x is interpolated by a 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 call is stopped; otherwise, the octree will recurse along the parent node to the eight child nodes; the topology creation actually generates the topological structure of the dual grid, which is topologically dual to the octree.

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

Images