Digital image-based coronary artery bypass graft planning method

Abstract

The application provides a coronary artery bypass grafting planning method based on digital images, which comprises preprocessing, three-dimensional segmentation, construction of a three-dimensional coronary artery anatomical model, multi-scale decomposition, extraction of a blood flow dynamics parameter vector, optimization estimation and fusion, functional evaluation of a stenosis lesion, determination of a set of preferred distal anastomosis candidate positions, planning of a virtual bridge blood vessel running path, post-implantation hemodynamic reevaluation of the virtual bridge blood vessel, and selection of an optimal bypass grafting surgery planning scheme.The application realizes multi-scale and refined characterization of the hemodynamic state of the coronary artery.

Description

Technical Field

[0001] This invention relates to the field of medical image processing technology, and in particular to a method for coronary artery bypass grafting planning based on digital images. Background Technology

[0002] Coronary artery bypass grafting (CABG) is the gold standard surgical procedure for treating multivessel coronary artery disease and left main coronary artery disease. This surgery uses a graft (such as the left internal mammary artery, great saphenous vein, or radial artery) to bypass the narrowed or occluded segment of the coronary artery, restoring blood flow to the distal myocardium. The success of the surgery largely depends on two key decisions: the location of the distal anastomosis, i.e., where the anastomosis point between the graft and the target vessel should be chosen to maximize the restoration of blood flow to the area distal to the stenosis; and the optimization of the graft's course, i.e., whether the spatial path of the graft from its origin to the anastomosis is smooth and free from the risk of torsion or excessive angulation.

[0003] Currently, preoperative planning for coronary artery bypass grafting (CABG) surgery primarily relies on the surgeon's visual interpretation of coronary computed tomography (CT) angiography images and intraoperative exploration. The surgeon observes two-dimensional or three-dimensional reconstructed coronary artery images, combines this with personal experience to determine the location and severity of the stenosis, and selects the location of the distal anastomosis based on experience. This planning method has the following inherent limitations.

[0004] Visual interpretation can only provide anatomical information about the stenosis, including its location, length, and percentage reduction in diameter, but it cannot provide information on the functional impact of the stenosis on distal blood flow perfusion. Fractional flow reserve (FFR) has been proven by numerous clinical studies to be the gold standard for assessing the functional significance of stenosis; however, current preoperative planning methods lack systematic calculation and application of FFR. Differences in surgeon experience lead to poor consistency in planning schemes, and there is a lack of repeatable quantitative decision-making criteria. The selection of the bypass graft path lacks three-dimensional spatial quantitative optimization, potentially leading to risks such as excessively long bypass grafts, excessively large bypass graft angles, or collisions between bypass grafts and surrounding tissues.

[0005] In recent years, computational fluid dynamics (CFD)-based coronary artery blood flow simulation technology has emerged. This technology simulates blood flow within the coronary arteries on a computer and calculates hemodynamic parameters such as the fractional flow reserve (FVR). For example, existing literature has disclosed a rapid hemodynamic planning method for personalized coronary artery bypass grafting (CABG) by calculating the FVR through the establishment of a personalized coronary artery lumped parameter model. However, existing CFD simulation methods still have the following shortcomings. First, the computational load is enormous; a single simulation may take several hours or even days, making it difficult to routinely apply in clinical practice. Second, existing methods are usually based on single anatomical image data, which is mostly derived from coronary computed tomography (CT) angiography images. The lack of multimodal image fusion leads to the loss of key information such as vessel wall thickness and plaque composition, affecting simulation accuracy. Third, existing methods lack a systematic multi-objective optimization mechanism for the selection of distal anastomosis locations and the planning of bypass graft pathways. There is a lack of data feedback and iterative optimization between different steps, and the generation of the planning scheme is a one-time, unidirectional process. Fourth, the accuracy of two-dimensional or three-dimensional segmentation of vascular trees in existing methods is limited by the slender and tortuous shape of coronary arteries, resulting in insufficient boundary integrity and trajectory continuity of the segmentation results. Summary of the Invention

[0006] This invention provides a coronary artery bypass graft planning method based on digital imaging to solve one or more of the problems mentioned above.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A coronary artery bypass graft (CABG) planning method based on digital imaging includes: acquiring multimodal coronary digital image data and preprocessing it to obtain preprocessed image data; performing three-dimensional segmentation of the coronary artery vascular tree on the preprocessed image data to obtain a three-dimensional segmentation mask; constructing a three-dimensional anatomical model of the coronary arteries based on the three-dimensional segmentation mask and the preprocessed image data, and calculating the blood flow velocity field distribution data and pressure field distribution data within the coronary artery vascular tree based on the three-dimensional anatomical model of the coronary arteries; performing multi-scale decomposition of the blood flow velocity field distribution data and pressure field distribution data using a symplectic geometric mode decomposition algorithm to obtain the symplectic geometric mode components of the velocity field and the pressure field; extracting hemodynamic parameter vectors for each segment of the coronary artery vascular tree based on the symplectic geometric mode components of the velocity field and the pressure field; optimizing and estimating the hemodynamic parameter vectors using a Gaussian process regression algorithm, and then applying the hemodynamic parameter vectors to the Gaussian process regression algorithm. The optimized estimates obtained from regression are fused to obtain a fused hemodynamic parameter vector. Based on this fused hemodynamic parameter vector, a functional assessment of the stenosis lesion is performed to obtain a set of functional assessment results. A set of candidate distal anastomosis locations is generated based on this set of results, and a preferred set of candidate distal anastomosis locations is determined from this set. A virtual bypass vessel path is planned based on the preferred set of candidate distal anastomosis locations. A hemodynamic reassessment is performed after virtual bypass vessel implantation based on this path, obtaining updated hemodynamic parameter vectors for each vessel segment and the virtual bypass vessel's own hemodynamic parameter vectors. Based on these updated and virtual bypass vessel hemodynamic parameter vectors, each virtual bypass vessel planning scheme is comprehensively scored, and the virtual bypass vessel planning scheme with the highest comprehensive score is selected as the optimal bypass surgery planning scheme.

[0008] In summary, the present invention has at least the following beneficial effects: This invention employs a symplectic geometric mode decomposition algorithm to perform multi-scale decomposition of the coronary artery blood flow velocity and pressure fields. This decomposes complex hemodynamic information into modal components with clear physical meaning, achieving a multi-scale, refined characterization of the coronary artery hemodynamic state. The decomposition results are not only used for subsequent hemodynamic parameter extraction but also interact with the Gaussian process regression algorithm through velocity-pressure mode pairs. This provides Gaussian process regression with hemodynamically rich prior features, enabling the model to achieve high estimation accuracy with a small number of training samples and effectively reducing the computational burden of computational fluid dynamics simulations.

[0009] This invention constructs a vascular network kernel function, consisting of a weighted sum of a squared exponential kernel function, a topological kernel function, and a hemodynamic kernel function. This enables Gaussian process regression models to simultaneously utilize anatomical similarity, topological proximity, and hemodynamic similarity between vascular segments for parameter estimation, significantly improving the estimation accuracy and computational efficiency of hemodynamic parameters. The introduction of the topological and hemodynamic kernel functions in the vascular network kernel function overcomes the limitation of traditional Gaussian process regression kernel functions that are based solely on the spatial distance of input features, achieving the encoding of the coronary artery tree topology and hemodynamic physical laws.

[0010] This invention establishes a closed-loop iterative optimization framework, feeding back the hemodynamic reassessment results after virtual bridging vessel implantation to the distal anastomosis candidate site generation stage. This achieves adaptive optimization of the planning scheme, overcoming the deficiency of existing methods where the unidirectional process cannot be adjusted based on prediction results. The deviation information feedback mechanism and weight coefficient adaptive adjustment mechanism in the closed-loop iterative optimization framework ensure that the optimization process has a clear technical direction and convergence guarantee.

[0011] This invention improves the boundary integrity and trajectory continuity of coronary artery tree segmentation by incorporating a boundary integrity enhancement module and a trajectory continuity enhancement module into a dynamic graph convolutional network, and employing a fuzzy C-means clustering algorithm to calculate fuzzy membership degrees. This provides a high-quality anatomical model foundation for subsequent hemodynamic simulations. The fuzzy membership degrees in the boundary integrity enhancement module can quantitatively characterize the uncertainty of voxel points belonging to the vessel boundary, thereby applying an adaptive penalty weight to the prediction error of the boundary region in the loss function, improving the segmentation network's sensitivity to vessel boundary detection. Attached Figure Description

[0012] Figure 1 This is a schematic diagram of a coronary artery bypass graft planning method based on digital imaging.

[0013] Figure 2 This is a flowchart of the segmentation process for a dynamic graph convolutional network.

[0014] Figure 3 Flowchart for generating candidate locations and planning paths for multi-target remote anastomosis joints.

[0015] Figure 4 A flowchart for closed-loop iterative optimization and verification. Detailed Implementation

[0016] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0017] like Figure 1 As shown, this embodiment provides a coronary artery bypass graft planning method based on digital images, including the following steps.

[0018] Step S1: Acquisition and Preprocessing of Multimodal Coronary Digital Images. Acquire the patient's multimodal coronary digital image data. This data includes coronary computed tomography (CT) angiography, intravascular ultrasound (IU) imaging, and optical coherence tomography (OCT) imaging. CT angiography data is acquired using a 256-slice or higher spiral CT scanner with a slice thickness of 0.5 mm to 0.625 mm and a reconstruction matrix of 512 × 512, containing coronary artery image sequences covering 70% to 80% of the diastolic RR interval. IU is acquired using an IU catheter within the coronary artery lumen at a retraction speed of 0.5 mm / s to 1.0 mm / s, with an ultrasound frequency of 20 MHz to 45 MHz. The acquisition of intravascular ultrasound (IU) and optical coherence tomography (OCT) images was performed synchronously with the patient's electrocardiogram (ECG) signal. During post-processing, image frames within 70% to 80% of the RR interval were extracted from the synchronously recorded ECG signals for registration with the coronary computed tomography (CT) angiography data. OCT images were acquired via an OCT catheter within the coronary artery lumen at a retraction rate of 20 to 40 mm / s, with a center wavelength of 1300 nm.

[0019] The preprocessing of multimodal coronary digital imaging data includes three sub-steps: image denoising, image registration, and image enhancement.

[0020] Image denoising was performed on each frame of image data using a three-dimensional block matching denoising algorithm based on nonlocal means. This denoising process was applied to each voxel in the coronary computed tomography angiography image data. Its denoised grayscale value Calculated using the following formula:

[0021] ; in, Voxel points in coronary computed tomography angiography (CT) images The original grayscale value, Represented by voxel points The three-dimensional search neighborhood centered on the center voxel point With voxel points The nonlocal weighting coefficients between them. The calculation formula is:

[0022] ; in, Represented by voxel points The image block vector centered on the image. Represented by voxel points The image block vector centered on the image. Represents image block vector With image block vector The square of the Euclidean distance between them These are the filter parameters. The value is set to 1.5 to 2.0 times the estimated noise standard deviation in coronary computed tomography angiography image data. This range achieves an optimal balance between sufficient denoising and structural preservation.

[0023] Image registration processing aligns the intravascular ultrasound (IU) and optical coherence tomography (OCT) images of the coronary arteries to the spatial coordinate system of the coronary computed tomography (CT) angiography images. The image registration process employs a mutual information-based three-dimensional rigid registration algorithm, using the CT angiography images as the fixed image and the IU or OCT images as the floating image, by maximizing the mutual information metric. To optimize spatial transformation parameters :

[0024] ; in, This represents the grayscale value distribution of a fixed image. This represents the grayscale value distribution of the floating image. Represents spatial transformation parameters. Represents information entropy. Represents joint information entropy, This indicates the composition of functions. Indicates the spatial transformation parameters Effect on floating images That is, the transformed image is obtained by rotating and translating the floating image.

[0025] After the 3D rigid registration algorithm based on mutual information is completed, the registration accuracy is evaluated using the target registration error. The method for measuring the target registration error is as follows: [The text abruptly shifts to a different topic] ...on a fixed image... and floating images At least six corresponding anatomical landmarks are manually or automatically labeled in the image. These landmarks include coronary artery bifurcation points, vessel bends, and calcified plaque center points. The average Euclidean distance between all corresponding landmarks is calculated as the target registration error. If an automatic labeling method is used, at least six corresponding anatomical landmarks are automatically identified based on the lumen boundary features and bifurcation point morphology features in the multimodal images. The method extracts the lumen centerline path and bifurcation point coordinates from both fixed and floating images, using the consistency of the bifurcation point topology as the primary matching criterion and the similarity of the lumen cross-sectional shape as a secondary matching criterion, automatically generating corresponding landmark pairs. If fewer than six bifurcation points are automatically identified, additional landmarks are extracted, including the curvature extrema of the lumen centerline and the edge points of calcified plaques in the vessel wall. If the spatial distance deviation between any pair of automatically labeled landmarks exceeds 2.0 mm, manual correction is performed. The registration accuracy is deemed to meet the requirements when the target registration error is less than 1.0 mm; when the target registration error is greater than or equal to 1.0 mm, the spatial transformation parameters are reinitialized. The registration optimization process is repeated until the target registration error value is less than 1.0 mm.

[0026] Image enhancement processing employs an enhancement algorithm based on the multi-scale retina-cortex theory to enhance the registered image data. For each voxel... Enhanced grayscale value Calculated using the following formula:

[0027] ; in, This indicates that the registered image data is located at voxel points. The gray values ​​are normalized to the interval [0,1]. Indicates the total number of scales. Indicates scale index. Indicates the first Weighting coefficients for each scale The standard deviation is expressed as The three-dimensional Gaussian kernel function, This represents the convolution operator. Represents the three-dimensional Gaussian kernel function With normalized gray values At voxel points The convolution result at the point.

[0028] After preprocessing, an enhanced multimodal coronary digital image dataset is obtained. .in, This represents preprocessed coronary computed tomography (CT) angiography image data. This represents preprocessed intravascular ultrasound imaging data of the coronary arteries. This represents the preprocessed coronary optical coherence tomography (OCT) image data.

[0029] Step S2: 3D segmentation of the coronary artery vascular tree based on a dynamic graph convolutional network. A dynamic graph convolutional network is used to process the enhanced multimodal coronary digital image dataset obtained in Step S1. A 3D segmentation of the coronary artery tree is performed to obtain a 3D segmentation mask. The dynamic graph convolutional network comprises two core components: a boundary integrity enhancement module and a trajectory continuity enhancement module.

[0030] The model construction process of dynamic graph convolutional networks is as follows.

[0031] Define the initial graph structure of the coronary artery vascular tree. .in, This represents the set of nodes in the initial graph structure. Indicates the first The spatial coordinate vectors of the initial nodes For node index and , The initial total number of nodes; Represents the set of edges of the initial graph structure. Indicates the first The initial node and the first Edges between the initial nodes For node index and , .

[0032] Initial graph structure The construction method is as follows: using coronary computed tomography (CT) angiography image data. The rough centerline of the coronary arteries serves as the framework, and samples are uniformly taken radially at each center point along this rough centerline. Each node is connected to an adjacent center point, and edges are also connected between adjacent radial nodes at the same center point. The value of is between 8 and 16, which achieves the optimal balance between segmentation accuracy and computational efficiency.

[0033] The training process for a dynamic graph convolutional network is as follows.

[0034] Construct the training dataset for the dynamic graph convolutional network. The training dataset contains... Group training samples, The total number of training samples is [number]. Each training sample set includes input multimodal coronary digital image data and its corresponding ground truth values ​​for 3D segmentation annotation of the coronary artery tree. The ground truth values ​​for 3D segmentation annotation of the coronary artery tree are obtained by taking the majority consensus result after independent annotation by three or more physicians with more than five years of experience in cardiovascular imaging diagnosis.

[0035] Define the loss function for dynamic graph convolutional networks. The composite loss function consists of the segmentation loss term. Boundary integrity loss term Trajectory continuity loss term The weighted sum of the three terms constitutes:

[0036] ; in, These are the weighting coefficients for the segmentation loss term. The weighting coefficients for the boundary integrity loss term. These are the weighting coefficients for the trajectory continuity loss term. The weighting coefficients for the segmentation loss term. The value is 0.6, which is the weighting coefficient of the boundary integrity loss term. The value is 0.2, which is the weighting coefficient of the trajectory continuity loss term. The value is 0.2. The specific values ​​of the above weight coefficients are obtained through hyperparameter optimization on the validation set using a grid search method, with a grid search range of 0.1 to 0.9 and a step size of 0.1. Segmentation loss term. A combination of Dice loss function and cross-entropy loss function is used:

[0037] ; in, The total number of voxels. voxel point index and , The first prediction for a dynamic graph convolutional network The probability value that an individual pixel belongs to a coronary artery. For the first The ground truth value of an individual pixel is used to label it. A ground truth value of 1 indicates that it belongs to a coronary artery, and a value of 0 indicates that it does not belong to a coronary artery. To prevent the denominator from being zero, a smoothing parameter is used. The value is .

[0038] Boundary integrity loss term Modeling based on fuzzy logic. Boundary integrity loss term. The calculation formula is:

[0039] ; in, For the first Fuzzy membership degree of individual points, fuzzy membership degree Based on the fuzzy C-means clustering algorithm, according to the... The grayscale value, gradient magnitude, and neighborhood texture features of individual pixels are calculated; the objective function of the fuzzy C-means clustering algorithm is... ,in In this embodiment, the number of clusters is [number]. The cluster centers correspond to three categories: lumen, vessel wall, and background. In this embodiment, the fuzzy index is used. ; For the first The feature vector of an individual pixel is composed of the gray value, gradient magnitude and neighborhood texture features of that pixel. For the first Cluster centers; For the first Individual points belong to the first The fuzzy membership degree of each cluster satisfies The iteration stopping condition is: ,in This is the preset iteration stopping threshold. The first prediction for a dynamic graph convolutional network The probability value of an individual pixel belonging to the blood vessel boundary; For the first True value of blood vessel boundary annotation for individual pixel points.

[0040] Trajectory continuity loss term Computation is based on dynamic graph convolutional networks. Trajectory continuity loss term. The calculation formula is:

[0041] ; in, This represents the set of edges in the current layer of a dynamic graph convolutional network. Represents the set of edges The total number of sides in the middle. Represents a node With nodes The edge between, This represents the first dynamic graph convolutional network. Layer nodes The hidden layer feature vectors This represents the first dynamic graph convolutional network. Layer nodes The hidden layer feature vectors Represents the hidden layer feature vector With hidden layer feature vectors The square of the Euclidean distance between them.

[0042] The encoder path of the dynamic graph convolutional network contains five convolutional modules. Each module consists of a 3D convolutional layer, a batch normalization layer, and a modified linear unit activation function layer connected in sequence. The first convolutional module's 3D convolutional layer contains 32 kernels of size 3×3×3; the second contains 64 kernels of size 3×3×3; the third contains 128 kernels of size 3×3×3; the fourth contains 256 kernels of size 3×3×3; and the fifth contains 512 kernels of size 3×3×3. Each convolutional module is followed by a 3×3×3 max-pooling layer with a stride of 2 to reduce the spatial resolution of the feature map and expand the receptive field.

[0043] The decoder path of the dynamic graph convolutional network contains five upsampling modules. Each upsampling module is connected to the corresponding layer in the encoder path for feature fusion via skip connections. Each upsampling module contains a 3D transposed convolutional layer, a batch normalization layer, and a modified linear unit activation function layer connected in sequence. The first upsampling module's 3D transposed convolutional layer contains 256 convolutional kernels, each with a size of 3×3×3 and a stride of 2. The second upsampling module's 3D transposed convolutional layer contains 128 convolutional kernels, each with a size of 3×3×3 and a stride of 2. The third upsampling module's 3D transposed convolutional layer contains 64 convolutional kernels, each with a size of 3×3×3 and a stride of 2. The fourth upsampling module's 3D transposed convolutional layer contains 32 convolutional kernels, each with a size of 3×3×3 and a stride of 2. The fifth upsampling module's 3D transposed convolutional layer contains 1 convolutional kernel, each with a size of 1×1×1 and a stride of 1. The output is a 3D segmentation probability map of the coronary artery vascular tree. The input to each upsampling module is the concatenation of the upsampling output of the previous layer and the feature map of the corresponding layer of the encoder path in the channel dimension. Skip connections pass the detailed features of the encoder path to the decoder path to restore the edge refinement of the segmentation result.

[0044] The dynamic graph convolutional network is trained using the backpropagation algorithm and an adaptive moment estimation optimizer. The specific parameters of the adaptive moment estimation optimizer are set as follows: the initial learning rate is... The exponential decay rate estimated by the first moment The exponential decay rate estimated by the second moment is 0.9. The numerical stability constant is 0.999. for The weight decay coefficient is The training process employs a learning rate cosine annealing scheduling strategy, meaning the learning rate is adjusted according to... The pattern is that the learning rate gradually decreases from the initial learning rate to... ,in Indicates the first The learning rate during each iteration The initial learning rate, This represents the total number of training epochs. Training continues until the loss function is reached. Converging to below a preset threshold. The preset threshold is set to... .

[0045] The application process of the dynamic graph convolutional network model is as follows.

[0046] The enhanced multimodal coronary digital imaging data set obtained in step S1 Input the trained dynamic graph convolutional network. The dynamic graph convolutional network first extracts multi-scale features from multimodal coronary digital image data through the encoder path, obtaining a set of multi-scale feature maps. .in, In this embodiment, the number of layers in the encoder path is... =5, Indicates the first Feature map of the layer For layer index and .

[0047] like Figure 2 As shown, the boundary integrity enhancement module receives a set of multi-scale feature maps. As input, the probability value of each voxel point belonging to the blood vessel boundary is calculated through fuzzy logic reasoning. and the probability value Feedback is sent to the corresponding layer in the decoder path.

[0048] The trajectory continuity enhancement module is based on a dynamic graph convolutional network for the initial graph structure. Perform iterative updates. In the... In the next iteration, node features The update formula is: ;in, Index of iteration count and , This represents the total number of iterations. Indicates the first Nodes in the next iteration eigenvectors, Indicates the first Nodes in the next iteration eigenvectors, Represents a node The set of neighboring nodes, Represents a node The degree, Represents a node The degree, Indicates the first The learnable weight matrix for the next iteration. This represents a non-linear activation function. (Non-linear activation function) A modified linear unit activation function is used. Total number of iterations. The value of is between 3 and 6, which achieves the optimal balance between feature aggregation sufficiency and branch discrimination preservation.

[0049] The trajectory continuity enhancement module updates the node features in each iteration. Feedback to the decoder path of the dynamic graph convolutional network.

[0050] The decoder path of the dynamic graph convolutional network integrates feedback information from the boundary integrity enhancement module and the trajectory continuity enhancement module, upsamples layer by layer, and fuses multi-scale features to output a 3D segmentation probability map of the coronary artery tree. 3D segmentation probability map Each voxel point takes a probability value between 0 and 1.

[0051] 3D segmentation probability map Thresholding is performed to obtain a 3D segmentation mask for the coronary artery tree. :

[0052] ; in, This represents the spatial coordinate vector of a voxel point. This represents the probability map of 3D segmentation at voxel points. The probability value at that location. The segmentation threshold is the threshold value. The value is 0.5.

[0053] Step S3: Multi-scale decomposition of coronary blood flow field based on symplectic geometric mode decomposition.

[0054] Three-dimensional segmentation mask based on the coronary artery tree obtained in step S2 Combined with the enhanced multimodal coronary digital image data set obtained in step S1 We constructed a patient-specific three-dimensional anatomical model of the coronary arteries and calculated the blood flow velocity field distribution data within the coronary artery vascular tree based on the three-dimensional anatomical model of the coronary arteries. and pressure field distribution data .

[0055] The method for constructing a three-dimensional anatomical model of the coronary arteries is as follows: using a three-dimensional segmentation mask. Using the moving cubes algorithm as input, a 3D surface mesh model of the coronary artery vascular tree was extracted. The surface mesh model underwent quality evaluation and optimization. Quality evaluation indicators included: aspect ratio less than 3.0, minimum interior angle greater than 20 degrees, and normal direction pointing outwards from the vessel. Triangular faces that did not meet these quality indicators were locally re-meshed, and faces with inconsistent normal directions underwent normal flipping correction. After quality evaluation and optimization, the total number of triangular faces in the surface mesh model was controlled between 100,000 and 500,000. Coronary intravascular ultrasound imaging data was used. Coronary optical coherence tomography (OCT) images For reference, the vessel wall thickness and lumen diameter of the three-dimensional surface mesh model were locally refined.

[0056] Preferably, before starting the computational fluid dynamics simulation, mesh independence verification is performed: three different density volume meshes (coarse mesh, medium mesh, and fine mesh, with the number of units in the coarse mesh being 0.5 times that of the medium mesh, and the number of units in the fine mesh being 2.0 times that of the medium mesh) are used to simulate the same three-dimensional anatomical model of the coronary artery. The blood flow reserve fraction of the target vessel segment is calculated under the three meshes. When the difference in blood flow reserve fraction between the coarse mesh and the medium mesh, and between the medium mesh and the fine mesh, is less than 5%, the medium mesh is determined to meet the mesh independence requirement.

[0057] Blood flow velocity field distribution data and pressure field distribution data The calculation method is as follows: using a three-dimensional anatomical model of the coronary artery as the computational domain, the incompressible Navier-Stokes equations are solved using computational fluid dynamics.

[0058] ;

[0059] ; in, This represents blood density, with values ​​ranging from 1050 kg / m³ to 1060 kg / m³. This represents the dynamic viscosity of blood, ranging from 0.0035 Pascal-second to 0.0045 Pascal-second. Represents the gradient operator. Represents the Laplace operator. Represents a time variable.

[0060] The boundary conditions for the computational fluid dynamics simulation are set as follows. Inlet boundary conditions: The left and right coronary artery ostia of the three-dimensional coronary anatomy model are used as inlets, and the patient-individualized physiological flow waveform is used as the inlet boundary condition, based on the coronary computed tomography angiography image data obtained in step S1. The CT time-density curve of the middle aortic root was used to estimate the total coronary blood flow using the patient's body surface area normalization method, and then allocated to each inlet according to the proportion of the lumen diameter of each coronary artery branch. The method for determining "each coronary artery branch" is based on the three-dimensional segmentation mask of the coronary vascular tree obtained in step S2. A skeletonization algorithm was used to extract the centerline of the coronary artery tree. Based on the bifurcation points of the centerline, the coronary artery tree was divided into branches of various levels. The main branches directly connecting the ostia of the left and right coronary arteries were identified as the left anterior descending artery, the left circumflex artery, and the right coronary artery, respectively. The diameter of each branch lumen was determined. The measurement method is as follows: along the centerline of the branch vessel, at positions 2.0 mm to 5.0 mm from the branch origin, the arithmetic mean of the lumen diameter on a cross-section perpendicular to the centerline is taken as the lumen diameter of that branch. Blood flow at each inlet... According to the diameter of the branch lumen where the inlet is located The distribution is based on the proportion of the square of the diameter to the sum of the squares of the diameters of all inlet branch lumens:

[0061] ; in The total number of entrances. For entry index and , For entry index and The normalization method for patient body surface area was the Dubose formula. ,in Indicates the patient's height (in centimeters). This indicates the patient's weight (in kilograms). The unit is square meters, total coronary blood flow. , The value ranges from 3.0 liters to 4.0 liters per minute per square meter of body surface area.

[0062] Export boundary conditions: Using all distal exits of the coronary artery 3D anatomical model as exits, coupled resistance boundary conditions are adopted. ,in Indicates the outlet boundary resistance. This indicates the average pressure at the outlet. This represents the average flow rate at the outlet. The empirical scaling factor is set to 0.5 to 1.0. Vessel wall boundary conditions: The vessel wall of the three-dimensional anatomical model of the coronary artery adopts a no-slip boundary condition, meaning the blood flow velocity at the vessel wall surface is zero, and the vessel wall is set as a rigid wall.

[0063] Time step of computational fluid dynamics simulation The setting method is as follows ,in This represents the minimum side length of all grid cells in the volume mesh. Indicates the maximum expected blood flow velocity in the simulation (via Estimate, Peak aortic flow velocity). The initial value of the Courant number is 0.3. The upper limit of the time step is set to the duration of the complete cardiac cycle. One-thousandth of the calculated time step. Less than At seconds, Set as Second.

[0064] The symplectic geometric mode decomposition algorithm was used to analyze the blood flow velocity field distribution data. Perform multi-scale decomposition. Define blood flow velocity field distribution data. The phase space in which it is located is an symplectic manifold ,in It is a smooth manifold. It is in the symptotic form. In the symptotic manifold Define Hamiltonian function :

[0065] ; in, The computational domain representing the three-dimensional anatomical model of the coronary arteries. Represents blood flow velocity field distribution data Conjugate momentum field distribution data. Blood flow velocity field distribution data. In symptotic manifolds The evolution on satisfies the Szymborsaurus-Lagrange equation:

[0066] ; in, Represents the Lagrange function, , Represents blood flow velocity field distribution data For time variables The first derivative.

[0067] The blood flow velocity field distribution data is analyzed using the symplectic geometric mode decomposition algorithm. Decomposed into Each symplectic geometric modal component and a residual component :

[0068] ; in The value ranges from 8 to 12. Each symplectic geometric modal component... The acquisition process is as follows: in the symptotic manifold The upper structure and the first The symptotic subspace corresponding to each mode Xinzi Space The basis vectors are obtained from the blood flow velocity field distribution data. covariance matrix It is obtained by performing symplectic eigenvalue decomposition. The mathematical definition of symplectic eigenvalue decomposition is solving the generalized symplectic eigenvalue problem. ,in For symplectic matrix , It is the identity matrix. It is a zero matrix. For the first Each symplectic eigenvalue, In order to be with the first The symplectic eigenvectors corresponding to each symplectic eigenvalue. Satisfying the symplectic orthogonality condition (when (Time). After the solution is completed, the symplectic eigenvectors are used as the basis for the calculation. Spanning a symplectic subspace of basis vectors Then, the blood flow velocity field distribution data Projected into symptotic space Above, obtain the first Each symplectic geometric modal component .

[0069] Symmetric geometric mode components Coronary computed tomography (CT) angiography data Correlation analysis was performed on the extracted vascular morphological features: energy distribution was calculated. ,in This represents the duration of a complete cardiac cycle, and then the energy distribution... Spatial correlation analysis was performed with vascular morphological features, and the Pearson correlation coefficient was calculated. ,when A significant correlation is determined when the absolute value is greater than 0.6. The symplectic geometric modal components are arranged from high to low frequency: high-frequency modal components correspond to local turbulence and vortex structures, mid-frequency modal components correspond to mainstream flow modes, and low-frequency modal components correspond to overall infusion characteristics.

[0070] Pressure field distribution data The symplectic geometric mode decomposition algorithm is also used for decomposition:

[0071] ; in Indicates the first Each pressure field has a symplectic geometric modal component. This represents the residual components of the pressure field. The symplectic geometric modal components of the velocity field are also represented. Symplectic geometric modal components of the pressure field According to the same modality index To pair up and obtain Group velocity-pressure mode pair.

[0072] Step S4: Extraction of hemodynamic parameters based on symplectic geometric mode decomposition results.

[0073] Based on the symplectic geometric modal components of the velocity field obtained in step S3 and the symplectic geometric modal components of the pressure field Hemodynamic parameters of each segment of the coronary artery vascular tree were extracted.

[0074] Based on 3D segmentation mask A skeletonization algorithm was used to extract the centerline of the coronary artery tree, and the coronary artery tree was divided into sections based on the bifurcation points of the centerline. Each blood vessel segment. From its proximal center point coordinates Coordinates of the far center point and centerline path The only certainty is that among them .

[0075] For each blood vessel segment The following hemodynamic parameters are calculated based on the symplectic geometric modal components of the velocity field.

[0076] The first hemodynamic parameter is the mean blood flow velocity. :

[0077] ; in Indicates the centerline path Length, Indicates the path along the centerline The arc length infinitesimal element.

[0078] The second hemodynamic parameter is the blood flow velocity pulsation intensity. :

[0079] ; in Indicates vascular segment internal voxel point The average blood flow velocity vector within a complete cardiac cycle.

[0080] The third hemodynamic parameter is the wall shear stress. :

[0081] ; in This is the unit vector of the direction of the external normal to the blood vessel wall. Using the nearest neighbor triangular facet Calculation: Let the coordinates of its three vertices be... , and Then the unit vector of the outward normal direction , Take as .

[0082] The following hemodynamic parameters are calculated based on the symplectic geometric modal components of the pressure field.

[0083] The fourth hemodynamic parameter is the pressure gradient across the vessel segment. :

[0084] ; in Indicates the first The pressure value of the symplectic geometric modal component of the pressure field at the near-end center point. Indicates the first The pressure value at the far-end center point of the symplectic geometric modal component of the pressure field.

[0085] The fifth hemodynamic parameter is the fractional flow reserve. :

[0086] ; Obtain each vascular segment hemodynamic parameter vector .

[0087] Step S5: Estimation of hemodynamic parameters based on Gaussian process regression of vascular network kernel function.

[0088] The hemodynamic parameters obtained in step S4 are processed using a Gaussian process regression algorithm based on the vascular network kernel function. Perform regression modeling and optimization estimation.

[0089] Define Gaussian process priors ,in This represents the input feature vector. Represents the mean function, Represents the kernel function. Mean function. The zero-mean prior is set as the zero function. The rationale for using a zero-mean prior is that, in the absence of prior information about the absolute values ​​of hemodynamic parameters, a zero-mean prior is a reasonable choice as the minimum assumption, ensuring that the Gaussian process regression model is dominated by the kernel function in its predictions, thus avoiding the introduction of human bias. Input feature vector From vascular segment Its anatomical and morphological features constitute:

[0090] ; in The length of the blood vessel segment. The average lumen diameter, Minimum lumen diameter, The maximum lumen diameter, For the mean curvature, For the branch angle, Percentage of narrow area.

[0091] Calculation method: along the blood vessel segment centerline The lumen cross-section was extracted with a sampling step size of 0.1 mm to calculate the lumen area. A "stable continuous region of lumen diameter" was identified in the lumen diameter sequence as a normal reference segment, with the coefficient of variation as the criterion. If the value is less than 0.05, calculate the average lumen diameter of the region as the diameter of the adjacent normal lumen. ,pass Obtain the reference normal lumen area. If or The area measured in square millimeters was determined to be a completely occlusive lesion. Set to 1.0; otherwise pass. calculate.

[0092] Define the vascular network kernel function :

[0093] ; in For the quadratic exponential kernel function:

[0094] ; in For signal variance parameters, It is a length-scaled diagonal matrix (7×7 dimensions). For the topological kernel function:

[0095] ; in Here is the variance parameter of the topological kernel function. For topological distance, is the length scale parameter of the topological kernel function. For hemodynamic kernel functions:

[0096] ; in The variance parameter of the hemodynamic kernel function. This refers to the parameter vector obtained through computational fluid dynamics simulation in step S4. is the length scale parameter of the hemodynamic kernel function.

[0097] Model training: Building the training dataset The hyperparameter set is optimized using the marginal likelihood maximization method. The conjugate gradient method is used for optimization. Initial values: , , The logarithmic form of the marginal likelihood function is:

[0098] ; in This is a column vector consisting of parameter vectors for all vessel segments. The matrix formed by all input feature vectors, It is the kernel matrix.

[0099] For the target vascular segment Its posterior distribution mean is:

[0100] ; The posterior distribution variance is:

[0101] ; Will and Obtain by weighted fusion :

[0102] ; Fusion weighting coefficient ,in Values to .

[0103] Step S6: Functional assessment of stenotic lesions based on fused hemodynamic parameters.

[0104] Based on the fused hemodynamic parameter vector obtained in step S5 For each segment of the coronary artery tree Perform functional assessment of the stenotic lesion. Functional assessment is based on the vector of hemodynamic parameters after fusion. Fractional blood flow reserve As the core indicator.

[0105] Define a grading standard for the functional severity of stenotic lesions. Based on fractional flow reserve. The numerical range, including the vascular segment Stenotic lesions are classified into the following four grades: Grade 1 is significant hemodynamic stenosis, when... Time to determine blood vessel segment There is significant hemodynamic stenosis; the second grade is hemodynamically critical stenosis, when Grade 1 is determined to have critical hemodynamic stenosis; Grade 3 is mild hemodynamic stenosis. The first level is considered to have mild hemodynamic stenosis; the fourth level is considered to have normal hemodynamics. The hemodynamic status was determined to be normal at that time.

[0106] For vascular segments identified as having significant hemodynamic stenosis, further extraction of the location information of the stenotic segment is required. Length information and minimum lumen diameter Location information of the narrow segment The method for determining the diameter is as follows: Within the identified continuous narrow segment region, search for the voxel point with the smallest lumen diameter, and use the spatial coordinates of this voxel point as... If multiple local minima exist with the same lumen diameter, then the arithmetic mean of the spatial coordinates of these local minima is taken as the minimum value. Length information of the narrow segment The method for determining this is: the diameter of adjacent normal lumens The value already calculated in step S5 can be used directly without recalculation; along the blood vessel segment centerline Extending from the narrowed region proximally and distally, select a continuous region with stable lumen diameter within a range of 3.0 mm to 5.0 mm from the edge of the narrowed region. The criteria for determining this "continuous region with stable lumen diameter" are the same as those defined in step S5. Obtain the diameter of adjacent normal lumens. Then, the length information of the narrow segment via the blood vessel segment centerline Measuring lumen diameter less than The length of the continuous region is determined.

[0107] Calculate the comprehensive hemodynamic severity score for each stenotic segment. :

[0108] ; in The weighting factor for the fractional blood flow reserve term is set to 0.5. The weighting coefficient for the narrow length term is set to 0.25. The weighting factor for the degree of narrowness of the lumen diameter is set to 0.25.

[0109] Obtain a collection of functional assessment results of stenosis lesions in all segments of the coronary artery tree. .

[0110] Step S7: Generation of multi-target distal anastomosis candidate locations based on hemodynamic and morphological features.

[0111] The set of functional assessment results of stenosis lesions obtained in step S6 For each target vessel segment requiring bypass grafting, a set of candidate locations for distal anastomosis is generated.

[0112] Define morphological constraints for the selection of distal anastomosis locations. For each vessel segment... Candidate locations of its distal anastomosis The following five morphological constraints must be met simultaneously.

[0113] The first morphological constraint is the lumen diameter constraint: candidate location of distal anastomosis. lumen diameter at It must be greater than or equal to 2.0 mm. This threshold is the minimum target vessel diameter recognized in coronary artery bypass grafting surgery to ensure long-term patency of the bypass graft.

[0114] The second morphological constraint is the landing area length constraint: based on the candidate location of the distal anastomosis. Centered on, along the blood vessel segment Centerline Within the landing area extending 5.0 mm forward and backward, the diameter of the tube is... It must always be greater than or equal to 1.5 mm. When the blood vessel segment... Total length When the diameter is less than 10 mm, adjust the position to the distal anastomosis candidate location. The total length of the vascular segment is the length extended from the centerline towards the anterior and distal ends. One-third of the diameter and not less than 1.0 mm. If If the diameter is less than 2 mm, skip that blood vessel segment.

[0115] The third morphological constraint is the calcification load constraint: candidate location of distal anastomosis. Calcification load index at the site The threshold must be less than or equal to 0.3. This threshold is determined based on clinical studies (referencing Kaiser et al., JACC 2022), which showed that when the coronary artery calcification burden index exceeds 0.3, the vessel wall compliance at the anastomosis suture site decreases by more than 40%, significantly increasing the risk of postoperative anastomotic leakage and restenosis. Calcification burden index Coronary computed tomography (CT) angiography image data obtained in step S1 mid voxel point CT value calculation at the location: ,in Represented by voxel points The number of voxels with a CT value greater than or equal to 130 Henlein units within a spherical neighborhood centered at a radius of 1.0 mm. This represents the total number of primes within the spherical neighborhood. 130 Henle units is an internationally recognized standard CT value threshold for identifying coronary artery calcification plaques.

[0116] The fourth morphological constraint is the bifurcation point distance constraint: candidate locations of distal anastomoses. The distance along the centerline from the nearest vascular bifurcation point should be greater than or equal to 5.0 mm to avoid the impact of blood flow disturbance in the bifurcation point area on the anastomosis.

[0117] The fifth morphological constraint is the vessel wall quality constraint: candidate location of distal anastomosis. Blood vessel wall thickness at the location The vessel wall thickness must be greater than or equal to 0.3 mm and less than or equal to 1.5 mm. Intravascular ultrasound imaging data of the coronary arteries obtained in step S1 Extract from.

[0118] For each voxel that satisfies the five morphological constraints described above, it is marked as a candidate location for distal anastomosis. The set of all candidate locations for distal anastomosis is then obtained. ,in For candidate position index and , This represents the total number of candidate locations for distal anastomosis. Indicates the first Spatial coordinates of candidate locations for distal anastomosis. If If the target vessel segment is not identified, the candidate location generation and subsequent scoring steps for that target vessel segment are skipped, and the target vessel segment is marked as "Unable to plan distal anastomosis". If all target vessel segments are marked as "Unable to plan distal anastomosis", the subsequent steps are terminated and a prompt message is output.

[0119] Define myocardial perfusion demand index The calculation is divided into two stages. The first stage: traversing all candidate positions. Calculate the CT value attenuation rate of the left ventricular myocardium in its distal blood supply area. The method for determining the boundary of the blood supply area is as follows: along the central line of the coronary artery tree from the voxel point... The position is traced point by point in the direction of distance, and for each centerline point on the tracing path... Calculate the centerline point The number of voxel points belonging to the left ventricular myocardial tissue within a radial neighborhood centered at a radius of 2.0 mm. The method for identifying "voxel points belonging to left ventricular myocardial tissue" is as follows: based on CT value determination, when the CT value of a voxel point is within the range of 20 to 60 HJ units, it is determined to belong to left ventricular myocardial tissue. Less than the threshold If the blood supply area boundary is exceeded, the index set of all voxel points belonging to the left ventricular myocardial tissue within the blood supply area boundary is compared with the current candidate position. Associated storage. Threshold. The individualized adjustment method is as follows: In the middle region of the left ventricular myocardium, three locations along the long axis of the left ventricle—the basal segment, the middle segment, and the apical segment—are selected respectively to measure the average number of myocardial tissue voxels in the radial neighborhood. ,Will Set as 50% and rounded down, when When the value is less than 10, the value is 3.

[0120] CT value attenuation rate The calculation method is as follows: using voxel points Starting from the centerline, trace along the centerline to the boundary of the blood supply area, recording each centerline point along the tracing path. The average CT value of the corresponding myocardial region The slope of the line fitted using the least squares method is... :

[0121] ; in Indicates from voxel point To the center line point Arc length along the centerline Indicates all The arithmetic mean, Indicates all The arithmetic mean.

[0122] After traversing all candidate positions, the determination is made. and .

[0123] Phase Two: For each candidate position ,pass Calculate its myocardial perfusion demand index. When At this point, the difference in CT value attenuation rate among all candidate locations is negligible, meaning that the myocardial perfusion demand in the distal blood supply area of ​​each candidate location is the same. Therefore, the myocardial perfusion demand index for all candidate locations is set to the median of the normalized interval, 0.5, to avoid bias towards any candidate location in multi-objective scoring. Post-definition and .

[0124] Constructing a multi-objective scoring function for candidate locations of distal anastomosis :

[0125] ; in Weighting coefficients for hemodynamic optimization objectives. Weighting coefficients for morphological adaptation targets. The weighting coefficients for the infusion optimization objective.

[0126] Hemodynamic optimization goals :

[0127] ; Morphological adaptation target :

[0128] ; Infusion optimization target :

[0129] ; Calculate the multi-objective score for each candidate position Sort by rating from highest to lowest, and select the top-rated... The candidate locations are selected as the preferred candidate locations for distal anastomosis. ,in The value ranges from 3 to 5.

[0130] Step S8: Virtual bridge vessel path planning based on preferred distal anastomosis candidate locations.

[0131] like Figure 3 As shown, the preferred distal anastomosis candidate location set obtained in step S7 Combined with the three-dimensional segmentation mask obtained in step S2 The three-dimensional anatomical model of the coronary arteries constructed in step S3 is used to identify each preferred candidate location for distal anastomosis. Plan the travel path of the virtual bridge blood vessels.

[0132] Define the set of origin locations of the bypass vessels. Including the origin of the left internal mammary artery Location of proximal anastomosis of the great saphenous vein and the origin of the radial artery Three sources. For each and the origin of each bridging vessel Planning a road from arrive Virtual bridge vascular pathway .

[0133] A path planning algorithm based on the fast traversal method is employed to search for the optimal path in the external space of the three-dimensional anatomical model of the coronary arteries. The computational domain boundary is defined with respect to the geometric center of the three-dimensional anatomical model of the coronary arteries. Centered on, along , , The bounding box region is extended outward by 20 mm in each axial direction. If the starting point or ending point is outside the bounding box region, the region is extended adaptively. After redefining the boundary of the computational domain, the equations of the process function need to be solved again to ensure that the arrival time field is updated within the expanded computational domain.

[0134] Cost function of path planning algorithm Cost based on path length Path curvature cost Path collision cost The weighted sum of the three terms constitutes:

[0135] ; in , , .

[0136] Path length cost ,in This represents the arc length infinitesimal element.

[0137] Path curvature cost ,in Indicates the position of the arc length on the path The curvature at that point.

[0138] Path collision cost ,in Represents points on the path The minimum Euclidean distance to the surface of the anatomical model. Millimeters.

[0139] Specific implementation of the fast travel method: defining the equation of the process. ,in Indicates the distance from the starting point to the spatial point. The shortest arrival time, The local velocity is given. The equations are solved numerically using an upwind difference scheme:

[0140] ; in , , Three-dimensional spatial meshes are respectively , , The step size along the axis is determined based on the voxel spatial resolution of the coronary computed tomography angiography image data. After solving, the optimal path is obtained by tracing backward from the endpoint along the descent direction of the arrival time gradient. For each The path with the minimum cost function is selected from the optimal paths corresponding to the origins of all bridging vessels. .

[0141] Obtain each preferred distal anastomosis candidate location The corresponding optimal virtual bridging vessel course Then, the morphological parameters of the virtual graft vessels were calculated: graft vessel length. For path Total arc length; average diameter of bridging vessels The determination is based on the type of graft vessel and the distal lumen diameter of the target vessel segment. When the graft vessel is the left internal mammary artery... When it is the great saphenous vein When it is the radial artery ; Bridge vessel anastomosis angle It is the angle between the tangent direction of the path at the endpoint and the tangent direction of the centerline of the target vessel segment at the endpoint.

[0142] Step S9: Hemodynamic reassessment and bypass surgery planning based on virtual bridging vessel implantation.

[0143] Based on the optimal virtual bridging vessel course obtained in step S8 The virtual bypass vessel model, along with its morphological parameters, is implanted into the three-dimensional coronary artery anatomical model constructed in step S3. The geometric parameters of the virtual bypass vessel model include: a circular cross-sectional shape, a wall thickness of 0.2 mm to 0.5 mm (determined according to the type of bypass vessel), and a length of... , diameter is .

[0144] Virtual bridging vessel model and coronary artery 3D anatomy model at the distal anastomosis location The geometric connection method at the location is: Centered along the centerline of the target vessel segment Smooth fusion is performed within a region extending 3.0 mm proximally and distally, using cubic spline interpolation to generate the data from... Diameter of the lumen to the target vessel segment The smooth diameter transition curve is obtained and the vessel wall surface mesh is regenerated.

[0145] The computational fluid dynamics simulation and symplectic geometric mode decomposition process in step S3 and the hemodynamic parameter extraction process in step S4 were repeatedly executed on the three-dimensional anatomical model of the coronary arteries after the virtual graft was implanted. During repeated execution, the boundary conditions were adjusted as follows: the inlet boundary conditions included the original left and right coronary artery openings, as well as the graft inlet at the origin of the graft. The graft inlet flow waveform was obtained through... (Reuse step S7 and) The corresponding blood supply area boundary results, due to ,and The corresponding set of voxel indices has been stored in step S7, based on spatial coordinates. (Matching can extract the voxel index set within the corresponding blood supply area boundary) and calculate the target average flow rate. , The value is set to 60 ml to 100 ml per minute per 100 grams of myocardial tissue, and then distributed to each phase of the cardiac cycle to generate physiological flow waveforms; the outlet boundary conditions include all existing distal outlets and the connection between the distal anastomosis of the bypass vessel and the target vessel segment (using flow coupling boundary conditions).

[0146] Obtain updated hemodynamic parameter vectors for each segment of the coronary artery vascular tree after virtual graft implantation. and the hemodynamic parameter vector of the virtual bridge vessel itself . Includes mean blood flow velocity Wall shear stress Cross-vascular pressure gradient and fractional blood flow reserve . The calculation method involves extracting the average pressure at the origin of the bypass vessel. Average pressure at the distal anastomosis site , .

[0147] Calculate the hemodynamic improvement performance index for each virtual bridging vessel planning scheme. :

[0148] ; in This is an index set of target vessel segments that require bypass grafting. The number of elements in it. The pre-implantation fusion blood flow reserve fraction. This is the updated value after implantation. If... If it is an empty set, then Set to 0 and terminate subsequent execution. If Then truncate to 1.0, if Then it is truncated to 0.

[0149] Simultaneously calculate the hemodynamic risk indicators of the virtual graft vessel itself. :

[0150] ; in Pascal, mmHg .

[0151] Comprehensive calculation of the overall score of the planning scheme ,in , The scheme with the highest overall score is selected as the optimal bypass graft surgical plan and output.

[0152] In some embodiments, after performing a hemodynamic reassessment after virtual graft implantation (i.e., step S9), the procedure further includes a step of “multi-scale blood flow complexity and energy stability assessment based on symplectic geometric modal components” for the distal anastomosis and the mid-segment of the graft. Suppose that after step S3, the symplectic geometric mode decomposition, in the virtual bridge vessel path The time-series signal of the velocity field was extracted from the distal anastomosis site to obtain the fluctuation function of the blood flow velocity amplitude over time at that observation point. :

[0153] ; in, For the first The velocity field has symplectic geometric modal components. To optimize the three-dimensional spatial coordinates of the distal anastomosis, This represents the Euclidean norm of a vector.

[0154] The first step is to calculate the multiscale sample entropy (MSE) of blood flow velocity to quantify the physiological complexity of blood flow at that location: For the original time series Perform discrete sampling to obtain a length of discrete sequence For the scale factor ( ,in Values ​​ranging from 10 to 20 represent different timescales of physiological rhythm levels, constructing coarse-grained sequences. :

[0155] ; For a given embedding dimension (This embodiment takes) ) and similarity tolerance (This embodiment takes) ,in (This represents the standard deviation of the calculated sequence), and the sample entropy of the coarse-grained sequence. :

[0156] ; in, This indicates that there are two sequences of length 1 in the sequence. The subsequences in similarity tolerance The probability of a match. Indicates length is The probability of a subsequence matching.

[0157] Define the blood flow complexity index at the anastomosis site. The sum of sample entropy under all scaling factors:

[0158] ; The second step is to calculate the Energy Gradient Dissipation (EGD) of the bypass vessel under pulsating blood flow, which is used to characterize the stability of energy loss when blood flow passes through the bypass vessel. Based on the pressure field symplectic geometric mode components obtained in step S3 The instantaneous pressure at the inlet (proximal end) and outlet (distal anastomosis) of the bypass graft is extracted. The instantaneous pressure difference between the two ends of the bypass graft is defined. and instantaneous flow :

[0159] ;

[0160] ; in, The spatial coordinates of the origin of the bridging vessel. This is the cross-sectional area at the inlet of the bridging vessel. Let be the outward normal direction vector of the area infinitesimal element.

[0161] Define the energy gradient dissipation index based on the pulsation period :

[0162] ; in, The duration of a complete heartbeat cycle. This index. The physical meaning of this value is: the degree of drastic change in instantaneous power (the product of pressure difference and flow rate) over time in pulsating blood flow. The larger the value, the more intense the energy oscillation of blood flow within the bypass graft, and the higher the risk of non-physiological turbulence.

[0163] The third step is to construct long-term patency prediction and assessment factors based on the aforementioned blood flow complexity index and energy dissipation index. :

[0164] ; in, The reference complexity index is the mean value obtained from statistics based on a large sample of normal coronary artery blood flow data (from public blood flow databases or normal control group CT-FFR data). An energy dissipation safety threshold set based on experience (in this embodiment, it is taken as...) ); and The weighting coefficients are (in this embodiment, they are respectively taken as) and ).

[0165] like The calculated value is lower than the preset safety lower limit threshold. (This embodiment takes) If the current virtual vascular graft scheme is deemed to have a high risk of long-term intimal hyperplasia or restenosis, then in the iterative optimization of step S10, the multi-objective scoring function of step S7 is forced. Add a blood flow stability penalty item Then, adaptive weight adjustment and path optimization are performed again until... Or it may reach the maximum number of iterations.

[0166] Step S10: Verification and optimization iteration of bypass surgery planning scheme.

[0167] The optimal bypass graft surgery plan generated in step S9 is visualized and verified on the three-dimensional anatomical model of the coronary artery. The virtual graft vessel's course is checked to see if it overlaps with the surface of the anatomical model, and whether the virtual graft vessel model is excessively distorted. The fractional flow reserve before and after implantation is displayed and compared using color mapping.

[0168] The optimal surgical plan is evaluated based on pre-defined evaluation criteria, including three criteria: the first evaluation criterion is the fractional flow reserve restoration criterion, which is applied after implantation of all target vessel segments. The second assessment criterion is the hemodynamic safety criterion of the bypass graft. Pascal The third evaluation criterion is the anatomical fit criterion, which is the minimum distance between the virtual bridging vessel path and the surface of the anatomical model. Millimeters and matching angle Spend.

[0169] like Figure 4 As shown, if the optimal solution simultaneously satisfies all three criteria, the output is the final solution. If any criterion is not met, the deviation information is fed back to the distal anastomosis candidate location generation process in step S7 to adjust the multi-objective scoring function. , , Adjust the weighting coefficients and re-execute steps S7, S8, and S9. Define an iteration counter. The initial value is 0, and it increases by 1 after each re-execution. achieve Stop iterating when the time comes.

[0170] Weighting coefficient adjustment method: If the first evaluation criterion is not met, increase the weighting coefficient. Adjust the formula to ,in For all target vessel segments that do not meet the criteria The average value, ; , Adjusted and normalized. If the second evaluation criterion is not met, increase... If the third evaluation criterion is not met, the increase will be... The adjustment methods are the same.

[0171] Iterate until all three criteria are met simultaneously or the maximum number of iterations is reached. If all criteria cannot be met after reaching the maximum number of iterations, output the solution with the highest overall score and include deviation information.

Claims

1. A method for coronary artery bypass graft planning based on digital imaging, characterized in that, include: Acquire multimodal coronary digital image data and perform preprocessing to obtain preprocessed image data; The preprocessed image data is subjected to three-dimensional segmentation of the coronary artery vascular tree to obtain a three-dimensional segmentation mask; A three-dimensional anatomical model of the coronary arteries is constructed based on the three-dimensional segmentation mask and preprocessed image data, and blood flow velocity field distribution data and pressure field distribution data within the coronary artery vascular tree are calculated based on the three-dimensional anatomical model of the coronary arteries. The symplectic geometric mode decomposition algorithm is used to perform multi-scale decomposition on the blood flow velocity field distribution data and pressure field distribution data to obtain the symplectic geometric mode components of the velocity field and the pressure field. Based on the symplectic geometric mode components of the velocity field and the pressure field, the hemodynamic parameter vectors of each segment of the coronary artery vascular tree are extracted. The hemodynamic parameter vector is optimized and estimated using the Gaussian process regression algorithm. The hemodynamic parameter vector is then fused with the optimized estimate obtained from the Gaussian process regression to obtain the fused hemodynamic parameter vector. Based on the fused hemodynamic parameter vector, a functional assessment of the stenosis lesion is performed to obtain a set of functional assessment results for the stenosis lesion. A set of candidate locations for distal anastomosis is generated based on the set of functional assessment results of the stenosis lesion, and a preferred set of candidate locations for distal anastomosis is determined from the set of candidate locations. The virtual bridge vessel course is planned based on the preferred set of distal anastomosis candidate locations. Based on the virtual bridge vessel's course, a re-evaluation of hemodynamics after virtual bridge vessel implantation is performed to obtain the updated hemodynamic parameter vectors of each vessel segment after virtual bridge vessel implantation and the hemodynamic parameter vectors of the virtual bridge vessel itself. Based on the updated hemodynamic parameter vector and the hemodynamic parameter vector of the virtual bypass vessel itself, each virtual bypass vessel planning scheme is comprehensively scored, and the virtual bypass vessel planning scheme with the highest comprehensive score is selected as the optimal bypass surgery planning scheme.

2. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The multimodal coronary digital imaging data includes coronary computed tomography angiography data, coronary intravascular ultrasound data, and coronary optical coherence tomography data.

3. The coronary artery bypass graft planning method based on digital imaging according to claim 2, characterized in that, The specific method for preprocessing the multimodal coronary digital imaging data is as follows: The multimodal coronary digital image data is denoised using a three-dimensional block matching denoising algorithm based on nonlocal means. A three-dimensional rigid registration algorithm based on mutual information is adopted, using the coronary computed tomography angiography image data as the fixed image and the coronary intravascular ultrasound image data or the coronary optical coherence tomography image data as the floating image, to register the coronary intravascular ultrasound image data and the coronary optical coherence tomography image data into the spatial coordinate system of the coronary computed tomography angiography image data; An enhancement algorithm based on the multi-scale retinal-cortical theory is used to enhance the registered image data.

4. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The three-dimensional segmentation of the coronary artery tree is achieved using a dynamic graph convolutional network, which includes a boundary integrity enhancement module and a trajectory continuity enhancement module. The boundary integrity enhancement module receives a set of multi-scale feature maps output by the encoder path of the dynamic graph convolutional network, calculates the probability value of each voxel point belonging to the blood vessel boundary through fuzzy logic reasoning, and feeds back the probability value to the corresponding layer of the decoder path of the dynamic graph convolutional network. The trajectory continuity enhancement module iteratively updates the initial graph structure based on a dynamic graph convolutional network. In each iteration, the node features are updated as follows: the feature vector of the current node is weighted and aggregated with the feature vectors of its neighboring nodes according to the reciprocal of the square root of the product of the degrees of the neighboring nodes. The aggregated result is added to the feature vector of the current node and then linearly transformed through a learnable weight matrix. After activation by a nonlinear activation function, the updated node features are output. The trajectory continuity enhancement module feeds back the updated node features from each iteration to the decoder path of the dynamic graph convolutional network.

5. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The training process of the dynamic graph convolutional network is as follows: A training dataset for a dynamic graph convolutional network is constructed. The training dataset contains multiple sets of training samples. Each set of training samples contains input multimodal coronary digital image data and its corresponding ground truth values ​​of three-dimensional segmentation annotations of the coronary artery vascular tree. The loss function of the dynamic graph convolutional network is defined as a composite loss function, which is composed of a weighted sum of a segmentation loss term, a boundary integrity loss term, and a trajectory continuity loss term. The segmentation loss term adopts a combination of the Dice loss function and the cross-entropy loss function; The boundary integrity loss term is modeled based on fuzzy logic, and the fuzzy membership degree of voxel points is calculated by the fuzzy C-means clustering algorithm. The dynamic graph convolutional network is trained using the backpropagation algorithm and the adaptive moment estimation optimizer until the loss function converges to below a preset threshold.

6. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The processing method of the symplectic geometric mode decomposition algorithm is as follows: The phase space where the blood flow velocity field distribution data is located is defined as a symplectic manifold, and a Hamiltonian function is defined on the symplectic manifold, wherein the Hamiltonian function represents the total energy of the blood flow velocity field distribution data. The evolution of the blood flow velocity field distribution data on the symplectic manifold satisfies the symplectic Euler-Lagrange equation; On the symplectic manifold, a symplectic subspace corresponding to each mode is constructed. The basis vectors of the symplectic subspace are obtained by performing symplectic eigenvalue decomposition on the covariance matrix of the blood flow velocity field distribution data. The symplectic eigenvalue decomposition is a problem of solving generalized symplectic eigenvalues. The blood flow velocity field distribution data is projected onto each symptotic subspace to obtain each symptotic geometric mode component; The symplectic geometric mode components of the velocity field and the symplectic geometric mode components of the pressure field are paired according to the same mode index to obtain multiple sets of velocity-pressure mode pairs.

7. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The Gaussian process regression algorithm uses a vascular network kernel function; The vascular network kernel function is composed of a weighted sum of three terms: a squared exponential kernel function, a topological kernel function, and a hemodynamic kernel function. The quadratic exponential kernel function calculates the similarity between different vascular segments based on the anatomical and morphological features of the vascular segments, including segment length, average lumen diameter, minimum lumen diameter, maximum lumen diameter, average curvature, branching angle, and percentage of stenosis area. The topological kernel function calculates topological correlation based on the topological distance between different vascular segments, where the topological distance is the number of nodes traversed along the central line of the coronary artery tree from one vascular segment to another. The hemodynamic kernel function calculates hemodynamic similarity based on the Euclidean distance between hemodynamic parameter vectors obtained through computational fluid dynamics simulation between different vascular segments.

8. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The specific method for determining the preferred set of distal anastomosis candidate locations from the candidate location set is as follows: Define morphological constraints for the selection of distal anastomosis location, including lumen diameter constraints, landing zone length constraints, calcification load constraints, bifurcation point distance constraints, and vessel wall quality constraints; A multi-objective scoring function for candidate locations of distal anastomosis is constructed. The multi-objective scoring function is composed of a weighted sum of hemodynamic optimization objective, morphological adaptation objective and perfusion optimization objective. The hemodynamic optimization objective is calculated based on the fused hemodynamic parameter vector; the morphological adaptation objective is calculated based on the local morphological features at the candidate location; and the perfusion optimization objective is calculated based on the myocardial perfusion demand index of the distal blood supply area at the candidate location. Calculate the multi-objective score for each candidate location of distal anastomosis, sort them from high to low score, and select the top N candidate locations with the highest scores as the preferred set of candidate locations for distal anastomosis.

9. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, The specific method for planning the virtual bridge blood vessel pathway is as follows: Define a set of origin locations for the bypass vessels, including the origin of the left internal mammary artery, the proximal anastomosis location of the great saphenous vein, and the origin of the radial artery. For each preferred distal anastomosis candidate location and each graft origin location, a path planning algorithm based on the fast travel method is used to plan a virtual graft travel path from the graft origin location to the preferred distal anastomosis candidate location. The cost function of the path planning algorithm is a weighted sum of path length cost, path curvature cost, and path collision cost. For each preferred distal anastomosis candidate location, the path with the minimum cost function is selected from all virtual graft paths corresponding to the origin of the graft vessels as the optimal virtual graft path for that preferred distal anastomosis candidate location.

10. The coronary artery bypass graft planning method based on digital imaging according to claim 1, characterized in that, After generating the optimal bypass surgery plan, the following verification and optimization iteration steps are also included: The optimal bypass surgery plan was visualized and validated on the three-dimensional anatomical model of the coronary arteries. The optimal bypass surgery plan is evaluated based on preset surgical planning evaluation criteria, which include the fractional flow reserve restoration criterion, the bypass vessel hemodynamic safety criterion, and the anatomical adaptation criterion. When the optimal bypass surgery planning scheme simultaneously satisfies the fractional flow reserve restoration criterion, the bypass vessel hemodynamic safety criterion, and the anatomical adaptation criterion, the scheme is output as the final bypass surgery planning scheme. When the optimal bypass surgery plan does not meet any of the fractional flow reserve recovery criterion, the bypass vessel hemodynamic safety criterion, and the anatomical adaptation criterion, the deviation information of the unmet evaluation criterion is fed back to the distal anastomosis candidate location generation step. The weight coefficients of each optimization objective in the multi-objective scoring function are adjusted, and the distal anastomosis candidate location generation, virtual bypass vessel path planning, and bypass surgery plan generation steps are re-executed to generate an updated bypass surgery plan. Repeat the above iterative optimization process until the generated bypass transplantation surgical plan simultaneously meets all evaluation criteria, or the number of iterations reaches the preset maximum number of iterations; If a bypass surgery plan that simultaneously meets all evaluation criteria cannot be generated after the number of iterations reaches the preset maximum number of iterations, the bypass surgery plan with the highest comprehensive score will be output, along with deviation information for the evaluation criteria that are not met.

Images