Automatic Vascular Extraction Method and System Based on CT Images

Abstract

This invention relates to the field of medical image processing technology, and particularly to an automatic vascular extraction method and system based on CT images. The method includes: Step 1: Construction of a local vascular coordinate system and dual-domain separation; Step 2: Dynamic dual-domain fusion, performing contrast-adaptive weighted fusion of the spatial domain enhanced image and the fundamental frequency vascular image; Step 3: Orientation-constrained morphological reconstruction: constructing structural elements with their major axes parallel to the vascular direction based on the vascular principal direction vector v; sequentially performing directional opening and closing operations along the vascular principal direction vector v; Step 4: Anatomy-driven vascular repair and correction: based on the generated calcified region mask Mcal marking region, tracing the contours of adjacent slices along the vascular principal direction vector v and fitting a repair path, correcting the topology at bifurcation points based on branch angles and hemodynamic constraints. By comprehensively utilizing repair techniques such as spatial domain enhancement and frequency domain separation, vascular structures can be accurately and stably extracted from CT images.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, and in particular to a method and system for automatic blood vessel extraction based on CT images. Background Technology

[0002] In the field of image post-processing for coronary CT angiography (CCTA), accurate extraction of vascular structures is a crucial prerequisite for diagnosing coronary artery disease and assessing the degree of stenosis. However, existing technologies face the following inherent limitations:

[0003] Coronary artery calcification (CAC) is a typical manifestation of atherosclerosis, with its CT value (usually >450 HU) significantly overlapping with the vascular lumen (80-120 HU). Traditional threshold segmentation methods (such as region growing and level sets) face a dilemma in calcified areas:

[0004] Setting a high threshold to exclude calcification will lead to the rupture of blood vessels distal to the calcification (false detection rate >40%).

[0005] If a low threshold is used, calcified plaques adhere to the blood vessel wall, causing false connectivity (false detection rate of up to 35%).

[0006] Furthermore, calcification can account for up to 70% of the cross-sectional area of ​​a blood vessel, and its distribution is irregular. This problem is particularly prominent in elderly patients and diabetic patients. A ruptured blood vessel centerline can lead to an error of more than 30% in the calculation of the stenosis rate, which directly affects the decision-making process for revascularization surgery.

[0007] Furthermore, the instantaneous displacement during cardiac systole can reach 15-20 mm, resulting in a step-like artifact at the vessel edges. The limitations of existing morphological methods (such as isotropic expansion corrosion) are:

[0008] Unable to distinguish between artifact-induced "spiculated" structures and real branch vessels; the morphological difference of the same vessel is >50% in diastolic and systolic images.

[0009] Coronary arteries move periodically with the heart, and their displacement trajectory is non-linear. Traditional methods rely on static image processing, which directly leads to artifacts being misidentified as vascular branches, increasing the misdiagnosis rate of bifurcation lesions by 25% and affecting the selection of stent size.

[0010] Therefore, there is an urgent need for an automatic blood vessel extraction method and system based on CT images to solve the above problems. Summary of the Invention

[0011] To achieve the above objectives, the present invention provides a method and system for automatic blood vessel extraction based on CT images, wherein the method for automatic blood vessel extraction based on CT images includes:

[0012] Step 1: Construction of the local coordinate system of the blood vessel and separation of the two domains:

[0013] Input CT angiography data and determine the principal direction vector v of the blood vessel through Hessian matrix feature analysis;

[0014] Construct a local cylindrical coordinate system with the principal direction vector v of the blood vessel as the axis, and simultaneously execute the following steps:

[0015] a: Spatial domain vessel enhancement filtering: A Frangi filter with dynamically adjusted radius is used, whose radius adaptively decreases with the depth of the vessel layer to generate a spatial domain enhancement map;

[0016] b: Frequency domain calcification feature separation: Standard vascular frequency band and calcification feature frequency band are separated by fast Fourier transform to generate a fundamental frequency vascular map;

[0017] Calcification region mask Mcal is generated based on calcification feature frequency bands;

[0018] Step 2: Dynamic dual-domain fusion, which involves adaptively weighting and fusing the spatial domain enhanced map with the fundamental frequency angiography.

[0019] Step 3: Orientation-constrained morphological reconstruction:

[0020] Construct a structural element whose major axis is parallel to the direction of the blood vessel based on the main direction vector v of the blood vessel.

[0021] Perform directional opening and directional closing operations sequentially along the principal direction vector v of the blood vessel;

[0022] Step 4: Anatomy-driven vascular repair and correction:

[0023] Based on the generated calcified region mask Mcal marking region, the contours of adjacent slices are traced along the main vascular direction vector v and the repair path is fitted. At the bifurcation point, the topology is corrected according to the branch angle and hemodynamic constraints.

[0024] Preferably, the process of dynamically adjusting the radius in step 1 includes:

[0025] Identifying vascular hierarchy: When a bifurcation feature is detected, a hierarchy increment coefficient is calculated based on the bifurcation angle, which is determined by the spatial angle between the direction vectors of adjacent branches.

[0026] The initial radius of the main trunk layer is set to the system resolution multiplied by a preset scaling factor, which is obtained by statistically analyzing the correlation between the diameter of the main trunk vessel and the resolution.

[0027] The reduction in the branch layer radius is calculated based on the hierarchical increment coefficient and the bifurcation angle: for each unit increase in the bifurcation angle, the radius is reduced according to the nonlinear decay rate corresponding to the hierarchical depth, and the nonlinear decay rate is determined by the branch vessel diameter decay model.

[0028] Preferably, the separation of the standard vascular frequency band and the calcification characteristic frequency band in step 1 includes:

[0029] Method for determining standard vascular frequency bands: In the frequency domain space, locate the energy concentration band corresponding to the preset vascular density range. Its boundary is determined by analyzing the cumulative probability of the frequency domain energy distribution of the vascular region. When the cumulative probability reaches the first critical value, the low-frequency boundary is defined, and when it reaches the second critical value, the high-frequency boundary is defined.

[0030] Method for separating calcification characteristic frequency bands: Extracting frequency domain components above a calcification density threshold, wherein the calcification density threshold is dynamically set through the following process:

[0031] a: Mark typical calcified sample areas in the airspace;

[0032] b: Calculate the ratio of the peak frequency domain energy of a typical calcified sample region to the average energy of the entire image;

[0033] c: When the ratio exceeds a preset multiple, the frequency corresponding to the peak energy in the frequency domain is used as the center of the reference frequency band.

[0034] Preferably, the contrast adaptive weighted fusion in step 2 includes:

[0035] Calculate the local grayscale standard deviation σ, and select the fusion weighting coefficient based on the interval where σ is located:

[0036] When σ is in the high noise region, the spatial enhancement map weight is set as the sum of the baseline weight and the noise compensation weight;

[0037] When σ is in the low noise range, the weight of the fundamental frequency angiography is increased proportionally to the signal-to-noise ratio.

[0038] The boundaries of the high-noise region and the low-noise region are determined in the following way:

[0039] a: Histogram of the standard deviation distribution of the background area of ​​the entire image;

[0040] b: Use the standard deviation corresponding to the peak value of the distribution as the upper limit of the low noise interval;

[0041] c: The standard deviation exceeding the upper limit of the low noise range by a preset multiple is taken as the lower limit of the high noise range.

[0042] Preferably, the construction of the structural elements in step 3 includes:

[0043] Major axis length setting: The minor axis diameter d is obtained by fitting the current vessel cross-section ellipse. The length L and d satisfy L = k × d, where k is the curvature coupling coefficient. The curvature coupling coefficient is adjusted inversely proportionally to the radius of curvature of the vessel centerline. The radius of curvature is calculated through the following process:

[0044] a: Extract the local centerline segment along the v direction;

[0045] b: Fit a cubic spline curve and calculate the curvature extrema;

[0046] c: Use the distance between extreme points as the estimated value of the radius of curvature;

[0047] The minor axis length is fixed as the product of the layer thickness and the pixel size, and its direction is determined by the secondary direction vector of the blood vessel cross section.

[0048] Preferably, step 4, which involves tracking the contours of adjacent slices and fitting the repair path, specifically includes:

[0049] Methods for tracing the contours of adjacent slices:

[0050] a: Extract the convex hull feature point set Pt from the calcified edge of the current slice;

[0051] b: Establish a search plane on adjacent slices with the projection direction of v as the normal vector;

[0052] c: Project Pt onto the search plane to generate the target region, and match the next slice contour points within the target region;

[0053] The fitting and repair path rules include:

[0054] a: Connect the matching point set to generate the initial spatial curve;

[0055] b: Calculate the rate of change of curvature of each segment of the curve. When the rate of change exceeds the reasonable anatomical threshold, insert control points for smoothing.

[0056] c: Gray-scale assignment adopts a neighborhood propagation model, with the repair path as the center line, extending to both sides with a preset width. The preset width is adaptively adjusted according to the gradient change rate of the blood vessel diameter.

[0057] Preferably, step 4, correcting the topology at the bifurcation point based on the branch angle and hemodynamic constraints, includes:

[0058] Abnormal branch angle detection:

[0059] a: Calculate the angles θ1 and θ2 between the direction vector of each branch and the parent blood vessel vector;

[0060] b: If min(θ1,θ2) < angle threshold or |θ1-θ2| > asymmetric threshold, it is judged as abnormal;

[0061] Hemodynamic reconstitution process:

[0062] a: Measure the cross-sectional area Am of the mother vessel and the cross-sectional areas A1 and A2 of its branches;

[0063] b: Calculate the cross-sectional area ratio η = (A1 + A2) / Am;

[0064] c: When η exceeds the physiologically reasonable range, the branch diameter is redistributed according to the law of blood flow conservation, specifically including:

[0065] An optimization equation for branch diameter is established, with the objective function being to minimize the variance between the cross-sectional area ratio and the physiological target value. The constraints include the range of branch angle variation, curvature continuity, and diameter attenuation rate.

[0066] d: Reconstruct the bifurcation topology using the optimized diameter.

[0067] Preferably, the angle threshold is obtained by statistically analyzing a database of healthy blood vessel bifurcations, and the specific process includes:

[0068] a: Collect bifurcation angle samples from different anatomical locations;

[0069] b: Calculate the variance of the angular distribution by grouping by blood vessel type;

[0070] c: Take the lower limit of the threshold as the value of the mean of each group minus three times the variance;

[0071] The method for setting the reasonable physiological range:

[0072] a: Analyze the relationship between the cross-sectional area ratio and pressure loss in blood flow simulation experiments;

[0073] b: When the rate of increase of pressure loss exceeds the critical slope, record the ratio boundary at this point;

[0074] c: Take the intersection of the boundaries of each anatomical location as the reasonable range for the whole.

[0075] Preferably, the fitting and repair path in step 4 further includes grayscale assignment, specifically including:

[0076] Selection of reference areas for normal blood vessels:

[0077] a: Tracing back along the v direction to the nearest uncalcified section;

[0078] b: Extract all pixels that are topologically connected to the repaired region within this section;

[0079] Weighted calculation rules:

[0080] a: Spatial weights: Distributed inversely proportional to the three-dimensional Euclidean distance;

[0081] b: Structural weights: adjusted positively based on the gradient intensity of the blood vessel wall;

[0082] c: Dynamic decay factor: It decays exponentially as the distance between the repair path and the reference section increases, and the decay rate is calibrated by the vascular grayscale propagation model.

[0083] Accordingly, embodiments of the present invention also provide an automatic blood vessel extraction system based on CT images, including a memory configured to store instructions, a processor configured to retrieve the instructions from the memory, and capable of implementing any of the automatic blood vessel extraction methods based on CT images described in the embodiments of the present invention when executing the instructions.

[0084] The beneficial effects of this invention are:

[0085] 1. The dynamic dual-domain fusion method proposed in this invention combines spatial domain enhanced maps and fundamental frequency vascular maps, and uses frequency domain calcification feature separation technology to accurately separate calcified regions from vascular lumen regions, thereby effectively avoiding missed and false detections. In particular, the frequency domain-based separation method makes the identification of calcified regions more accurate, avoiding errors caused by improper threshold selection in traditional methods, thus greatly improving the accuracy of vessel extraction.

[0086] 2. This invention utilizes vascular repair and correction technology. Based on calcified region masks and the principal direction vector of the blood vessel, it can trace the contours of adjacent slices along the principal direction of the blood vessel and fit a repair path. Particularly at bifurcation points, through hemodynamic and branch angle corrections, the topological structure of the blood vessel is further repaired. This technology effectively avoids the problem of vascular centerline breakage, ensuring the accuracy of vascular stenosis rate calculation and reducing calculation errors caused by vascular breakage.

[0087] 3. The orientation-constrained morphological reconstruction method proposed in this invention can construct directional opening and closing operations based on the principal direction vector of blood vessels, eliminating spur structures caused by artifacts and accurately extracting real blood vessel branches. Furthermore, by combining a dynamically adjusted radius Frangi filter with the principal direction vector of blood vessels, the extraction stability of blood vessels in diastolic and systolic images is effectively improved, reducing artifact misidentification and significantly lowering the misdiagnosis rate of bifurcation lesions caused by artifacts.

[0088] 4. This invention utilizes vascular repair and correction techniques, combined with hemodynamic constraints, to correct the topology at vascular bifurcation points. At the bifurcation point, the branch diameter is corrected based on the branch angle and hemodynamic constraints, avoiding an increase in the misdiagnosis rate of bifurcation lesions. By reallocating vascular branch diameters and optimizing the topology, the accuracy of branch vessel extraction is further improved, reducing the misdiagnosis rate, especially in stent size selection. Attached Figure Description

[0089] To more clearly illustrate the technical solutions in this invention or the prior art, the accompanying drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, those skilled in the art can obtain other drawings based on these drawings without creative effort.

[0090] Figure 1 This is a flowchart of the steps of the method of the present invention;

[0091] Figure 2 This is a flowchart illustrating the process of dynamically adjusting the radius as described in step 1 of the method of the present invention.

[0092] Figure 3 A flowchart illustrating the steps of setting the physiologically reasonable range for the method of the present invention. Detailed Implementation

[0093] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. It should also be noted that, to make the embodiments more comprehensive, the following embodiments are the best and preferred embodiments, and those skilled in the art can use other alternative methods to implement some well-known technologies; moreover, the accompanying drawings are only for more specific description of the embodiments and are not intended to specifically limit the present invention.

[0094] Please see Figures 1-3 This invention provides an automatic blood vessel extraction method based on CT images. In step 1, three-dimensional image data obtained through CT angiography (CCTA) is first input into the system for extraction and analysis of blood vessel structures. Hessian matrix feature analysis is a method for calculating the orientation and curvature of local structures. In this step, the blood vessel region in the CT image is analyzed using the Hessian matrix to determine the principal direction vector v of the blood vessel, which represents the principal axis direction of the blood vessel.

[0095] Once the principal direction vector v of the blood vessel is obtained, a local cylindrical coordinate system is constructed using this vector as the axis. This local coordinate system helps to better adapt to the local geometry of the blood vessel structure, providing a spatial reference for subsequent blood vessel extraction.

[0096] After establishing the local cylindrical coordinate system, a Frangi filter with dynamically adjusted radius is used for spatial enhancement. The Frangi filter's radius is adaptively adjusted, allowing for precise enhancement of vascular structural information at different vascular depth levels while suppressing noise.

[0097] The image was decomposed in the frequency domain using Fast Fourier Transform (FFT) to separate the standard vascular frequency band from the calcification feature frequency band. This frequency domain analysis effectively distinguished calcified regions (regions with high CT values) from vascular structures (regions with low CT values). The generated fundamental frequency vascular map displayed the main features of the blood vessels, while the calcification region mask Mcal precisely identified the calcified regions.

[0098] In step 2, the spatial domain enhancement map and the fundamental frequency vascular map are subjected to contrast-adaptive weighted fusion. This weighted fusion effectively combines the vascular structure obtained from spatial domain enhancement filtering with the main vascular features separated from the frequency domain, thereby improving the accuracy and robustness of vascular extraction. The adaptive weighted fusion method automatically adjusts the weights based on the image features of each region, ensuring that artifacts are avoided while enhancing the vascular structure.

[0099] In step 3, a structuring element is constructed along the major axis of the blood vessel based on the principal direction vector v. Then, orientation opening and closing operations are performed sequentially along the principal direction of the blood vessel. This operation can remove noise and artifacts from the edges of the blood vessel and smooth and repair its shape, ensuring that the true shape of the blood vessel is correctly extracted.

[0100] In step 4, the system identifies and marks calcified regions based on the generated calcified region mask Mcal. During the vessel repair phase, a continuous vessel path is fitted by tracing the contours of adjacent slices along the principal vessel direction vector v. This process ensures the consistency and accuracy of the vessel extraction results. Particularly at bifurcation points, the system further corrects the vessel topology based on the vessel branch angle and hemodynamic constraints, avoiding vessel breakage and misidentification caused by interference from calcified regions.

[0101] In one possible implementation, the vascular hierarchy is first identified by analyzing bifurcation points in the vascular structure. In angiographic images, blood vessels typically exhibit a bifurcation pattern, and these bifurcation points are crucial for identifying vascular hierarchy. By calculating the bifurcation angle at each bifurcation point, the hierarchy to which the bifurcation belongs can be determined. The bifurcation angle is calculated using the direction vectors of adjacent vascular branches; specifically, it involves calculating the spatial angle between the direction vectors of adjacent branches.

[0102] Once a bifurcation feature is detected, the system calculates a hierarchical increment coefficient based on the bifurcation angle. Generally, the larger the angle of the branching vessel, the greater the difference in diameter and morphology between the new branch vessel and the original main vessel. Therefore, the hierarchical increment coefficient increases accordingly, thereby adjusting the vessel radius. This coefficient is calculated by setting rules based on the physiological characteristics of vessel bifurcation.

[0103] For the main blood vessel section, its initial radius is set to the system resolution multiplied by a preset scaling factor. This scaling factor is obtained through statistical analysis of the correlation between the main blood vessel diameter and the system resolution. In other words, this scaling factor is set based on historical data or the diameter distribution of the main blood vessel during actual scans, thereby ensuring that the radius of the main blood vessel section is set reasonably and accurately during blood vessel extraction.

[0104] As the number of branching levels in a blood vessel increases, its radius gradually decreases. The amount of radius reduction during this process is determined by a combined calculation of the level increase coefficient and the bifurcation angle. Specifically, for every unit increase in the bifurcation angle, the vessel radius decreases according to a non-linear decay rate. This non-linear decay rate is determined by a branch vessel diameter decay model. Based on vascular physiological characteristics, this decay model provides an accurate radius reduction ratio by modeling the diameter decay patterns of different branching vessel types.

[0105] In one possible implementation, the standard vascular frequency band is defined by performing a frequency domain transformation on CT images to find the energy concentration band corresponding to the density distribution of the vascular region. In the frequency domain, a vascular region typically manifests as a specific energy band, the density distribution of which is related to the morphology, thickness, and other characteristics of the blood vessel. To accurately determine this frequency band, the boundaries of the energy concentration band need to be set based on the frequency domain energy distribution of the vascular region.

[0106] Boundary determination is achieved by analyzing the cumulative probability of energy distribution in the vascular region in the frequency domain. When the cumulative probability reaches the first critical value, the low-frequency boundary is defined; when the cumulative probability reaches the second critical value, the high-frequency boundary is defined. This ensures that the frequency band extraction of the vascular region is not affected by noise, and that the features of the vascular region are extracted as accurately as possible, based on the frequency characteristics of the vascular region.

[0107] Calcification features typically manifest as high-density areas, resulting in prominent calcification signals in CT images. To effectively extract calcification feature bands, it is first necessary to label typical calcification regions in the spatial domain image. These regions represent areas with high calcification density, usually consisting of relatively firm tissue or abnormal areas.

[0108] By transforming the frequency domain of the labeled calcified sample region, the peak value of its frequency domain energy is calculated and compared with the average frequency domain energy of the entire image. This ratio reflects the prominence of the calcified region in the frequency domain. If the ratio exceeds a preset multiple, it indicates that the frequency domain energy of that region is relatively prominent and can be used as the core frequency band of the calcification feature band.

[0109] Based on the ratio of the frequency domain energy peaks, when this ratio exceeds a preset multiple, the system will use the frequency corresponding to the frequency domain energy peak as the center of the reference frequency band. This frequency band is the calcification characteristic frequency band, and its frequency domain range includes regions with high calcification density. The extraction of the calcification characteristic frequency band plays an important role in identifying calcified areas in blood vessels, helping to avoid misclassifying calcified areas as normal vascular structures.

[0110] In one possible implementation, grayscale values ​​are statistically analyzed for each local region, and the grayscale standard deviation σ for that region is calculated. The standard deviation reflects the degree of variation in grayscale values ​​within a local region, i.e., the contrast and noise characteristics of the image. The larger the standard deviation of a local region, the higher the contrast and the stronger the noise in that region. Conversely, regions with smaller standard deviations have lower contrast and less noise.

[0111] Based on the calculated local grayscale standard deviation σ, different fusion weight coefficients are selected to adapt to different noise and contrast conditions.

[0112] If the grayscale standard deviation σ of a local region falls within the high-noise range, it indicates that the region has high noise levels, which may affect the accuracy of vessel extraction. In this case, the spatial domain enhancement map (by enhancing the details and edge information of the image) needs to rely more on noise compensation to mitigate the interference of noise on vessel extraction. Therefore, the fusion weight is set as the sum of the baseline weight and the noise compensation weight, ensuring that more noise compensation components are retained in areas with higher noise levels, thereby improving the robustness of the final image.

[0113] If the grayscale standard deviation σ of a local area is small, it indicates that the area has less noise and lower contrast. To ensure accurate extraction of vascular information, the weight of the fundamental frequency angiogram (which typically contains less noise and more accurate vascular boundary information) should be increased proportionally to the signal-to-noise ratio (SNR). Regions with high SNR usually have clearer vascular extraction results, so these regions should be given higher fusion weights to retain richer vascular information.

[0114] First, the standard deviation distribution of the background region in the entire image is statistically analyzed, generating a histogram of the standard deviations. This effectively reveals the noise level in different regions of the entire image. Background regions typically have high noise levels and lack obvious vascular structures; therefore, the distribution of their standard deviations is of significant reference value.

[0115] By analyzing the histogram of the standard deviation distribution of the background region in the entire image, the peak location of the distribution can be determined. The standard deviation value corresponding to this peak can be regarded as the upper limit of the low-noise interval. The low-noise interval refers to areas with less noise and lower contrast, which typically contain relatively clear vascular features.

[0116] A preset multiple of the upper standard deviation of the low-noise interval is used as the lower limit of the high-noise interval. The high-noise interval represents areas with strong noise, which may contain artifacts or other irrelevant image features, and therefore require noise compensation.

[0117] In one possible implementation, to extract the morphological information of the blood vessel more accurately, it is necessary to fit the blood vessel cross-section. By fitting an ellipse to the current cross-section of the blood vessel, the main features of the blood vessel cross-section, including the minor axis and the major axis, can be obtained.

[0118] Specifically, the minor axis diameter (d): The minor axis diameter d is one of the parameters after ellipse fitting. It represents the shortest radius of the blood vessel cross section and reflects the lateral dimension of the blood vessel cross section.

[0119] Major axis length (L): The major axis L is the length of the semi-major axis after ellipse fitting, reflecting the longitudinal dimension of the blood vessel in cross-section. There is a specific relationship between the major axis L and the minor axis diameter d, namely L = k × d, where k is the curvature coupling coefficient. This relationship indicates that the major axis length and the minor axis diameter are related through the curvature coupling coefficient k. The curvature coupling coefficient is inversely proportional to the radius of curvature of the blood vessel, meaning that the greater the curvature of the blood vessel, the more the ratio of the major axis to the minor axis will be adjusted.

[0120] The curvature coupling coefficient k is inversely proportional to the radius of curvature of the vessel's centerline. The radius of curvature reflects the degree of curvature of the vessel's centerline; that is, the greater the curvature of the vessel, the smaller the radius of curvature. To ensure an accurate description of the vascular structure, the ratio of the major axis to the minor axis is adjusted according to the radius of curvature at the curved sections of the vessel, resulting in a more precise vascular geometry.

[0121] Extracting the local centerline segment along the v-direction: When extracting the curvature of a blood vessel, the first step is to extract the local centerline segment along the v-direction. This step is to determine the shape of the vessel's centerline, as the curvature of the vessel is closely related to the shape of its centerline; therefore, accurate extraction of the centerline segment is crucial.

[0122] The extracted centerline segment is fitted as a cubic spline curve. Spline curve fitting can smooth the centerline of the blood vessel and reduce the impact of noise on curvature calculation. In this way, we can obtain more accurate curvature information.

[0123] Based on spline curve fitting, the extreme points of curvature are further calculated. These extreme points represent the points where the degree of vascular tortuosity changes, helping us to identify areas where vascular tortuosity is more pronounced.

[0124] The distance between the extreme points of curvature is used as an estimate of the radius of curvature. The radius of curvature is an important parameter describing the degree of tortuosity of blood vessels; a smaller radius of curvature indicates greater tortuosity, while a larger radius of curvature indicates a straighter blood vessel.

[0125] The minor axis length is fixed as the product of the CT image slice thickness and the pixel size. Slice thickness is typically a physical parameter of the image, while pixel size is a measure of image resolution. The minor axis length is determined to ensure the consistency and accuracy of blood vessels across different image slices.

[0126] The direction of the minor axis is determined by the secondary direction vector of the blood vessel cross-section. The secondary direction vector is perpendicular to the plane of the blood vessel cross-section and is usually related to the physiological orientation of the blood vessel. Therefore, the shape and orientation of the blood vessel cross-section directly affect the definition of the minor axis.

[0127] In one possible implementation, the first step is to extract the edge features of the blood vessels at the calcified edges of the current CT slice. Calcification typically occurs in the hardened areas of blood vessels, and analyzing the calcified edges can help extract the shape information of the vessels. By calculating the convex hull of these calcified edges, the feature point set Pt of the blood vessel edges can be determined. The convex hull is the smallest convex polygon formed around the blood vessel contour based on these feature points. Extracting these feature points provides an accurate starting reference for subsequent contour tracking.

[0128] For adjacent slices, a search plane needs to be established along the v-direction, using this direction as the normal vector for the search. The v-direction refers to the direction perpendicular to the current slice plane, typically corresponding to the longitudinal extension direction of the blood vessel. By establishing such a search plane on adjacent slices, the matching process can be effectively guided, making contour tracking from one slice to the next more accurate.

[0129] The set of feature points Pt at the edge of the blood vessel is projected onto the search plane to form the target region. Then, within the target region, contour points of adjacent slices are matched. This process mainly involves comparing the contour features of two slices to find the optimal matching contour points. In this way, the contour of the blood vessel can be accurately traced, thereby establishing the continuity of the blood vessel morphology between different slices.

[0130] Once the contour point sets of adjacent slices are matched, these matching points can be connected to generate a preliminary spatial curve. This spatial curve represents the continuous morphology of blood vessels across different slices and forms the basis for subsequent path repair.

[0131] The generated initial spatial curve needs further analysis to determine its smoothness. By calculating the rate of change of curvature in each segment of the curve, abrupt changes or unnatural bends in the vascular morphology can be detected. When the rate of change of curvature exceeds a preset anatomically reasonable threshold, it indicates that the vascular morphology in that part does not conform to the natural physiological structure and needs to be repaired. At this point, control points can be inserted into the curve to ensure a smooth transition within that region. The insertion of control points helps to correct the vascular morphology, making it conform to reasonable anatomical requirements.

[0132] Grayscale assignment is a crucial step in the restoration process. To restore the vascular morphology of the path region, a neighborhood propagation model is used for grayscale assignment. The restoration path is used as the center line, extending a certain width to both sides for restoration. This width is adaptively adjusted based on the rate of change of the vessel diameter gradient. That is, when the vessel diameter changes significantly, the width of the restoration path adaptively expands or shrinks to ensure that the restoration process matches the actual structure of the vessel. In this way, the complete morphology of the vessel can be restored more accurately.

[0133] In one possible implementation, at the bifurcation point of the blood vessel, it is necessary to analyze the direction of the branching vessels. First, the direction vectors of the parent vessel and each branch are extracted. The parent vessel vector refers to the direction of extension of the parent vessel from the bifurcation point, while the branch vectors are along the direction of the branch vessels. By calculating the angle between the direction vector of each branch and the direction vector of the parent vessel, two angles are obtained: θ1 and θ2, representing the angles between the parent vessel and branch 1 and branch 2, respectively.

[0134] If the angles θ1 and θ2 of the branch do not conform to normal physiological patterns with the angle of the parent vessel—that is, if the included angle θ1 or θ2 is less than a preset angle threshold, or the difference between the two angles |θ1-θ2| is greater than a preset asymmetry threshold—then the branch is considered abnormal. In this case, the abnormal branch angle may be due to image noise, abnormal changes in vessel morphology, or a vessel bifurcation structure that does not conform to physiological norms. This detection process is a crucial step in ensuring the accuracy and reliability of the vessel extraction results.

[0135] At the bifurcation point of a blood vessel, the cross-sectional area of ​​the parent vessel and its branches must first be measured. This process can estimate the cross-sectional area using the diameter of the vessel in CT images or curve fitting methods. The cross-sectional area Am of the parent vessel refers to the cross-sectional area of ​​the parent vessel at the bifurcation point, and the branch cross-sectional areas A1 and A2 are the cross-sectional areas of branch 1 and branch 2, respectively. Measuring the cross-sectional area is fundamental to understanding hemodynamic characteristics and performing subsequent optimizations.

[0136] The ratio (η) of the total cross-sectional area of ​​each branch to the cross-sectional area of ​​the parent vessel is calculated. This ratio reflects the blood flow load of the vascular branches. If the ratio exceeds a reasonable physiological range, it may lead to hemodynamic imbalance, further affecting normal blood flow. Therefore, calculating this ratio is crucial for assessing the blood flow load of vascular branches.

[0137] If the calculated cross-sectional area ratio η exceeds the preset physiologically reasonable range, it indicates an abnormality in the bifurcation structure of the blood vessel, requiring adjustment according to the law of blood flow conservation. The law of blood flow conservation requires maintaining the flow rate within the vascular system; therefore, adjusting the branch diameters allows for a more rational distribution of blood flow. Specifically, an optimization equation for branch diameter is established, aiming to minimize the variance between the cross-sectional area ratio and the physiological target value. The constraints of the optimization equation include:

[0138] Branch angle variation range: Ensure that the optimized branch angle is within a reasonable range to avoid excessively sharp bends.

[0139] Curvature continuity: Ensure the continuity and smoothness of the vascular path, and avoid unnatural breaks or sharp turns.

[0140] Diameter decay rate: The change in blood vessel diameter with distance should conform to the normal decay trend, avoiding unnatural vasodilation or vasoconstriction.

[0141] After optimizing the branch diameter, the optimized diameter is used to reconstruct the bifurcation topology. By adjusting the branch diameter and angle, a vascular branching structure that conforms to physiological requirements is obtained. This process helps restore the normal hemodynamic state of the blood vessels, making the vascular system more in line with physiological laws and improving the accuracy and reliability of the automatic vascular extraction results.

[0142] By accurately calculating the angle between the branch direction vector and the parent vessel vector, abnormal vascular branch angles can be effectively detected. This process eliminates errors caused by scanning noise or variations in vascular morphology, ensuring that the branching structure conforms to actual anatomical patterns and reducing the risk of misdiagnosis.

[0143] By calculating the cross-sectional area ratio of vascular branches and optimizing the branch diameters based on the law of blood flow conservation, a more rational blood flow distribution within the vascular system can be ensured. The optimized vessel diameter can better simulate actual blood flow characteristics, thereby improving the match between vessel morphology and hemodynamics. This is of great significance for vascular modeling in medical image processing and analysis.

[0144] By accurately modeling and optimizing hemodynamics, the topology of vascular bifurcation points can be reconstructed, ensuring the continuity and stability of the vascular system in its branching structure. This optimization process effectively avoids problems such as blood flow imbalance and irregular morphology in the vascular system, improving the robustness of the automatic vascular extraction algorithm.

[0145] By combining hemodynamic constraints and anatomical principles for branch topology correction, the actual physiological structure can be better simulated during vessel extraction. This results in vessel extraction that not only conforms to anatomical standards but also better meets actual blood flow requirements. This has significant clinical application value for medical imaging diagnosis and surgical planning.

[0146] In one possible implementation, to obtain statistical data on vascular bifurcation angles, it is first necessary to scan blood vessels at different anatomical locations and collect bifurcation angle data at each location. This bifurcation angle data comes from a large number of healthy blood vessel samples, covering bifurcation points in different anatomical locations, such as the brain, heart, and limbs. Using high-precision CT imaging technology, the angle of each bifurcation point is accurately measured, forming a large sample database to provide data support for subsequent angle statistical analysis.

[0147] The collected angle samples were categorized according to vessel type (e.g., large arteries, microvessels, etc.). Then, statistical analysis was performed on the bifurcation angle data within each vessel type group, calculating the variance of the angle distribution for each group. The variance reflects the fluctuation of the bifurcation angles in each vessel group, thus providing a mathematical basis for setting subsequent angle thresholds. Vessel types with larger variances may indicate more dispersed or variable angle distributions, while those with smaller variances tend to have more concentrated angles.

[0148] After calculating the mean and variance of the angles for each blood vessel type, a lower limit for the angle threshold is set. Specifically, for each blood vessel type, the mean of the bifurcation angles for that group is calculated, and three times the variance is subtracted from this value to obtain the lower limit for the angle threshold of that group. This lower limit represents the value that, under normal healthy conditions, the bifurcation angle will not be less than, and is thus used for subsequent detection and correction of abnormal angles. The setting of three times the variance ensures a relatively strict screening condition, effectively removing abnormal blood vessel bifurcation angles.

[0149] In blood flow simulation, the dynamic characteristics of blood flow in blood vessels are numerically simulated to analyze the relationship between the cross-sectional area ratio and blood pressure loss. As the cross-sectional area ratio of blood vessel branches changes, the blood pressure loss also changes; an excessively large cross-sectional area ratio may lead to blood flow obstruction or even hemodynamic imbalance. Through simulation experiments, a mathematical model of the relationship between the cross-sectional area ratio and pressure loss is constructed to provide a reasonable reference.

[0150] Continuing with blood flow simulation experiments, when the change in the cross-sectional area ratio causes the rate of increase in blood pressure loss to exceed the critical slope, the boundary value of this ratio is recorded. The critical slope refers to the point in hemodynamics where pressure loss begins to increase sharply. This critical point reflects excessively high blood flow resistance in the vascular branch, exceeding the normal range. The recorded boundary value of the cross-sectional area ratio at this point is the pressure loss threshold for this type of blood vessel.

[0151] For blood vessels at different anatomical locations, the critical pressure loss points for each location are analyzed, and the intersection of the ratio boundaries at each location is used as the definition of a globally reasonable range. This means that the globally reasonable range is a comprehensive result considering the hemodynamic characteristics of each anatomical site, ensuring that the cross-sectional area ratio of the blood vessel does not exceed the physiologically permissible range regardless of the anatomical location. This ensures that the blood vessel extraction model has consistent physiological rationality across different anatomical regions.

[0152] In one possible implementation, during automated vessel extraction, a normal vascular region must first be identified and selected as a reference. To ensure the reliability of the selected region, the process begins by backtracking along the v-direction of the vessel (i.e., the longitudinal direction) until the nearest uncalcified vascular section is found. Calcified vessels typically cause distortion or interference in the imaging signal; therefore, sections containing calcified areas must be avoided to ensure the healthy vascular characteristics of the reference region. This step helps ensure that the selected reference region is a normal vessel that conforms to physiological characteristics, avoiding misfitting due to calcification or diseased areas.

[0153] Once an uncalcified section is located, all pixels that are topologically connected to the area requiring repair will be extracted from that section. Topologically connected pixels are those that have a spatial or adjacency relationship with the repair area; the connectivity of these pixels helps preserve the morphology and structure of the original blood vessel during the repair process. Pixels extracted in this way can serve as grayscale reference values ​​for normal blood vessels, providing high-quality input data for subsequent repair.

[0154] During the grayscale assignment process of the restoration path, to ensure a reasonable distribution of pixel values ​​in the restoration area, three-dimensional Euclidean distance is used to evaluate the spatial relationship between the reference area and the restoration area. Based on this relationship, the pixel values ​​in the restoration area are weighted inversely proportional to their spatial distance from the reference area. Pixels that are closer are assigned higher weights, while pixels that are farther away are assigned lower weights. This weighting calculation helps ensure that the grayscale values ​​of the restoration area are closer to the actual situation of the reference area, avoiding excessive influence from pixels far from the reference area on the restoration result.

[0155] The structural weights are related to the gradient intensity of the vessel wall, which reflects the degree of change at the vessel edge, i.e., the detailed information of the vessel morphology. Regions with high gradient intensity (usually the vessel wall edges) should be assigned higher structural weights to ensure these details are better preserved during the restoration process. Regions with low gradient intensity (such as the central part of the vessel) are assigned lower structural weights. In this way, the restoration algorithm can take into account the morphological features of the vessel during reconstruction, especially the details of the vessel wall, ensuring that the restoration result is more realistic and conforms to physiological structure.

[0156] The dynamic attenuation factor is set to account for the distance variation between the repair path and the reference region. As the repair path moves further away from the reference section, the attenuation factor decays exponentially, thus gradually reducing the influence of the repaired area. The attenuation rate is calibrated using a vascular grayscale propagation model, which optimizes the attenuation factor setting based on the vascular grayscale distribution and hemodynamic information. This attenuation strategy helps ensure that the repair process does not overly rely on data far from the reference region, thereby improving the accuracy and physiological rationality of the repair.

[0157] Accordingly, the present invention also provides an automated blood vessel extraction system based on CT images, including a memory configured to store instructions, a processor configured to retrieve the instructions from the memory, and capable of implementing the automated blood vessel extraction method based on CT images according to any embodiment of the present invention when executing the instructions.

[0158] This invention encompasses any substitutions, modifications, equivalent methods, and solutions made within the spirit and scope of this invention. To provide the public with a thorough understanding of this invention, specific details are described in detail in the following preferred embodiments; however, those skilled in the art will fully understand the invention even without these details. Furthermore, to avoid unnecessary misunderstanding of the essence of this invention, well-known methods, processes, procedures, components, and circuits are not described in detail.

[0159] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements and modifications can be made without departing from the principle of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. An automatic blood vessel extraction method based on CT images, characterized in that, include: Step 1: Construction of the local coordinate system of the blood vessel and separation of the two domains: Input CT angiography data and determine the principal direction vector v of the blood vessel through Hessian matrix feature analysis; Construct a local cylindrical coordinate system with the principal direction vector v of the blood vessel as the axis, and simultaneously execute the following steps: a: Spatial domain vessel enhancement filtering: A Frangi filter with dynamically adjusted radius is used, whose radius adaptively decreases with the depth of the vessel layer to generate a spatial domain enhancement map; b: Frequency domain calcification feature separation: The standard vascular frequency band and the calcification feature frequency band are separated by fast Fourier transform to generate a fundamental frequency vascular map, and a calcification region mask M is generated based on the calcification feature frequency band. cal ; Step 2: Dynamic dual-domain fusion, which involves adaptively weighting and fusing the spatial domain enhanced map with the fundamental frequency angiography. Step 3: Orientation-constrained morphological reconstruction: Construct a structural element whose major axis is parallel to the direction of the blood vessel based on the main direction vector v of the blood vessel. Perform directional opening and directional closing operations sequentially along the principal direction vector v of the blood vessel; Step 4: Anatomy-driven vascular repair and correction: Based on the generated calcification region mask M cal The region is marked, and the contours of adjacent slices are traced along the main vascular direction vector v and the repair path is fitted. At the bifurcation point, the topology is corrected according to the branch angle and hemodynamic constraints.

2. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, The process of dynamically adjusting the radius described in step 1 includes: Identifying vascular hierarchy: When a bifurcation feature is detected, a hierarchy increment coefficient is calculated based on the bifurcation angle, which is determined by the spatial angle between the direction vectors of adjacent branches. The initial radius of the main trunk layer is set to the system resolution multiplied by a preset scaling factor, which is obtained by statistically analyzing the correlation between the diameter of the main trunk vessel and the resolution. The reduction in the branch layer radius is calculated based on the hierarchical increment coefficient and the bifurcation angle: for each unit increase in the bifurcation angle, the radius is reduced according to the nonlinear decay rate corresponding to the hierarchical depth, and the nonlinear decay rate is determined by the branch vessel diameter decay model.

3. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, Step 1, which involves separating the standard vascular frequency band from the calcification characteristic frequency band, includes: Method for determining standard vascular frequency bands: In the frequency domain space, locate the energy concentration band corresponding to the preset vascular density range. Its boundary is determined by analyzing the cumulative probability of the frequency domain energy distribution of the vascular region. When the cumulative probability reaches the first critical value, the low-frequency boundary is defined, and when it reaches the second critical value, the high-frequency boundary is defined. Method for separating calcification characteristic frequency bands: Extracting frequency domain components above a calcification density threshold, wherein the calcification density threshold is dynamically set through the following process: a: Mark typical calcified sample areas in the airspace; b: Calculate the ratio of the peak frequency domain energy of a typical calcified sample region to the average energy of the entire image; c: When the ratio exceeds a preset multiple, the frequency corresponding to the peak energy in the frequency domain is used as the center of the reference frequency band.

4. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, The contrast adaptive weighted fusion in step 2 includes: Calculate the local grayscale standard deviation σ, and select the fusion weighting coefficient based on the interval where σ is located: When σ is in the high noise region, the spatial enhancement map weight is set as the sum of the baseline weight and the noise compensation weight; When σ is in the low noise range, the weight of the fundamental frequency angiography is increased proportionally to the signal-to-noise ratio. The boundaries of the high-noise region and the low-noise region are determined in the following way: a: Histogram of the standard deviation distribution of the background area of ​​the entire image; b: Use the standard deviation corresponding to the peak value of the distribution as the upper limit of the low noise interval; c: The standard deviation exceeding the upper limit of the low noise range by a preset multiple is taken as the lower limit of the high noise range.

5. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, The construction of structural elements in step 3 includes: Major axis length setting: The minor axis diameter d is obtained by fitting the current vessel cross-section ellipse. The length L and d satisfy L = k × d, where k is the curvature coupling coefficient. The curvature coupling coefficient is adjusted inversely proportionally to the radius of curvature of the vessel centerline. The radius of curvature is calculated through the following process: a: Extract the local centerline segment along the v direction; b: Fit a cubic spline curve and calculate the curvature extrema; c: Use the distance between extreme points as the estimated value of the radius of curvature; The minor axis length is fixed as the product of the layer thickness and the pixel size, and its direction is determined by the secondary direction vector of the blood vessel cross section.

6. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, Step 4, which involves tracking the contours of adjacent slices and fitting the repair path, specifically includes: Methods for tracing the contours of adjacent slices: a: Extract the convex hull feature point set P from the calcification edge of the current slice. t ; b: Establish a search plane on adjacent slices with the projection direction of v as the normal vector; c: P t Project the target region onto the search plane, and match the next slice contour points within the target region; The fitting and repair path rules include: a: Connect the matching point set to generate the initial spatial curve; b: Calculate the rate of change of curvature of each segment of the curve. When the rate of change exceeds the reasonable anatomical threshold, insert control points for smoothing. c: Gray-scale assignment adopts a neighborhood propagation model, with the repair path as the center line, extending to both sides with a preset width. The preset width is adaptively adjusted according to the gradient change rate of the blood vessel diameter.

7. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, Step 4, correcting the topology at the bifurcation point based on the branch angle and hemodynamic constraints, includes: Abnormal branch angle detection: a: Calculate the angles θ1 and θ2 between the direction vector of each branch and the parent blood vessel vector; b: If min(θ1,θ2) < angle threshold or |θ1-θ2| > asymmetric threshold, it is judged as abnormal; Hemodynamic reconstitution process: a: Measure the cross-sectional area A of the mother blood vessel m and the branch cross-sectional areas A1, A2; b: Calculate the cross-sectional area ratio η = (A1 + A2) / A m ; c: When η exceeds the physiologically reasonable range, the branch diameter is redistributed according to the law of blood flow conservation, specifically including: An optimization equation for branch diameter is established, with the objective function being to minimize the variance between the cross-sectional area ratio and the physiological target value. The constraints include the range of branch angle variation, curvature continuity, and diameter attenuation rate. d: Reconstruct the bifurcation topology using the optimized diameter.

8. The automatic blood vessel extraction method based on CT images according to claim 7, characterized in that: The angle threshold is obtained by statistically analyzing a database of healthy blood vessel bifurcations. The specific process includes: a: Collect bifurcation angle samples from different anatomical locations; b: Calculate the variance of the angular distribution by grouping by blood vessel type; c: Take the lower limit of the threshold as the value of the mean of each group minus three times the variance; The method for setting the reasonable physiological range: a: Analyze the relationship between the cross-sectional area ratio and pressure loss in blood flow simulation experiments; b: When the rate of increase of pressure loss exceeds the critical slope, record the ratio boundary at this point; c: Take the intersection of the boundaries of each anatomical location as the reasonable range for the whole.

9. The automatic blood vessel extraction method based on CT images according to claim 1, characterized in that, The fitting and repair path in step 4 also includes grayscale assignment, specifically including: Selection of reference areas for normal blood vessels: a: Tracing back along the v direction to the nearest uncalcified section; b: Extract all pixels that are topologically connected to the repaired region within this section; Weighted calculation rules: a: Spatial weights: Distributed inversely proportional to the three-dimensional Euclidean distance; b: Structural weights: adjusted positively based on the gradient intensity of the blood vessel wall; c: Dynamic decay factor: It decays exponentially as the distance between the repair path and the reference section increases, and the decay rate is calibrated by the vascular grayscale propagation model.

10. An automated blood vessel extraction system based on CT images, characterized in that, The system includes a memory configured to store instructions, a processor configured to retrieve the instructions from the memory, and, when executing the instructions, to implement the automatic blood vessel extraction method based on CT images as described in any one of claims 1-9.

Images