A fully automatic segmentation method for tracheal tree based on differential geometry

Abstract

This invention discloses a fully automated tracheal tree segmentation method based on differential geometry, comprising the following steps: 1) surface modeling of lung CT scans; 2) calculation of principal curvature and principal direction; 3) filtering of non-tracheal regions; 4) removal of non-tracheal regions; 5) extraction of connected regions; 6) removal of connected regions with planar shape features; 7) fusion of tracheal tree sets S generated by different CT values; 8) obtaining the organ tree surface model. This invention belongs to the field of fully automated tracheal tree segmentation technology, specifically through a fully automated tracheal tree segmentation method based on differential geometry, which is fully automated, easy to implement, insensitive to noisy images, and produces segmented tracheal trees without leakage or blockage issues.

Description

Technical Field

[0001] The present invention belongs to the technical field of fully automatic segmentation of tracheal trees, and specifically refers to a fully automatic segmentation method of tracheal trees based on differential geometry. Background Art

[0002] Accurate segmentation and quantitative analysis of bronchial tree morphology play a very important role in clinical evaluation of the presence and severity of specific diseases (e.g., reduced bronchial lumen area, thinning of wall thickness, etc.). Currently, computed tomography (CT) imaging technology can realize 3D visualization of lung structure, but cannot provide separate display and quantitative description of the bronchial tree. Due to the complexity of the bronchial tree structure, manual segmentation is very time-consuming and labor-intensive. At the same time, the subjective judgment of the operator usually introduces large manual errors. Therefore, an automatic and reliable bronchial tree segmentation method is of great significance for the diagnosis of clinical lung diseases.

[0003] Leakage and obstruction are the two main challenges of automatic segmentation of the tracheal tree. The main reason for leakage and obstruction is the partial volume effect in computed tomography (CT) images, which reduces the contrast between the wall of the tracheal tree and the air in the airway. Leakage will cause the segmented tracheal tree to merge with the surrounding lung parenchyma; while obstruction will cause the tracheal tree to break or the trachea to be discontinuous. In addition, image noise and artifacts will also increase the difficulty of tracheal tree segmentation. Especially for patients with lung diseases, such as chronic obstructive pulmonary disease or interstitial lung disease, the segmentation of the tracheal tree will be more difficult. It cannot be fully automatic, is sensitive to image noise, the segmented tracheal tree is prone to leakage and obstruction, and the algorithm is complex to implement.

[0004] The patent application with application number CN201510009239.X discloses a method for obtaining a three-dimensional tracheal tree from a lung CT image, including extracting a number of image features for constructing an energy term in a cost function; using a multi-kernel learning method, using the combined kernel function of the features to embed into a three-dimensional bronchial seed point extraction algorithm to obtain the seed points that constitute the bronchial segment; and continuously obtaining each independent bronchial segment based on the continuity of the seed points that constitute the bronchial segment. The invention uses a variant of the regional growing method, first segmenting the main bronchus and the bronchial segment, and then "stitching" each bronchial segment with the main bronchus, but there are still a large number of tiny terminal bronchi that have not been segmented, and the stitching operation is more sensitive to noise images, which may cause false bronchi to be incorrectly connected. At the same time, a large number of complex "stitching" algorithms will inevitably lead to an increase in computing costs. Summary of the invention

[0005] To solve the above existing problems, the present invention provides a fully automatic tracheal tree segmentation method based on differential geometry, which is fully automatic, easy to implement, insensitive to noisy images, and the segmented tracheal tree has no leakage and obstruction problems.

[0006] The technical solution adopted by the present invention is as follows: a fully automatic tracheal tree segmentation method based on differential geometry, comprising the following steps:

[0007] 1) Surface modeling of lung CT: The grayscale value range [-900HU, -450HU] is sampled at an equal distance with a sampling interval of 10HU. The CT values ​​of the tracheal tree at different positions vary greatly, and the grayscale value range [-900HU, -450HU] can cover all CT values ​​of the tracheal tree. According to the sampled CT values, the MarchingCube algorithm is used iteratively to model the isosurface of the lung CT to obtain triangular patches.

[0008] 2) Calculation of principal curvature and principal direction: Laplacian smoothing is performed on the triangular facets modeled in step 1) using an area smoothing operator. The specific smoothing formula is as follows:

[0009]

[0010] Where Aj represents the area of ​​the jth adjacent triangle with adjacent vertices, where v j Represents the vertex v i The adjacent j-th vertex;

[0011] 3) Filtering of non-tracheal areas: Establish a classification table of lung tissue principal curvature features: Classify the soft tissue shape type features of the human lung area according to the principal curvature values ​​calculated in step 2). The lung tissue in the human body is mainly divided into the following four types: spherical, such as nodules; flat, such as lung boundaries; convex cylinders, such as vascular trees and concave cylinders, such as tracheal trees. The relationship between the four soft tissue shapes and the principal curvature values ​​is shown in Table 1:

[0012]

[0013] Table 1: Classification table of the main curvature characteristics of lung tissue, where r is the radius of the cylinder and the sphere

[0014] Eliminate non-tracheal areas: According to the above results, after eliminating the sphere, plane, and convex cylinder, what remains are the candidate tracheal area vertices. In order to ensure that no tracheal vertices are missed, this step has three conditions: Cmax ≥ -0.05, unit: 1 / mm, |Cmin| ≥ 0.2, unit: 1 / mm and Amax > 30°. If one of the above three conditions is met, it is identified as a non-tracheal area vertex, where Amax represents the angle between the target vertex and the adjacent vertex in the direction of the minimum principal curvature; the elimination range formed by the three conditions will be smaller than the actual The set of vertices in the non-tracheal area will bring in some lung boundary fragments with planar features; extract connected areas: divide the vertex set determined to be the tracheal area into several connected areas according to the connectivity relationship between the vertices; eliminate connected areas with planar shape features: rotate the mean direction d of the minimum principal curvature of the connected area to the Z axis, synchronously rotate the normal vector n of the vertex in the connected area and project it to the XOY plane perpendicular to the Z axis to obtain the projected normal vector n', calculate the distribution angle θ of n' in the XOY plane, and θ<3π / 2 is the lung boundary fragment with planar features;

[0015] 4) Fusion of tracheal tree set S generated by different CT values: Map the tracheal tree from geometric space to volume data space, calculate the bounding box of S, calculate the minimum and maximum values ​​of the vertices in S on the XYZ axis [Xmin, Xmax, Ymin, Ymax, Zmin, Zmax], for each element in S, sample a single tracheal tree at equal distances along the Z axis [Zmin, Zmax], calculate the intersection contour of each plane perpendicular to the Z axis and the tracheal tree surface model, use the raster scanning algorithm to fill the intersection contour, fill the voxels within the contour with 1, fill the voxels outside the contour with 0, and perform a AND operation on the voxels within the contour filled with 1;

[0016] 5) Generating an organ tree surface model: The tracheal tree data generated in step 4) is subjected to isosurface modeling to obtain a final tracheal tree surface model.

[0017] Furthermore, the isosurface modeling described in 8) adopts MC algorithm.

[0018] The present invention achieves the following beneficial effects by adopting the above scheme: the present invention provides a fully automatic tracheal tree segmentation method based on differential geometry, which is fully automatic, easy to implement, insensitive to noise images, and the segmented tracheal tree has no leakage and obstruction problems.

[0019] Instruction Manual

[0020] Figure 1 A schematic diagram of a fully automatic segmentation process of a tracheal tree based on differential geometry in the present invention;

[0021] Figure 2The minimum curvature direction distribution map of the cylindrical blood vessel connected area and the plane connected area of ​​the tracheal tree fully automatic segmentation method based on differential geometry of the present invention;

[0022] Figure 3 This is a diagram showing the effect of tracheal tree segmentation according to a fully automatic tracheal tree segmentation method based on differential geometry in the present invention. DETAILED DESCRIPTION

[0023] The technical solutions in the embodiments of the present invention are described clearly and completely below. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0024] The technical solution adopted by the present invention is as follows: a fully automatic tracheal tree segmentation method based on differential geometry, comprising the following steps:

[0025] 1) Surface modeling of lung CT: The grayscale value range [-900HU, -450HU] is sampled at an equal distance with a sampling interval of 10HU. The CT values ​​of the tracheal tree at different positions vary greatly, and the grayscale value range [-900HU, -450HU] can cover all CT values ​​of the tracheal tree. According to the sampled CT values, the MarchingCube algorithm is used iteratively to model the isosurface of the lung CT to obtain triangular patches.

[0026] 2) Calculation of principal curvature and principal direction: Laplacian smoothing is performed on the triangular facets modeled in step 1) using an area smoothing operator. The specific smoothing formula is as follows:

[0027]

[0028] Where Aj represents the area of ​​the jth adjacent triangle with adjacent vertices, where v j Represents the vertex v i The adjacent j-th vertex;

[0029] 3) Filtering of non-tracheal areas: Establish a classification table of lung tissue principal curvature features: Classify the soft tissue shape type features of the human lung area according to the principal curvature values ​​calculated in step 2). The lung tissue in the human body is mainly divided into the following four types: spherical, such as nodules; flat, such as lung boundaries; convex cylinders, such as vascular trees and concave cylinders, such as tracheal trees. The relationship between the four soft tissue shapes and the principal curvature values ​​is shown in Table 1:

[0030]

[0031] Table 1: Classification table of the main curvature characteristics of lung tissue, where r is the radius of the cylinder and the sphere

[0032] Eliminate non-tracheal areas: According to the above results, after eliminating these vertices, the remaining ones are candidate tracheal area vertices. In order to ensure that no tracheal vertices are missed, the elimination range formed by the following three conditions will be smaller than the set of actual non-tracheal area vertices, and some lung boundary fragments with planar features will be introduced;

[0033] Extracting connected areas: dividing the vertex set determined as the trachea area into several connected areas according to the connectivity relationship between the vertices;

[0034] Eliminate connected areas with planar shape features: align the mean direction d of the minimum principal curvature of the connected area to the Z axis, rotate the normal vector n of the vertex of the connected area synchronously and project it to the XOY plane perpendicular to the Z axis to obtain the projected normal vector n', calculate the distribution angle θ of n' in the XOY plane, and θ<3π / 2 is the lung boundary fragment with planar features, such as Figure 2 As shown;

[0035] 4) Fusion of tracheal tree set S generated by different CT values: Map the tracheal tree from geometric space to volume data space, calculate the bounding box of S, calculate the minimum and maximum values ​​of the vertices in S on the XYZ axis [Xmin, Xmax, Ymin, Ymax, Zmin, Zmax], for each element in S, sample a single tracheal tree at equal distances along the Z axis [Zmin, Zmax], calculate the intersection contour of each plane perpendicular to the Z axis and the tracheal tree surface model, use the raster scanning algorithm to fill the intersection contour, fill the voxels within the contour with 1, fill the voxels outside the contour with 0, and perform a AND operation on the voxels within the contour filled with 1;

[0036] 5) Obtaining an organ tree surface model: performing isosurface modeling on the tracheal tree data generated in step 4) to obtain a final tracheal tree surface model;

[0037] Furthermore, the CT values ​​described in 1) are iteratively modeled on isosurfaces of lung CT using a MarchingCube algorithm.

[0038] Furthermore, the Laplacian smoothing described in 2) can effectively filter the problem of rough triangular facets after modeling caused by image noise, while not changing the topological structure of the triangular facets, thereby ensuring the accurate calculation of the subsequent principal curvature and principal direction.

[0039] Furthermore, 3) the maximum principal curvature and the minimum principal curvature are Cmax and Cmin respectively, and Dmax and Dmin are the curvature directions of Cmax and Cmin respectively.

[0040] Furthermore, the three conditions described in 3) are Cmax≥-0.05, unit: 1 / mm, |Cmin|≥0.2, unit: 1 / mm and Amax>30°. If one of the above three conditions is met, it is identified as a vertex in the non-tracheal area, where Amax represents the angle between the target vertex and the adjacent vertex in the direction of the minimum principal curvature.

[0041] Furthermore, the isosurface modeling described in 8) adopts MC algorithm.

[0042] The invention discloses a fully automatic tracheal tree segmentation method based on differential geometry, which is fully automatic, easy to implement, insensitive to noise images, and has no leakage and obstruction problems in the segmented tracheal tree.

[0043] The above descriptions are merely embodiments of the present invention and are not intended to limit the patent scope of the present invention. Any equivalent structure or equivalent process transformation made using the contents of the present invention specification, or directly or indirectly applied in other related technical fields, are also included in the patent protection scope of the present invention.

Claims

1. A fully automatic tracheal tree segmentation method based on differential geometry, characterized in that: The following steps are involved: 1) Surface modeling of lung CT: The grayscale value range [-900HU, -450HU] is sampled at an equal distance with a sampling interval of 10HU. The CT values ​​of the tracheal tree at different positions vary greatly, and the grayscale value range [-900HU, -450HU] can cover all CT values ​​of the tracheal tree. According to the sampled CT values, the MarchingCube algorithm is used iteratively to model the isosurface of the lung CT to obtain triangular patches. 2) Calculation of principal curvature and principal direction: Laplacian smoothing is performed on the triangular facets modeled in step 1) using an area smoothing operator. The specific smoothing formula is as follows: Among them A j represents the area of ​​the jth adjacent triangle with adjacent vertices, where v j Represents the vertex v i The adjacent j-th vertex, C max and C min Represent the maximum principal curvature and the minimum principal curvature, D max and D min Corresponding to C max and C min The main direction of 3) Filtering of non-tracheal areas: Establishing a classification table of lung tissue principal curvature features: Classifying the soft tissue shape type features of the human lung area according to the principal curvature values ​​calculated in step 2). The lung tissue in the human body is mainly divided into the following four types: sphere, plane, convex cylinder and concave cylinder; Eliminating non-tracheal areas: According to the above results, after eliminating the sphere, plane and convex cylinder, what remains are the candidate tracheal area vertices. In order to ensure that no tracheal vertices are missed, this step has three conditions: Cmax ≥ -0.05, unit: 1 / mm, |Cmin| ≥ 0.2, unit: 1 / mm and Amax>30°. If one of the above three conditions is met, it is identified as a non-tracheal area vertex, where Amax represents the angle between the target vertex and the adjacent vertex in the direction of the minimum principal curvature; the elimination range formed by the three conditions will be smaller than the set of real non-tracheal area vertices, and will bring in some lung boundary fragments with planar features; extract connected areas: divide the vertex set determined as the tracheal area into several connected areas according to the connectivity relationship between the vertices; eliminate connected areas with planar shape features: rotate the mean direction d of the minimum principal curvature of the connected area to the Z axis, synchronously rotate the normal vector n of the vertex in the connected area and project it to the XOY plane perpendicular to the Z axis to obtain the projection normal vector n', calculate the distribution angle θ of n' in the XOY plane, θ < 3π / 2 is the lung boundary fragment with planar features; 4) Fusion of tracheal tree set S generated by different CT values: Map the tracheal tree from geometric space to volume data space, calculate the bounding box of S, calculate the minimum and maximum values ​​of the vertices in S on the XYZ axis [Xmin, Xmax, Ymin, Ymax, Zmin, Zmax], for each element in S, sample a single tracheal tree at equal distances along the Z axis [Zmin, Zmax], calculate the intersection contour of each plane perpendicular to the Z axis and the tracheal tree surface model, use the raster scanning algorithm to fill the intersection contour, fill the voxels inside the contour with 1, and fill the voxels outside the contour with 0, and finally generate a final volume data based on the volume data of set S. The voxel value assignment rule in the final volume data is as follows: if there is an individual data at this position in set S with a voxel value of 1 at this position, then the voxel value at this position is assigned to 1, otherwise it is assigned to 0; 5) Generating an organ tree surface model: The tracheal tree data generated in step 4) is subjected to isosurface modeling to obtain a final tracheal tree surface model.

2. The method for fully automatic segmentation of tracheal tree based on differential geometry according to claim 1, characterized in that: The isosurface modeling adopts the MC algorithm.

Images