Method for extracting pulmonary artery and vein based on CT image

Abstract

The application discloses a lung artery and vein extraction method based on a CT image and belongs to the technical field of medical image processing. The method comprises the following steps: performing anisotropic diffusion filtering on a three-dimensional CT image to obtain a smooth image; segmenting a lung parenchyma region from the smooth image to obtain a lung parenchyma mask; obtaining a lung parenchyma image based on the lung parenchyma mask and the smooth image mask; calculating multi-scale Hessian matrix eigenvalues on the lung parenchyma image; calculating a tubular structure enhancement response value according to the eigenvalues to obtain a multi-scale tubular structure enhancement image; performing connected domain analysis on the enhancement image, and screening out lung artery and vein candidate regions according to size and shape characteristics; and distinguishing the lung artery and vein in the candidate regions by using local gray distribution and morphological characteristics to obtain an extraction result.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, specifically a method for extracting pulmonary arteries and veins based on CT images. Background Technology

[0002] Accurate extraction of pulmonary arteries and veins from chest CT images is a crucial step in the diagnosis, surgical planning, and pathological analysis of pulmonary vascular diseases. Current techniques primarily rely on tubular structure enhancement filtering to segment pulmonary vessels as a whole, followed by attempts to separate arteries and veins. This approach first enhances the tubular structures in the image using tubular filters such as multi-scale Hessian matrices, then extracts the complete pulmonary vascular tree through thresholding or region growing. However, this approach has significant drawbacks. In the filtering preprocessing stage, conventional Gaussian smoothing or multi-scale filters, while suppressing noise, blur vessel boundaries, especially for small vessels and branch tips, where edge information is easily destroyed by smoothing operations. This leads to offsets in subsequent Hessian matrix calculations, resulting in vessel response fragmentation or diameter estimation errors. In the final arteriovenous classification stage, traditional methods often rely solely on the overall geometric morphology of the vessels, such as segmenting them by analyzing the hierarchical structure or course of the vascular tree. While the pulmonary arteries and pulmonary veins exhibit clear anatomical differences near the hilum of the lung, their courses and branching morphologies are highly similar in the terminal regions of the lung parenchyma. Relying solely on coarse global geometric analysis makes accurate differentiation difficult, leading to a high misclassification rate of peripheral vessels and affecting the integrity and accuracy of arteriovenous separation. Therefore, achieving precise differentiation between pulmonary arteries and pulmonary veins through preprocessing that effectively smooths noise and maintains clear vascular boundaries has become an urgent problem to be solved. Summary of the Invention

[0003] The purpose of this invention is to provide a method for extracting pulmonary arteries and veins based on CT images. In the preprocessing stage, the method uses anisotropic diffusion filtering to suppress noise while maintaining the sharpness of the blood vessel boundaries. In the subsequent classification stage, it comprehensively utilizes the local gray-level distribution and skeleton morphology features of the candidate region to improve the completeness and accuracy of pulmonary artery and vein extraction.

[0004] To achieve the above objectives, the present invention provides the following technical solution: The present invention provides a method for extracting pulmonary arteries and veins based on CT images, the method comprising the following steps:

[0005] Anisotropic diffusion filtering is applied to 3D CT images to obtain smooth images. By applying anisotropic diffusion filtering to the images, noise in the images can be effectively suppressed while maintaining the clarity of blood vessel edges, providing a high-quality image foundation for subsequent accurate segmentation and blood vessel extraction.

[0006] The lung parenchyma region is segmented from the smooth image to obtain a lung parenchyma mask. This step can accurately separate the lung parenchyma from surrounding bone, muscle and other tissues, limiting the processing scope to the interior of the lung parenchyma, thereby eliminating interference from extrapulmonary structures on blood vessel extraction.

[0007] A lung parenchyma image is obtained by performing a masking operation on the smoothed image. By preserving only the image information inside the lung parenchyma through the masking operation, the computational area for blood vessel extraction is further reduced, thus improving processing efficiency.

[0008] The Hessian matrix was calculated on lung parenchyma images to obtain Hessian matrix eigenvalues ​​at multiple scales. By calculating the Hessian matrix at multiple scales, tubular structures of different diameters can be detected with greater sensitivity, enhancing the response capability to pulmonary vascular branches at various levels.

[0009] The tubular structure enhancement response value for each voxel is calculated based on the eigenvalues ​​of the Hessian matrix at each scale, resulting in a multi-scale tubular structure enhancement image. Constructing a tubular response function based on eigenvalues ​​significantly improves the contrast of tubular vascular structures in the image, suppresses sheet-like or dot-like non-vascular structures, and makes the pulmonary arteries and veins more prominent in the image.

[0010] Connectivity analysis is performed on enhanced images of multi-scale tubular structures to screen candidate regions for pulmonary arteries and pulmonary veins based on the size and shape characteristics of the connected regions. This step leverages the spatially continuous, elongated tubular shape of vascular trees to eliminate isolated noise points and non-target structures whose morphology does not conform to vascular characteristics, thereby rapidly locating candidate vascular regions.

[0011] Within the candidate regions of pulmonary arteries and pulmonary veins, local grayscale distribution and morphological features are used to distinguish between pulmonary arteries and pulmonary veins, resulting in pulmonary artery and pulmonary vein extraction. By integrating grayscale statistical information and the topological morphological features of the skeleton lines, it is possible to effectively distinguish between arteries and veins that exhibit differences due to variations in blood oxygen content and course, thus achieving arteriovenous separation and extraction of pulmonary vessels.

[0012] As a preferred technical solution of the present invention, the specific process of obtaining a smooth image by anisotropic diffusion filtering of a three-dimensional CT image includes: using the three-dimensional CT image as the initial image, setting the diffusion time step and the number of iterations; in each iteration, calculating the diffusion coefficient of each voxel point based on the gray-level gradient of the current image; updating the gray-level value of each voxel point in the current image based on the diffusion coefficient to obtain the updated image; repeating the iteration process until the number of iterations is reached, and using the finally obtained updated image as the smooth image. Specifically, when calculating the diffusion coefficient, for each voxel point in the current image, the gray-level gradient vector of that voxel point in three-dimensional space is calculated, and the diffusion coefficient of that voxel point is calculated based on the magnitude of the gray-level gradient vector. When the magnitude of the gray-level gradient vector is less than a preset gradient threshold, the diffusion coefficient value tends towards a first preset value; when the magnitude of the gray-level gradient vector is greater than the preset gradient threshold, the diffusion coefficient value tends towards a second preset value. This allows the filtering process to perform strong smoothing in flat areas with smaller gradients, while weakening smoothing in edge areas with larger gradients, achieving the beneficial effect of edge preservation and noise reduction.

[0013] As another preferred technical solution of the present invention, the specific process of segmenting the lung parenchyma region from a smooth image to obtain a lung parenchyma mask includes: performing grayscale thresholding on the smooth image to obtain an initial lung parenchyma region; performing morphological opening operation on the initial lung parenchyma region to separate the lung parenchyma from the surrounding tissue; performing connected component analysis on the region after morphological opening operation, and selecting the two connected components with the largest volume as the left and right lung parenchyma regions; performing morphological closing operation on the left and right lung parenchyma regions to fill the cavities inside the lung parenchyma, thereby obtaining a lung parenchyma mask. This process can automatically and completely extract the left and right lung fields and repair internal cavities caused by structures such as blood vessels and bronchi, resulting in a continuous and complete lung parenchyma mask.

[0014] As a further preferred technical solution of the present invention, the process of calculating the Hessian matrix on a lung parenchyma image to obtain the eigenvalues ​​of the Hessian matrix at multiple scales includes: setting a set of Gaussian kernel scale parameters and traversing each scale parameter; for the current scale parameter, using the Gaussian kernel corresponding to that scale to perform convolution filtering on the lung parenchyma image to obtain the filtered image; calculating the second-order partial derivative at each voxel point in the filtered image to construct the Hessian matrix of that voxel point; and solving for the eigenvalues ​​of the Hessian matrix to obtain the three eigenvalues ​​of that voxel point at the current scale. By introducing multi-scale Gaussian filtering, it is possible to match blood vessels of different thicknesses, so that both the main blood vessels near the hilum and the peripheral blood vessels of the lung can be effectively enhanced.

[0015] Accordingly, the process of calculating the tubular structure enhancement response value of each voxel based on the Hessian matrix eigenvalues ​​at each scale includes: for each voxel, obtaining the three Hessian matrix eigenvalues ​​of that voxel at each scale; sorting the three eigenvalues ​​at each scale according to their absolute values; calculating the tubular response value of that voxel at that scale based on the ratio of the first eigenvalue to the second eigenvalue and the ratio of the first eigenvalue to the third eigenvalue among the sorted eigenvalues; traversing all scales and taking the maximum value of the tubular response value across all scales as the final tubular structure enhancement response value of that voxel; and traversing all voxels in the lung parenchyma image to obtain a multi-scale tubular structure enhancement image. This method utilizes the characteristic that the eigenvalues ​​of tubular structures have a specific magnitude relationship under ideal conditions to construct an enhancement function that enhances blood vessels while suppressing other morphological structures.

[0016] As another preferred technical solution of the present invention, the process of performing connected component analysis on multi-scale tubular structure enhanced images and screening candidate regions for pulmonary arteries and pulmonary veins includes: binarizing the multi-scale tubular structure enhanced images to obtain binarized tubular structure images; labeling the binarized tubular structure images with connected components to obtain multiple connected components; calculating the number of voxels, major axis length, and minor axis length of each connected component; selecting connected components with a voxel count greater than a preset voxel count threshold as initial connected components; and selecting connected components with a ratio of major axis length to minor axis length greater than a preset aspect ratio threshold as candidate regions for pulmonary arteries and pulmonary veins based on the combined screening of volume and aspect ratio. Through the combined screening of volume and aspect ratio, tiny noise particles and clump-like non-vascular structures can be reliably excluded, while retaining slender tubular vascular candidate targets.

[0017] In another preferred embodiment of the present invention, the process of distinguishing pulmonary arteries and pulmonary veins within candidate regions of pulmonary arteries and pulmonary veins using local grayscale distribution and morphological features includes: extracting the CT values ​​of all voxel points within each candidate region of the original 3D CT image; calculating the mean and standard deviation of the CT values ​​of voxel points within each candidate region to obtain the grayscale distribution features of the candidate region; for each candidate region, calculating the skeleton line of the candidate region, obtaining the number of bifurcation points of the skeleton line and the average distance between the bifurcation points to obtain the morphological features of the candidate region; inputting the grayscale distribution features and morphological features of each candidate region into a preset classifier to obtain the category label of the candidate region, which is pulmonary artery or pulmonary vein; and classifying the candidate regions into the pulmonary artery set and the pulmonary vein set respectively according to the category label to obtain the pulmonary artery extraction result and the pulmonary vein extraction result. By integrating grayscale statistical features with the morphological topological features of the vascular tree and using a classifier for automatic discrimination, the accuracy and robustness of arteriovenous differentiation are significantly improved.

[0018] In calculating the morphological features of candidate regions, a morphological thinning operation is performed on the binarized image of the candidate regions to obtain skeleton line images. Skeleton line intersections are identified in the skeleton line images, and the number of intersections is counted as the number of bifurcation points. The Euclidean distance from each skeleton line intersection to its adjacent intersections is calculated, and all Euclidean distances are averaged to obtain the average distance between bifurcation points. Skeleton line features can reflect the differences in bifurcation patterns between arteries and veins, providing an effective basis for the classifier's discrimination.

[0019] Preferably, the classifier construction and usage process is as follows: A training sample set labeled with pulmonary artery and pulmonary vein tags is pre-collected, with each training sample containing grayscale distribution features and morphological features; a support vector machine classifier is trained using the training sample set to obtain the trained classifier; the grayscale distribution features and morphological features of each candidate region are combined into a feature vector and input into the trained classifier; the probability value output by the classifier is obtained, and if the probability value is greater than a preset classification threshold, it is determined to be a pulmonary artery; otherwise, it is determined to be a pulmonary vein. Using a pre-trained supervised classification model, accurate separation of arteries and veins can be achieved based on automatic extraction, avoiding the limitations of manually setting rules.

[0020] The technical effects and advantages provided by the present invention in the above technical solution are as follows:

[0021] By employing anisotropic diffusion filtering to process 3D CT images, this scheme adaptively controls the diffusion level using local gray-level gradient information. Strong diffusion is applied in areas with gentle gray-level changes to effectively suppress noise, while diffusion is restricted in vascular boundary areas with large gray-level gradients. This effectively smooths internal image noise, improves image quality, and preserves the edge details of the vessel walls completely and sharply. This processing provides a more accurate gradient information basis for subsequent calculation of the Hessian matrix, making the response of the Hessian matrix eigenvalues ​​at vessel boundaries more precise. It avoids the weakening, breakage, and artificial increase in vessel diameter caused by blurred boundaries in conventional Gaussian smoothing, ensuring continuous response and true size representation of small vascular structures at the terminal pulmonary parenchyma. By simultaneously extracting gray-level distribution features based on CT values ​​and bifurcation morphology features based on skeleton line analysis within candidate regions of the pulmonary artery and pulmonary vein for classification, this scheme integrates tissue density information with branching structure characteristics. The varying oxygen content of blood flow within the pulmonary artery leads to a systematic difference in its CT value distribution compared to the pulmonary vein. Extracting the mean and standard deviation of voxel CT values ​​within the region captures this inherent density distribution attribute. Furthermore, the dendritic branching of the pulmonary artery is typically denser and more complex than that of the accompanying pulmonary vein. Calculating the number of bifurcation points on the skeleton lines of the candidate region and the average distance between these bifurcation points allows for a quantitative description of the topological complexity of different vessels. Integrating local grayscale features and skeleton morphology features into the classifier for comprehensive decision-making effectively addresses the challenge of highly similar geometric shapes between arteries and veins in the terminal regions of the lung parenchyma. This significantly improves the accuracy of classifying small vessel branches, reduces classification confusion caused by relying solely on global shape analysis, and ultimately yields spatially coherent and reliable overall extraction results for vessel attribution. Attached Figure Description

[0022] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0023] Figure 1 This is a flowchart of a method for extracting pulmonary arteries and veins based on CT images;

[0024] Figure 2 This is a flowchart of the anisotropic diffusion filtering process for 3D CT images;

[0025] Figure 3 This is a flowchart of lung parenchyma mask extraction;

[0026] Figure 4 This is a flowchart of the multi-scale tubular structure reinforcement process;

[0027] Figure 5This is a flowchart of the pulmonary artery and pulmonary vein candidate region extraction process;

[0028] Figure 6 This is a flowchart of the classification of candidate regions for pulmonary arteries and pulmonary veins;

[0029] Figure 7 This is a schematic diagram of pulmonary artery and pulmonary vein support vector machine classification based on grayscale distribution features. Detailed Implementation

[0030] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0031] See Figure 1 This invention provides a method for extracting pulmonary arteries and veins based on CT images, comprising the following steps: performing anisotropic diffusion filtering on a three-dimensional CT image to obtain a smooth image; segmenting the lung parenchyma region from the smooth image to obtain a lung parenchyma mask; performing a masking operation on the smooth image based on the lung parenchyma mask to obtain a lung parenchyma image; calculating the Hessian matrix on the lung parenchyma image to obtain Hessian matrix eigenvalues ​​at multiple scales; calculating the tubular structure enhancement response value of each voxel point based on the Hessian matrix eigenvalues ​​at each scale to obtain a multi-scale tubular structure enhancement image; performing connected component analysis on the multi-scale tubular structure enhancement image, and selecting candidate regions for pulmonary arteries and pulmonary veins based on the size and shape features of the connected components; and distinguishing pulmonary arteries and pulmonary veins within the candidate regions for pulmonary arteries and pulmonary veins using local grayscale distribution and morphological features to obtain pulmonary artery extraction results and pulmonary vein extraction results.

[0032] Example 1:

[0033] In specific implementation, please refer to Figure 2 The 3D CT images are derived from CT scanning equipment, and the voxel gray values ​​correspond to the X-ray attenuation coefficients. The implementation method of anisotropic diffusion filtering is explained in detail.

[0034] First, define the diffusion time step and the number of iterations. The diffusion time step is denoted as... Its value affects the magnitude of grayscale value updates in each iteration. Diffusion time step The value ranges from 0.1 to 0.25, and in one embodiment it is set to... The number of iterations is denoted as... Determine based on the desired smoothness. The value ranges from 5 to 20, and in one implementation, it is set to... The 3D CT image to be processed is used as the initial image. .

[0035] The process proceeds iteratively, with the current image as input and the updated image as output for each iteration. For the... iteration Increasing from 1 to The current image is denoted as Based on the current image Calculate the diffusion coefficient for each voxel, then update the gray value of that voxel based on the diffusion coefficient, and generate the updated image. . No. After the next iteration is completed, As a smoothed image.

[0036] The diffusion coefficient of each voxel is calculated as follows: In three-dimensional space, for the current image... The median coordinate is Given a voxel point, calculate the gray-level gradient vector of that voxel point. The gray-level gradient vector is a three-dimensional vector, whose three components are along... First-order partial derivative of grayscale along the axis ,along First-order partial derivative of grayscale along the axis and along First-order partial derivative of grayscale along the axis The partial derivatives are approximated using the central difference method:

[0037]

[0038]

[0039]

[0040] Calculate the magnitude of the gray-level gradient vector Its expression is .

[0041] Based on the magnitude of the gray-level gradient vector Calculate the diffusion coefficient of voxel points Diffusion coefficient The calculation formula is:

[0042]

[0043] in, To preset the gradient threshold, The magnitude of the aforementioned gray-level gradient vector, The diffusion coefficient is the preset gradient threshold. Used to distinguish between flat areas and edge areas, when much smaller When, diffusion coefficient It tends towards the first preset value, which is 1; when Much larger When, diffusion coefficient It tends towards the second preset value, which is 0. (Preset gradient threshold) The gradient threshold is set based on the grayscale gradient characteristics of the vessel wall and lumen in 3D CT images. In contrast-enhanced CT images of lung vessels, the gradient modulus at the vessel wall is typically distributed between 50 and 100, therefore a preset gradient threshold is used. The value is set to 70. Using the above diffusion coefficient formula, a diffusion coefficient close to 1 can be obtained inside the blood vessel lumen where the grayscale changes are gradual, achieving strong smoothness; while a diffusion coefficient close to 0 can be obtained at the blood vessel wall where the grayscale changes drastically, preserving edge information.

[0044] Diffusion coefficient at voxel point Then, the grayscale value of the voxel is updated. The update method is: based on the coordinates... Centered on the voxel point, consider the orientations of its six nearest neighbor voxels, i.e., positive and negative. Axis, positive and negative Axis and positive and negative Axial directions. For each direction, calculate the average of the diffusion coefficients of adjacent voxels and the central voxel, multiply it by the grayscale difference between adjacent voxels and the central voxel, then sum the results for all six directions and multiply by the diffusion time step. Finally, the resulting increment is added to the grayscale value of the current voxel. Specifically, the updated grayscale value... The expression is:

[0045]

[0046] in, This indicates that the coordinates in the current image are The grayscale value of the voxel point. This represents the grayscale value of that voxel in the updated image; Indicates the diffusion coefficient of the central voxel; For directional indication, Indicates along direction The coordinates of voxels adjacent to the central voxel. This represents the gray value of the adjacent voxel. This represents the diffusion coefficient of the adjacent voxel point; This is the preset diffusion time step. Traverse all voxels in the current image and perform the above update operation to obtain the updated image. .

[0047] Repeat the above iterative process until the number of iterations reaches... Stop the iteration when the time is right, and output the final image. The output is a smoothed image. This smoothed image maintains clear blood vessel boundaries while suppressing grayscale fluctuations caused by noise within the lumen and lung parenchyma.

[0048] Example 2:

[0049] In specific implementation, please refer to Figure 3 The process of segmenting the lung parenchyma region from a smooth image to obtain a lung parenchyma mask includes four processing steps: grayscale thresholding, morphological opening operation, connected component analysis, and morphological closing operation.

[0050] The smoothed image is segmented using grayscale thresholding to obtain the initial lung parenchyma region. The CT value of the lung parenchyma region falls within a specific range in the smoothed image, determined by the attenuation characteristics of lung tissue with air. A lower limit for the grayscale threshold is set. and grayscale threshold upper limit Lower limit of grayscale threshold The value range is -1024 to -800 Henness units, which is the upper limit of the grayscale threshold. The value range is -500 to -200 Henle units. In one embodiment, the lower limit of the grayscale threshold... Set to -950 Henle units, upper limit of grayscale threshold. Set to -400 Henle units. Perform a voxel-by-voxel scan on the smoothed image, selecting images with gray values ​​greater than the lower limit of the gray value threshold. And less than the upper limit of the grayscale threshold The voxel points are labeled as foreground and the remaining voxel points are labeled as background, thereby generating a binarized initial lung parenchyma region image, where the foreground voxel value is 1 and the background voxel value is 0.

[0051] A morphological opening operation is performed on the initial lung parenchyma region to separate the lung parenchyma from the surrounding tissue. The morphological opening operation consists of an erosion followed by an expansion operation. The erosion operation uses a three-dimensional spherical structuring element, the radius of which is denoted as [missing information]. ,radius The value is 3 voxels. During the erosion process, a 3D spherical structuring element is used as a template to traverse each foreground voxel in the initial lung parenchyma image. When all voxels within the area covered by the 3D spherical structuring element are foreground, that voxel is retained as foreground; otherwise, it is set as background, resulting in the eroded image. Subsequently, a dilation operation is performed on the eroded image. The dilation operation uses the same 3D spherical structuring element to traverse each voxel in the eroded image. When at least one foreground voxel exists within the area covered by the 3D spherical structuring element, the voxel corresponding to the center of the structuring element is set as foreground; otherwise, it is set as background, resulting in a binary image after morphological opening. Morphological opening can break the fine connections between the lung parenchyma and surrounding tissues such as the trachea and blood vessel roots, making the lung parenchyma region a relatively independent connected body.

[0052] Connectivity analysis was performed on the regions after morphological opening, and the two largest connected regions were selected as the left and right lung parenchyma regions. The connectivity analysis adopted the three-dimensional six-adjacent connectivity criterion, meaning that two foreground voxels are considered connected if they share a face in three-dimensional space. The binary image after morphological opening was traversed, and connected region labels were assigned to all foreground voxels. The number of voxels contained in each connected region label was counted, and this number was used as the volume of the connected region. The connected regions were sorted from largest to smallest volume, and the two largest connected regions were extracted. The labels corresponding to these two connected regions were retained as foreground, and the voxels corresponding to the labels of other connected regions were set as background, resulting in the binary images of the left and right lung parenchyma regions.

[0053] A morphological closure operation is performed on the left and right lung parenchyma regions to fill cavities within the lung parenchyma, obtaining a lung parenchyma mask. The morphological closure operation consists of a dilation followed by an erosion operation. The dilation operation uses a three-dimensional spherical structuring element, the radius of which is denoted as [missing information]. ,radius The value is taken from 5 voxels. Radius The rationale for setting the radius to 5 voxels is that cavities with diameters ranging from 2 to 8 voxels may form within the lung parenchyma due to the passage of blood vessels and bronchi. Choosing a radius of 5 voxels allows for the coverage of most small cavities while avoiding overfilling that could cause the lung parenchyma boundary to expand outwards. During the expansion process, this 3D spherical structuring element is used to traverse the binary images of the left and right lung parenchyma regions. The center of the location containing a foreground voxel within the area covered by the structuring element is placed as the foreground, filling small cavities while slightly expanding the outer boundary of the lung parenchyma. Subsequently, an erosion operation is performed on the expanded image, using the same radius. A three-dimensional spherical structuring element is used to traverse the dilated image. Only the central voxel is retained as foreground if all voxels within the structuring element's coverage area are foreground; otherwise, it is set as background. After erosion, the outer boundary of the lung parenchyma is restored to a position close to its pre-dilation state, while the internal cavities remain filled. The foreground region obtained after morphological closing is used as a lung parenchyma mask. The lung parenchyma mask is a binary three-dimensional image, with foreground voxels labeled 1 to represent lung parenchyma regions and background voxels labeled 0 to represent non-lung parenchyma regions.

[0054] Example 3:

[0055] In specific implementation, please refer to Figure 4 The lung parenchyma image is obtained by masking the smoothed image using a lung parenchyma mask, setting the gray values ​​of voxels outside the lung parenchyma mask in the smoothed image to 0. The lung parenchyma image retains only the gray-level information within the lung parenchyma region, and subsequent Hessian matrix calculations are limited to non-zero voxel points in the lung parenchyma image.

[0056] A set of Gaussian kernel scaling parameters is defined to construct a multi-scale space. The Gaussian kernel scaling parameters are denoted as […]. , Controls the width of the Gaussian smoothing, in units of voxel side length. (Settings...) Each scale parameter is denoted as... The selection of the scale parameter set needs to cover the diameter range of pulmonary vessels, which ranges from 2 to 10 voxels. Therefore, the scale parameter set takes values ​​of {1.0, 1.5, 2.0, 2.5, 3.0}, for a total of 5 scales. Iterate through each scale parameter and perform the following operations for each scale parameter.

[0057] For the current scale parameters The lung parenchyma image is then convolved and filtered using a Gaussian kernel corresponding to this scale to obtain the filtered image. A three-dimensional Gaussian kernel is constructed. Its expression is:

[0058]

[0059] in, , , This represents the voxel offset within the Gaussian kernel relative to the center point. This is the current scale parameter. The radius of the Gaussian kernel is taken as... To ensure that the Gaussian kernel covers more than 99% of the energy, the constructed Gaussian kernel is convolved with the lung parenchyma image in three dimensions to obtain the filtered image. The specific operation of 3D convolution is as follows: for each coordinate in the lung parenchyma image... For each voxel point, align the center of the Gaussian kernel with that voxel point, calculate the weighted sum of the voxel values ​​and Gaussian kernel coefficients within the Gaussian kernel coverage area, and use this weighted sum as the filtered image. The grayscale value of this voxel. For voxel points near the boundary in the lung parenchyma image, a filling method using copied boundary pixels is adopted.

[0060] Filtered image Calculate the second-order partial derivative at each voxel point to construct the Hessian matrix for that voxel point. For coordinates... At the voxel point, six second-order partial derivatives need to be calculated: , , , , , All second-order partial derivatives are calculated using the central difference approximation. For example, its approximate calculation formula is: Cross-partial derivatives The approximate calculation formula is:

[0061] The other four second-order partial derivatives were calculated in the same way. The six second-order partial derivatives are then organized into one... The real symmetric matrix is ​​the value of the current voxel at the scale. Hessian matrix :

[0062]

[0063] Solving the Hessian matrix The eigenvalues ​​of the matrix are given. Since the Hessian matrix is ​​a real symmetric matrix, all its eigenvalues ​​are real numbers. The Jacobi iteration method is used to find three eigenvalues, which are denoted as . , , Store each voxel at the scale. The three characteristic values ​​are used for subsequent calculation of tubular response values.

[0064] After calculating the eigenvalues ​​of the Hessian matrix at all scales, for the coordinates in the lung parenchyma image For each voxel point, obtain the three Hessian matrix eigenvalues ​​of that voxel point at each scale. For each scale... , to obtain eigenvalues , , Sort these three feature values ​​in ascending order according to their absolute values ​​to obtain the sorted feature values, denoted as . , , ,satisfy .in The first feature value after sorting. The second feature value after sorting. This is the third characteristic value after sorting.

[0065] Based on the first characteristic value after sorting With the second eigenvalue The ratio, and the first eigenvalue With the third eigenvalue The ratio is used to calculate the scale of the voxel point. tubular response value Tubular response value The calculation formula is:

[0066]

[0067] in, Representing voxel points at scale The tubular response value ranges from 0 to 1. The closer the value is to 1, the more tubular the structure of the voxel point is at that scale. This represents the absolute value of the first feature after sorting. This represents the absolute value of the second feature after sorting. This represents the absolute value of the third feature after sorting. This is the cross-sectional roundness sensitivity coefficient, which controls the degree of suppression for non-circular cross-sections. The value ranges from 0.3 to 0.7. In one embodiment... Take a fixed value of 0.5. The basis for choosing 0.5 is: for an ideal tubular cross-section, the eigenvalue satisfies When the ratio is close to 1, substituting it into the formula results in an exponential term close to 1, maintaining a high response. When the cross-section deviates from a circular shape, the ratio decreases, the exponential term decays, and the decay rate changes from... When the control value is set to 0.5, it can provide a balanced assessment of tubular structures across a wide range of vascular morphologies. The spherical structure suppression coefficient controls the suppression intensity on the spherical structure. The value ranges from 0.3 to 0.7. In one embodiment... Take a fixed value of 0.5. The rationale for choosing 0.5 is that, for spherical structures, the three eigenvalues ​​are of similar magnitude. Approximately 1, in the formula When the term approaches 0, the overall response is suppressed; for tubular structures As the value approaches 0, this term approaches 1. Taking 0.5 allows the filter to achieve a better balance in distinguishing between spherical and tubular structures. Indicated by An exponential function with base 0; the first term of the product This term measures the degree to which the cross-sectional shape deviates from a circle; for tubular structures, this term tends to 1, while for planar or flat structures, it decreases rapidly; the second term of the product... This measure of grayscale change along the principal axis is used to distinguish between tubular and spherical structures; this term tends to 1 for tubular structures and to 0 for spherical structures. Multiplying these two terms ensures that the structure only satisfies the condition of having both a cylindrical cross-section and low curvature along the principal axis. Achieve a higher value.

[0068] Traversing the tubular response values ​​at all scales of this voxel point The maximum value is selected as the final tubular structure enhancement response value for that voxel point, denoted as . The operation of taking the maximum value at multiple scales allows each voxel point to be automatically matched to the scale that best suits its local blood vessel diameter.

[0069] Iterate through all non-zero voxels in the lung parenchyma image, performing the above operation for each voxel to generate a 3D image of the same size as the lung parenchyma image. The value of each voxel is... This image is a multi-scale tubular structure enhanced image. In the multi-scale tubular structure enhanced image, the high-response areas correspond to the pulmonary artery and pulmonary vein regions with tubular morphology, while the background and other tissues of the lung parenchyma have lower response values.

[0070] Example 4:

[0071] In specific implementation, please refer to Figure 5 The multi-scale tubular structure enhancement image is binarized to obtain a binarized tubular structure image. The grayscale value of each voxel in the multi-scale tubular structure enhancement image is the final tubular structure enhancement response value for that voxel; regions with higher response values ​​are more likely to exhibit a tubular morphology. A binarization threshold is set. Binarization threshold The value ranges from 0.05 to 0.25. In one implementation, the binarization threshold... Set to 0.1. Binarization threshold. The rationale for setting the threshold to 0.1 is as follows: the distribution range of the tubular structure enhancement response value in multi-scale tubular structure enhanced images is 0 to 1. The response value of the lung vascular region is usually concentrated above 0.3, while the response value of the background lung parenchyma is usually concentrated below 0.05. Selecting 0.1 can preserve most of the vascular regions while eliminating background noise from the lung parenchyma, avoiding the introduction of low-response non-tubular structures into subsequent processing. For each voxel in the multi-scale tubular structure enhanced image, the tubular structure enhancement response value of that voxel is compared with a binarization threshold. Comparison: When the response value is greater than the binarization threshold If the value of the voxel point is 1 in the binary tubular structure image, then the value is 0 otherwise.

[0072] Connected component labeling was performed on the binarized tubular structure image to obtain multiple connected components. The 3D 26-adjacency connectivity criterion was used for labeling, meaning that two voxels with a value of 1 are considered connected if they share a face, edge, or vertex in 3D space. The rationale for using 3D 26-adjacency is that pulmonary vessels have continuity in 3D space, and 26-adjacency can maintain the integrity of connected components at vessel bifurcation or changes in direction, avoiding the segmentation of the same vessel into multiple isolated connected component fragments due to overly strict adjacency definitions. A two-pass scanning method was used for connected component labeling: the first scan traversed the binarized tubular structure image voxel by voxel, assigning a temporary label to each voxel with a value of 1, and recording the temporary labels of adjacent voxels as equivalence relations; the second scan replaced all temporary labels in the equivalence relations with a unified minimum temporary label, forming the final connected component labels. A 3D image with connected component labels was output, with each connected component identified by a unique label value.

[0073] For each connected component in a 3D image with connected component labels, calculate the number of voxels, the major axis length, and the minor axis length of that connected component. (Number of voxels) This is obtained by counting the number of voxels within the connected component. For a connected component... The number of its voxels is denoted as .

[0074] The calculation of the major and minor axis lengths is based on the covariance matrix of the connected component in three-dimensional space. For a connected component... Extract the 3D coordinates of all voxel points within the connected component and construct a coordinate set. ,in Represents connected components The Middle Individual primal points in The coordinate values ​​along the axis. Represents connected components The Middle Individual primal points in The coordinate values ​​along the axis. Represents connected components The Middle Individual primal points in The coordinate values ​​along the axis. Represents connected components The number of voxels. Calculate the connected components. centroid coordinates ,in:

[0075] , , .

[0076] Construct connected components The three-dimensional covariance matrix :

[0077]

[0078] in, Represents connected components The three-dimensional covariance matrix is A real symmetric matrix; Represents connected components The Middle Individual primal points in Coordinate values ​​along the axis; Represents connected components The Middle Individual primal points in Coordinate values ​​along the axis; Represents connected components The Middle Individual primal points in Coordinate values ​​along the axis; Represents connected components All voxel points The mean value of the coordinates along the axis; Represents connected components All voxel points The mean value of the coordinates along the axis; Represents connected components All voxel points The mean value of the coordinates along the axis; Represents connected components The number of voxels.

[0079] For the three-dimensional covariance matrix Perform eigenvalue decomposition to obtain the three-dimensional covariance matrix. Let be a real symmetric matrix with three non-negative real eigenvalues. Arrange the three eigenvalues ​​in descending order, denoted as . ,in Represents connected components The largest eigenvalue of the covariance matrix, Represents connected components The intermediate eigenvalues ​​of the covariance matrix, Represents connected components The smallest eigenvalue of the covariance matrix. Connected region. Major axis length and Proportional, connected components minor axis length and Proportional to the length of the major axis. Set as minor axis length Set as The basis for using a coefficient of 4 is that, under the assumption of a multivariate normal distribution, an interval of 4 times the standard deviation covers approximately 95% of the distribution range, which can better reflect the extent of the spread of the connected domain along each principal axis.

[0080] Iterate through all connected components to obtain the number of voxels in each component. Major axis length and minor axis length .

[0081] Connected components with a number of voxels greater than a preset voxel count threshold are selected as initial connected components. The preset voxel count threshold is denoted as... , The value ranges from 100 to 500. In one embodiment... Set to 200. Preset voxel count threshold. The threshold of 200 is set as follows: the number of connected voxels in pulmonary vessels in 3D CT images is typically in the hundreds to thousands, while the number of connected voxels caused by noise is usually less than 100. Selecting 200 as the threshold effectively filters out noise-induced connected regions while retaining meaningful vascular connected regions. For connected regions... ,when When, the connected component is included in the initial connected component set; when When this happens, the connected component is discarded.

[0082] Connected regions whose ratio of major axis length to minor axis length is greater than a preset aspect ratio threshold are selected as candidate regions for pulmonary arteries and pulmonary veins. For each connected region in the initial screening set... Calculate the length of the major axis With minor axis length ratio The preset aspect ratio threshold is denoted as... , The value ranges from 2.5 to 5.0. In one embodiment... Set to 3.0. Preset aspect ratio threshold. The rationale for setting it to 3.0 is that pulmonary vessels present as slender tubular structures with a length-to-width ratio typically greater than 3, while non-tubular structures such as pulmonary nodules or lymph nodes have a length-to-width ratio close to 1. Selecting 3.0 effectively distinguishes between tubular vascular structures and non-tubular mass structures. When, mark the connected region as a candidate region for pulmonary artery and pulmonary vein; when If the connection is not found, the connected region is discarded. All the final candidate regions for pulmonary arteries and pulmonary veins constitute a candidate region set, which is used in subsequent arteriovenous differentiation steps.

[0083] Example 5:

[0084] In specific implementation, please refer to Figure 6 After obtaining the candidate regions of pulmonary artery and pulmonary vein, each candidate region is classified using local grayscale distribution features and morphological features within the candidate regions of pulmonary artery and pulmonary vein to distinguish between pulmonary artery and pulmonary vein.

[0085] For each pulmonary artery and pulmonary vein candidate region, the CT values ​​of all voxel points within that region are extracted from the original 3D CT image. The original 3D CT image contains complete grayscale information and has not been altered by any smoothing or enhancement processing. Using a binary mask of the pulmonary artery and pulmonary vein candidate regions, the voxel CT values ​​in the original 3D CT image corresponding to the mask coverage area are read as a set of values. The mean of this set of values ​​is calculated. and standard deviation , and The grayscale distribution characteristics constituting the candidate regions of the pulmonary artery and pulmonary vein. Mean The standard deviation reflects the average attenuation level of blood against the vessel wall within the candidate region of the pulmonary artery and pulmonary vein. This reflects the degree of dispersion of CT values ​​within the candidate regions of the pulmonary artery and pulmonary vein.

[0086] In practice, morphological features are obtained by calculating the skeleton lines of candidate pulmonary artery and pulmonary vein regions. Morphological thinning is performed on the binarized images of the candidate pulmonary artery and pulmonary vein regions to obtain skeleton line images. The morphological thinning operation employs a three-dimensional thinning algorithm, iteratively removing boundary voxel points while maintaining the topological structure until skeleton lines of single-voxel width are obtained. Intersection points of the skeleton lines are identified in the skeleton line image. The method for determining skeleton line intersection points is as follows: for coordinates on the skeleton line... For a voxel point at a given location, examine the number of skeleton line voxels connected to that voxel point within its 26-neighborhood in 3D. If the number of connected voxels is greater than or equal to 3, the voxel point is identified as a skeleton line intersection; if the number is 2, it is a non-intersection skeleton line voxel point; and if the number is 1, it is a skeleton line endpoint. Count the number of voxels identified as skeleton line intersections; this count is recorded as the number of bifurcation points. .

[0087] In the skeleton line image, calculate the Euclidean distance from each skeleton line intersection to its adjacent skeleton line intersections. Construct an undirected graph with all skeleton line intersections as nodes and skeleton line paths as edges. For each skeleton line intersection, its adjacent skeleton line intersections are defined as other skeleton line intersections directly connected to that intersection in the undirected graph. Intersection with its adjacent skeleton line Calculate the Euclidean distance between two points. Euclidean distance The calculation formula is:

[0088]

[0089] in, Indicates the intersection of skeleton lines Intersection with skeleton line The Euclidean distance between them; Indicates the intersection of skeleton lines In the coordinate system of a three-dimensional CT image Axis coordinate values; Indicates the intersection of skeleton lines exist Axis coordinate values; Indicates the intersection of skeleton lines exist Axis coordinate values; Indicates the intersection of skeleton lines exist Axis coordinate values; Indicates the intersection of skeleton lines exist Axis coordinate values; Indicates the intersection of skeleton lines exist Axis coordinate values. Traverse all pairs of intersections of adjacent skeleton lines in the undirected graph, calculate all Euclidean distances, and take the arithmetic mean of all calculated Euclidean distances to obtain the average distance between bifurcation points, denoted as . Number of branching points Average distance between the bifurcation point and the bifurcation point Morphological features constituting the candidate regions of the pulmonary artery and pulmonary vein.

[0090] In practical implementation, the grayscale distribution characteristics are determined by the mean. and standard deviation Composition, morphological characteristics are determined by the number of bifurcation points Average distance between the bifurcation point and the bifurcation point The four features are combined into a four-dimensional feature vector, which is then used as the input to the preset classifier.

[0091] The default classifier is a Support Vector Machine (SVM) classifier, which uses a Gaussian radial basis function as the kernel function. The specific steps for training the SVM classifier are as follows: a training sample set labeled with pulmonary artery and pulmonary vein tags is pre-collected. Each training sample corresponds to a candidate region for either a pulmonary artery or pulmonary vein, and the category is determined manually by a physician on a 3D CT image. Each training sample contains four feature values: mean... Standard deviation Number of branching points Average distance between the bifurcation point and the bifurcation point The support vector machine (SVM) classifier is trained by using the feature vectors of the training sample set as input and the class labels as output. During training, the SVM classifier finds the optimal classification hyperplane by solving a quadratic programming problem, with the kernel function being a Gaussian radial basis function. ,in and Representing two eigenvectors, These are kernel function parameters that control the influence range of a single sample. Set to 0.125, penalty coefficient. Set to 1.0. Kernel function parameters. The value is set to 0.125 because: after standardization, the numerical values ​​of each dimension of the eigenvectors are distributed between 0 and 1; the expected value of the Euclidean distance between eigenvectors in the four-dimensional feature space is approximately 1.0; and after substituting this into the Gaussian radial basis function... Taking 0.125 will make The value is approximately -0.125, corresponding to a kernel function value of approximately 0.882, which falls within the sensitive range. This avoids premature saturation of the kernel function value, preventing underfitting or overfitting of the model. (Penalty coefficient) The value is set to 1.0 because it maintains a balance between training error and model complexity, avoiding overfitting or underfitting. After training, a trained support vector machine classifier is obtained, which stores the support vectors and their corresponding Lagrange multipliers.

[0092] For each candidate region of pulmonary artery and pulmonary vein, the mean of the gray-level distribution features of that candidate region is calculated. Standard deviation and the number of morphological bifurcation points Average distance between bifurcation points These features are combined into a single feature vector. This feature vector is then input into a trained Support Vector Machine (SVM) classifier. The SVM classifier calculates the signed distance from this feature vector to the classification hyperplane and obtains the probability value through a sigmoid function. probability value The value ranges from 0 to 1. Set the preset classification threshold. Preset classification threshold The value is 0.5. Preset classification threshold. The basis for setting it to 0.5 is: when the probability value When the probability value is greater than 0.5, the support vector machine classifier is more inclined to classify the sample as positive, and the classification boundary is an unbiased decision boundary. When the probability value is [not specified], the category label of the candidate regions for pulmonary artery and pulmonary vein is determined to be pulmonary artery; when [not specified], the category label of the candidate regions for pulmonary artery and pulmonary vein is determined to be pulmonary artery. When this occurs, the category label of the candidate regions for pulmonary artery and pulmonary vein is determined to be pulmonary vein.

[0093] Based on category labels, candidate regions for pulmonary arteries and pulmonary veins are categorized into pulmonary artery sets and pulmonary vein sets, respectively. After traversing all candidate regions for pulmonary arteries and pulmonary veins, the pulmonary artery set contains all candidate regions identified as pulmonary arteries, and the pulmonary vein set contains all candidate regions identified as pulmonary veins. The pulmonary artery set is used as the pulmonary artery extraction result, and the pulmonary vein set is used as the pulmonary vein extraction result.

[0094] See Figure 7 This figure illustrates the classification distribution of candidate regions for pulmonary arteries and pulmonary veins in the grayscale feature space. The horizontal axis represents the mean CT value M (in HU, ranging from approximately 150 to 600) of the candidate regions, and the vertical axis represents the standard deviation S (in HU, ranging from approximately 20 to 140). In the figure, circles represent pulmonary artery samples, squares represent pulmonary vein samples, and the dashed line represents the decision boundary of the trained support vector machine classifier.

[0095] As shown in the figure, pulmonary artery samples mainly clustered in the range of mean (M) approximately 350 to 550 HU, with standard deviation (S) concentrated in the range of 40 to 100 HU, exhibiting moderate mean and low to moderate dispersion. Pulmonary vein samples were mainly distributed in the range of mean (M) approximately 200 to 400 HU, with standard deviation (S) concentrated in the range of 70 to 130 HU, exhibiting lower mean and higher dispersion. The two types of samples showed clear distinction in mean and standard deviation space; pulmonary artery samples had a higher overall mean and lower standard deviation, while pulmonary vein samples had a lower mean and higher standard deviation.

[0096] The support vector machine (SVM) classification decision boundary exhibits a slightly non-linear trend, reaching a minimum standard deviation (S) of approximately 70 HU around the mean M (around 350-400 HU). This boundary effectively separates the two classes of samples in the image. This decision boundary setting effectively utilizes the combined features of the mean M and the standard deviation S, achieving the distinction between pulmonary artery and pulmonary vein candidate regions, and providing a basis for subsequent arteriovenous determination of candidate regions based on this classifier.

[0097] Overall, the data in this figure intuitively reflects the technical solution of the support vector machine classifier in Example 5, which uses gray-level distribution features to effectively distinguish between pulmonary arteries and pulmonary veins. It verifies the importance of gray-level mean and standard deviation as classification features and the rationality of the classification boundary.

[0098] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for extracting pulmonary arteries and veins based on CT images, characterized in that, Includes the following steps: Anisotropic diffusion filtering is applied to the 3D CT image to obtain a smooth image. The lung parenchyma region is then segmented from the smooth image to obtain a lung parenchyma mask. A smooth image is masked using a lung parenchyma mask to obtain a lung parenchyma image. The Hessian matrix is ​​then calculated on the lung parenchyma image to obtain the eigenvalues ​​of the Hessian matrix at multiple scales. The tubular structure enhancement response value of each voxel is calculated based on the eigenvalues ​​of the Hessian matrix at each scale to obtain a multi-scale tubular structure enhancement image; Connectivity analysis was performed on the enhanced images of multi-scale tubular structures. Candidate regions for pulmonary arteries and pulmonary veins were selected based on the size and shape features of the connected regions. The pulmonary arteries and pulmonary veins were distinguished by local gray-scale distribution and morphological features, and the extraction results of pulmonary arteries and pulmonary veins were obtained.

2. The method for extracting pulmonary arteries and veins based on CT images according to claim 1, characterized in that, The specific steps for performing anisotropic diffusion filtering on the 3D CT image to obtain a smooth image are as follows: Using 3D CT images as initial images, the diffusion time step and number of iterations are set. In each iteration, the diffusion coefficient of each voxel is calculated based on the gray-level gradient of the current image; The grayscale value of each voxel in the current image is updated based on the diffusion coefficient to obtain the updated image; The iteration process is repeated until the number of iterations is reached, and the finally obtained updated image is used as the smoothed image.

3. The method for extracting pulmonary arteries and veins based on CT images according to claim 2, characterized in that, The specific steps for calculating the diffusion coefficient of each voxel based on the gray-level gradient of the current image are as follows: For each voxel in the current image, calculate the gray-level gradient vector of that voxel in three-dimensional space; The diffusion coefficient of the voxel is calculated based on the magnitude of the gray-level gradient vector. Specifically, when the magnitude of the gray-level gradient vector is less than the preset gradient threshold, the diffusion coefficient tends to the first preset value; when the magnitude of the gray-level gradient vector is greater than the preset gradient threshold, the diffusion coefficient tends to the second preset value.

4. The method for extracting pulmonary arteries and veins based on CT images according to claim 1, characterized in that, The specific steps for segmenting the lung parenchyma region from the smooth image to obtain the lung parenchyma mask are as follows: The smoothed image is segmented by grayscale thresholding to obtain the initial lung parenchyma region; Morphological opening operations were performed on the initial lung parenchyma region to separate the lung parenchyma from the surrounding tissues; Connectivity analysis was performed on the regions after morphological opening, and the two largest connected regions were selected as the left and right lung parenchyma regions. Morphological closure operations are performed on the left and right lung parenchyma regions to fill the cavities inside the lung parenchyma and obtain a lung parenchyma mask.

5. The method for extracting pulmonary arteries and veins based on CT images according to claim 1, characterized in that, The specific steps for calculating the Hessian matrix on lung parenchyma images and obtaining the eigenvalues ​​of the Hessian matrix at multiple scales are as follows: Define a set of Gaussian kernel scaling parameters and iterate through each scaling parameter; For the current scale parameter, the lung parenchyma image is convolved and filtered using the Gaussian kernel corresponding to that scale to obtain the filtered image. Calculate the second-order partial derivative at each voxel point in the filtered image and construct the Hessian matrix for that voxel point; Solve for the eigenvalues ​​of the Hessian matrix to obtain the three eigenvalues ​​of the voxel at the current scale.

6. The method for extracting pulmonary arteries and veins based on CT images according to claim 5, characterized in that, The specific steps for calculating the tubular structure enhancement response value of each voxel point based on the Hessian matrix eigenvalues ​​at each scale to obtain a multi-scale tubular structure enhancement image are as follows: For each voxel, obtain the three Hessian matrix eigenvalues ​​of that voxel at each scale; The three feature values ​​at each scale are sorted according to their absolute values ​​to obtain the sorted feature values; The tubular response value of the voxel at this scale is calculated based on the ratio of the first eigenvalue to the second eigenvalue and the ratio of the first eigenvalue to the third eigenvalue in the sorted eigenvalues. Traverse all scales to obtain the tubular response value of the voxel at each scale, and take the maximum value of the tubular response value among all scales as the final tubular structure enhancement response value of the voxel. By traversing all voxel points in the lung parenchyma image, a multi-scale tubular structure enhancement image is obtained.

7. The method for extracting pulmonary arteries and veins based on CT images according to claim 1, characterized in that, The specific steps for performing connected component analysis on the enhanced multi-scale tubular structure image, and selecting candidate regions for pulmonary arteries and pulmonary veins based on the size and shape features of the connected components, are as follows: Binarization processing is performed on the enhanced image of the multi-scale tubular structure to obtain a binarized tubular structure image; Connected component labeling is performed on the binarized tubular structure image to obtain multiple connected components; Calculate the number of voxels, major axis length, and minor axis length for each connected component; Connected components with a number greater than a preset voxel number threshold are selected based on the number of voxels and used as initial connected components. Connected regions whose ratio of major axis length to minor axis length is greater than a preset aspect ratio threshold are selected as candidate regions for pulmonary arteries and pulmonary veins.

8. The method for extracting pulmonary arteries and veins based on CT images according to claim 1, characterized in that, The specific steps for distinguishing pulmonary arteries and pulmonary veins within the candidate regions of pulmonary arteries and pulmonary veins using local grayscale distribution and morphological features to obtain pulmonary artery extraction results and pulmonary vein extraction results are as follows: On the original 3D CT image, extract the CT values ​​of all voxel points in each candidate region of pulmonary artery and pulmonary vein; Calculate the mean and standard deviation of the voxel CT values ​​within each candidate region to obtain the grayscale distribution characteristics of that candidate region; For each candidate region, calculate the skeleton line of the candidate region, obtain the number of bifurcation points of the skeleton line and the average distance between the bifurcation points, and obtain the morphological features of the candidate region. The grayscale distribution features and morphological features of each candidate region are input into a preset classifier to obtain the category label of the candidate region, wherein the category label is pulmonary artery or pulmonary vein; Candidate regions are categorized into pulmonary artery and pulmonary vein collections based on category labels, resulting in pulmonary artery extraction and pulmonary vein extraction results.

9. The method for extracting pulmonary arteries and veins based on CT images according to claim 8, characterized in that, The specific steps for calculating the skeleton line of each candidate region, obtaining the number of bifurcation points of the skeleton line and the average distance between the bifurcation points, and obtaining the morphological features of the candidate region are as follows: For each candidate region, perform morphological thinning on the binarized image of the candidate region to obtain the skeleton line image of the candidate region; In the skeleton line image, identify the intersections of the skeleton lines and count the number of skeleton line intersections as the number of bifurcation points. In the skeleton line image, calculate the Euclidean distance between each skeleton line intersection point and its adjacent skeleton line intersection point, and average all the calculated Euclidean distances to obtain the average distance between the bifurcation points.

10. The method for extracting pulmonary arteries and veins based on CT images according to claim 8, characterized in that, The specific steps for inputting the grayscale distribution features and morphological features of each candidate region into a preset classifier to obtain the category label of the candidate region are as follows: A training sample set labeled with pulmonary artery and pulmonary vein tags is collected in advance. Each training sample contains grayscale distribution features and morphological features. A support vector machine classifier is trained using the training sample set to obtain the trained classifier; The grayscale distribution features and morphological features of each candidate region are combined into a feature vector, and the feature vector is input into the trained classifier. Obtain the probability value output by the classifier. When the probability value is greater than the preset classification threshold, the category label of the candidate region is determined to be pulmonary artery; otherwise, it is determined to be pulmonary vein.

Images